G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsBarrierCore.h
Go to the documentation of this file.
1 
19 #pragma once
20 
21 
22 #include <gsUtils/gsStopwatch.h>
23 #include <gsCore/gsLinearAlgebra.h>
24 #include <gsCore/gsBasis.h>
25 
26 #include <gsAssembler/gsExprHelper.h>
27 #include <gsAssembler/gsExprEvaluator.h>
28 #include <gsAssembler/gsExprAssembler.h>
29 
30 
31 #ifdef gsHLBFGS_ENABLED
32 #include <gsHLBFGS/gsHLBFGS.h>
33 #endif
34 
35 //#include <gsLBFGSpp/gsLBFGSpp.h>
36 //#include <gsPreAA/gsPreAA.h>
37 
38 namespace gismo {
48 template<typename T = real_t>
50 gsVector<T> convertMultiPatchToFreeVector(const gsMultiPatch<T> &mp,
51  const gsDofMapper &mapper) {
52  gsVector<T> freeVec(mapper.freeSize());
53 
54  for (size_t iptch = 0; iptch != mp.nPatches(); iptch++) {
55  for (index_t idim = 0; idim != mp.targetDim(); idim++) {
56  for (size_t idof = 0; idof != mapper.patchSize(iptch, idim); idof++)
57  // if it is possible to just loop over the free index
58  // since this function is called very often during optimization
59  if (mapper.is_free(idof, iptch, idim)) {
60  index_t idx = mapper.index(idof, iptch, idim);
61  freeVec(idx) = mp.patch(iptch).coefs()(idof, idim);
62  }
63  }
64  }
65  return freeVec;
66 }
67 
69 template<typename T = real_t>
70 void convertFreeVectorToMultiPatch(const gsVector<T> &gsFreeVec,
71  const gsDofMapper &mapper,
72  gsMultiPatch<T> &mp) {
73  for (size_t iptch = 0; iptch != mp.nPatches(); iptch++) {
74  for (index_t idim = 0; idim != mp.targetDim(); idim++) {
75  for (size_t idof = 0; idof != mapper.patchSize(iptch, idim); idof++)
76  // if it is possible to just loop over the free index
77  // since this function is called very often during optimization
78  if (mapper.is_free(idof, iptch, idim)) {
79  index_t idx = mapper.index(idof, iptch, idim);
80  mp.patch(iptch).coefs()(idof, idim) = gsFreeVec(idx);
81  }
82  }
83  }
84 }
85 
87 #ifdef gsHLBFGS_ENABLED
88 template<typename T>
89 void setOptimizerOptions(gsHLBFGS<T> &optimizer, const gsOptionList &options) {
90  optimizer.options().setInt("MaxIterations",
91  options.askInt("qi_MaxIterations", 1e4));
92  optimizer.options().setReal("MinGradientLength",
93  options.askReal("qi_MinGradientLength", 1e-4));
94  optimizer.options().setReal("MinStepLength",
95  options.askReal("qi_MinStepLength", 1e-4));
96  optimizer.options().setInt("Verbose", options.askInt("Verbose", 0));
97 }
98 #endif
99 
101 inline void verboseLog(const std::string &message, const index_t &verbose) {
102  if (verbose > 0) { gsInfo << message << "\n"; }
103 }
104 
106 template<short_t d, typename T= real_t>
107 class gsBarrierCore
108 {
109  private:
110  typedef typename gsExprAssembler<T>::geometryMap geometryMap;
111  typedef typename gsExprAssembler<T>::space space;
112  typedef typename gsExprAssembler<T>::solution solution;
113 
114  public:
116  static gsMultiPatch<T> compute(const gsMultiPatch<T> &mp,
117  const gsDofMapper &mapper,
118  const gsOptionList &options);
119 
121  static T computeArea(const gsMultiPatch<T> &mp);
122 
124  static gsOptionList defaultOptions();
125 
126  private:
128  static gsMultiPatch<T> computeBarrierPatch(const gsMultiPatch<T> &mp,
129  const gsDofMapper &mapper,
130  const gsOptionList &options);
131 
133  static gsMultiPatch<T> computePenaltyPatch(const gsMultiPatch<T> &mp,
134  const gsDofMapper &mapper,
135  const gsOptionList &options);
136 
138  static gsMultiPatch<T> computePenaltyPatch2(const gsMultiPatch<T> &mp,
139  const gsDofMapper &mapper,
140  const gsOptionList &options);
141 
143  static gsMultiPatch<T> computePDEPatch(const gsMultiPatch<T> &mp,
144  const gsDofMapper &mapper,
145  const gsOptionList &options);
146 
148  static gsMultiPatch<T> computeVHPatch(const gsMultiPatch<T> &mp,
149  const gsDofMapper &mapper,
150  const gsOptionList &options);
151 
152  protected:
160  static T computeAreaInterior(const gsMultiPatch<T> &mp);
161 
169  static T computeAreaBoundary(const gsMultiPatch<T> &mp);
170 
171  private:
172  static void foldoverElimination(const gsMultiPatch<T> &mp,
173  const gsDofMapper &mapper,
174  gsVector<T> &initialGuessVector,
175  const T &scaledArea,
176  const gsOptionList &options);
177 
178  static void qualityImprovement(const gsMultiPatch<T> &mp,
179  const gsDofMapper &mapper,
180  gsVector<T> &initialGuessVector,
181  const T &scaledArea,
182  const gsOptionList &options);
183 };
184 
185 #ifdef gsHLBFGS_ENABLED
186 template<short_t d, typename T = real_t>
187 class gsObjFoldoverFree : public gsOptProblem<T> {
188 
189  private:
190  typedef typename gsExprAssembler<T>::geometryMap geometryMap;
191  typedef typename gsExprAssembler<T>::space space;
192  typedef typename gsExprAssembler<T>::solution solution;
193 
194  public:
195  gsObjFoldoverFree(const gsMultiPatch<T> &patches, gsDofMapper mapper);
196 
198  T evalObj(const gsAsConstVector<T> &u) const override;
199 
201  // and stores it in result.
202  void gradObj_into(const gsAsConstVector<T> &u,
203  gsAsVector<T> &result) const override;
204 
205  void setDelta(const T &delta) { m_eps = delta; };
206 
208  gsOptionList options() { return m_options; }
209 
211  void defaultOptions();
212 
214  void addOptions(const gsOptionList &options);
215 
217  void applyOptions(const gsOptionList &options);
218 
219  protected:
220  mutable gsMultiPatch<T> m_mp;
221  const gsDofMapper m_mapper;
222  const gsMultiBasis<T> m_mb;
223 
224  mutable gsExprEvaluator<T> m_evaluator;
225  mutable gsExprAssembler<T> m_assembler;
226 
227  gsOptionList m_options;
228 
229  T m_eps = T (1e-2);
230 };
231 
232 template<short_t d, typename T>
233 class gsObjQualityImprovePt : public gsOptProblem<T> {
234  private:
235  typedef typename gsExprAssembler<T>::geometryMap geometryMap;
236  typedef typename gsExprAssembler<T>::space space;
237  typedef typename gsExprAssembler<T>::solution solution;
238 
239  public:
240  gsObjQualityImprovePt(const gsMultiPatch<T> &patches, gsDofMapper mapper);
241 
243  T evalObj(const gsAsConstVector<T> &u) const override;
244 
246  // and stores it in result.
247  void gradObj_into(const gsAsConstVector<T> &u,
248  gsAsVector<T> &result) const override;
249 
251  gsOptionList options() { return m_options; }
252 
254  void defaultOptions();
255 
257  void addOptions(const gsOptionList &options);
258 
260  void applyOptions(const gsOptionList &options);
261 
262  private:
263  template<short_t _d>
264  typename std::enable_if<_d == 2, T>::type
265  evalObj_impl(const gsAsConstVector<T> &u) const;
266 
267  template<short_t _d>
268  typename std::enable_if<_d == 3, T>::type
269  evalObj_impl(const gsAsConstVector<T> &u) const;
270 
271  template<short_t _d>
272  typename std::enable_if<_d != 2 && _d != 3, T>::type
273  evalObj_impl(const gsAsConstVector<T> &u) const {
274  GISMO_ERROR("The dimension of target domain should be 2 or 3.");
275  }
276 
277  template<short_t _d>
278  typename std::enable_if<_d == 2, T>::type
279  gradObj_into_impl(const gsAsConstVector<T> &u,
280  gsAsVector<T> &result) const;
281 
282  template<short_t _d>
283  typename std::enable_if<_d == 3, T>::type
284  gradObj_into_impl(const gsAsConstVector<T> &u,
285  gsAsVector<T> &result) const;
286 
287  template<short_t _d>
288  typename std::enable_if<_d != 2 && _d != 3, T>::type
289  gradObj_into_impl(const gsAsConstVector<T> &u,
290  gsAsVector<T> &result) const {
291  GISMO_ERROR("The dimension of target domain should be 2 or 3.");
292  }
293 
294  protected:
295  mutable gsMultiPatch<T> m_mp;
296  const gsDofMapper m_mapper;
297  const gsMultiBasis<T> m_mb;
298 
299  mutable gsExprEvaluator<T> m_evaluator;
300  mutable gsExprAssembler<T> m_assembler;
301 
302  gsOptionList m_options;
303 
304  T m_lambda1 = 1.0, m_lambda2 = 1.0;
305 };
306 
307 template<short_t d, typename T>
308 class gsObjVHPt : public gsOptProblem<T> {
309  private:
310  typedef typename gsExprAssembler<T>::geometryMap geometryMap;
311  typedef typename gsExprAssembler<T>::space space;
312  typedef typename gsExprAssembler<T>::solution solution;
313 
314  public:
315  gsObjVHPt(const gsMultiPatch<T> &patches,
316  gsDofMapper mapper);
317 
318  T evalObj(const gsAsConstVector<T> &u) const final;
319 
320  void gradObj2_into(const gsAsConstVector<T> &u,
321  gsAsVector<T> &result) const;
322 
323  gsOptionList &options() { return m_options; }
324 
325  void setEps(T tol) { m_eps = tol; }
326 
327  void defaultOptions();
328 
329  void addOptions(const gsOptionList &options);
330 
331  void applyOptions(const gsOptionList &options);
332 
333  private:
334  template<short_t _d>
335  typename std::enable_if<_d == 2, T>::type
336  evalObj_impl(const gsAsConstVector<T> &u) const;
337 
338  template<short_t _d>
339  typename std::enable_if<_d == 3, T>::type
340  evalObj_impl(const gsAsConstVector<T> &u) const;
341 
342  template<short_t _d>
343  typename std::enable_if<_d != 2 && _d != 3, T>::type
344  evalObj_impl(
345  const gsAsConstVector<T> &u) const {GISMO_NO_IMPLEMENTATION; }
346 
347  template<short_t _d>
348  typename std::enable_if<_d == 2, T>::type
349  gradObj_into_impl(const gsAsConstVector<T> &u,
350  gsAsVector<T> &result) const;
351 
352  template<short_t _d>
353  typename std::enable_if<_d == 3, T>::type
354  gradObj_into_impl(const gsAsConstVector<T> &u,
355  gsAsVector<T> &result) const;
356 
357  template<short_t _d>
358  typename std::enable_if<_d != 2 && _d != 3, T>::type
359  gradObj_into_impl(const gsAsConstVector<T> &u,
360  gsAsVector<T> &result) const {GISMO_NO_IMPLEMENTATION; }
361 
362  protected:
363  mutable gsMultiPatch<T> m_mp;
364  const gsDofMapper m_mapper;
365  const gsMultiBasis<T> m_mb;
366 
367  mutable gsExprEvaluator<T> m_evaluator;
368  mutable gsExprAssembler<T> m_assembler;
369 
370  gsOptionList m_options;
371 
372  T m_lambda1, m_lambda2, m_eps = 4e-4; // 4e-4
373 };
374 
375 template<short_t d, typename T>
376 class gsObjPenaltyPt : public gsOptProblem<T> {
377  private:
378  typedef typename gsExprAssembler<T>::geometryMap geometryMap;
379  typedef typename gsExprAssembler<T>::space space;
380  typedef typename gsExprAssembler<T>::solution solution;
381 
382  public:
383  gsObjPenaltyPt(const gsMultiPatch<T> &patches,
384  gsDofMapper mapper);
385 
386  T evalObj(const gsAsConstVector<T> &u) const final;
387 
388  void gradObj_into(const gsAsConstVector<T> &u,
389  gsAsVector<T> &result) const final;
390 
391  gsOptionList &options() { return m_options; }
392 
393  void setEps(T tol) { m_eps = tol; }
394 
395  void defaultOptions();
396 
397  void addOptions(const gsOptionList &options);
398 
399  void applyOptions(const gsOptionList &options);
400 
401  private:
402 
403  template<short_t _d>
404  typename std::enable_if<_d == 2, T>::type
405  gradObj_into_impl(const gsAsConstVector<T> &u,
406  gsAsVector<T> &result) const;
407 
408  template<short_t _d>
409  typename std::enable_if<_d == 3, T>::type
410  gradObj_into_impl(const gsAsConstVector<T> &u,
411  gsAsVector<T> &result) const;
412 
413  template<short_t _d>
414  typename std::enable_if<_d != 2 && _d != 3, T>::type
415  gradObj_into_impl(const gsAsConstVector<T> &u,
416  gsAsVector<T> &result) const {GISMO_NO_IMPLEMENTATION; }
417 
418  protected:
419  mutable gsMultiPatch<T> m_mp;
420  const gsDofMapper m_mapper;
421  const gsMultiBasis<T> m_mb;
422 
423  mutable gsExprEvaluator<T> m_evaluator;
424  mutable gsExprAssembler<T> m_assembler;
425 
426  gsOptionList m_options;
427 
428  T m_lambda1, m_lambda2, m_eps = 1e-4; // 4e-4
429 };
430 
431 template<short_t d, typename T>
432 class gsObjPenaltyPt2 : public gsOptProblem<T> {
433  private:
434  typedef typename gsExprAssembler<T>::geometryMap geometryMap;
435  typedef typename gsExprAssembler<T>::space space;
436  typedef typename gsExprAssembler<T>::solution solution;
437 
438  public:
439  gsObjPenaltyPt2(const gsMultiPatch<T> &patches,
440  gsDofMapper mapper);
441 
442  T evalObj(const gsAsConstVector<T> &u) const final;
443 
444  void gradObj_into(const gsAsConstVector<T> &u,
445  gsAsVector<T> &result) const;
446 
447  void setEps(T tol) { m_eps = tol; }
448 
449  gsOptionList &options() { return m_options; }
450 
451  void defaultOptions();
452 
453  void addOptions(const gsOptionList &options);
454 
455  void applyOptions(const gsOptionList &options);
456 
457  private:
458  template<short_t _d>
459  typename std::enable_if<_d == 2, T>::type
460  evalObj_impl(const gsAsConstVector<T> &u) const;
461 
462  template<short_t _d>
463  typename std::enable_if<_d == 3, T>::type
464  evalObj_impl(const gsAsConstVector<T> &u) const;
465 
466  template<short_t _d>
467  typename std::enable_if<_d != 2 && _d != 3, T>::type
468  evalObj_impl(
469  const gsAsConstVector<T> &u) const {GISMO_NO_IMPLEMENTATION; }
470 
471  template<short_t _d>
472  typename std::enable_if<_d == 2, T>::type
473  gradObj_into_impl(const gsAsConstVector<T> &u,
474  gsAsVector<T> &result) const;
475 
476  template<short_t _d>
477  typename std::enable_if<_d == 3, T>::type
478  gradObj_into_impl(const gsAsConstVector<T> &u,
479  gsAsVector<T> &result) const;
480 
481  template<short_t _d>
482  typename std::enable_if<_d != 2 && _d != 3, T>::type
483  gradObj_into_impl(const gsAsConstVector<T> &u,
484  gsAsVector<T> &result) const {GISMO_NO_IMPLEMENTATION; }
485 
486  protected:
487  mutable gsMultiPatch<T> m_mp;
488  const gsDofMapper m_mapper;
489  const gsMultiBasis<T> m_mb;
490 
491  mutable gsExprEvaluator<T> m_evaluator;
492  mutable gsExprAssembler<T> m_assembler;
493 
494  gsOptionList m_options;
495 
496  T m_lambda1 = 1.0, m_lambda2 = 1.0, m_eps = 1e-3;
497 };
498 #endif
499 
503 template <typename Scalar = double>
504 class preAAParam {
505  public:
513  int m;
521  Scalar epsilon;
530  Scalar epsilon_rel;
537  int max_iterations;
542  int updatePreconditionerStep;
548 // Scalar beta;
552  bool usePreconditioning;
553 
554  public:
559  preAAParam() {
560  // clang-format off
561  m = 5;
562  epsilon = Scalar(1e-5);
563  epsilon_rel = Scalar(1e-5);
564  max_iterations = 100;
565  updatePreconditionerStep = 10;
566 // beta = Scalar(1.0);
567  usePreconditioning = false;
568  // clang-format on
569  }
570 
576  inline void check_param() const {
577  if (m < 0)
578  throw std::invalid_argument("'m' must be non-negative");
579  if (epsilon < 0)
580  throw std::invalid_argument("'epsilon' must be non-negative");
581  if (epsilon_rel < 0)
582  throw std::invalid_argument("'epsilon_rel' must be non-negative");
583  if (max_iterations < 0)
584  throw std::invalid_argument("'max_iterations' must be non-negative");
585 // if (beta < 0 || beta > 1.0)
586 // throw std::invalid_argument("'beta' must be between 0 and 1");
587  }
588 };
589 
590 namespace preAApp {
591 
592 template<typename Func>
593 inline void measureTime(Func func) {
594  auto beforeTime = std::chrono::high_resolution_clock::now();
595 
596  func(); // Call the function
597 
598  auto afterTime = std::chrono::high_resolution_clock::now();
599  auto duration = std::chrono::duration<double>(afterTime - beforeTime).count();
600  printf("\nTime passed by %.8f sec. \n", duration);
601 }
602 
606 template<typename Scalar = double>
607 class AndersonAcceleration {
608  public:
609  typedef std::function<gsSparseMatrix<Scalar>(gsVector<Scalar> const &)> Jacobian_t;
610  typedef std::function<gsVector<Scalar>(gsVector<Scalar> const &)> Residual_t;
611 
612  explicit AndersonAcceleration(const preAAParam<Scalar> &param) :
613  m_param(param) {
614  checkParam();
615  }
616 
618  const gsVector<Scalar> &compute(const gsVector<Scalar> &u0, const Residual_t &F,
619  const Jacobian_t &Jacobian) {
620  // print the current solver information
621  if (m_printInfo) { printSolverInfo();}
622 
623  // initialize the solver and preallocate memory
624  initialize(u0);
625 
626  // Iteration 0
627  updateG(F, Jacobian);
628 
629  // check if the initial guess is the solution
630  if (m_currResidualNorm < m_param.epsilon) {
631  printf("You are lucky, the initial guess is exactly the solution.\n\n\n");
632  return u0;
633  } else {
634  m_solution = m_currentG;
635  }
636 
637  // print iteration information
638  if (m_printInfo) {
639  printf("Iter. ||F(x)|| \n");
640  printIterationInfo();
641  }
642 
643  // Start iteration
644  startIteration(F, Jacobian);
645 
646  return m_solution;
647  }
648 
649  inline void enableIterationInfoPrinting() { m_printInfo = true; }
650  inline void disableIterationInfoPrinting() { m_printInfo = false; }
651 
652  private:
653  inline void initialize(const gsVector<Scalar> &u0) {
654  m_solution = u0;
655  m_dim = m_solution.size();
656  m_iter = 0;
657  m_columnIndex = 0;
658 
659  m_currentF.resize(m_dim);
660  m_prevdG.resize(m_dim, m_param.m);
661  m_prevdF.resize(m_dim, m_param.m);
662 
663  m_normalEquationMatrix.resize(m_param.m, m_param.m);
664  m_alpha.resize(m_param.m);
665  m_scaledF.resize(m_param.m);
666 
667  if (m_param.usePreconditioning) { m_preconditioner.resize(m_dim, m_dim); }
668  }
669 
670  inline void checkParam() {
671  m_param.check_param();
672  }
673 
674  inline void printSolverInfo() {
675  printf("\nAnderson Acceleration SOLVER: parameter settings... \n");
676  printf("depth = %d\n", m_param.m);
677  printf("use preconditioner = %d\n",
678  m_param.usePreconditioning);
679  printf("update preconditioner step = %d\n\n",
680  m_param.updatePreconditionerStep);
681  }
682 
683  inline void startIteration(const Residual_t &F, const Jacobian_t &Jacobian) {
684  m_iter = 1;
685  if (m_param.m == 0) {
686  performPicardIteration(F, Jacobian);
687  } else {
688  performAAIteration(F, Jacobian);
689  }
690  }
691 
692  inline void performPicardIteration(const Residual_t &F,
693  const Jacobian_t &Jacobian) {
694  while (m_iter < m_param.max_iterations
695  && m_currResidualNorm > m_param.epsilon) {
696  updateG(F, Jacobian);
697 
698  // update the solution
699  m_solution = m_currentG;
700 
701  printIterationInfo();
702 
703  ++m_iter;
704  trackIterationInfo();
705  }
706  }
707 
708  inline void performAAIteration(const Residual_t &F,
709  const Jacobian_t &Jacobian) {
710  m_prevdF.col(0) = -m_currentF;
711  m_prevdG.col(0) = -m_currentG;
712 
713  while (m_iter < m_param.max_iterations &&
714  m_currResidualNorm > m_param.epsilon) {
715  updateG(F, Jacobian);
716 
717  printIterationInfo();
718 
719  updatePrevdFAndPrevdG();
720 
721  // update alpha and solution (compute normal equation)
722  updateAlpha();
723  updateSolution();
724 
725  // update the column indices
726  updateColumnIndices();
727 
728  ++m_iter;
729  trackIterationInfo();
730  }
731  }
732 
734  inline void updateG(const Residual_t &F, const Jacobian_t &Jacobian) {
735  if (m_param.usePreconditioning) {
736  updateGWithPreconditioning(F, Jacobian);
737  } else {
738  updateGWithoutPreconditioning(F);
739  }
740  m_currentG = m_currentF + m_solution;
741  }
742 
744  inline void updateGWithPreconditioning(const Residual_t &F,
745  const Jacobian_t &Jacobian) {
746  if (!(m_iter % m_param.updatePreconditionerStep)) {
747  m_preconditioner = Jacobian(m_solution);
748  m_linearSolverPreconditioning.compute(m_preconditioner);
749  }
750  gsVector<Scalar> residual = F(m_solution);
751  m_currResidualNorm = residual.norm();
752 
753  m_currentF = -m_linearSolverPreconditioning.solve(residual);
754  }
755 
757  inline void updateGWithoutPreconditioning(const Residual_t &F) {
758  m_currentF = F(m_solution);
759  m_currResidualNorm = m_currentF.norm();
760  }
761 
762  inline void printIterationInfo() {
763  if (m_printInfo) {
764  printf(" %d %.4e\n", m_iter, cast<Scalar,double>( m_currResidualNorm ));
765  }
766  }
767 
768  inline void updatePrevdFAndPrevdG() {
769  m_prevdF.col(m_columnIndex) += m_currentF;
770  m_prevdG.col(m_columnIndex) += m_currentG;
771 
772  // scale previous dF for better numerical stability
773  Scalar scale = std::max(EPSILON, m_prevdF.col(m_columnIndex).norm());
774  m_scaledF(m_columnIndex) = scale;
775  m_prevdF.col(m_columnIndex) /= scale;
776  }
777 
779  inline void updateAlpha() {
780  // compute m_mk
781  m_mk = std::min(m_param.m, m_iter);
782 
783  if (m_mk == 1) {
784  m_alpha(0) = 0;
785  Scalar dF_squaredNorm = m_prevdF.col(m_columnIndex).squaredNorm();
786  m_normalEquationMatrix(0, 0) = dF_squaredNorm;
787  Scalar dF_norm = math::sqrt(dF_squaredNorm);
788 
789  // For better numerical stability
790  if (dF_norm > EPSILON) {
791  // compute alpha = (dF * F) / (dF * dF)
792  m_alpha(0) = (m_prevdF.col(m_columnIndex) / dF_norm).dot(
793  m_currentF / dF_norm);
794  }
795  } else {
796  // Update the normal equation matrix
797  // for the column and row corresponding to the new dF column.
798  // note: only one column and one row are needed to be updated.
799  gsVector<Scalar> new_inner_prod = (m_prevdF.col(m_columnIndex).transpose()
800  * m_prevdF.block(0, 0, m_dim, m_mk)).transpose();
801  m_normalEquationMatrix.block(m_columnIndex, 0, 1, m_mk) =
802  new_inner_prod.transpose();
803  m_normalEquationMatrix.block(0, m_columnIndex, m_mk, 1) =
804  new_inner_prod;
805 
806  // Solve normal equation: A^{T} A x = A^{T} b
807  m_linearSolver.compute(m_normalEquationMatrix.block(0, 0, m_mk, m_mk));
808  m_alpha.head(m_mk) = m_linearSolver.solve(
809  m_prevdF.block(0, 0, m_dim, m_mk).transpose() * m_currentF);
810  }
811  }
812 
814  inline void updateSolution() {
815  // Update the current solution (x) using the rescaled alpha
816  m_solution = m_currentG - m_prevdG.block(0, 0, m_dim, m_mk) *
817  ((m_alpha.head(m_mk).array()
818  / m_scaledF.head(m_mk).array()).matrix());
819  }
820 
821  void updateColumnIndices() {
822  m_columnIndex = (m_columnIndex + 1) % m_param.m;
823  m_prevdF.col(m_columnIndex) = -m_currentF;
824  m_prevdG.col(m_columnIndex) = -m_currentG;
825  }
826 
827  inline void trackIterationInfo() {
828  if (m_trackResidualNorm) { m_residualList.push_back(m_currResidualNorm); }
829  }
830 
831  // Linear solver for the Least-Square problem
832  gsEigen::FullPivLU<gsMatrix<Scalar>> m_linearSolver;
833  // Parameters
834  const preAAParam<Scalar> m_param;
835  int m_dim = -1;
836  // Iteration
837  int m_iter = 0;
838  int m_mk = 0;
839  int m_columnIndex = 0;
840  // Residual
841  Scalar m_currResidualNorm;
842  // Solution
843  gsVector<Scalar> m_solution;
844  gsVector<Scalar> m_currentF;
845  gsVector<Scalar> m_currentG;
846  // Anderson extrapolation
847  gsMatrix<Scalar> m_prevdG;
848  gsMatrix<Scalar> m_prevdF;
849  gsMatrix<Scalar> m_normalEquationMatrix;
850  gsVector<Scalar> m_alpha;
851  gsVector<Scalar> m_scaledF;
852  // Preconditioning
853  gsSparseMatrix<Scalar> m_preconditioner;
854  gsEigen::SparseLU<gsSparseMatrix<Scalar>> m_linearSolverPreconditioning;
855  // Print information
856  bool m_printInfo = true;
857 
858  // track residual norm history
859  bool m_trackResidualNorm = false;
860  std::vector<Scalar> m_residualList;
861 
862  const Scalar EPSILON = 1e-14;
863 };
864 
865 }
866 
867 }// namespace gismo
868 
869 #ifndef GISMO_BUILD_LIB
870 #include GISMO_HPP_HEADER(gsBarrierCore.hpp)
871 #endif
gismo::gsBarrierCore
gsBarrierCore
Definition: gsBarrierCore.h:107
gismo::gsBarrierCore::computeAreaBoundary
static T computeAreaBoundary(const gsMultiPatch< T > &mp)
Computes the area of a multipatch representing boundary curves.
Definition: gsBarrierCore.hpp:114
gismo::preAAParam::max_iterations
int max_iterations
Definition: gsBarrierCore.h:537
gismo::preAAParam::check_param
void check_param() const
Definition: gsBarrierCore.h:576
gismo::gsOptProblem
Class defining an optimization problem.
Definition: gsOptProblem.h:24
gismo::preAAParam::epsilon
Scalar epsilon
Definition: gsBarrierCore.h:521
GISMO_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition: gsDebug.h:129
gismo::gsExprAssembler::geometryMap
gsExprHelper< T >::geometryMap geometryMap
Geometry map type.
Definition: gsExprAssembler.h:58
gismo::preAAParam
Definition: gsBarrierCore.h:504
gismo::preAApp::AndersonAcceleration::updateAlpha
void updateAlpha()
update the coefficients by solving a Least-Square problem
Definition: gsBarrierCore.h:779
gismo::gsBarrierCore::computePDEPatch
static gsMultiPatch< T > computePDEPatch(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
PDE-based methods, including H2 and H1.
Definition: gsBarrierCore.hpp:1046
gismo::gsDofMapper::is_free
bool is_free(index_t i, index_t k=0, index_t c=0) const
Returns true if local dof i of patch k is not eliminated.
Definition: gsDofMapper.h:382
gismo::preAApp::AndersonAcceleration::compute
const gsVector< Scalar > & compute(const gsVector< Scalar > &u0, const Residual_t &F, const Jacobian_t &Jacobian)
perform Anderson acceleration iteration
Definition: gsBarrierCore.h:618
gsBasis.h
Provides declaration of Basis abstract interface.
gismo::preAApp::AndersonAcceleration
Anderson acceleration solver and its (preconditioned) variants.
Definition: gsBarrierCore.h:607
index_t
#define index_t
Definition: gsConfig.h:32
gsExprHelper.h
Generic expressions helper.
gismo::preAApp::AndersonAcceleration::updateGWithPreconditioning
void updateGWithPreconditioning(const Residual_t &F, const Jacobian_t &Jacobian)
update fixed point function with preconditioning
Definition: gsBarrierCore.h:744
gismo::preAApp::AndersonAcceleration::updateGWithoutPreconditioning
void updateGWithoutPreconditioning(const Residual_t &F)
update fixed point function without preconditioning
Definition: gsBarrierCore.h:757
gismo::preAApp::AndersonAcceleration::updateG
void updateG(const Residual_t &F, const Jacobian_t &Jacobian)
update fixed point function
Definition: gsBarrierCore.h:734
gsStopwatch.h
Timing functions.
gismo::gsDofMapper
Maintains a mapping from patch-local dofs to global dof indices and allows the elimination of individ...
Definition: gsDofMapper.h:68
gismo::verboseLog
void verboseLog(const std::string &message, const index_t &verbose)
helper function to set optimizer options
Definition: gsBarrierCore.h:101
gismo::gsBarrierCore::computeVHPatch
static gsMultiPatch< T > computeVHPatch(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
variational harmonic method
gismo::preAAParam::epsilon_rel
Scalar epsilon_rel
Definition: gsBarrierCore.h:530
gsExprEvaluator.h
Generic expressions evaluator.
gismo::gsDofMapper::index
index_t index(index_t i, index_t k=0, index_t c=0) const
Returns the global dof index associated to local dof i of patch k.
Definition: gsDofMapper.h:325
gismo::gsAsVector
Creates a mapped object or data pointer to a vector without copying data.
Definition: gsLinearAlgebra.h:129
gismo::convertFreeVectorToMultiPatch
void convertFreeVectorToMultiPatch(const gsVector< T > &gsFreeVec, const gsDofMapper &mapper, gsMultiPatch< T > &mp)
Convert free control points from a vector into a multi-patch.
Definition: gsBarrierCore.h:70
gismo::gsBarrierCore::computeAreaInterior
static T computeAreaInterior(const gsMultiPatch< T > &mp)
Computes the area of a multipatch representing a full domain.
Definition: gsBarrierCore.hpp:94
gismo::gsMultiPatch::targetDim
short_t targetDim() const
Dimension of the target space.
Definition: gsMultiPatch.h:192
gismo::gsBarrierCore::computeBarrierPatch
static gsMultiPatch< T > computeBarrierPatch(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
Barrier function-based method.
gismo::gsMultiPatch::patch
gsGeometry< T > & patch(size_t i) const
Return the i-th patch.
Definition: gsMultiPatch.h:226
gismo::gsMultiPatch::nPatches
size_t nPatches() const
Number of patches.
Definition: gsMultiPatch.h:208
gismo::preAApp::AndersonAcceleration::updateSolution
void updateSolution()
update the solution
Definition: gsBarrierCore.h:814
gismo::gsBarrierCore::compute
static gsMultiPatch< T > compute(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
construct analysis-suitable parameterization
Definition: gsBarrierCore.hpp:129
gsInfo
#define gsInfo
Definition: gsDebug.h:43
gismo::preAAParam::updatePreconditionerStep
int updatePreconditionerStep
Definition: gsBarrierCore.h:542
gismo::expr::gsFeSpace
Definition: gsDirichletValues.h:23
gismo::gsMultiPatch
Container class for a set of geometry patches and their topology, that is, the interface connections ...
Definition: gsMultiPatch.h:33
gismo::gsDofMapper::patchSize
size_t patchSize(const index_t k, const index_t c=0) const
Returns the total number of patch-local DoFs that live on patch k for component c.
Definition: gsDofMapper.h:480
gismo::gsBarrierCore::defaultOptions
static gsOptionList defaultOptions()
Default options.
Definition: gsBarrierCore.hpp:24
GISMO_ERROR
#define GISMO_ERROR(message)
Definition: gsDebug.h:118
gismo::gsAsConstVector
Creates a mapped object or data pointer to a const vector without copying data.
Definition: gsLinearAlgebra.h:130
gsLinearAlgebra.h
This is the main header file that collects wrappers of Eigen for linear algebra.
gismo::gsBarrierCore::computePenaltyPatch2
static gsMultiPatch< T > computePenaltyPatch2(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
Penalty function-based method (2)
gismo::preAAParam::preAAParam
preAAParam()
Definition: gsBarrierCore.h:559
gismo::gsOptionList
Class which holds a list of parameters/options, and provides easy access to them. ...
Definition: gsOptionList.h:32
gismo::gsDofMapper::freeSize
index_t freeSize() const
Returns the number of free (not eliminated) dofs.
Definition: gsDofMapper.h:436
gismo::convertMultiPatchToFreeVector
gsVector< T > convertMultiPatchToFreeVector(const gsMultiPatch< T > &mp, const gsDofMapper &mapper)
Computes a patch parametrization given a set of boundary geometries. Parametrization is not guarantee...
Definition: gsBarrierCore.h:50
gismo::preAAParam::usePreconditioning
bool usePreconditioning
Definition: gsBarrierCore.h:552
gismo::gsExprAssembler::solution
expr::gsFeSolution< T > solution
Solution type.
Definition: gsExprAssembler.h:61
gismo::gsBarrierCore::computeArea
static T computeArea(const gsMultiPatch< T > &mp)
Compute the area of the computational domain.
Definition: gsBarrierCore.hpp:88
gsExprAssembler.h
Generic expressions matrix assembly.
gismo::preAAParam::m
int m
Definition: gsBarrierCore.h:513
gismo::gsBarrierCore::computePenaltyPatch
static gsMultiPatch< T > computePenaltyPatch(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
Penalty function-based method (1)