G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
gsXBraidMultigrid.h
1#include <gismo.h>
2#include <string>
3
4namespace gismo {
5
11 template<class T>
12 struct gsXBraidMultigridBase
13 {
14 protected:
15 int maxIter;
16 int numLevels;
17 int numSmoothing;
18 int typeBCHandling;
19 int typeCycle_h;
20 int typeCycle_p;
21 int typeLumping;
22 int typeProjection;
23 int typeSmoother;
24 gsMatrix<> hp;
25 T tol;
26
27 public:
29 gsXBraidMultigridBase()
30 : maxIter(100000),
31 numLevels(1),
32 numSmoothing(1),
33 typeBCHandling(1),
34 typeCycle_h(2),
35 typeCycle_p(1),
36 typeLumping(1),
37 typeProjection(1),
38 typeSmoother(1),
39 tol(1e-8)
40 {}
41
42 void setMaxIter(int maxIter)
43 { this->maxIter = maxIter; }
44
45 void setTolerance(T tol)
46 { this->tol = tol; }
47
48 void setNumLevels(int numLevels, int typeProjection, int numDegree)
49 {
50 if(typeProjection == 1)
51 {
52 this->numLevels = numLevels - numDegree + 2;
53 }
54 else
55 {
56 this->numLevels = numLevels;
57 }
58 }
59
60 void setNumSmoothing(int numSmoothing)
61 { this->numSmoothing = numSmoothing; }
62
63 void setTypeBCHandling(int typeBCHandling)
64 { this->typeBCHandling = typeBCHandling; }
65
66 void setTypeCycle_h(int typeCycle_h)
67 { this->typeCycle_h = typeCycle_h; }
68
69 void setTypeCycle_p(int typeCycle_p)
70 { this->typeCycle_p = typeCycle_p; }
71
72 void setTypeLumping(int typeLumping)
73 { this->typeLumping = typeLumping; }
74
75 void setTypeProjection(int typeProjection)
76 { this->typeProjection = typeProjection; }
77
78 void setTypeSmoother(int typeSmoother)
79 { this->typeSmoother = typeSmoother; }
80
81 void setCoarsening(gsMatrix<> hp)
82 { this->hp = hp; }
83
84 virtual gsXBraidMultigridBase& compute(const gsSparseMatrix<T>& mat, const T tstep, const int& numDegree, index_t typeMethod)
85 {
86 // Get arguments explicitly
87 gsMatrix<T> x = gsMatrix<>::Zero(mat.rows(),1);
88 gsMatrix<T> b = gsMatrix<>::Zero(mat.rows(),1);
89 gsFunctionExpr<> rhs("1",2);
90 int iterTot = 1;
91 int typeMultigrid = 2;
92 int typeCoarseOperator = 1;
93
95 setup(rhs,
96 x,
97 b,
98 iterTot,
99 numLevels,
100 numDegree,
101 typeMultigrid,
102 hp,
103 typeCoarseOperator,
104 tstep,
105 typeMethod);
106
107
108 return *this; }
109
110 virtual gsMatrix<T> solveWithGuess(const gsMatrix<T>& b,
111 const gsMatrix<T>& x0)
112 {
113 // Get arguments explicitly
114 gsMatrix<T> x(x0);
115 x = x0;
116
117 gsFunctionExpr<> rhs("1",2);
118 int iterTot = 1;
119 int typeMultigrid = 2;
120 int typeCoarseOperator = 1;
121
123 solve(rhs,
124 x,
125 b,
126 iterTot,
127 numLevels,
128 typeMultigrid,
129 hp,
130 typeCoarseOperator);
131 return x;
132 }
133
135 virtual void solveMG(const gsMatrix<T> & rhs,
136 std::vector<memory::shared_ptr<gsMultiBasis<T> > > m_bases,
137 gsMatrix<T>& x,
138 const int& numLevels,
139 gsBoundaryConditions<T> bcInfo,
140 gsMultiPatch<T> mp,
141 std::vector<gsSparseMatrix<T> >& m_prolongation_P,
142 std::vector<gsSparseMatrix<T> >& m_restriction_P,
143 std::vector<gsMatrix<T> >& m_prolongation_M,
144 std::vector<gsMatrix<T> >& m_restriction_M,
145 std::vector<gsSparseMatrix<T> >& m_prolongation_H,
146 std::vector<gsSparseMatrix<T> >& m_restriction_H,
147 const gsMatrix<T>& hp)
148 {
149 if ( numLevels == 1)
150 {
151 solvecoarse(rhs, x, numLevels);
152 return;
153 }
154
155 if (hp(std::max(numLevels-2,0),0) == 0 )
156 {
157 gsMatrix<T> fineRes, coarseRes, fineCorr, coarseCorr, postRes;
158 presmoothing(rhs, x, numLevels, fineRes, hp);
159 restriction(fineRes, coarseRes, numLevels, m_bases,
160 bcInfo, mp,
161 m_prolongation_P, m_restriction_P,
162 m_prolongation_M, m_restriction_M,
163 m_prolongation_H, m_restriction_H, hp);
164 //coarseRes.setZero(coarseRes.rows(),1);
165 coarseCorr.setZero(coarseRes.rows(),1);
166 for( int j = 0 ; j < (typeCycle_p == 2 ? 2 : 1) ; j++)
167 {
168 solveMG(coarseRes, m_bases, coarseCorr, numLevels-1,
169 bcInfo, mp,
170 m_prolongation_P, m_restriction_P,
171 m_prolongation_M, m_restriction_M,
172 m_prolongation_H, m_restriction_H, hp);
173 }
174 prolongation(coarseCorr, fineCorr, numLevels, m_bases,
175 bcInfo, mp,
176 m_prolongation_P, m_restriction_P,
177 m_prolongation_M, m_restriction_M,
178 m_prolongation_H, m_restriction_H, hp);
179 postsmoothing(rhs, x, numLevels, fineCorr, postRes,
180 hp);
181 }
182
183 if (hp(std::max(numLevels-2,0),0) == 1 )
184 {
185 gsMatrix<T> fineRes, coarseRes, fineCorr, coarseCorr, postRes;
186 presmoothing(rhs, x, numLevels, fineRes, hp);
187 restriction(fineRes, coarseRes, numLevels, m_bases,
188 bcInfo, mp,
189 m_prolongation_P, m_restriction_P,
190 m_prolongation_M, m_restriction_M,
191 m_prolongation_H, m_restriction_H, hp);
192 //coarseRes.setZero(coarseRes.rows(),1);
193 coarseCorr.setZero(coarseRes.rows(),1);
194 for( int i = 0 ; i < (typeCycle_h == 2 ? 2 : 1) ; i++)
195 {
196 solveMG(coarseRes, m_bases, coarseCorr, numLevels-1,
197 bcInfo, mp,
198 m_prolongation_P, m_restriction_P,
199 m_prolongation_M, m_restriction_M,
200 m_prolongation_H, m_restriction_H, hp);
201 }
202 prolongation(coarseCorr, fineCorr, numLevels, m_bases,
203 bcInfo, mp,
204 m_prolongation_P, m_restriction_P,
205 m_prolongation_M, m_restriction_M,
206 m_prolongation_H, m_restriction_H, hp);
207 postsmoothing(rhs,x, numLevels, fineCorr, postRes,
208 hp);
209 }
210 }
211
212 virtual void setup(const gsFunctionExpr<T> & rhs,
213 gsMatrix<T>& x,
214 gsMatrix<T> f,
215 const int& iterTot,
216 const int& numLevels,
217 const int& numDegree,
218 const int& typeMultigrid,
219 const gsMatrix<T>& hp,
220 const int& typeCoarseOperator,
221 T tstep,
222 index_t typeMethod){}
223
224 virtual void solve(const gsFunctionExpr<T> & rhs,
225 gsMatrix<T>& x,
226 gsMatrix<T> f,
227 const int& iterTot,
228 const int& numLevels,
229 const int& typeMultigrid,
230 const gsMatrix<T>& hp,
231 const int& typeCoarseOperator){}
232
234 virtual void presmoothing(const gsMatrix<T>& rhs,
235 gsMatrix<T>& x,
236 const int& numLevels,
237 gsMatrix<T> & fineRes ,
238 const gsMatrix<T>& hp) = 0;
239
241 virtual void postsmoothing(const gsMatrix<T>& rhs,
242 gsMatrix<T>& x,
243 const int& numLevels,
244 gsMatrix<T> & fineCorr,
245 gsMatrix<T> & postRes,
246 const gsMatrix<T>& hp) = 0;
247
249 virtual void solvecoarse(const gsMatrix<T>& rhs,
250 gsMatrix<T>& x,
251 const int& numLevels) = 0;
252
254 virtual gsSparseMatrix<T> prolongation_P(const int& numLevels,
255 std::vector<memory::shared_ptr<gsMultiBasis<T> > > m_bases) = 0;
256
258 virtual gsSparseMatrix<T> restriction_P(const int& numLevels,
259 std::vector<memory::shared_ptr<gsMultiBasis<T> > > m_bases) = 0;
260
262 virtual gsMatrix<T> prolongation_M(const int& numLevels,
263 std::vector<memory::shared_ptr<gsMultiBasis<T> > > m_bases) = 0;
264
266 virtual gsMatrix<T> restriction_M(const int& numLevels,
267 std::vector<memory::shared_ptr<gsMultiBasis<T> > > m_bases) = 0;
268
270 virtual void prolongation(const gsMatrix<T>& Xcoarse,
271 gsMatrix<T>& Xfine,
272 const int& numLevels,
273 std::vector<memory::shared_ptr<gsMultiBasis<T> > > m_bases,
274 gsBoundaryConditions<T> bcInfo,
275 gsMultiPatch<T> mp,
276 std::vector<gsSparseMatrix<T> >& m_prolongation_P,
277 std::vector<gsSparseMatrix<T> >& m_restriction_P,
278 std::vector<gsMatrix<T> >& m_prolongation_M,
279 std::vector<gsMatrix<T> >& m_restriction_M,
280 std::vector<gsSparseMatrix<T> >& m_prolongation_H,
281 std::vector<gsSparseMatrix<T> >& m_restriction_H,
282 const gsMatrix<T>& hp)
283 {
284 if (hp(numLevels-2,0) == 1)
285 {
286 Xfine = m_prolongation_H[numLevels-2]*Xcoarse;
287 }
288 else
289 {
290 if (typeLumping == 1)
291 {
292 gsMatrix<T> temp = m_prolongation_P[numLevels-2]*Xcoarse;
293 gsMatrix<T> M_L_inv = (m_prolongation_M[numLevels-2]).array().inverse();
294 Xfine = (M_L_inv).cwiseProduct(temp);
295 }
296 else
297 {
298 // Define the low and high order bases
299 gsMultiBasis<T> basesL = *m_bases[numLevels-2];
300 gsMultiBasis<T> basesH = *m_bases[numLevels-1];
301 typedef gsExprAssembler<real_t>::geometryMap geometryMap;
302 typedef gsExprAssembler<real_t>::variable variable;
303 typedef gsExprAssembler<real_t>::space space;
304
305 // Determine matrix M (high_order * high_order)
306 gsExprAssembler<real_t> ex2(1,1);
307 geometryMap G2 = ex2.getMap(mp);
308 space w_n = ex2.getSpace(basesH ,1, 0);
309 w_n.setInterfaceCont(0);
310 if (typeBCHandling == 1)
311 {
312 w_n.setup(bcInfo, dirichlet::l2Projection, 0);
313 //#w_n.addBc(bcInfo.get("Dirichlet"));
314 }
315 ex2.setIntegrationElements(basesH);
316 ex2.initSystem();
317 ex2.assemble(w_n * meas(G2) * w_n.tr());
318
319 // Prolongate Xcoarse to Xfine
320 gsMatrix<T> temp = m_prolongation_P[numLevels-2]*Xcoarse;
321 gsSparseMatrix<T> M = ex2.matrix();
322 gsConjugateGradient<T> CGSolver(M);
323 CGSolver.setTolerance(1e-12);
324 CGSolver.solve(temp,Xfine);
325 }
326 }
327 }
328
330 virtual void restriction(const gsMatrix<T>& Xfine,
331 gsMatrix<T>& Xcoarse,
332 const int& numLevels,
333 std::vector<memory::shared_ptr<gsMultiBasis<T> > > m_bases,
334 gsBoundaryConditions<T> bcInfo,
335 gsMultiPatch<T> mp,
336 std::vector<gsSparseMatrix<T> >& m_prolongation_P,
337 std::vector<gsSparseMatrix<T> >& m_restriction_P,
338 std::vector<gsMatrix<T> >& m_prolongation_M,
339 std::vector<gsMatrix<T> >& m_restriction_M,
340 std::vector<gsSparseMatrix<T> >& m_prolongation_H,
341 std::vector<gsSparseMatrix<T> >& m_restriction_H,
342 const gsMatrix<T>& hp)
343 {
344 if (hp(numLevels-2,0) == 1)
345 {
346 Xcoarse = m_restriction_H[numLevels-2]*Xfine;
347 }
348 else
349 {
350 if (typeLumping == 1)
351 {
352 // Standard way
353 gsMatrix<T> temp = m_restriction_P[numLevels-2]*Xfine;
354 gsMatrix<T> M_L_inv = (m_restriction_M[numLevels-2]).array().inverse();
355 Xcoarse = (M_L_inv).cwiseProduct(temp);
356 }
357 else
358 {
359 // Define the low and high order bases
360 gsMultiBasis<T> basesL = *m_bases[numLevels-2];
361 gsMultiBasis<T> basesH = *m_bases[numLevels-1];
362 typedef gsExprAssembler<real_t>::geometryMap geometryMap;
363 typedef gsExprAssembler<real_t>::variable variable;
364 typedef gsExprAssembler<real_t>::space space;
365
366 // Determine matrix M (low_order * low_order)
367 gsExprAssembler<real_t> ex2(1,1);
368 geometryMap G2 = ex2.getMap(mp);
369 space w_n = ex2.getSpace(basesL, 1, 0);
370 w_n.setInterfaceCont(0);
371 if (typeBCHandling == 1)
372 {
373 w_n.setup(bcInfo, dirichlet::l2Projection, 0);
374 //#w_n.addBc(bcInfo.get("Dirichlet"));
375 }
376 ex2.setIntegrationElements(basesL);
377 ex2.initSystem();
378 ex2.assemble(w_n * meas(G2) * w_n.tr());
379
380 // Restrict Xfine to Xcoarse
381 gsMatrix<T> temp = m_restriction_P[numLevels-2]*Xfine;
382 gsSparseMatrix<T> M = ex2.matrix();
383 gsConjugateGradient<T> CGSolver(M);
384 CGSolver.setTolerance(1e-12);
385 CGSolver.solve(temp, Xcoarse);
386 }
387 }
388 }
389 };
390
399 template<class T, class CoarseSolver>
400 struct gsXBraidMultigrid : public gsXBraidMultigridBase<T>
401 {
402 private:
403
405 typedef gsXBraidMultigridBase<T> Base;
406
408 memory::shared_ptr<gsMultiPatch<T> > m_mp_ptr;
409
411 memory::shared_ptr<gsBoundaryConditions<T> > m_bcInfo_ptr;
412
414 std::vector<memory::shared_ptr<gsMultiBasis<T> > > m_bases;
415
417 std::vector< gsSparseMatrix<T> > m_prolongation_P;
418
420 std::vector< gsSparseMatrix<T> > m_restriction_P;
421
423 std::vector< gsMatrix<T> > m_prolongation_M;
424
426 std::vector< gsMatrix<T> > m_restriction_M;
427
429 std::vector< gsSparseMatrix<T> > m_prolongation_H;
430
432 std::vector< gsSparseMatrix<T> > m_restriction_H;
433
435 std::vector< std::vector< gsSparseMatrix<T> > > m_ILUT;
436
438 std::vector< std::vector < gsEigen::PermutationMatrix<Dynamic,Dynamic,index_t> > > m_P;
439
441 std::vector < std::vector < gsEigen::PermutationMatrix<Dynamic,Dynamic,index_t> > > m_Pinv;
442
444 std::vector< typename gsPreconditionerOp<T>::Ptr > m_SCMS;
445
447 std::vector< gsSparseMatrix<T> > m_operator;
448
450 std::vector < std::vector< gsSparseMatrix<T> > > m_block_operator;
451
453 std::vector < std::vector < gsSparseMatrix<T> > > m_ddB;
454
456 std::vector < std::vector < gsSparseMatrix<T> > > m_ddC;
457
459 std::vector < std::vector < gsMatrix<T> > > m_ddBtilde;
460
462 std::vector < std::vector < gsMatrix<T> > > m_ddCtilde;
463
465 std::vector < gsMatrix<T> > m_A_aprox;
466
468 std::vector < gsSparseMatrix<T> > m_S;
469
471 std::vector < std::vector< int > > m_shift;
472
473 public:
474
475 // Constructor
476 gsXBraidMultigrid(const gsMultiPatch<T> & mp,
477 const gsMultiBasis<T> & bases,
478 const gsBoundaryConditions<T> & bcInfo)
479 {
480 m_mp_ptr = memory::make_shared_not_owned(&mp);
481 m_bcInfo_ptr = memory::make_shared_not_owned(&bcInfo);
482 m_bases.push_back(memory::make_shared_not_owned(&bases));
483 }
484
485 virtual ~gsXBraidMultigrid() {}
486
487 public:
488
490 void setup(const gsFunctionExpr<T> & rhs,
491 gsMatrix<T>& x,
492 gsMatrix<T> f,
493 const int& iterTot,
494 const int& numLevels,
495 const int& numDegree,
496 const int& typeMultigrid,
497 const gsMatrix<T>& hp,
498 const int& typeCoarseOperator,
499 T tstep,
500 index_t typeMethod)
501 {
502 for (int i = 1; i < numLevels; i++)
503 {
504 m_bases.push_back(give(m_bases.back()->clone()));
505 switch((int) hp(i-1,0) )
506 {
507 case 0 : (Base::typeProjection == 1 ?
508 m_bases.back()->degreeIncrease(numDegree-1) :
509 m_bases.back()->degreeIncrease()); break;
510
511 case 1 : m_bases.back()->uniformRefine(); break;
512
513 case 2: m_bases.back()->uniformRefine();
514 m_bases.back()->degreeIncrease(); break;
515 }
516 }
517
518 // Generate sequence of matrix K and M
519 m_operator.resize(numLevels);
520 gsStopwatch clock;
521 //gsInfo << "|| Multigrid hierarchy ||" <<std::endl;
522
523 for (int i = 0; i < numLevels; i++)
524 {
525 //gsInfo << "Level " << i+1 << " " ;
526 //gsInfo << "Degree: " << m_bases[i]->degree() << ", Ndof: " << m_bases[i]->totalSize() <<std::endl;
527
528 typedef typename gsExprAssembler<T>::geometryMap geometryMap;
529 typedef typename gsExprAssembler<T>::variable variable;
530 typedef typename gsExprAssembler<T>::space space;
531 typedef typename gsExprAssembler<T>::solution solution;
532
533 gsExprAssembler<T> K, M;
534
535 // Set the bases
536 K.setIntegrationElements(*m_bases[i]);
537 M.setIntegrationElements(*m_bases[i]);
538
539 // Set the geometry map
540 geometryMap G_K = K.getMap(*m_mp_ptr);
541 geometryMap G_M = M.getMap(*m_mp_ptr);
542
543 // Set the discretization space
544 space u_K = K.getSpace(*m_bases[i]);
545 space u_M = M.getSpace(*m_bases[i]);
546 u_K.setInterfaceCont(0);
547 u_M.setInterfaceCont(0);
548 u_K.setup(*m_bcInfo_ptr, dirichlet::l2Projection, 0);
549 u_M.setup(*m_bcInfo_ptr, dirichlet::l2Projection, 0);
550 //#u_K.addBc( m_bcInfo_ptr->get("Dirichlet") );
551 //#u_M.addBc( m_bcInfo_ptr->get("Dirichlet") );
552
553 // Set the source term
554 auto ff_K = K.getCoeff(rhs, G_K);
555 auto ff_M = M.getCoeff(rhs, G_M);
556
557 // Initialize and assemble the system matrix
558 K.initSystem();
559 K.assemble( igrad(u_K, G_K) * igrad(u_K, G_K).tr() * meas(G_K), u_K * ff_K * meas(G_K) );
560
561 // Initialize and assemble the mass matrix
562 M.initSystem();
563 M.assemble( u_M * u_M.tr() * meas(G_M), u_M * ff_M * meas(G_M) );
564
565
566 m_operator[i] = M.matrix() + tstep*K.matrix();
567 switch(typeMethod)
568 {
569 case 0: m_operator[i] = M.matrix(); break;
570 case 1: m_operator[i] = M.matrix() + tstep*K.matrix(); break;
571 case 2: m_operator[i] = M.matrix() + 0.5*tstep*K.matrix();
572 }
573
574 }
575 real_t Time_Assembly = clock.stop();
576 GISMO_UNUSED(Time_Assembly);
577
578
579 // Resize vector of operators
580 m_prolongation_P.resize(numLevels-1);
581 m_prolongation_M.resize(numLevels-1);
582 m_prolongation_H.resize(numLevels-1);
583 m_restriction_P.resize(numLevels-1);
584 m_restriction_M.resize(numLevels-1);
585 m_restriction_H.resize(numLevels-1);
586
587 // Determine prolongation/restriction operators in p
588 clock.restart();
589 for (int i = 1; i < numLevels; i++)
590 {
591 if (hp(i-1,0) == 0)
592 {
593 m_prolongation_P[i-1] = prolongation_P(i+1, m_bases);
594 m_restriction_P[i-1] = m_prolongation_P[i-1].transpose(); //restriction_P(i+1, m_bases);
595 m_prolongation_M[i-1] = prolongation_M(i+1, m_bases);
596 m_restriction_M[i-1] = restriction_M(i+1, m_bases);
597 }
598 }
599
600 // Determine prolongation/restriction operators in h
601 gsSparseMatrix<real_t, RowMajor> transferMatrix;
602 gsOptionList options;
603 Base::typeBCHandling == 1 ? options.addInt("DirichletStrategy","",dirichlet::elimination) : options.addInt("DirichletStrategy","",dirichlet::nitsche);
604 for(int i = 1; i < numLevels; i++)
605 {
606 if (hp(i-1,0) == 1)
607 {
608 gsMultiBasis<T> m_bases_copy = *m_bases[i];
609 m_bases_copy.uniformCoarsen_withTransfer(transferMatrix,*m_bcInfo_ptr,options);
610 m_prolongation_H[i-1] = transferMatrix;
611 m_restriction_H[i-1] = m_prolongation_H[i-1].transpose();
612 }
613 }
614 real_t Time_Transfer = clock.stop();
615 GISMO_UNUSED(Time_Transfer);
616
617 // Obtain operators with Galerkin projection (TO DO)
618 clock.restart();
619 if (typeCoarseOperator == 2)
620 {
621 for (int i = numLevels-1; i > -1; i--)
622 {
623 if (hp(hp.rows()-1,0) == 0)
624 {
625 if (hp(std::min(i,hp.rows()-1),0) == 1)
626 {
627 m_operator[i] = m_restriction_H[i]*m_operator[i+1]*m_prolongation_H[i];
628 }
629 }
630 else
631 {
632 if (hp(std::min(i,hp.rows()-1),0) == 1 && i > 0)
633 {
634 m_operator[i-1] = m_restriction_H[i-1]*m_operator[i]*m_prolongation_H[i-1];
635 }
636 }
637 }
638 }
639 real_t Time_Assembly_Galerkin = clock.stop();
640 GISMO_UNUSED(Time_Assembly_Galerkin);
641
642 // Setting up the subspace corrected mass smoother
643 clock.restart();
644 if (Base::typeSmoother == 3)
645 {
646 // Generate sequence of SCM smoothers
647 m_SCMS.resize(numLevels);
648 gsOptionList opt;
649 opt.addReal("Scaling","",0.12);
650 for(int i = 0 ; i < numLevels ; i++)
651 {
652 m_SCMS[i] = setupSubspaceCorrectedMassSmoother(m_operator[i], *m_bases[i], *m_bcInfo_ptr, opt, Base::typeBCHandling);
653 }
654 }
655 real_t Time_SCMS = clock.stop();
656 GISMO_UNUSED(Time_SCMS);
657
658 // Determine ILUT factorizations at each level
659 clock.restart();
660 int numPatch = m_mp_ptr->nPatches();
661
662 if (Base::typeSmoother == 1)
663 {
664 // Generate factorizations (ILUT)
665 m_ILUT.resize(numLevels);
666 m_P.resize(numLevels);
667 m_Pinv.resize(numLevels);
668 for(int i = 0; i < numLevels; i++)
669 {
670 m_ILUT[i].resize(1);
671 m_P[i].resize(1);
672 m_Pinv[i].resize(1);
673 if (Base::typeProjection == 2)
674 {
675 gsEigen::IncompleteLUT<real_t> ilu;
676 ilu.setFillfactor(1);
677 ilu.compute(m_operator[i]);
678 m_ILUT[i][0] = ilu.factors();
679 m_P[i][0] = ilu.fillReducingPermutation();
680 m_Pinv[i][0] = ilu.inversePermutation();
681 }
682 else
683 {
684 if (i == numLevels-1) // Only at finest level
685 {
686 gsEigen::IncompleteLUT<real_t> ilu;
687 ilu.setFillfactor(1);
688 ilu.compute(m_operator[i]);
689 m_ILUT[i][0] = ilu.factors();
690 m_P[i][0] = ilu.fillReducingPermutation();
691 m_Pinv[i][0] = ilu.inversePermutation();
692 }
693 }
694 }
695 }
696 real_t Time_ILUT_Factorization = clock.stop();
697 GISMO_UNUSED(Time_ILUT_Factorization);
698
699 clock.restart();
700 if (Base::typeSmoother == 5)
701 {
702 int shift0 = 0;
703 m_ddB.resize(numLevels);
704 m_ddC.resize(numLevels);
705 m_ddBtilde.resize(numLevels);
706 m_ddCtilde.resize(numLevels);
707
708 m_ILUT.resize(numLevels);
709 m_P.resize(numLevels);
710 m_Pinv.resize(numLevels);
711 m_shift.resize(numLevels);
712 m_S.resize(numLevels);
713
714 for(int i = 0 ; i < numLevels ; i++)
715 {
716 m_shift[i].resize(numPatch+1);
717 m_ILUT[i].resize(numPatch+1);
718 m_P[i].resize(numPatch+1);
719 m_Pinv[i].resize(numPatch+1);
720
721 // Use of partition functions
722 std::vector<gsVector<index_t> > interior, boundary;
723 std::vector<std::vector<gsVector<index_t> > > interface;
724 std::vector<gsMatrix<index_t> > global_interior, global_boundary;
725 std::vector<std::vector<gsMatrix<index_t> > > global_interface;
726 //m_bases[i]->partition(interior,boundary,interface,global_interior,global_boundary,global_interface);
727 for(int l=0; l< numPatch; l++)
728 {
729 m_shift[i][l] = global_interior[l].rows();
730 }
731 m_shift[i][numPatch] = 0;
732 m_shift[i][numPatch] = m_operator[i].rows() - accumulate(m_shift[i].begin(),m_shift[i].end(),0);
733
734 // Put shift on zero
735 shift0 = 0;
736 for(int j = 0 ; j < numPatch ; j++)
737 {
738 const gsSparseMatrix<T> block = m_operator[i].block(shift0,shift0,m_shift[i][j],m_shift[i][j]);
739 gsEigen::IncompleteLUT<real_t> ilu;
740 ilu.setFillfactor(1);
741 ilu.compute(block);
742 m_ILUT[i][j] = ilu.factors();
743
744 m_P[i][j] = ilu.fillReducingPermutation();
745 m_Pinv[i][j] = ilu.inversePermutation();
746 shift0 = shift0 + m_shift[i][j];
747
748 }
749
750 shift0 = 0;
751 // Obtain the blocks of the matrix
752 m_ddB[i].resize(numPatch+1);
753 m_ddC[i].resize(numPatch+1);
754
755 for(int j = 0 ; j < numPatch+1 ; j++)
756 {
757 m_ddB[i][j] = m_operator[i].block(m_operator[i].rows()-m_shift[i][numPatch],shift0,m_shift[i][numPatch],m_shift[i][j]);
758 m_ddC[i][j] = m_operator[i].block(shift0,m_operator[i].cols()-m_shift[i][numPatch],m_shift[i][j],m_shift[i][numPatch]);
759 shift0 = shift0 + m_shift[i][j];
760 }
761 shift0 = 0;
762 }
763
764 m_A_aprox.resize(numLevels);
765 for(int i = 0 ; i < numLevels ; i++)
766 {
767 // Define the A_aprox matrix
768 m_A_aprox[i] = gsSparseMatrix<T>(m_operator[i].rows(),m_operator[i].cols());
769
770 // Retrieve a block of each patch
771 for(int k=0; k< numPatch; k++)
772 {
773 m_A_aprox[i].block(shift0,shift0,m_shift[i][k],m_shift[i][k]) = m_ILUT[i][k];
774 shift0 = shift0 + m_shift[i][k];
775 }
776 shift0 = 0;
777 m_ddBtilde[i].resize(numPatch);
778 m_ddCtilde[i].resize(numPatch);
779
780 for(int j=0 ; j < numPatch ; j ++)
781 {
782 m_ddBtilde[i][j] = gsSparseMatrix<T>(m_shift[i][j],m_shift[i][numPatch]);
783 m_ddCtilde[i][j] = gsSparseMatrix<T>(m_shift[i][j],m_shift[i][numPatch]);
784 for(int k=0 ; k < m_shift[i][numPatch]; k++)
785 {
786 gsMatrix<T> Brhs = m_ddC[i][j].col(k);
787 gsMatrix<T> Crhs = m_ddC[i][j].col(k);
788 m_ddBtilde[i][j].col(k) = m_ILUT[i][j].template triangularView<gsEigen::Upper>().transpose().solve(Brhs);
789 m_ddCtilde[i][j].col(k) = m_ILUT[i][j].template triangularView<gsEigen::UnitLower>().solve(Crhs);
790 }
791 }
792
793 // Define matrix S
794 m_S[i] = m_ddC[i][numPatch];
795 for(int l = 0 ; l < numPatch ; l++)
796 {
797 m_S[i] = m_S[i] - m_ddBtilde[i][l].transpose()*m_ddCtilde[i][l];
798 }
799
800 // Fill matrix A_aprox
801 for(int m = 0 ; m < numPatch ; m++)
802 {
803 m_A_aprox[i].block(shift0,m_A_aprox[i].rows() - m_shift[i][numPatch],m_shift[i][m],m_shift[i][numPatch]) = m_ddCtilde[i][m];
804 m_A_aprox[i].block(m_A_aprox[i].rows() - m_shift[i][numPatch],shift0,m_shift[i][numPatch],m_shift[i][m]) = m_ddBtilde[i][m].transpose();
805 shift0 = shift0 + m_shift[i][m];
806 }
807 shift0 = 0;
808
809 // Perform ILUT on the S-matrix!
810 gsEigen::IncompleteLUT<real_t> ilu;
811 ilu.setFillfactor(1);
812 gsSparseMatrix<T> II = m_S[i];
813 ilu.compute(II);
814 m_A_aprox[i].block(m_A_aprox[i].rows() - m_shift[i][numPatch],m_A_aprox[i].rows() - m_shift[i][numPatch],m_shift[i][numPatch],m_shift[i][numPatch]) = ilu.factors();
815 }
816 }
817
818 real_t Time_Block_ILUT_Factorization = clock.stop();
819 GISMO_UNUSED(Time_Block_ILUT_Factorization);
820
821 // gsInfo << "\n|| Setup Timings || " <<std::endl;
822 // gsInfo << "Total Assembly time: " << Time_Assembly <<std::endl;
823 // gsInfo << "Total ILUT factorization time: " << Time_ILUT_Factorization <<std::endl;
824 // gsInfo << "Total block ILUT factorization time: " << Time_Block_ILUT_Factorization <<std::endl;
825 // gsInfo << "Total SCMS time: " << Time_SCMS <<std::endl;
826 // gsInfo << "Total setup time: " << Time_Assembly_Galerkin + Time_Assembly + Time_Transfer + Time_ILUT_Factorization + Time_SCMS <<std::endl;
827 }
828
830 void solve(const gsFunctionExpr<T> & rhs,
831 gsMatrix<T>& x,
832 gsMatrix<T> f,
833 const int& iterTot,
834 const int& numLevels,
835 const int& typeMultigrid,
836 const gsMatrix<T>& hp,
837 const int& typeCoarseOperator)
838 {
839 gsStopwatch clock;
840 gsMatrix<T> b = f;
841
842 // Determine residual and L2 error
843 real_t r0 = (m_operator[numLevels-1]*x - b).norm();
844 real_t r = r0;
845 int iter = 1;
846
847 // Solve with p-multigrid method
848 real_t r_old = r0;
849 clock.restart();
850 // Adjusted stopping criterion!!
851 while( r/b.norm() > Base::tol && iter < Base::maxIter )
852 {
853 // Call solver from base class
854 Base::solveMG(b, m_bases, x, numLevels,
855 *m_bcInfo_ptr, *m_mp_ptr,
856 m_prolongation_P, m_restriction_P,
857 m_prolongation_M, m_restriction_M,
858 m_prolongation_H, m_restriction_H, hp);
859
860 r = (m_operator[numLevels-1]*x - b).norm();
861 if ( r_old < r)
862 {
863 gsInfo << "Residual increased during solving!!! " <<std::endl;
864 }
865 r_old = r;
866 //gsInfo << "Residual after cycle " << iter << " equals: " << r <<std::endl;
867 iter++;
868 }
869 real_t Time_Solve = clock.stop();
870 GISMO_UNUSED(Time_Solve);
871
872 // gsInfo << "\n|| Solver information || " <<std::endl;
873 // gsInfo << "Solver converged in " << Time_Solve << " seconds!" <<std::endl;
874 // gsInfo << "Solver converged in: " << iter-1 << " iterations!" <<std::endl;
875 }
876
877 private:
878
880 virtual void solvecoarse(const gsMatrix<T>& rhs,
881 gsMatrix<T>& x,
882 const int& numLevels)
883 {
884 //Direct solver (LU factorization)
885 CoarseSolver solver;
886 solver.analyzePattern(m_operator[numLevels-1]);
887 solver.factorize(m_operator[numLevels-1]);
888 x = solver.solve(rhs);
889 }
890
892 virtual gsMatrix<T> prolongation_M(const int& numLevels,
893 std::vector<memory::shared_ptr<gsMultiBasis<T> > > m_bases)
894 {
895 // Define the low and high order bases
896 gsMultiBasis<T> basesL = *m_bases[numLevels-2];
897 gsMultiBasis<T> basesH = *m_bases[numLevels-1];
898
899 // Determine matrix M (high_order * high_order)
900 typedef gsExprAssembler<real_t>::geometryMap geometryMap;
901 typedef gsExprAssembler<real_t>::variable variable;
902 typedef gsExprAssembler<real_t>::space space;
903 gsExprAssembler<real_t> ex2(1,1);
904 geometryMap G2 = ex2.getMap(*m_mp_ptr);
905 space w_n = ex2.getSpace(basesH ,1, 0);
906 w_n.setInterfaceCont(0);
907 if (Base::typeBCHandling == 1)
908 {
909 w_n.setup(*m_bcInfo_ptr, dirichlet::l2Projection, 0);
910 //#w_n.addBc(m_bcInfo_ptr->get("Dirichlet"));
911 }
912 ex2.setIntegrationElements(basesH);
913 ex2.initSystem();
914 ex2.assemble(w_n * meas(G2) );
915 return ex2.rhs();
916 }
917
919 virtual gsSparseMatrix<T> prolongation_P(const int& numLevels,
920 std::vector<memory::shared_ptr<gsMultiBasis<T> > > m_bases)
921 {
922 // Define the low and high order bases
923 gsMultiBasis<T> basesL = *m_bases[numLevels-2];
924 gsMultiBasis<T> basesH = *m_bases[numLevels-1];
925
926 // Determine matrix P (high_order * low_order)
927 typedef gsExprAssembler<real_t>::geometryMap geometryMap;
928 gsExprAssembler<real_t> ex(1,1);
929 geometryMap G = ex.getMap(*m_mp_ptr);
930 typedef gsExprAssembler<real_t>::variable variable;
931 typedef gsExprAssembler<real_t>::space space;
932 space v_n = ex.getSpace(basesH ,1, 0);
933 v_n.setInterfaceCont(0);
934 space u_n = ex.getTestSpace(v_n , basesL);
935 u_n.setInterfaceCont(0);
936 if (Base::typeBCHandling == 1)
937 {
938 v_n.setup(*m_bcInfo_ptr, dirichlet::l2Projection, 0);
939 u_n.setup(*m_bcInfo_ptr, dirichlet::l2Projection, 0);
940 //#v_n.addBc(m_bcInfo_ptr->get("Dirichlet"));
941 //#u_n.addBc(m_bcInfo_ptr->get("Dirichlet"));
942 }
943 ex.setIntegrationElements(basesH);
944 ex.initSystem();
945 ex.assemble(u_n*meas(G) * v_n.tr());
946 gsSparseMatrix<T> P = ex.matrix().transpose();
947 return P;
948 }
949
951 virtual gsMatrix<T> restriction_M(const int& numLevels,
952 std::vector<memory::shared_ptr<gsMultiBasis<T> > > m_bases)
953 {
954 // Define the low and high order bases
955 gsMultiBasis<T> basesL = *m_bases[numLevels-2];
956 gsMultiBasis<T> basesH = *m_bases[numLevels-1];
957
958 // Determine matrix M (low_order * low_order)
959 typedef gsExprAssembler<real_t>::geometryMap geometryMap;
960 typedef gsExprAssembler<real_t>::variable variable;
961 typedef gsExprAssembler<real_t>::space space;
962 gsExprAssembler<real_t> ex2(1,1);
963 geometryMap G2 = ex2.getMap(*m_mp_ptr);
964 space w_n = ex2.getSpace(basesL ,1, 0);
965 w_n.setInterfaceCont(0);
966 if (Base::typeBCHandling == 1)
967 {
968 w_n.setup(*m_bcInfo_ptr, dirichlet::l2Projection, 0);
969 //#w_n.addBc(m_bcInfo_ptr->get("Dirichlet"));
970 }
971 ex2.setIntegrationElements(basesL);
972 ex2.initSystem();
973 ex2.assemble(w_n * meas(G2) );
974 return ex2.rhs();
975 }
976
978 virtual gsSparseMatrix<T> restriction_P(const int& numLevels,
979 std::vector<memory::shared_ptr<gsMultiBasis<T> > > m_bases)
980 {
981 // Define the low and high order bases
982 gsMultiBasis<T> basesL = *m_bases[numLevels-2];
983 gsMultiBasis<T> basesH = *m_bases[numLevels-1];
984
985 // Determine matrix P (high_order * low_order)
986 gsExprAssembler<real_t> ex(1,1);
987 typedef gsExprAssembler<real_t>::geometryMap geometryMap;
988 geometryMap G = ex.getMap(*m_mp_ptr);
989
990 typedef gsExprAssembler<real_t>::variable variable;
991 typedef gsExprAssembler<real_t>::space space;
992 space v_n = ex.getSpace(basesH ,1, 0);
993 v_n.setInterfaceCont(0);
994 space u_n = ex.getTestSpace(v_n , basesL);
995 u_n.setInterfaceCont(0);
996 if (Base::typeBCHandling == 1)
997 {
998 u_n.setup(*m_bcInfo_ptr, dirichlet::l2Projection, 0);
999 v_n.setup(*m_bcInfo_ptr, dirichlet::l2Projection, 0);
1000 //#u_n.addBc(m_bcInfo_ptr->get("Dirichlet"));
1001 //#v_n.addBc(m_bcInfo_ptr->get("Dirichlet"));
1002 }
1003 ex.setIntegrationElements(basesH);
1004 ex.initSystem();
1005 ex.assemble(u_n * meas(G)* v_n.tr());
1006 gsSparseMatrix<T> P = ex.matrix();
1007 return P;
1008 }
1009
1011 virtual void presmoothing(const gsMatrix<T>& rhs,
1012 gsMatrix<T>& x,
1013 const int& numLevels,
1014 gsMatrix<T> & fineRes,
1015 const gsMatrix<T>& hp)
1016 {
1017 //gsInfo << "Residual before presmoothing: " << (rhs-m_operator[numLevels-1]*x).norm() << " at level " << numLevels <<std::endl;
1018 for(int i = 0 ; i < Base::numSmoothing ; i++)
1019 {
1020 if (Base::typeSmoother == 1)
1021 {
1022 if (hp(numLevels-2,0) == 1 && hp(hp.rows()-1,0) == 0)
1023 {
1024 internal::gaussSeidelSweep(m_operator[numLevels-1],x,rhs);
1025 }
1026 else
1027 {
1028 gsMatrix<T> e;
1029 gsMatrix<T> d = rhs-m_operator[numLevels-1]*x;
1030 e = m_Pinv[numLevels-1][0]*d;
1031 e = m_ILUT[numLevels-1][0].template triangularView<gsEigen::UnitLower>().solve(e);
1032 e = m_ILUT[numLevels-1][0].template triangularView<gsEigen::Upper>().solve(e);
1033 e = m_P[numLevels-1][0]*e;
1034 x = x + e;
1035 }
1036 }
1037 if (Base::typeSmoother == 2)
1038 {
1039 internal::gaussSeidelSweep(m_operator[numLevels-1],x,rhs);
1040 }
1041 if (Base::typeSmoother == 3)
1042 {
1043 m_SCMS[numLevels-1]->step(rhs,x);
1044 }
1045 if (Base::typeSmoother == 5)
1046 {
1047 if (hp(numLevels-2,0) == 1 && hp(hp.rows()-1,0) == 0)
1048 {
1049 internal::gaussSeidelSweep(m_operator[numLevels-1],x,rhs);
1050 }
1051 else
1052 {
1053 gsMatrix<T> e;
1054 gsMatrix<T> d = rhs-m_operator[numLevels-1]*x;
1055 e = m_A_aprox[numLevels-1].template triangularView<gsEigen::UnitLower>().solve(d);
1056 e = m_A_aprox[numLevels-1].template triangularView<gsEigen::Upper>().solve(e);
1057 x = x + e;
1058 }
1059 }
1060 }
1061 // gsInfo << "Residual after presmoothing: " << (rhs-m_operator[numLevels-1]*x).norm() << " at level " << numLevels <<std::endl;
1062 fineRes = m_operator[numLevels-1]*x - rhs;
1063 }
1064
1066 virtual void postsmoothing(const gsMatrix<T>& rhs,
1067 gsMatrix<T>& x,
1068 const int& numLevels,
1069 gsMatrix<T> & fineCorr,
1070 gsMatrix<T> & postRes,
1071 const gsMatrix<T>& hp)
1072 {
1073 real_t alpha = 1;
1074 x = x - alpha*fineCorr;
1075 //gsInfo << "Residual before postsmoothing: " << (rhs-m_operator[numLevels-1]*x).norm() << " at level " << numLevels <<std::endl;
1076
1077 for(int i = 0 ; i < Base::numSmoothing ; i++)
1078 {
1079 if (Base::typeSmoother == 1)
1080 {
1081 if (hp(numLevels-2,0) == 1 && hp(hp.rows()-1,0) == 0)
1082 {
1083 internal::gaussSeidelSweep(m_operator[numLevels-1],x,rhs);
1084 }
1085 else
1086 {
1087 gsMatrix<T> e;
1088 gsMatrix<T> d = rhs-m_operator[numLevels-1]*x;
1089 e = m_Pinv[numLevels-1][0]*d;
1090 e = m_ILUT[numLevels-1][0].template triangularView<gsEigen::UnitLower>().solve(e);
1091 e = m_ILUT[numLevels-1][0].template triangularView<gsEigen::Upper>().solve(e);
1092 e = m_P[numLevels-1][0]*e;
1093 x = x + e;
1094 }
1095 }
1096 if (Base::typeSmoother == 2)
1097 {
1098 internal::gaussSeidelSweep(m_operator[numLevels-1],x,rhs);
1099 }
1100 if (Base::typeSmoother == 3)
1101 {
1102 m_SCMS[numLevels-1]->step(rhs,x);
1103 }
1104 if (Base::typeSmoother == 5)
1105 {
1106 if (hp(numLevels-2,0) == 1 && hp(hp.rows()-1,0) == 0)
1107 {
1108 internal::gaussSeidelSweep(m_operator[numLevels-1],x,rhs);
1109 }
1110 else
1111 {
1112 gsMatrix<T> e;
1113 gsMatrix<T> d = rhs-m_operator[numLevels-1]*x;
1114 e = m_A_aprox[numLevels-1].template triangularView<gsEigen::UnitLower>().solve(d);
1115 e = m_A_aprox[numLevels-1].template triangularView<gsEigen::Upper>().solve(e);
1116 x = x + e;
1117 }
1118 }
1119 postRes = rhs - m_operator[numLevels-1]*x;
1120 // gsInfo << "Residual after postsmoothing: " << (rhs-m_operator[numLevels-1]*x).norm() << " at level " << numLevels <<std::endl;
1121 }
1122 }
1123 };
1124
1125// Create the subspace corrected mass smoother
1126gsPreconditionerOp<>::Ptr setupSubspaceCorrectedMassSmoother(const gsSparseMatrix<>& matrix, const gsMultiBasis<>& mb, const gsBoundaryConditions<>& bc, const gsOptionList& opt, const int &typeBCHandling)
1127{
1128 const short_t dim = mb.topology().dim();
1129
1130 // Setup dof mapper
1131 gsDofMapper dm;
1132 mb.getMapper(
1133 typeBCHandling == 1 ? (dirichlet::strategy)opt.askInt("DirichletStrategy",11) : (dirichlet::strategy)opt.askInt("DirichletStrategy",14),
1134 (iFace ::strategy)opt.askInt("InterfaceStrategy", 1),
1135 bc,
1136 dm,
1137 0
1138 );
1139 const index_t nTotalDofs = dm.freeSize();
1140
1141 // Decompose the whole domain into components
1142 std::vector< std::vector<patchComponent> > components = mb.topology().allComponents(true);
1143 const index_t nr_components = components.size();
1144
1145 // Setup Dirichlet boundary conditions
1146 gsBoundaryConditions<> dir_bc;
1147 for( index_t ps=0; ps < 2*dim; ++ps )
1148 dir_bc.addCondition( 0, 1+ps, condition_type::dirichlet, NULL );
1149
1150 // Setup transfer matrices and local preconditioners
1151 std::vector< gsSparseMatrix<real_t,RowMajor> > transfers;
1152 transfers.reserve(nr_components);
1153 std::vector< gsLinearOperator<>::Ptr > ops;
1154 ops.reserve(nr_components);
1155
1156 for (index_t i=0; i<nr_components; ++i)
1157 {
1158 gsMatrix<index_t> indices;
1159 std::vector<gsBasis<>::uPtr> bases = mb.componentBasis_withIndices(components[i],dm,indices,true);
1160 index_t sz = indices.rows();
1161 gsSparseEntries<> se;
1162 se.reserve(sz);
1163 for (index_t i=0; i<sz; ++i)
1164 se.add(indices(i,0),i,real_t(1));
1165 gsSparseMatrix<real_t,RowMajor> transfer(nTotalDofs,sz);
1166 transfer.setFrom(se);
1167 if (sz>0)
1168 {
1169 if (bases[0]->dim() == dim)
1170 {
1171 GISMO_ASSERT ( bases.size() == 1, "Only one basis is expected for each patch." );
1172 ops.push_back(
1173 gsPatchPreconditionersCreator<>::subspaceCorrectedMassSmootherOp(
1174 *(bases[0]),
1175 dir_bc,
1176 gsOptionList(),
1177 opt.getReal("Scaling")
1178 )
1179 );
1180 }
1181 else
1182 {
1183 gsSparseMatrix<> mat = transfer.transpose() * matrix * transfer;
1184 ops.push_back( makeSparseCholeskySolver(mat) );
1185 }
1186 transfers.push_back(give(transfer));
1187 }
1188 }
1189 return gsPreconditionerFromOp<>::make(makeMatrixOp(matrix), gsAdditiveOp<>::make(transfers, ops));
1190}
1191
1192
1193} // namespace gismo
gismo::expr::gsFeSpace
Definition gsExpressions.h:973
gismo::expr::gsFeVariable
Definition gsExpressions.h:928
gismo::gsAdditiveOp
Generic preconditioner which applies an arbitrary linear operator to the residual.
Definition gsAdditiveOp.h:51
gismo::gsBoundaryConditions
Class containing a set of boundary conditions.
Definition gsBoundaryConditions.h:342
gismo::gsBoundaryConditions::addCondition
void addCondition(int p, boxSide s, condition_type::type t, gsFunctionSet< T > *f, short_t unknown=0, bool parametric=false, int comp=-1)
Adds another boundary condition.
Definition gsBoundaryConditions.h:650
gismo::gsConjugateGradient
The conjugate gradient method.
Definition gsConjugateGradient.h:30
gismo::gsDofMapper
Maintains a mapping from patch-local dofs to global dof indices and allows the elimination of individ...
Definition gsDofMapper.h:69
gismo::gsDofMapper::freeSize
index_t freeSize() const
Returns the number of free (not eliminated) dofs.
Definition gsDofMapper.h:436
gismo::gsExprAssembler
Definition gsExprAssembler.h:33
gismo::gsExprAssembler::setIntegrationElements
void setIntegrationElements(const gsMultiBasis< T > &mesh)
Sets the domain of integration.
Definition gsExprAssembler.h:161
gismo::gsExprAssembler::geometryMap
gsExprHelper< T >::geometryMap geometryMap
Geometry map type.
Definition gsExprAssembler.h:65
gismo::gsExprAssembler::getTestSpace
space getTestSpace(const gsFunctionSet< T > &mp, index_t dim=1, index_t id=0)
Definition gsExprAssembler.h:217
gismo::gsExprAssembler::getSpace
space getSpace(const gsFunctionSet< T > &mp, index_t dim=1, index_t id=0)
Definition gsExprAssembler.h:191
gismo::gsExprAssembler::rhs
const gsMatrix< T > & rhs() const
Returns the right-hand side vector(s)
Definition gsExprAssembler.h:154
gismo::gsExprAssembler::getMap
geometryMap getMap(const gsFunctionSet< T > &g)
Registers g as an isogeometric geometry map and return a handle to it.
Definition gsExprAssembler.h:186
gismo::gsExprAssembler::initSystem
void initSystem(const index_t numRhs=1)
Initializes the sparse system (sparse matrix and rhs)
Definition gsExprAssembler.h:315
gismo::gsExprAssembler::getCoeff
variable getCoeff(const gsFunctionSet< T > &func)
Definition gsExprAssembler.h:285
gismo::gsExprAssembler::assemble
void assemble(const expr &... args)
Adds the expressions args to the system matrix/rhs The arguments are considered as integrals over the...
Definition gsExprAssembler.h:1077
gismo::gsExprAssembler::solution
expr::gsFeSolution< T > solution
Solution type.
Definition gsExprAssembler.h:68
gismo::gsExprAssembler::matrix
const gsSparseMatrix< T > & matrix() const
Returns the left-hand global matrix.
Definition gsExprAssembler.h:127
gismo::gsFunctionExpr
Class defining a multivariate (real or vector) function given by a string mathematical expression.
Definition gsFunctionExpr.h:52
gismo::gsGenericStopwatch::stop
double stop()
Return elapsed time in seconds.
Definition gsStopwatch.h:83
gismo::gsGenericStopwatch::restart
void restart()
Start taking the time.
Definition gsStopwatch.h:80
gismo::gsIterativeSolver::setTolerance
void setTolerance(T tol)
Set the tolerance for the error criteria on the relative residual error (default: 1e-10)
Definition gsIterativeSolver.h:223
gismo::gsIterativeSolver::solve
void solve(const VectorType &rhs, VectorType &x)
Solves the linear system and stores the solution in x.
Definition gsIterativeSolver.h:114
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsMultiBasis
Holds a set of patch-wise bases and their topology information.
Definition gsMultiBasis.h:37
gismo::gsMultiBasis::uniformCoarsen_withTransfer
void uniformCoarsen_withTransfer(gsSparseMatrix< T, RowMajor > &transfer, const gsBoundaryConditions< T > &boundaryConditions, const gsOptionList &assemblerOptions, int numKnots=1, index_t unk=0)
Coarsen every basis uniformly.
Definition gsMultiBasis.hpp:222
gismo::gsMultiPatch
Container class for a set of geometry patches and their topology, that is, the interface connections ...
Definition gsMultiPatch.h:100
gismo::gsOptionList
Class which holds a list of parameters/options, and provides easy access to them.
Definition gsOptionList.h:33
gismo::gsOptionList::addInt
void addInt(const std::string &label, const std::string &desc, const index_t &value)
Adds a option named label, with description desc and value value.
Definition gsOptionList.cpp:201
gismo::gsOptionList::getReal
Real getReal(const std::string &label) const
Reads value for option label from options.
Definition gsOptionList.cpp:44
gismo::gsOptionList::addReal
void addReal(const std::string &label, const std::string &desc, const Real &value)
Adds a option named label, with description desc and value value.
Definition gsOptionList.cpp:211
gismo::gsOptionList::askInt
index_t askInt(const std::string &label, const index_t &value=0) const
Reads value for option label from options.
Definition gsOptionList.cpp:117
gismo::gsPatchPreconditionersCreator
Provides robust preconditioners for single patch geometries.
Definition gsPatchPreconditionersCreator.h:30
gismo::gsPreconditionerFromOp::make
static uPtr make(BasePtr underlying, BasePtr preconditioner, T tau=(T) 1)
Make function returning a smart pointer.
Definition gsPreconditioner.h:188
gismo::gsPreconditionerOp::Ptr
memory::shared_ptr< gsPreconditionerOp > Ptr
Shared pointer for gsLinearOperator.
Definition gsPreconditioner.h:47
gismo::gsSparseEntries
Class that provides a container for triplets (i,j,value) to be filled in a sparse matrix.
Definition gsSparseMatrix.h:34
gismo::gsSparseMatrix
Sparse matrix class, based on gsEigen::SparseMatrix.
Definition gsSparseMatrix.h:139
gismo.h
Main header to be included by clients using the G+Smo library.
short_t
#define short_t
Definition gsConfig.h:35
index_t
#define index_t
Definition gsConfig.h:32
GISMO_UNUSED
#define GISMO_UNUSED(x)
Definition gsDebug.h:112
gsInfo
#define gsInfo
Definition gsDebug.h:43
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gismo::memory::make_shared_not_owned
shared_ptr< T > make_shared_not_owned(const T *x)
Creates a shared pointer which does not eventually delete the underlying raw pointer....
Definition gsMemory.h:189
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::give
S give(S &x)
Definition gsMemory.h:266
gismo::boundary
Struct that defines the boundary sides and corners and types of a geometric object.
Definition gsBoundary.h:56
gismo::condition_type::dirichlet
@ dirichlet
Dirichlet type.
Definition gsBoundaryConditions.h:31
gismo::gsXBraidMultigridBase
The p-multigrid base class provides the basic methods (smoothing, prolongation, restriction) for impl...
Definition gsXBraidMultigrid.h:13
gismo::gsXBraidMultigridBase::restriction
virtual void restriction(const gsMatrix< T > &Xfine, gsMatrix< T > &Xcoarse, const int &numLevels, std::vector< memory::shared_ptr< gsMultiBasis< T > > > m_bases, gsBoundaryConditions< T > bcInfo, gsMultiPatch< T > mp, std::vector< gsSparseMatrix< T > > &m_prolongation_P, std::vector< gsSparseMatrix< T > > &m_restriction_P, std::vector< gsMatrix< T > > &m_prolongation_M, std::vector< gsMatrix< T > > &m_restriction_M, std::vector< gsSparseMatrix< T > > &m_prolongation_H, std::vector< gsSparseMatrix< T > > &m_restriction_H, const gsMatrix< T > &hp)
Restrict fine space function to coarse space.
Definition gsXBraidMultigrid.h:330
gismo::gsXBraidMultigridBase::solvecoarse
virtual void solvecoarse(const gsMatrix< T > &rhs, gsMatrix< T > &x, const int &numLevels)=0
Apply coarse solver (pure virtual method)
gismo::gsXBraidMultigridBase::prolongation
virtual void prolongation(const gsMatrix< T > &Xcoarse, gsMatrix< T > &Xfine, const int &numLevels, std::vector< memory::shared_ptr< gsMultiBasis< T > > > m_bases, gsBoundaryConditions< T > bcInfo, gsMultiPatch< T > mp, std::vector< gsSparseMatrix< T > > &m_prolongation_P, std::vector< gsSparseMatrix< T > > &m_restriction_P, std::vector< gsMatrix< T > > &m_prolongation_M, std::vector< gsMatrix< T > > &m_restriction_M, std::vector< gsSparseMatrix< T > > &m_prolongation_H, std::vector< gsSparseMatrix< T > > &m_restriction_H, const gsMatrix< T > &hp)
Prolongate coarse space function to fine space.
Definition gsXBraidMultigrid.h:270
gismo::gsXBraidMultigridBase::solveMG
virtual void solveMG(const gsMatrix< T > &rhs, std::vector< memory::shared_ptr< gsMultiBasis< T > > > m_bases, gsMatrix< T > &x, const int &numLevels, gsBoundaryConditions< T > bcInfo, gsMultiPatch< T > mp, std::vector< gsSparseMatrix< T > > &m_prolongation_P, std::vector< gsSparseMatrix< T > > &m_restriction_P, std::vector< gsMatrix< T > > &m_prolongation_M, std::vector< gsMatrix< T > > &m_restriction_M, std::vector< gsSparseMatrix< T > > &m_prolongation_H, std::vector< gsSparseMatrix< T > > &m_restriction_H, const gsMatrix< T > &hp)
Apply p-multigrid solver to given right-hand side on level l.
Definition gsXBraidMultigrid.h:135
gismo::gsXBraidMultigridBase::restriction_M
virtual gsMatrix< T > restriction_M(const int &numLevels, std::vector< memory::shared_ptr< gsMultiBasis< T > > > m_bases)=0
Prolongate coarse space function to fine space.
gismo::gsXBraidMultigridBase::postsmoothing
virtual void postsmoothing(const gsMatrix< T > &rhs, gsMatrix< T > &x, const int &numLevels, gsMatrix< T > &fineCorr, gsMatrix< T > &postRes, const gsMatrix< T > &hp)=0
Apply fixed number of smoothing steps (pure virtual method)
gismo::gsXBraidMultigridBase::restriction_P
virtual gsSparseMatrix< T > restriction_P(const int &numLevels, std::vector< memory::shared_ptr< gsMultiBasis< T > > > m_bases)=0
Prolongate coarse space function to fine space.
gismo::gsXBraidMultigridBase::presmoothing
virtual void presmoothing(const gsMatrix< T > &rhs, gsMatrix< T > &x, const int &numLevels, gsMatrix< T > &fineRes, const gsMatrix< T > &hp)=0
Apply fixed number of smoothing steps (pure virtual method)
gismo::gsXBraidMultigridBase::compute
virtual gsXBraidMultigridBase & compute(const gsSparseMatrix< T > &mat, const T tstep, const int &numDegree, index_t typeMethod)
Definition gsXBraidMultigrid.h:84
gismo::gsXBraidMultigridBase::prolongation_P
virtual gsSparseMatrix< T > prolongation_P(const int &numLevels, std::vector< memory::shared_ptr< gsMultiBasis< T > > > m_bases)=0
Prolongate coarse space function to fine space.
gismo::gsXBraidMultigridBase::solveWithGuess
virtual gsMatrix< T > solveWithGuess(const gsMatrix< T > &b, const gsMatrix< T > &x0)
Definition gsXBraidMultigrid.h:110
gismo::gsXBraidMultigridBase::gsXBraidMultigridBase
gsXBraidMultigridBase()
Constructor.
Definition gsXBraidMultigrid.h:29
gismo::gsXBraidMultigridBase::prolongation_M
virtual gsMatrix< T > prolongation_M(const int &numLevels, std::vector< memory::shared_ptr< gsMultiBasis< T > > > m_bases)=0
Prolongate coarse space function to fine space.
gismo::gsXBraidMultigrid
The p-multigrid class implements a generic p-multigrid solver that can be customized by passing assem...
Definition gsXBraidMultigrid.h:401
gismo::gsXBraidMultigrid::prolongation_P
virtual gsSparseMatrix< T > prolongation_P(const int &numLevels, std::vector< memory::shared_ptr< gsMultiBasis< T > > > m_bases)
Construct prolongation operator at level numLevels.
Definition gsXBraidMultigrid.h:919
gismo::gsXBraidMultigrid::restriction_P
virtual gsSparseMatrix< T > restriction_P(const int &numLevels, std::vector< memory::shared_ptr< gsMultiBasis< T > > > m_bases)
Construct restriction operator at level numLevels.
Definition gsXBraidMultigrid.h:978
gismo::gsXBraidMultigrid::m_ILUT
std::vector< std::vector< gsSparseMatrix< T > > > m_ILUT
std::vector of factorized operators
Definition gsXBraidMultigrid.h:435
gismo::gsXBraidMultigrid::m_SCMS
std::vector< typename gsPreconditionerOp< T >::Ptr > m_SCMS
std::vector of SCM smoother object
Definition gsXBraidMultigrid.h:444
gismo::gsXBraidMultigrid::m_Pinv
std::vector< std::vector< gsEigen::PermutationMatrix< Dynamic, Dynamic, index_t > > > m_Pinv
std::vector of factorized operators
Definition gsXBraidMultigrid.h:441
gismo::gsXBraidMultigrid::m_ddCtilde
std::vector< std::vector< gsMatrix< T > > > m_ddCtilde
std::vector of std::vector of block operator objects
Definition gsXBraidMultigrid.h:462
gismo::gsXBraidMultigrid::m_ddBtilde
std::vector< std::vector< gsMatrix< T > > > m_ddBtilde
std::vector of std::vector of block operator objects
Definition gsXBraidMultigrid.h:459
gismo::gsXBraidMultigrid::setup
void setup(const gsFunctionExpr< T > &rhs, gsMatrix< T > &x, gsMatrix< T > f, const int &iterTot, const int &numLevels, const int &numDegree, const int &typeMultigrid, const gsMatrix< T > &hp, const int &typeCoarseOperator, T tstep, index_t typeMethod)
Set-up p-multigrid solver.
Definition gsXBraidMultigrid.h:490
gismo::gsXBraidMultigrid::m_restriction_M
std::vector< gsMatrix< T > > m_restriction_M
std::vector of restriction operators
Definition gsXBraidMultigrid.h:426
gismo::gsXBraidMultigrid::presmoothing
virtual void presmoothing(const gsMatrix< T > &rhs, gsMatrix< T > &x, const int &numLevels, gsMatrix< T > &fineRes, const gsMatrix< T > &hp)
Apply fixed number of presmoothing steps.
Definition gsXBraidMultigrid.h:1011
gismo::gsXBraidMultigrid::m_prolongation_H
std::vector< gsSparseMatrix< T > > m_prolongation_H
std::vector of prolongation operators
Definition gsXBraidMultigrid.h:429
gismo::gsXBraidMultigrid::m_operator
std::vector< gsSparseMatrix< T > > m_operator
std::vector of operator objects
Definition gsXBraidMultigrid.h:447
gismo::gsXBraidMultigrid::prolongation_M
virtual gsMatrix< T > prolongation_M(const int &numLevels, std::vector< memory::shared_ptr< gsMultiBasis< T > > > m_bases)
Construct prolongation operator at level numLevels.
Definition gsXBraidMultigrid.h:892
gismo::gsXBraidMultigrid::m_bases
std::vector< memory::shared_ptr< gsMultiBasis< T > > > m_bases
std::vector of multi-basis objects
Definition gsXBraidMultigrid.h:414
gismo::gsXBraidMultigrid::solve
void solve(const gsFunctionExpr< T > &rhs, gsMatrix< T > &x, gsMatrix< T > f, const int &iterTot, const int &numLevels, const int &typeMultigrid, const gsMatrix< T > &hp, const int &typeCoarseOperator)
Apply p-multigrid solver to given right-hand side on level l.
Definition gsXBraidMultigrid.h:830
gismo::gsXBraidMultigrid::m_A_aprox
std::vector< gsMatrix< T > > m_A_aprox
std::vector of std::vector of block operator objects
Definition gsXBraidMultigrid.h:465
gismo::gsXBraidMultigrid::m_restriction_H
std::vector< gsSparseMatrix< T > > m_restriction_H
std::vector of restriction operators
Definition gsXBraidMultigrid.h:432
gismo::gsXBraidMultigrid::m_bcInfo_ptr
memory::shared_ptr< gsBoundaryConditions< T > > m_bcInfo_ptr
Shared pointer to boundary conditions.
Definition gsXBraidMultigrid.h:411
gismo::gsXBraidMultigrid::m_mp_ptr
memory::shared_ptr< gsMultiPatch< T > > m_mp_ptr
Shared pointer to multi-patch geometry.
Definition gsXBraidMultigrid.h:408
gismo::gsXBraidMultigrid::m_P
std::vector< std::vector< gsEigen::PermutationMatrix< Dynamic, Dynamic, index_t > > > m_P
std::vector of factorized operators
Definition gsXBraidMultigrid.h:438
gismo::gsXBraidMultigrid::m_ddC
std::vector< std::vector< gsSparseMatrix< T > > > m_ddC
std::vector of std::vector of block operator objects
Definition gsXBraidMultigrid.h:456
gismo::gsXBraidMultigrid::m_shift
std::vector< std::vector< int > > m_shift
std::vector of std::vector of shift objects
Definition gsXBraidMultigrid.h:471
gismo::gsXBraidMultigrid::restriction_M
virtual gsMatrix< T > restriction_M(const int &numLevels, std::vector< memory::shared_ptr< gsMultiBasis< T > > > m_bases)
Construct restriction operator at level numLevels.
Definition gsXBraidMultigrid.h:951
gismo::gsXBraidMultigrid::m_S
std::vector< gsSparseMatrix< T > > m_S
std::vector of std::vector of block operator objects
Definition gsXBraidMultigrid.h:468
gismo::gsXBraidMultigrid::m_ddB
std::vector< std::vector< gsSparseMatrix< T > > > m_ddB
std::vector of std::vector of block operator objects
Definition gsXBraidMultigrid.h:453
gismo::gsXBraidMultigrid::m_block_operator
std::vector< std::vector< gsSparseMatrix< T > > > m_block_operator
std::vector of std::vector of block operator objects
Definition gsXBraidMultigrid.h:450
gismo::gsXBraidMultigrid::m_restriction_P
std::vector< gsSparseMatrix< T > > m_restriction_P
std::vector of restriction operators
Definition gsXBraidMultigrid.h:420
gismo::gsXBraidMultigrid::postsmoothing
virtual void postsmoothing(const gsMatrix< T > &rhs, gsMatrix< T > &x, const int &numLevels, gsMatrix< T > &fineCorr, gsMatrix< T > &postRes, const gsMatrix< T > &hp)
Apply fixed number of postsmoothing steps.
Definition gsXBraidMultigrid.h:1066
gismo::gsXBraidMultigrid::m_prolongation_P
std::vector< gsSparseMatrix< T > > m_prolongation_P
std::vector of prolongation operators
Definition gsXBraidMultigrid.h:417
gismo::gsXBraidMultigrid::m_prolongation_M
std::vector< gsMatrix< T > > m_prolongation_M
std::vector of prolongation operators
Definition gsXBraidMultigrid.h:423
gismo::gsXBraidMultigrid::solvecoarse
virtual void solvecoarse(const gsMatrix< T > &rhs, gsMatrix< T > &x, const int &numLevels)
Apply coarse solver.
Definition gsXBraidMultigrid.h:880
gismo::gsXBraidMultigrid::Base
gsXBraidMultigridBase< T > Base
Base class type.
Definition gsXBraidMultigrid.h:405