G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsXBraidMultigrid.h
1 #include <gismo.h>
2 #include <string>
3 
4 namespace 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
1126 gsPreconditionerOp<>::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::gsXBraidMultigrid::m_restriction_P
std::vector< gsSparseMatrix< T > > m_restriction_P
std::vector of restriction operators
Definition: gsXBraidMultigrid.h:420
gismo::gsMultiBasis::componentBasis_withIndices
gsBasis< T >::uPtr componentBasis_withIndices(patchComponent pc, const gsDofMapper &dm, gsMatrix< index_t > &indices, bool no_lower=true) const
Returns the basis that corresponds to the component.
Definition: gsMultiBasis.hpp:262
gismo::gsSparseEntries
Class that provides a container for triplets (i,j,value) to be filled in a sparse matrix...
Definition: gsSparseMatrix.h:33
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::gsExprAssembler
Definition: gsExprAssembler.h:30
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::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::gsXBraidMultigridBase::gsXBraidMultigridBase
gsXBraidMultigridBase()
Constructor.
Definition: gsXBraidMultigrid.h:29
gismo::gsXBraidMultigridBase
The p-multigrid base class provides the basic methods (smoothing, prolongation, restriction) for impl...
Definition: gsXBraidMultigrid.h:12
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. Usefull to refer to objects which should not be destroyed.
Definition: gsMemory.h:189
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::gsExprAssembler::geometryMap
gsExprHelper< T >::geometryMap geometryMap
Geometry map type.
Definition: gsExprAssembler.h:58
gismo::condition_type::dirichlet
Dirichlet type.
Definition: gsBoundaryConditions.h:31
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::gsXBraidMultigrid::m_shift
std::vector< std::vector< int > > m_shift
std::vector of std::vector of shift objects
Definition: gsXBraidMultigrid.h:471
gismo::gsBoxTopology::dim
short_t dim() const
Dimension of the boxes.
Definition: gsBoxTopology.h:98
gismo.h
Main header to be included by clients using the G+Smo library.
gismo::gsXBraidMultigrid::m_restriction_M
std::vector< gsMatrix< T > > m_restriction_M
std::vector of restriction operators
Definition: gsXBraidMultigrid.h:426
gismo::gsXBraidMultigridBase::compute
virtual gsXBraidMultigridBase & compute(const gsSparseMatrix< T > &mat, const T tstep, const int &numDegree, index_t typeMethod)
Definition: gsXBraidMultigrid.h:84
short_t
#define short_t
Definition: gsConfig.h:35
gismo::gsExprAssembler::getTestSpace
space getTestSpace(const gsFunctionSet< T > &mp, index_t dim=1, index_t id=0)
Definition: gsExprAssembler.h:192
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::boundary
Struct that defines the boundary sides and corners and types of a geometric object.
Definition: gsBoundary.h:55
gismo::gsXBraidMultigrid::m_prolongation_M
std::vector< gsMatrix< T > > m_prolongation_M
std::vector of prolongation operators
Definition: gsXBraidMultigrid.h:423
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::gsConjugateGradient
The conjugate gradient method.
Definition: gsConjugateGradient.h:29
gismo::gsAdditiveOp
Generic preconditioner which applies an arbitrary linear operator to the residual.
Definition: gsAdditiveOp.h:50
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::give
S give(S &x)
Definition: gsMemory.h:266
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:756
gismo::gsXBraidMultigrid::m_prolongation_H
std::vector< gsSparseMatrix< T > > m_prolongation_H
std::vector of prolongation operators
Definition: gsXBraidMultigrid.h:429
index_t
#define index_t
Definition: gsConfig.h:32
gismo::gsXBraidMultigrid::m_bcInfo_ptr
memory::shared_ptr< gsBoundaryConditions< T > > m_bcInfo_ptr
Shared pointer to boundary conditions.
Definition: gsXBraidMultigrid.h:411
gismo::gsXBraidMultigrid::Base
gsXBraidMultigridBase< T > Base
Base class type.
Definition: gsXBraidMultigrid.h:405
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition: gsMatrix.h:38
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::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_ASSERT
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
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::gsDofMapper
Maintains a mapping from patch-local dofs to global dof indices and allows the elimination of individ...
Definition: gsDofMapper.h:68
gismo::gsExprAssembler::rhs
const gsMatrix< T > & rhs() const
Returns the right-hand side vector(s)
Definition: gsExprAssembler.h:129
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::gsExprAssembler::matrix
const gsSparseMatrix< T > & matrix() const
Returns the left-hand global matrix.
Definition: gsExprAssembler.h:116
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::solvecoarse
virtual void solvecoarse(const gsMatrix< T > &rhs, gsMatrix< T > &x, const int &numLevels)
Apply coarse solver.
Definition: gsXBraidMultigrid.h:880
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_ddCtilde
std::vector< std::vector< gsMatrix< T > > > m_ddCtilde
std::vector of std::vector of block operator objects
Definition: gsXBraidMultigrid.h:462
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::gsMultiBasis
Holds a set of patch-wise bases and their topology information.
Definition: gsMultiBasis.h:36
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
The p-multigrid class implements a generic p-multigrid solver that can be customized by passing assem...
Definition: gsXBraidMultigrid.h:400
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_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::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::gsGenericStopwatch::restart
void restart()
Start taking the time.
Definition: gsStopwatch.h:80
gismo::gsXBraidMultigrid::m_prolongation_P
std::vector< gsSparseMatrix< T > > m_prolongation_P
std::vector of prolongation operators
Definition: gsXBraidMultigrid.h:417
gismo::gsXBraidMultigrid::m_operator
std::vector< gsSparseMatrix< T > > m_operator
std::vector of operator objects
Definition: gsXBraidMultigrid.h:447
gsInfo
#define gsInfo
Definition: gsDebug.h:43
gismo::gsGenericStopwatch::stop
double stop()
Return elapsed time in seconds.
Definition: gsStopwatch.h:83
gismo::expr::gsFeSpace
Definition: gsDirichletValues.h:23
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::gsXBraidMultigrid::m_restriction_H
std::vector< gsSparseMatrix< T > > m_restriction_H
std::vector of restriction operators
Definition: gsXBraidMultigrid.h:432
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:161
gismo::gsMultiPatch
Container class for a set of geometry patches and their topology, that is, the interface connections ...
Definition: gsMultiPatch.h:33
gismo::gsExprAssembler::getSpace
space getSpace(const gsFunctionSet< T > &mp, index_t dim=1, index_t id=0)
Definition: gsExprAssembler.h:166
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::gsGenericStopwatch
A Stopwatch object can be used to measure execution time of code, algorithms, etc.
Definition: gsStopwatch.h:72
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_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_Pinv
std::vector< std::vector< gsEigen::PermutationMatrix< Dynamic, Dynamic, index_t > > > m_Pinv
std::vector of factorized operators
Definition: gsXBraidMultigrid.h:441
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::gsBoundaryConditions
Class containing a set of boundary conditions.
Definition: gsBoundaryConditions.h:341
gismo::gsFunctionExpr
Class defining a multivariate (real or vector) function given by a string mathematical expression...
Definition: gsFunctionExpr.h:51
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::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_UNUSED
#define GISMO_UNUSED(x)
Definition: gsDebug.h:112
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::gsPatchPreconditionersCreator
Provides robust preconditioners for single patch geometries.
Definition: gsPatchPreconditionersCreator.h:29
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::gsBoxTopology::allComponents
std::vector< std::vector< patchComponent > > allComponents(bool combineCorners=false) const
Returns all components representing the topology.
Definition: gsBoxTopology.cpp:294
gismo::expr::meas
EIGEN_STRONG_INLINE meas_expr< T > meas(const gsGeometryMap< T > &G)
The measure of a geometry map.
Definition: gsExpressions.h:4557
gismo::gsXBraidMultigridBase::solveWithGuess
virtual gsMatrix< T > solveWithGuess(const gsMatrix< T > &b, const gsMatrix< T > &x0)
Definition: gsXBraidMultigrid.h:110
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::gsXBraidMultigridBase::solvecoarse
virtual void solvecoarse(const gsMatrix< T > &rhs, gsMatrix< T > &x, const int &numLevels)=0
Apply coarse solver (pure virtual method)
gismo::gsOptionList
Class which holds a list of parameters/options, and provides easy access to them. ...
Definition: gsOptionList.h:32
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_ddBtilde
std::vector< std::vector< gsMatrix< T > > > m_ddBtilde
std::vector of std::vector of block operator objects
Definition: gsXBraidMultigrid.h:459
gismo::gsXBraidMultigrid::m_ILUT
std::vector< std::vector< gsSparseMatrix< T > > > m_ILUT
std::vector of factorized operators
Definition: gsXBraidMultigrid.h:435
gismo::gsExprAssembler::setIntegrationElements
void setIntegrationElements(const gsMultiBasis< T > &mesh)
Sets the domain of integration.
Definition: gsExprAssembler.h:136
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::gsExprAssembler::solution
expr::gsFeSolution< T > solution
Solution type.
Definition: gsExprAssembler.h:61
gismo::gsExprAssembler::initSystem
void initSystem(const index_t numRhs=1)
Initializes the sparse system (sparse matrix and rhs)
Definition: gsExprAssembler.h:290
gismo::gsExprAssembler::getCoeff
variable getCoeff(const gsFunctionSet< T > &func)
Definition: gsExprAssembler.h:260
gismo::expr::gsFeVariable
Definition: gsExpressions.h:114
gismo::gsOptionList::getReal
Real getReal(const std::string &label) const
Reads value for option label from options.
Definition: gsOptionList.cpp:44
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::gsPreconditionerOp::Ptr
memory::shared_ptr< gsPreconditionerOp > Ptr
Shared pointer for gsLinearOperator.
Definition: gsPreconditioner.h:47