G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
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
38namespace gismo {
49template<typename T = real_t>
50gsVector<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
69template<typename T = real_t>
70void 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
88template<typename T>
89void 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
101inline void verboseLog(const std::string &message, const index_t &verbose) {
102 if (verbose > 0) { gsInfo << message << "\n"; }
103}
104
106template<short_t d, typename T= real_t>
107class 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
186template<short_t d, typename T = real_t>
187class 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
232template<short_t d, typename T>
233class 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
307template<short_t d, typename T>
308class 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
375template<short_t d, typename T>
376class 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
431template<short_t d, typename T>
432class 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
503template <typename Scalar = double>
504class 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
590namespace preAApp {
591
592template<typename Func>
593inline 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
606template<typename Scalar = double>
607class 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::expr::gsFeSpace
Definition gsExpressions.h:973
gismo::gsAsConstVector
Creates a mapped object or data pointer to a const vector without copying data.
Definition gsAsMatrix.h:285
gismo::gsAsVector
Creates a mapped object or data pointer to a vector without copying data.
Definition gsAsMatrix.h:239
gismo::gsBarrierCore
gsBarrierCore
Definition gsBarrierCore.h:108
gismo::gsBarrierCore::computePenaltyPatch
static gsMultiPatch< T > computePenaltyPatch(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
Penalty function-based method (1)
gismo::gsBarrierCore::computePenaltyPatch2
static gsMultiPatch< T > computePenaltyPatch2(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
Penalty function-based method (2)
gismo::gsBarrierCore::compute
static gsMultiPatch< T > compute(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
construct analysis-suitable parameterization
Definition gsBarrierCore.hpp:131
gismo::gsBarrierCore::computeAreaInterior
static T computeAreaInterior(const gsMultiPatch< T > &mp)
Computes the area of a multipatch representing a full domain.
Definition gsBarrierCore.hpp:96
gismo::gsBarrierCore::computeArea
static T computeArea(const gsMultiPatch< T > &mp)
Compute the area of the computational domain.
Definition gsBarrierCore.hpp:90
gismo::gsBarrierCore::computeVHPatch
static gsMultiPatch< T > computeVHPatch(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
variational harmonic method
gismo::gsBarrierCore::computeAreaBoundary
static T computeAreaBoundary(const gsMultiPatch< T > &mp)
Computes the area of a multipatch representing boundary curves.
Definition gsBarrierCore.hpp:116
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:1048
gismo::gsBarrierCore::computeBarrierPatch
static gsMultiPatch< T > computeBarrierPatch(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
Barrier function-based method.
gismo::gsBarrierCore::defaultOptions
static gsOptionList defaultOptions()
Default options.
Definition gsBarrierCore.hpp:26
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::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::gsDofMapper::freeSize
index_t freeSize() const
Returns the number of free (not eliminated) dofs.
Definition gsDofMapper.h:436
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::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::gsExprAssembler::geometryMap
gsExprHelper< T >::geometryMap geometryMap
Geometry map type.
Definition gsExprAssembler.h:65
gismo::gsExprAssembler::solution
expr::gsFeSolution< T > solution
Solution type.
Definition gsExprAssembler.h:68
gismo::gsMultiPatch
Container class for a set of geometry patches and their topology, that is, the interface connections ...
Definition gsMultiPatch.h:100
gismo::gsOptProblem
Class defining an optimization problem.
Definition gsOptProblem.h:25
gismo::gsOptionList
Class which holds a list of parameters/options, and provides easy access to them.
Definition gsOptionList.h:33
gismo::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
gismo::preAAParam
Definition gsBarrierCore.h:504
gismo::preAAParam::m
int m
Definition gsBarrierCore.h:513
gismo::preAAParam::updatePreconditionerStep
int updatePreconditionerStep
Definition gsBarrierCore.h:542
gismo::preAAParam::max_iterations
int max_iterations
Definition gsBarrierCore.h:537
gismo::preAAParam::epsilon
Scalar epsilon
Definition gsBarrierCore.h:521
gismo::preAAParam::epsilon_rel
Scalar epsilon_rel
Definition gsBarrierCore.h:530
gismo::preAAParam::check_param
void check_param() const
Definition gsBarrierCore.h:576
gismo::preAAParam::preAAParam
preAAParam()
Definition gsBarrierCore.h:559
gismo::preAAParam::usePreconditioning
bool usePreconditioning
Definition gsBarrierCore.h:552
gismo::preAApp::AndersonAcceleration
Anderson acceleration solver and its (preconditioned) variants.
Definition gsBarrierCore.h:607
gismo::preAApp::AndersonAcceleration::updateAlpha
void updateAlpha()
update the coefficients by solving a Least-Square problem
Definition gsBarrierCore.h:779
gismo::preAApp::AndersonAcceleration::updateGWithoutPreconditioning
void updateGWithoutPreconditioning(const Residual_t &F)
update fixed point function without preconditioning
Definition gsBarrierCore.h:757
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
gismo::preAApp::AndersonAcceleration::updateSolution
void updateSolution()
update the solution
Definition gsBarrierCore.h:814
gismo::preAApp::AndersonAcceleration::updateG
void updateG(const Residual_t &F, const Jacobian_t &Jacobian)
update fixed point function
Definition gsBarrierCore.h:734
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::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
gsBasis.h
Provides declaration of Basis abstract interface.
index_t
#define index_t
Definition gsConfig.h:32
GISMO_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition gsDebug.h:129
GISMO_ERROR
#define GISMO_ERROR(message)
Definition gsDebug.h:118
gsInfo
#define gsInfo
Definition gsDebug.h:43
gsExprAssembler.h
Generic expressions matrix assembly.
gsExprEvaluator.h
Generic expressions evaluator.
gsExprHelper.h
Generic expressions helper.
gsLinearAlgebra.h
This is the main header file that collects wrappers of Eigen for linear algebra.
gsStopwatch.h
Timing functions.
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::verboseLog
void verboseLog(const std::string &message, const index_t &verbose)
helper function to set optimizer options
Definition gsBarrierCore.h:101
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