G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsBarrierCore.hpp
Go to the documentation of this file.
1 
19 #pragma once
20 
21 namespace gismo {
22 
23 template<short_t d, typename T>
25  gsOptionList options;
26 
27  // Print output 0: no print, 1: summary, 2: iterations and summary
28  options.addInt("Verbose", "Print output level", 0);
29 
30  // Initialization Method: 0 Coons' patch (default), 1 Spring patch, 2: Cross-Ap. patch
31  options.addInt("InitialMethod", "Initialization method", 0);
32 
33  // Parameterization Method: 0 Barrier patch (default), 1 Penalty patch, 2: PDE patch
34  options.addInt("ParamMethod", "Parameterization method", 0);
35 
36  // Parameter and stopping criteria for foldover elimination step
37  options.addReal("ff_Delta", "Delta for foldover-free optimization", 0.05);
38  options.addInt("ff_MaxIterations",
39  "Max iterations for quality improvement",
40  1e4);
41  options.addReal("ff_MinGradientLength",
42  "Min gradient length for foldover-free optimization",
43  1e-20);
44  options.addReal("ff_MinStepLength",
45  "Min step length for foldover-free optimization",
46  1e-20);
47 
48  // Parameters and stopping criteria for quality improvement step
49  options.addInt("qi_MaxIterations",
50  "Max iterations for quality improvement",
51  1e4);
52  options.addReal("qi_MinGradientLength",
53  "Min gradient length for quality improvement",
54  1e-4);
55  options.addReal("qi_MinStepLength",
56  "Min step length for quality improvement",
57  1e-4);
58 
59  // Set quadrature rules for objective function and gradient evaluation
60  options.addReal("quA",
61  "Quadrature points: quA*deg + quB; For patchRule: Order of target space",
62  1.0);
63  options.addInt("quB",
64  "Quadrature points: quA*deg + quB; For patchRule: Regularity of target space",
65  1);
66  options.addInt("quRule",
67  "Quadrature rule [1:GaussLegendre,2:GaussLobatto,3:PatchRule]",
68  1);
69  options.addInt("overInt", "Apply over-integration?", 0);
70 
71  // Preconditioner type for AA: 0: NO preconditioning, 1: Full Jacobian preconditioner, 2: Diagonal Jacobian preconditioner, 3: Diagonal Block Jacobian preconditioner
72  options.addInt("AAPreconditionType", "Preconditioner type for AA", 0);
73 
74  // Update the preconditioner every N_update steps
75  options.addInt("N_update", "update preconditioner every N_update steps", 10);
76 
77  // window size m for our preconditioned AA solver
78  options.addInt("AAwindowsize", "window size for preconditioned AA solver", 5);
79 
80  // need improve the parameterization quality for PDE-based method?
81  options.addSwitch("needPDEH1", "improve quality by H1 discrezation?", true);
82 
83  return options;
84 }
85 
87 template<short_t d, typename T>
89  return computeAreaInterior(mp);
90 }
91 
93 template<short_t d, typename T>
95  // Creating a multi-basis from the provided multi-patch
96  gsMultiBasis<T> multiBasis(multiPatch);
97 
98  // Initializing an expression evaluator
99  gsExprEvaluator<T> evaluator;
100 
101  // Setting integration elements
102  evaluator.setIntegrationElements(multiBasis);
103 
104  // Getting the geometry map
105  geometryMap geometry = evaluator.getMap(multiPatch);
106 
107  // Computing the area by integrating over the determinant of the Jacobian of the geometry map
108  T area = evaluator.integral(jac(geometry).det());
109 
110  return area;
111 }
112 
113 template<short_t d, typename T>
115 {
116  GISMO_ERROR("NO_IMPLEMENTATION");
117 }
118 
119 enum class ParamMethod {
120  PDEPatch,
121  PenaltyPatch,
122  PenaltyPatch2,
123  BarrierPatch,
124  VariationalHarmonicPatch,
125  // Add new methods here...
126 };
127 
128 template<short_t d, typename T>
130  const gsDofMapper &mapper,
131  const gsOptionList &options) {
132  gsMultiPatch<T> result;
133 
134  ParamMethod
135  method = static_cast<ParamMethod>(options.askInt("ParamMethod", 0));
136 
137  switch (method) {
138  case ParamMethod::PenaltyPatch: {
139 #ifdef gsHLBFGS_ENABLED
140  result = computePenaltyPatch(mp, mapper, options);
141 #else
142  GISMO_ERROR("PenaltyPatch not available without gsHLBFGS. Please "
143  "compile with gsHLBGFS enabled");
144 #endif
145  break;
146  }
147  case ParamMethod::PenaltyPatch2: {
148 #ifdef gsHLBFGS_ENABLED
149  result = computePenaltyPatch2(mp, mapper, options);
150 #else
151  GISMO_ERROR("PenaltyPatch2 not available without gsHLBFGS. Please "
152  "compile with gsHLBGFS enabled");
153 #endif
154  break;
155  }
156  case ParamMethod::BarrierPatch: {
157 #ifdef gsHLBFGS_ENABLED
158  result = computeBarrierPatch(mp, mapper, options);
159 #else
160  GISMO_ERROR("BarrierPatch not available without gsHLBFGS. Please "
161  "compile with gsHLBGFS enabled");
162 #endif
163  break;
164  }
165  case ParamMethod::VariationalHarmonicPatch: {
166 #ifdef gsHLBFGS_ENABLED
167  result = computeVHPatch(mp, mapper, options);
168 #else
169  GISMO_ERROR("VHPatch not available without gsHLBFGS. Please "
170  "compile with gsHLBGFS enabled");
171 #endif
172  break;
173  }
174  case ParamMethod::PDEPatch:
175  default: {
176  if (method != ParamMethod::PDEPatch) {
177  gsWarn << "Invalid ParamMethod value. Defaulting to PDEPatch.\n";
178  }
179  result = computePDEPatch(mp, mapper, options);
180  break;
181  }
182  }
183  return result;
184 }
185 
186 #ifdef gsHLBFGS_ENABLED
187 // Modified Variational Harmonic Method
188 template<short_t d, typename T>
191  const gismo::gsDofMapper &mapper,
192  const gismo::gsOptionList &options) {
193 
194  verboseLog("Start variational harmonic parameterization construction...",
195  options.askInt("Verbose", 0));
196 
197  // get initial guess vector
198  gsVector<T> initialGuessVector = convertMultiPatchToFreeVector(mp, mapper);
199 
200  gsObjVHPt<d, T> objVHPt(mp, mapper);
201  objVHPt.applyOptions(options);
202 
203  gsHLBFGS<T> optimizer(&objVHPt);
204  setOptimizerOptions<T>(optimizer, options);
205 
206  optimizer.solve(initialGuessVector);
207  gsMultiPatch<T> result = mp;
208  convertFreeVectorToMultiPatch<T>(optimizer.currentDesign(), mapper, result);
209 
210  verboseLog("Finished!", options.askInt("Verbose", 0));
211 
212  return result;
213 }
214 
215 template<short_t d, typename T>
216 gsMultiPatch<T>
217 gsBarrierCore<d, T>::computePenaltyPatch(const gsMultiPatch<T> &mp,
218  const gsDofMapper &mapper,
219  const gsOptionList &options) {
220 
221  verboseLog("Start penalty function-based parameterization construction...",
222  options.askInt("Verbose", 0));
223 
224  // Compute scaledArea and initial guess vector
225  T scaledArea = computeAreaInterior(mp);
226 
227  // get initial guess vector
228  gsVector<T> initialGuessVector = convertMultiPatchToFreeVector(mp, mapper);
229 
230  gsObjPenaltyPt<d, T> objPenaltyPt(mp, mapper);
231  gsOptionList thisOptions = options;
232  thisOptions.addReal("qi_lambda1", "Sets the lambda_1 value for Emips", 1.0);
233  thisOptions.addReal("qi_lambda2",
234  "Sets the lambda 2 value for Eunif",
235  1.0 / pow(scaledArea, 2));
236  objPenaltyPt.applyOptions(thisOptions);
237 
238  gsHLBFGS<T> optimizer(&objPenaltyPt);
239  setOptimizerOptions<T>(optimizer, options);
240 
241  optimizer.solve(initialGuessVector);
242  gsMultiPatch<T> result = mp;
243  convertFreeVectorToMultiPatch<T>(optimizer.currentDesign(), mapper, result);
244 
245  verboseLog("Finished!", options.askInt("Verbose", 0));
246 
247  return result;
248 }
249 
250 template<short_t d, typename T>
251 gsMultiPatch<T>
252 gsBarrierCore<d, T>::computePenaltyPatch2(const gsMultiPatch<T> &mp,
253  const gsDofMapper &mapper,
254  const gsOptionList &options) {
255 
256  verboseLog("Penalty function-based (2) parameterization construction ...",
257  options.askInt("Verbose", 0));
258 
259  // Compute scaledArea and initial guess vector
260  T scaledArea = computeAreaInterior(mp);
261 
262  // get initial guess vector
263  gsVector<T> initialGuessVector = convertMultiPatchToFreeVector(mp, mapper);
264 
265  gsObjPenaltyPt2<d, T> objPenaltyPt(mp, mapper);
266  objPenaltyPt.setEps(options.askReal("penaltyFactor", 0.01) * scaledArea);
267 
268  gsOptionList thisOptions = options;
269  thisOptions.addReal("qi_lambda1", "Sets the lambda_1 value for Emips", 1.0);
270  thisOptions.addReal("qi_lambda2",
271  "Sets the lambda 2 value for Eunif",
272  1.0 / pow(scaledArea, 2));
273  objPenaltyPt.applyOptions(thisOptions);
274 
275  gsHLBFGS<T> optimizer(&objPenaltyPt);
276  setOptimizerOptions<T>(optimizer, options);
277 
278  optimizer.solve(initialGuessVector);
279  gsMultiPatch<T> result = mp;
280  convertFreeVectorToMultiPatch<T>(optimizer.currentDesign(), mapper, result);
281 
282  verboseLog("Finished!", options.askInt("Verbose", 0));
283 
284  return result;
285 }
286 
287 template<short_t d, typename T>
288 gsMultiPatch<T>
289  gsBarrierCore<d,T>::computeBarrierPatch(const gsMultiPatch<T> &mp,
290  const gsDofMapper &mapper,
291  const gsOptionList &options) {
292 
293  // Compute scaledArea and initial guess vector
294  T scaledArea = computeAreaInterior(mp);
295 
296  // get initial guess vector
297  gsVector<T> initialGuessVector = convertMultiPatchToFreeVector(mp, mapper);
298 
299  // STEP 2: foldover elimination step
300  foldoverElimination(mp, mapper, initialGuessVector, scaledArea, options);
301 
302  // STEP 3: parameterization quality improvement
303  qualityImprovement(mp, mapper, initialGuessVector, scaledArea, options);
304 
305  // Update the result with optimized parameterization
306  gsMultiPatch<T> result = mp;
307  convertFreeVectorToMultiPatch(initialGuessVector, mapper, result);
308 
309  return result;
310 }
311 
312 template<short_t d, typename T>
313 void gsBarrierCore<d, T>::foldoverElimination(const gsMultiPatch<T> &mp,
314  const gsDofMapper &mapper,
315  gsVector<T> &initialGuessVector,
316  const T &scaledArea,
317  const gsOptionList &options) {
318 
319  verboseLog("Start foldover elimination step...",
320  options.askInt("Verbose", 0));
321  constexpr T EPSILON = 1e-20;
322  constexpr int MAX_ITER = 10;
323  gsObjFoldoverFree<d, T> objFoldoverFree(mp, mapper);
324  objFoldoverFree.addOptions(options);
325 
326  gsHLBFGS<T> optFoldoverFree(&objFoldoverFree);
327  optFoldoverFree.options().setInt("MaxIterations",
328  options.askInt("ff_MaxIterations", 1e4));
329  optFoldoverFree.options().setReal("MinGradientLength",
330  options.askReal("ff_MinGradientLength",
331  1e-12));
332  optFoldoverFree.options().setReal("MinStepLength",
333  options.askReal("ff_MinStepLength", 1e-12));
334  optFoldoverFree.options().setInt("Verbose", options.askInt("Verbose", 0));
335 
336  T Efoldover = std::numeric_limits<T>::max();
337  for (index_t it = 0; it < MAX_ITER; ++it) {
338  T delta = pow(0.1, it) * 5e-2 * scaledArea; // parameter delta
339  objFoldoverFree.setDelta(delta);
340 
341  optFoldoverFree.solve(initialGuessVector);
342 
343  Efoldover = optFoldoverFree.objective();
344  initialGuessVector = optFoldoverFree.currentDesign();
345 
346  if (Efoldover <= EPSILON) { break; }
347  }
348 
349  if (Efoldover > EPSILON) {
350  throw std::runtime_error(
351  "Maximum iterations reached. The foldover-energy value is " +
352  std::to_string(Efoldover) +
353  ". This suggests there may be issues with the input data.");
354  }
355 }
356 
357 template<short_t d, typename T>
358 void gsBarrierCore<d, T>::qualityImprovement(const gsMultiPatch<T> &mp,
359  const gsDofMapper &mapper,
360  gsVector<T> &initialGuessVector,
361  const T &scaledArea,
362  const gsOptionList &options) {
363  verboseLog("Start parameterization quality improvement step...",
364  options.askInt("Verbose", 0));
365  gsObjQualityImprovePt<d, T> objQualityImprovePt(mp, mapper);
366 
367  gsOptionList thisOptions = options;
368  thisOptions.addReal("qi_lambda1", "Sets the lambda_1 value for Emips", 1.0);
369  thisOptions.addReal("qi_lambda2",
370  "Sets the lambda 2 value for Eunif",
371  1.0 / pow(scaledArea, 2));
372  objQualityImprovePt.applyOptions(thisOptions);
373 
374  gsHLBFGS<T> optQualityImprovePt(&objQualityImprovePt);
375  setOptimizerOptions<T>(optQualityImprovePt, options);
376 
377  optQualityImprovePt.solve(initialGuessVector);
378  initialGuessVector = optQualityImprovePt.currentDesign();
379 }
380 
381 template<short_t d, typename T>
382 gsObjFoldoverFree<d, T>::gsObjFoldoverFree(const gsMultiPatch<T> &patches,
383  gsDofMapper mapper)
384  :
385  m_mp(patches),
386  m_mapper(std::move(mapper)),
387  m_mb(m_mp) {
388  defaultOptions();
389  m_assembler.setIntegrationElements(m_mb);
390  m_evaluator = gsExprEvaluator<T>(m_assembler);
391 }
392 
393 template<short_t d, typename T>
394 void gsObjFoldoverFree<d, T>::defaultOptions() {
395  // @Ye, make this reasonable default options
396  m_options.addReal("ff_Delta", "Sets the delta value", 1e-2);
397 }
398 
399 template<short_t d, typename T>
400 void gsObjFoldoverFree<d, T>::addOptions(const gsOptionList &options) {
401  m_options.update(options, gsOptionList::addIfUnknown);
402 }
403 
404 template<short_t d, typename T>
405 void gsObjFoldoverFree<d, T>::applyOptions(const gsOptionList &options) {
406  m_eps = m_options.getReal("ff_Delta");
407  m_evaluator.options().update(m_options, gsOptionList::addIfUnknown);
408 }
409 
410 template<short_t d, typename T>
411 T gsObjFoldoverFree<d, T>::evalObj(const gsAsConstVector<T> &u) const {
412  convertFreeVectorToMultiPatch<T>(u, m_mapper, m_mp);
413  geometryMap G = m_evaluator.getMap(m_mp);
414 
415  auto EfoldoverFree = (m_eps - jac(G).det()).ppartval();
416  return m_evaluator.integral(EfoldoverFree);
417 }
418 
419 template<short_t d, typename T>
420 void gsObjFoldoverFree<d, T>::gradObj_into(const gsAsConstVector<T> &u,
421  gsAsVector<T> &result) const {
422  convertFreeVectorToMultiPatch<T>(u, m_mapper, m_mp);
423  geometryMap G = m_assembler.getMap(m_mp);
424 
425  // Only call these once if their results don't change.
426  space space1 = m_assembler.getSpace(m_mb, d);
427  space1.setupMapper(m_mapper);
428 
429  // |J|' w.r.t. physical coordinates x and y
430  auto derJacDet = frprod2(space1, jac(G).tr().adj());
431 
432  gsConstantFunction<T> zeroFunc(gsVector<T>::Zero(d), d);
433  auto zeroVar = m_evaluator.getVariable(zeroFunc);
434 
435  auto Eder = ternary(m_eps - jac(G).det(), -derJacDet, space1 * zeroVar);
436 
437  m_assembler.initSystem();
438  m_assembler.assemble(Eder);
439 
440  result.resize(m_assembler.rhs().rows());
441  std::copy(m_assembler.rhs().data(),
442  m_assembler.rhs().data() + m_assembler.rhs().rows(),
443  result.data());
444 }
445 
446 template<short_t d, typename T>
447 gsObjQualityImprovePt<d, T>::gsObjQualityImprovePt(
448  const gsMultiPatch<T> &patches,
449  gsDofMapper mapper)
450  :
451  m_mp(patches),
452  m_mapper(std::move(mapper)),
453  m_mb(m_mp) {
454  m_assembler.setIntegrationElements(m_mb);
455  m_evaluator = gsExprEvaluator<T>(m_assembler);
456 // defaultOptions();
457 }
458 
459 template<short_t d, typename T>
460 void gsObjQualityImprovePt<d, T>::defaultOptions() {
461  // @Ye, make this reasonable default options
462  m_options.addReal("qi_lambda1", "Sets the lambda 1 value", 1.0);
463  m_options.addReal("qi_lambda2", "Sets the lambda 2 value", 1.0);
464 }
465 
466 template<short_t d, typename T>
467 void gsObjQualityImprovePt<d, T>::addOptions(const gsOptionList &options) {
468  m_options.update(options, gsOptionList::addIfUnknown);
469 }
470 
471 template<short_t d, typename T>
472 void gsObjQualityImprovePt<d, T>::applyOptions(const gsOptionList &options) {
473  m_options.update(options, gsOptionList::addIfUnknown);
474  m_lambda1 = m_options.getReal("qi_lambda1");
475  m_lambda2 = m_options.getReal("qi_lambda2");
476  m_evaluator.options().update(m_options, gsOptionList::addIfUnknown);
477 }
478 
479 template<short_t d, typename T>
480 T gsObjQualityImprovePt<d, T>::evalObj(const gsAsConstVector<T> &u) const {
481  return evalObj_impl<d>(u);
482 }
483 
484 template<short_t d, typename T>
485 template<short_t _d>
486 typename std::enable_if<_d == 2, T>::type
487 gsObjQualityImprovePt<d, T>::evalObj_impl(const gsAsConstVector<T> &u) const {
488  convertFreeVectorToMultiPatch<T>(u, m_mapper, m_mp);
489  geometryMap G = m_evaluator.getMap(m_mp);
490 
491  if (m_evaluator.min(jac(G).det()) < 0) {
492  return std::numeric_limits<T>::max();
493  } else {
494  auto Ewinslow = jac(G).sqNorm() / jac(G).det();
495  auto Euniform = pow(jac(G).det(), 2);
496 
497  return m_evaluator.integral(m_lambda1 * Ewinslow + m_lambda2 * Euniform);
498  }
499 }
500 
501 template<short_t d, typename T>
502 template<short_t _d>
503 typename std::enable_if<_d == 3, T>::type
504 gsObjQualityImprovePt<d, T>::evalObj_impl(const gsAsConstVector<T> &u) const {
505 // gsMultiPatch<T> mp = m_mp;
506  convertFreeVectorToMultiPatch<T>(u, m_mapper, m_mp);
507 
508  geometryMap G = m_evaluator.getMap(m_mp);
509 
510  if (m_evaluator.min(jac(G).det()) < 0) {
511  return std::numeric_limits<T>::max();
512  } else {
513  auto Euniform = pow(jac(G).det(), 2);
514  auto Ewinslow = 0.5 * (jac(G).sqNorm() * jac(G).inv().sqNorm());
515 
516 // // another objective function term - Jiao et al. 2011
517 // auto Ewinslow = jac(G).sqNorm() / pow(jac(G).det(), 2.0 / 3.0);
518 
519  return m_evaluator.integral(m_lambda1 * Ewinslow + m_lambda2 * Euniform);
520  }
521 }
522 
523 template<short_t d, typename T>
524 void gsObjQualityImprovePt<d, T>::gradObj_into(const gsAsConstVector<T> &u,
525  gsAsVector<T> &result) const {
526  gradObj_into_impl<d>(u, result);
527 }
528 
529 template<short_t d, typename T>
530 template<short_t _d>
531 typename std::enable_if<_d == 2, T>::type
532 gsObjQualityImprovePt<d, T>::gradObj_into_impl(const gsAsConstVector<T> &u,
533  gsAsVector<T> &result) const {
534  convertFreeVectorToMultiPatch<T>(u, m_mapper, m_mp);
535 
536  geometryMap G = m_assembler.getMap(m_mp);
537 
538  space space1 = m_assembler.getSpace(m_mb, d); // 1D space!!
539  space1.setupMapper(m_mapper);
540 
541  // |J|' w.r.t. physical coordinates x and y
542  auto derJacDet = frprod2(space1, jac(G).tr().adj());
543 
544  auto Ewinslow = jac(G).sqNorm() / jac(G).det();
545  auto derEwinslow = 2.0 / jac(G).det() * (frprod2(space1, jac(G))) -
546  Ewinslow.val() / jac(G).det() * derJacDet;
547  auto derEuniform = 2.0 * jac(G).det() * derJacDet;
548 
549  m_assembler.initSystem();
550  m_assembler.assemble(m_lambda1 * derEwinslow + m_lambda2 * derEuniform);
551  result = gsAsVector<T>(const_cast<T *>(m_assembler.rhs().data()),
552  m_assembler.rhs().rows());
553  return EXIT_SUCCESS;
554 }
555 
556 template<short_t d, typename T>
557 template<short_t _d>
558 typename std::enable_if<_d == 3, T>::type
559 gsObjQualityImprovePt<d, T>::gradObj_into_impl(const gsAsConstVector<T> &u,
560  gsAsVector<T> &result) const {
561  convertFreeVectorToMultiPatch<T>(u, m_mapper, m_mp);
562  geometryMap G = m_assembler.getMap(m_mp);
563 
564  space space1 = m_assembler.getSpace(m_mb, d); // 1D space!!
565  space1.setupMapper(m_mapper);
566 
567  // |J|' w.r.t. physical coordinates x and y
568  auto derJacDet = frprod2(space1, jac(G).tr().adj());
569 
570  auto derEwinslow = frprod2(space1, (jac(G).inv().sqNorm() * jac(G) - jac(G)
571  .sqNorm() * (jac(G).tr() * jac(G) * jac(G).tr()).inv()));
572 
573 // // gradient of another objective function term
574 // auto Ewinslow = jac(G).sqNorm() / pow(jac(G).det(), 2.0 / 3.0);
575 // auto derEwinslow =
576 // 2.0 * frprod2(space1, jac(G)) / pow(jac(G).det(), 2.0 / 3.0) -
577 // 2.0 * Ewinslow.val() * derJacDet / jac(G).det() / 3.0;
578  auto derEuniform = 2.0 * jac(G).det() * derJacDet;
579 
580  m_assembler.initSystem();
581  m_assembler.assemble(m_lambda1 * derEwinslow + m_lambda2 * derEuniform);
582  result = gsAsVector<T>(const_cast<T *>(m_assembler.rhs().data()),
583  m_assembler.rhs().rows());
584  return EXIT_SUCCESS;
585 }
586 
587 template<short_t d, typename T>
588 gsObjVHPt<d, T>::gsObjVHPt(const gsMultiPatch<T> &patches,
589  gsDofMapper mapper)
590  :
591  m_mp(patches),
592  m_mapper(std::move(mapper)),
593  m_mb(m_mp) {
594  defaultOptions();
595  m_assembler.setIntegrationElements(m_mb);
596  m_evaluator = gsExprEvaluator<T>(m_assembler);
597 }
598 
599 template<short_t d, typename T>
600 void gsObjVHPt<d, T>::defaultOptions() {
601  // @Ye, make this reasonable default options
602  m_options.addReal("qi_lambda1", "Sets the lambda 1 value", 1.0);
603  m_options.addReal("qi_lambda2", "Sets the lambda 2 value", 1.0);
604 }
605 
606 template<short_t d, typename T>
607 void gsObjVHPt<d, T>::addOptions(const gsOptionList &options) {
608  m_options.update(options, gsOptionList::addIfUnknown);
609 }
610 
611 template<short_t d, typename T>
612 void gsObjVHPt<d, T>::applyOptions(const gsOptionList &options) {
613  m_options.update(options, gsOptionList::addIfUnknown);
614  m_evaluator.options().update(m_options, gsOptionList::addIfUnknown);
615 }
616 
617 template<short_t d, typename T>
618 T gsObjVHPt<d, T>::evalObj(const gsAsConstVector<T> &u) const {
619  return evalObj_impl<d>(u);
620 }
621 
622 template<short_t d, typename T>
623 template<short_t _d>
624 typename std::enable_if<_d == 2, T>::type
625 gsObjVHPt<d, T>::evalObj_impl(const gsAsConstVector<T> &u) const {
626  // Convert the free vector to multi-patch
627  convertFreeVectorToMultiPatch<T>(u, m_mapper, m_mp);
628 
629  // Get the map for the geometry
630  geometryMap G = m_evaluator.getMap(m_mp);
631 
632  // Get and set up the 1D space
633  space space1 = m_assembler.getSpace(m_mb, d);
634  space1.setupMapper(m_mapper);
635 
636  // Compute the metric matrix
637  auto metricMat = jac(G).tr() * jac(G);
638 
639  // Compute the Hessian
640  auto hessG = hess(G);
641 
642  // Calculate the scale factor
643  auto scale = metricMat.trace();
644 
645  // Calculate the Lx value
646  auto Lx = hessG % metricMat.adj() / scale.val();
647 
648  // Build the nonlinear system
649  auto nonlinearSystem = frprod3(space1, Lx);
650 
651  // Initialize and assemble the system
652  m_assembler.initSystem();
653  m_assembler.assemble(nonlinearSystem);
654 
655  // Create the result vector
656  gsVector<T> result = gsAsVector<T>(
657  const_cast<T *>(m_assembler.rhs().data()),
658  m_assembler.rhs().rows());
659 
660  return result.squaredNorm();
661 }
662 
663 template<short_t d, typename T>
664 template<short_t _d>
665 typename std::enable_if<_d == 3, T>::type
666 gsObjVHPt<d, T>::evalObj_impl(const gsAsConstVector<T> &u) const {
668 }
669 
670 template<short_t d, typename T>
671 void gsObjVHPt<d, T>::gradObj2_into(const gsAsConstVector<T> &u,
672  gsAsVector<T> &result) const {
673  gradObj_into_impl<d>(u, result);
674 }
675 
676 template<short_t d, typename T>
677 gsObjPenaltyPt<d, T>::gsObjPenaltyPt(const gsMultiPatch<T> &patches,
678  gsDofMapper mapper)
679  :
680  m_mp(patches),
681  m_mapper(std::move(mapper)),
682  m_mb(m_mp) {
683  defaultOptions();
684  m_assembler.setIntegrationElements(m_mb);
685  m_evaluator = gsExprEvaluator<T>(m_assembler);
686 }
687 
688 template<short_t d, typename T>
689 void gsObjPenaltyPt<d, T>::defaultOptions() {
690  // @Ye, make this reasonable default options
691  m_options.addReal("qi_lambda1", "Sets the lambda 1 value", 1.0);
692  m_options.addReal("qi_lambda2", "Sets the lambda 2 value", 1.0);
693 }
694 
695 template<short_t d, typename T>
696 void gsObjPenaltyPt<d, T>::addOptions(const gsOptionList &options) {
697  m_options.update(options, gsOptionList::addIfUnknown);
698 }
699 
700 template<short_t d, typename T>
701 void gsObjPenaltyPt<d, T>::applyOptions(const gsOptionList &options) {
702  m_options.update(options, gsOptionList::addIfUnknown);
703  m_lambda1 = m_options.getReal("qi_lambda1");
704  m_lambda2 = m_options.getReal("qi_lambda2");
705  m_evaluator.options().update(m_options, gsOptionList::addIfUnknown);
706 }
707 
708 template<short_t d, typename T>
709 T gsObjPenaltyPt<d, T>::evalObj(const gsAsConstVector<T> &u) const {
710  // Convert the free vector to multipatch
711  convertFreeVectorToMultiPatch<T>(u, m_mapper, m_mp);
712 
713  // Generate map for the geometry
714  geometryMap G = m_evaluator.getMap(m_mp);
715 
716  // Generate epsilon value
717  gsConstantFunction<T> eps1(pow(m_eps, 2.0), d);
718  auto eps = m_evaluator.getVariable(eps1);
719 
720  // Calculate chi value - penalty function
721  auto
722  chi = 0.5 * (jac(G).det() + pow(eps.val() + pow(jac(G).det(), 2.0), 0.5));
723 
724  // Calculate Ewinslow and Euniform
725  auto Ewinslow = jac(G).sqNorm() / pow(chi, (2.0 / static_cast<T>(d)));
726  auto Euniform = pow(jac(G).det(), 2.0);
727 
728  return m_evaluator.integral(m_lambda1 * Ewinslow + m_lambda2 * Euniform);
729 // return evalObj_impl<d>(u);
730 }
731 
732 template<short_t d, typename T>
733 void gsObjPenaltyPt<d, T>::gradObj_into(const gsAsConstVector<T> &u,
734  gsAsVector<T> &result) const {
735  gradObj_into_impl<d>(u, result);
736 }
737 
738 template<short_t d, typename T>
739 template<short_t _d>
740 typename std::enable_if<_d == 2, T>::type
741 gsObjPenaltyPt<d, T>::gradObj_into_impl(const gsAsConstVector<T> &u,
742  gsAsVector<T> &result) const {
743 
744  // Convert the free vector to multipatch
745  convertFreeVectorToMultiPatch<T>(u, m_mapper, m_mp);
746 
747  // Get map and space
748  geometryMap G = m_assembler.getMap(m_mp);
749  space space1 = m_assembler.getSpace(m_mb, d);
750  space1.setupMapper(m_mapper);
751 
752  // Derivative of Jacobian determinant
753  auto derJacDet = jac(space1) % jac(G).tr().adj();
754 
755  // Define constants for evaluation
756  gsConstantFunction<T> eps1(pow(m_eps, 2.0), d);
757  gsConstantFunction<T> unit1(1.0, d);
758  auto eps = m_evaluator.getVariable(eps1);
759  auto unit = m_evaluator.getVariable(unit1);
760 
761  // Compute common term, chi and chip
762  auto commonTerm = pow(eps.val() + pow(jac(G).det(), 2.0), 0.5);
763  auto chi = 0.5 * (jac(G).det() + commonTerm);
764  auto chip = 0.5 * (unit.val() + jac(G).det() / commonTerm);
765 
766  // Compute Ewinslow
767  auto Ewinslow = jac(G).sqNorm() / chi;
768 
769  // Compute derivatives of Ewinslow and Euniform
770  auto derEwinslow = 1.0 / chi * (2.0 * frprod2(space1, jac(G)) -
771  Ewinslow.val() * chip * derJacDet);
772  auto derEuniform = 2.0 * jac(G).det() * derJacDet;
773 
774  // Assemble the system
775  m_assembler.initSystem();
776  m_assembler.assemble(m_lambda1 * derEwinslow + m_lambda2 * derEuniform);
777 
778  // Update the result
779  result = gsAsVector<T>(const_cast<T *>(m_assembler.rhs().data()),
780  m_assembler.rhs().rows());
781 
782  return EXIT_SUCCESS;
783 }
784 
785 template<short_t d, typename T>
786 template<short_t _d>
787 typename std::enable_if<_d == 3, T>::type
788 gsObjPenaltyPt<d, T>::gradObj_into_impl(const gsAsConstVector<T> &u,
789  gsAsVector<T> &result) const {
790  const T twoThirds = 2.0 / 3.0;
791 
792  // Convert the free vector to multipatch
793  convertFreeVectorToMultiPatch<T>(u, m_mapper, m_mp);
794 
795  // Get map and space
796  geometryMap G = m_assembler.getMap(m_mp);
797  space space1 = m_assembler.getSpace(m_mb, d); // 1D space!!
798  space1.setupMapper(m_mapper);
799 
800  // Derivative of Jacobian determinant
801  auto derJacDet = frprod2(space1, jac(G).tr().adj());
802 
803  // Define constants for evaluation
804  gsConstantFunction<T> eps1(pow(m_eps, 2.0), d);
805  gsConstantFunction<T> unit1(1.0, d);
806  auto eps = m_evaluator.getVariable(eps1);
807  auto unit = m_evaluator.getVariable(unit1);
808 
809  // Compute common term, chi and chip
810  auto commonTerm = pow(eps.val() + pow(jac(G).det(), 2.0), 0.5);
811  auto chi = 0.5 * (jac(G).det() + commonTerm);
812  auto chip = 0.5 * (unit.val() + jac(G).det() / commonTerm);
813 
814  // Compute Ewinslow
815  auto Ewinslow = jac(G).sqNorm() / pow(chi, twoThirds);
816 
817  // Compute derivatives of Ewinslow and Euniform
818  auto derEwinslow = 2.0 * frprod2(space1, jac(G)) / pow(chi, twoThirds) -
819  twoThirds * Ewinslow / chi * chip * derJacDet;
820  auto derEuniform = 2.0 * jac(G).det() * derJacDet;
821 
822  // Assemble the system
823  m_assembler.initSystem();
824  m_assembler.assemble(m_lambda1 * derEwinslow + m_lambda2 * derEuniform);
825 
826  // Update the result
827  result = gsAsVector<T>(const_cast<T *>(m_assembler.rhs().data()),
828  m_assembler.rhs().rows());
829 
830  return EXIT_SUCCESS;
831 }
832 
833 // gsObjPenaltyPt2: my penalty function
834 template<short_t d, typename T>
835 gsObjPenaltyPt2<d, T>::gsObjPenaltyPt2(const gsMultiPatch<T> &patches,
836  gsDofMapper mapper)
837  :
838  m_mp(patches),
839  m_mapper(std::move(mapper)),
840  m_mb(m_mp) {
841  defaultOptions();
842  m_assembler.setIntegrationElements(m_mb);
843  m_evaluator = gsExprEvaluator<T>(m_assembler);
844 }
845 
846 template<short_t d, typename T>
847 void gsObjPenaltyPt2<d, T>::defaultOptions() {
848  // @Ye, make this reasonable default options
849  m_options.addReal("qi_lambda1", "Sets the lambda 1 value", 1.0);
850  m_options.addReal("qi_lambda2", "Sets the lambda 2 value", 1.0);
851 }
852 
853 template<short_t d, typename T>
854 void gsObjPenaltyPt2<d, T>::addOptions(const gsOptionList &options) {
855  m_options.update(options, gsOptionList::addIfUnknown);
856 }
857 
858 template<short_t d, typename T>
859 void gsObjPenaltyPt2<d, T>::applyOptions(const gsOptionList &options) {
860  m_options.update(options, gsOptionList::addIfUnknown);
861  m_lambda1 = m_options.getReal("qi_lambda1");
862  m_lambda2 = m_options.getReal("qi_lambda2");
863  m_evaluator.options().update(m_options, gsOptionList::addIfUnknown);
864 }
865 
866 template<short_t d, typename T>
867 T gsObjPenaltyPt2<d, T>::evalObj(const gsAsConstVector<T> &u) const {
868  return evalObj_impl<d>(u);
869 }
870 
871 template<short_t d, typename T>
872 template<short_t _d>
873 typename std::enable_if<_d == 2, T>::type
874 gsObjPenaltyPt2<d, T>::evalObj_impl(const gsAsConstVector<T> &u) const {
875  // Convert the free vector to multipatch
876  convertFreeVectorToMultiPatch<T>(u, m_mapper, m_mp);
877 
878  // Get the map G
879  geometryMap G = m_evaluator.getMap(m_mp);
880 
881  // Constant function for evaluation
882  gsConstantFunction<T> epsilon(m_eps, d);
883  auto eps = m_evaluator.getVariable(epsilon);
884 
885  // Calculation of chi with ternary operation
886  auto chiPPart = eps * ((jac(G).det() - eps.val()).exp());
887  auto chi =
888  ternary(eps.val() - jac(G).det(), chiPPart.val(), jac(G).det().val());
889 
890  // Calculation of Ewinslow and Euniform
891  auto Ewinslow = jac(G).sqNorm() / chi;
892 // auto Euniform = chi + 1.0 / chi;
893  auto Euniform = pow(jac(G).det(), 2.0);
894 
895  // Evaluate the integral and return
896  return m_evaluator.integral(m_lambda1 * Ewinslow + m_lambda2 * Euniform);
897 }
898 
899 template<short_t d, typename T>
900 template<short_t _d>
901 typename std::enable_if<_d == 3, T>::type
902 gsObjPenaltyPt2<d, T>::evalObj_impl(const gsAsConstVector<T> &u) const {
903  // Convert the free vector to multipatch
904  convertFreeVectorToMultiPatch<T>(u, m_mapper, m_mp);
905 
906  // Get the map G
907  geometryMap G = m_evaluator.getMap(m_mp);
908 
909  // Define epsilon constant function
910  gsConstantFunction<T> epsilon(m_eps, d);
911  auto eps = m_evaluator.getVariable(epsilon);
912 
913  // Compute chi part and chi
914  auto chiPPart = eps * ((jac(G).det() - eps.val()).exp());
915  auto chi =
916  ternary(eps.val() - jac(G).det(), chiPPart.val(), jac(G).det().val());
917 
918  // Compute Ewinslow
919  auto Ewinslow = 0.5 * (jac(G).sqNorm() * jac(G).inv().sqNorm()) *
920  pow(jac(G).det(), 2.0) / pow(chi, 2.0);
921 
922  // Compute Euniform
923 // auto Euniform = chi + 1.0 / chi;
924  auto Euniform = pow(jac(G).det(), 2.0);
925 
926  // Compute and return the integral
927  return m_evaluator.integral(m_lambda1 * Ewinslow + m_lambda2 * Euniform);
928 }
929 
930 template<short_t d, typename T>
931 void gsObjPenaltyPt2<d, T>::gradObj_into(const gsAsConstVector<T> &u,
932  gsAsVector<T> &result) const {
933  gradObj_into_impl<d>(u, result);
934 }
935 
936 template<short_t d, typename T>
937 template<short_t _d>
938 typename std::enable_if<_d == 2, T>::type
939 gsObjPenaltyPt2<d, T>::gradObj_into_impl(const gsAsConstVector<T> &u,
940  gsAsVector<T> &result) const {
941  // Convert the free vector to multipatch
942  convertFreeVectorToMultiPatch<T>(u, m_mapper, m_mp);
943 
944  // Get the map G
945  geometryMap G = m_assembler.getMap(m_mp);
946 
947  // Set up the space mapper
948  space spaceMapper = m_assembler.getSpace(m_mb, d);
949  spaceMapper.setupMapper(m_mapper);
950 
951  // Compute the derivative of |J| with respect to physical coordinates x and y
952  auto derJacDet = frprod2(spaceMapper, jac(G).tr().adj());
953 
954  // Define constant functions for epsilon and unity
955  gsConstantFunction<T> epsilon(m_eps, d);
956  gsConstantFunction<T> unity(1.0, d);
957 
958  // Get the variables for epsilon and unity
959  auto eps = m_evaluator.getVariable(epsilon);
960  auto unit = m_evaluator.getVariable(unity);
961 
962  // Define chi and its derivative chip
963  auto chiPPart = eps * ((jac(G).det() - eps.val()).exp());
964  auto chi =
965  ternary(eps.val() - jac(G).det(), chiPPart.val(), jac(G).det().val());
966  auto chip = ternary(eps.val() - jac(G).det(), chiPPart.val(), unit.val());
967 
968  // Define Ewinslow and its derivative
969  auto Ewinslow = jac(G).sqNorm() / chi;
970  auto derEwinslow = 1.0 / chi * (2.0 * frprod2(spaceMapper, jac(G)) -
971  Ewinslow.val() * chip * derJacDet);
972 
973  // Define the derivative of Euniform
974 // auto derEuniform = (chip - chip / pow(chi, 2)) * derJacDet;
975  auto derEuniform = 2.0 * jac(G).det() * derJacDet;
976 
977  // Assemble the system and set the result
978  m_assembler.initSystem();
979  m_assembler.assemble(m_lambda1 * derEwinslow + m_lambda2 * derEuniform);
980  result = gsAsVector<T>(const_cast<T *>(m_assembler.rhs().data()),
981  m_assembler.rhs().rows());
982 
983  // Return success
984  return EXIT_SUCCESS;
985 }
986 
987 template<short_t d, typename T>
988 template<short_t _d>
989 typename std::enable_if<_d == 3, T>::type
990 gsObjPenaltyPt2<d, T>::gradObj_into_impl(const gsAsConstVector<T> &u,
991  gsAsVector<T> &result) const {
992  // Convert the free vector to multipatch
993  convertFreeVectorToMultiPatch<T>(u, m_mapper, m_mp);
994 
995  // Initialize the map G
996  geometryMap G = m_assembler.getMap(m_mp);
997 
998  // Initialize and setup the space
999  space space1 = m_assembler.getSpace(m_mb, d);
1000  space1.setupMapper(m_mapper);
1001 
1002  // Compute the derivative of |J| w.r.t. physical coordinates x and y
1003  auto derJacDet = frprod2(space1, jac(G).tr().adj());
1004 
1005  // Initialize constant functions
1006  gsConstantFunction<T> epsilonFunction(m_eps, d);
1007  gsConstantFunction<T> unityFunction(1.0, d);
1008 
1009  // Get the corresponding variables
1010  auto eps = m_evaluator.getVariable(epsilonFunction);
1011  auto unit = m_evaluator.getVariable(unityFunction);
1012 
1013  // Compute the chi part
1014  auto chiPPart = eps * ((jac(G).det() - eps.val()).exp());
1015 
1016  // Ternary operation to compute chi and chip
1017  auto chi =
1018  ternary(eps.val() - jac(G).det(), chiPPart.val(), jac(G).det().val());
1019  auto chip = ternary(eps.val() - jac(G).det(), chiPPart.val(), unit.val());
1020 
1021  // Define additional computations
1022  auto jacFrobNorm2 = jac(G).sqNorm() * jac(G).inv().sqNorm();
1023  auto derJacFrob2 = frprod2(space1,
1024  (jac(G).inv().sqNorm() * jac(G) - jac(G).sqNorm()
1025  * (jac(G).tr() * jac(G) * jac(G).tr()).inv()));
1026 
1027  auto derEwinslow = derJacFrob2 * pow(jac(G).det(), 2) / pow(chi, 2)
1028  + jacFrobNorm2 / pow(chi, 2)
1029  * (jac(G).det() - chip / chi * pow(jac(G).det(), 2)) * derJacDet;
1030 // auto derEuniform = (chip - chip / pow(chi, 2)) * derJacDet;
1031  auto derEuniform = 2.0 * jac(G).det() * derJacDet;
1032 
1033  // Initialize and assemble the system
1034  m_assembler.initSystem();
1035  m_assembler.assemble(m_lambda1 * derEwinslow + m_lambda2 * derEuniform);
1036 
1037  // Compute and return the result
1038  result = gsAsVector<T>(const_cast<T *>(m_assembler.rhs().data()),
1039  m_assembler.rhs().rows());
1040  return EXIT_SUCCESS;
1041 }
1042 #endif
1043 
1044 template<short_t d, typename T>
1045 gsMultiPatch<T>
1047  const gismo::gsDofMapper &mapper,
1048  const gismo::gsOptionList &options) {
1049 
1050  // get initial guess vector
1051  gsVector<T> initialGuessVector = convertMultiPatchToFreeVector(mp, mapper);
1052 
1053 // std::string solverSetting = "Parameters setting:\n";
1054 // solverSetting += "\t\t\t Window size = " +
1055 // std::to_string(options.getInt("AAwindowsize")) + "\n";
1056 // solverSetting += "\t\t\t Update preconditioner every " +
1057 // std::to_string(options.getInt("N_update")) + " steps \n";
1058 // solverSetting += "\t\t\t Preconditioner type: " +
1059 // std::to_string(options.getInt("AAPreconditionType")) + "\n";
1060 // solverSetting += "\t\t\t Need improve by H1? " +
1061 // std::to_string(options.getSwitch("needPDEH1")) + "\n";
1062 // verboseLog(solverSetting, options.askInt("Verbose", 0));
1063 
1064  verboseLog("PDE-based parameterization construction ...\n",
1065  options.askInt("Verbose", 0));
1066 
1067  gsExprEvaluator<T> evaluator;
1068  gsExprAssembler<T> assembler;
1069 
1070  gsMultiBasis<T> mb(mp);
1071  assembler.setIntegrationElements(mb);
1072 
1074  bc.setGeoMap(mp);
1075  for (auto bit = mp.bBegin(); bit != mp.bEnd(); ++bit) {
1076  bc.addCondition(*bit, condition_type::dirichlet, nullptr);
1077  }
1078 
1079  space space1 = assembler.getSpace(mb, d); // 1D space!!
1080  space1.setup(bc, dirichlet::homogeneous, 0);
1081 
1082  // Function for the Residual
1083  typedef std::function<gsSparseMatrix<T>(gsVector<T> const &)> Jacobian_t;
1084  typedef std::function<gsVector<T>(gsVector<T> const &)> Residual_t;
1085 
1086  gsMultiPatch<T> mpSubstitute = mp;
1087  Residual_t Residual = [&assembler, &mapper, &mpSubstitute, &space1](
1088  gsVector<T> const &x) {
1089  // Convert free vector to MultiPatch
1090  convertFreeVectorToMultiPatch<T>(x, mapper, mpSubstitute);
1091 
1092  // Calculate geometry map
1093  geometryMap G = assembler.getMap(mpSubstitute);
1094 
1095  // Compute metric matrix and Hessian of G
1096  auto metricMat = jac(G).tr() * jac(G);
1097  auto hessG = hess(G);
1098 
1099  // Compute scale factor
1100  auto scale = metricMat.trace();
1101 
1102  // Compute Lx
1103  auto Lx = hessG % metricMat.adj() / scale.val();
1104 
1105  // Compute nonlinear system
1106  auto nonlinearSystem = frprod3(space1, Lx);
1107 
1108  // Assemble the system
1109  assembler.initSystem();
1110  assembler.assemble(nonlinearSystem);
1111 
1112  // Return the assembled system
1113  return assembler.rhs();
1114  };
1115 
1116  Jacobian_t Jacobian;
1117  int preconditionerType = options.askInt("AAPreconditionType", 0);
1118  switch (preconditionerType) {
1119  case 1:
1120  Jacobian = [&assembler, &mapper, &mpSubstitute, &space1](gsVector<T> const &x)
1121  {
1122  // diagonal Jacobian matrix as a preconditioner
1123  convertFreeVectorToMultiPatch<T>(x, mapper, mpSubstitute);
1124  geometryMap G = assembler.getMap(mpSubstitute);
1125 
1126  auto jacMat = jacScaledLxDiag(space1, G);
1127 
1128  assembler.initSystem();
1129  assembler.assemble(jacMat);
1130  return assembler.matrix();
1131  };
1132  break;
1133 
1134  case 2:
1135  Jacobian = [&assembler, &mapper, &mpSubstitute, &space1](gsVector<T> const &x)
1136  {
1137  // diagonal-block Jacobian matrix as a preconditioner
1138  convertFreeVectorToMultiPatch<T>(x, mapper, mpSubstitute);
1139  geometryMap G = assembler.getMap(mpSubstitute);
1140 
1141  auto jacMat = jacScaledLxDiagBlock(space1, G);
1142 
1143  assembler.initSystem();
1144  assembler.assemble(jacMat);
1145  return assembler.matrix();
1146  };
1147  break;
1148 
1149  default:
1150  // analytical Jacobian matrix as a preconditioner
1151  Jacobian = [&assembler, &mapper, &mpSubstitute, &space1](gsVector<T> const &x)
1152  {
1153  convertFreeVectorToMultiPatch<T>(x, mapper, mpSubstitute);
1154  geometryMap G = assembler.getMap(mpSubstitute);
1155 
1156  auto jacMat = jacScaledLx(space1, G);
1157 
1158  assembler.initSystem();
1159  assembler.assemble(jacMat);
1160  return assembler.matrix();
1161  };
1162  }
1163 
1164  preAAParam<T> param;
1165  param.m = options.askInt("AAwindowsize", 5);
1166  param.usePreconditioning = true;
1167  param.updatePreconditionerStep = options.askInt("N_update", 5);
1168  param.epsilon = 1e-5;
1169 
1170  preAApp::AndersonAcceleration<T> AASolver(param);
1171  gsVector<T> solVector = AASolver.compute(initialGuessVector,Residual, Jacobian);
1172 
1173  if (options.askSwitch("needPDEH1", true))
1174  {
1175  verboseLog("\nStart parameterization improvement by H1 discrezation...",options.askInt("Verbose", 0));
1176  Residual = [&assembler, &mapper, &mpSubstitute, &space1](gsVector<T> const &x)
1177  {
1178  // H1 discretization
1179  convertFreeVectorToMultiPatch<T>(x, mapper, mpSubstitute);
1180  geometryMap G = assembler.getMap(mpSubstitute);
1181 
1182  auto invJacMat = jac(G).inv();
1183  auto jacU = jac(space1) * invJacMat;
1184  auto nonlinearSystem = jacU % invJacMat;
1185 
1186  assembler.initSystem();
1187  assembler.assemble(nonlinearSystem);
1188  return assembler.rhs();
1189  };
1190 
1191  preconditionerType = options.askInt("AAPreconditionType", 0);
1192  switch (preconditionerType) {
1193  case 1:
1194  Jacobian = [&assembler, &mapper, &mpSubstitute, &space1](gsVector<T> const &x)
1195  {
1196  // diagonal block Jacobian matrix as a preconditioner
1197  convertFreeVectorToMultiPatch<T>(x, mapper, mpSubstitute);
1198  geometryMap G = assembler.getMap(mpSubstitute);
1199 
1200  auto jacMat = jacScaledLxH1DiagBlock(space1, G);
1201 
1202  assembler.initSystem();
1203  assembler.assemble(jacMat);
1204 
1205  return assembler.matrix();
1206  };
1207  break;
1208 
1209  default:
1210  Jacobian = [&assembler, &mapper, &mpSubstitute,&space1](gsVector<T> const &x)
1211  {
1212  // analytical Jacobian matrix as a preconditioner
1213  convertFreeVectorToMultiPatch<T>(x, mapper, mpSubstitute);
1214  geometryMap G = assembler.getMap(mpSubstitute);
1215 
1216  auto jacMat = jacScaledLxH1(space1, G);
1217 
1218  assembler.initSystem();
1219  assembler.assemble(jacMat);
1220 
1221  return assembler.matrix();
1222  };
1223  }
1224 
1225  preAApp::AndersonAcceleration<T> AASolver2(param);
1226  solVector = AASolver2.compute(solVector, Residual, Jacobian);
1227  }
1228 
1229  gsMultiPatch<T> result = mp;
1230  convertFreeVectorToMultiPatch<T>(solVector, mapper, result);
1231 
1232  verboseLog("Finished!", options.askInt("Verbose", 0));
1233 
1234  return result;
1235 }
1236 
1237 namespace expr {
1238 
1239 template<typename E1, typename E2>
1240 class frprod2_expr;
1241 template<typename E1, typename E2>
1242 class frprod3_expr;
1243 template<class E0, class E1, class E2>
1244 class ternary_expr; // ternary expression
1245 template<class E>
1246 class jacScaledLx_expr;
1247 template<class E>
1248 class jacScaledLxDiag_expr;
1249 template<class E>
1250 class jacScaledLxDiagBlock_expr;
1251 template<class E>
1252 class jacScaledLxH1_expr;
1253 template<class E>
1254 class jacScaledLxH1DiagBlock_expr;
1255 
1256 /*
1257 
1258  Expression for the product of the jacob of space and matrix A. Return term
1259  is a column vector (d*n x 1). The result should be the same as the result
1260  from jac(space) % A, whereas avoids redundant multiply with zero-value components.
1261 
1262  For 2D case:
1263  frprod2(space,A) = [a_{11}*dN_1/dxi_1, a_{12}*dN_1/dxi_2, a_{11}*dN_2/dxi_1,
1264  a_{12}*dN_2/dxi_2,...,a_{11}*dN_n/dxi_1, a_{12}*dN_n/dxi_2, a_{21}*dN_1/dxi_1,
1265  a_{22}*dN_1/dxi_2, a_{21}*dN_2/dxi_1, a_{22}*dN_2/dxi_2,...,a_{21}*dN_n/dxi_1,
1266  a_{22}*dN_n/dxi_2]
1267  For 3D case:
1268  frprod2(space,A) = [a_{11}*dN_1/dxi_1, a_{12}*dN_1/dxi_2, a_{13}*dN_1/dxi_3,...
1269  a_{11}*dN_2/dxi_1, a_{12}*dN_2/dxi_1, a_{13}*dN_2/dxi_3,...,
1270  a_{11}*dN_n/dxi_1, a_{12}*dN_n/dxi_1, a_{13}*dN_n/dxi_3,
1271  a_{21}*dN_1/dxi_1, a_{22}*dN_1/dxi_2, a_{23}*dN_1/dxi_3,...
1272  a_{21}*dN_2/dxi_1, a_{22}*dN_2/dxi_1, a_{23}*dN_2/dxi_3,...,
1273  a_{21}*dN_n/dxi_1, a_{22}*dN_n/dxi_1, a_{23}*dN_n/dxi_3,
1274  a_{31}*dN_1/dxi_1, a_{32}*dN_1/dxi_2, a_{33}*dN_1/dxi_3,...
1275  a_{31}*dN_2/dxi_1, a_{32}*dN_2/dxi_1, a_{33}*dN_2/dxi_3,...,
1276  a_{31}*dN_n/dxi_1, a_{32}*dN_n/dxi_1, a_{33}*dN_n/dxi_3,]
1277 
1278 NOTE: _u should be a space, _v should NOT be a space (fix with assert)
1279 
1280  */
1281 template<typename E1, typename E2>
1282 class frprod3_expr : public _expr<frprod3_expr<E1, E2> > {
1283  public:
1284  typedef typename E2::Scalar Scalar;
1285  enum { ScalarValued = 0, Space = E1::Space, ColBlocks = 0 };
1286 
1287  private:
1288  typename E1::Nested_t _u;
1289  typename E2::Nested_t _v;
1290 
1291  mutable gsMatrix<Scalar> res, bGrads, b;
1292 
1293  public:
1294 
1295  frprod3_expr(_expr<E1> const &u, _expr<E2> const &v)
1296  : _u(u), _v(v) {
1297  // gsInfo << "expression is space ? "<<E1::Space <<"\n"; _u.print(gsInfo);
1298  // GISMO_ASSERT(_u.rows() == _v.rows(),
1299  // "Wrong dimensions "<<_u.rows()<<"!="<<_v.rows()<<" in % operation");
1300  // GISMO_ASSERT(_u.cols() == _v.cols(),
1301  // "Wrong dimensions "<<_u.cols()<<"!="<<_v.cols()<<" in % operation");
1302  }
1303 
1304  const gsMatrix<Scalar> &eval(const index_t k) const //todo: specialize for nb==1
1305  {
1306  auto A = _v.eval(k);
1307  b = _u.eval(k);
1308 // gsDebugVar(b);
1309 
1310  res.noalias() = b * A;
1311  return res;
1312  }
1313 
1314  index_t rows() const { return 1; }
1315  index_t cols() const { return 1; }
1316 
1317  void parse(gsExprHelper<Scalar> &evList) const {
1318  _v.parse(evList);
1319 // evList.add(_u);
1320 // _u.data().flags |= NEED_GRAD;
1321  _u.parse(evList);
1322  }
1323 
1324  const gsFeSpace<Scalar> &rowVar() const { return _u.rowVar(); }
1325  const gsFeSpace<Scalar> &colVar() const { return _v.rowVar(); }
1326 
1327  void print(std::ostream &os) const {
1328  os << "(";
1329  _u.print(os);
1330  os << " % ";
1331  _v.print(os);
1332  os << ")";
1333  }
1334 };
1335 
1337 template<typename E1, typename E2>
1338 EIGEN_STRONG_INLINE
1339 frprod3_expr<E1, E2> const frprod3(E1 const &u,
1340  E2 const &M) {
1341  return frprod3_expr<E1, E2>(u, M);
1342 }
1343 
1344 template<class E0, class E1, class E2>
1345 class ternary_expr : public _expr<ternary_expr<E0, E1, E2> > {
1346  typename E0::Nested_t _u;
1347  typename E1::Nested_t _v;
1348  typename E2::Nested_t _w;
1349  public:
1350  typedef typename E1::Scalar Scalar;
1351 
1352  explicit ternary_expr(_expr<E0> const &u,
1353  _expr<E1> const &v,
1354  _expr<E2> const &w)
1355  :
1356  _u(u),
1357  _v(v),
1358  _w(w) {
1359  GISMO_ASSERT(E0::ScalarValued, "Condition must be scalar valued");
1360  GISMO_ASSERT((int) E1::ScalarValued == (int) E2::ScalarValued,
1361  "Both v and w must be scalar valued (or not).");
1362  GISMO_ASSERT((int) E1::ColBlocks == (int) E2::ColBlocks,
1363  "Both v and w must be colblocks (or not).");
1364  GISMO_ASSERT((int) E1::Space == (int) E2::Space,
1365  "Both v and w must be space (or not), but E1::Space = "
1366  << E1::Space << " and E2::Space = " << E2::Space);
1367  GISMO_ASSERT(_v.rows() == _w.rows(),
1368  "Rows of v and w differ. _v.rows() = " << _v.rows()
1369  << ", _w.rows() = "
1370  << _w.rows());
1371  GISMO_ASSERT(_v.cols() == _w.cols(),
1372  "Columns of v and w differ. _v.cols() = " << _v.cols()
1373  << ", _w.cols() = "
1374  << _w.cols());
1375  GISMO_ASSERT(_v.rowVar() == _w.rowVar(), "rowVar of v and w differ.");
1376  GISMO_ASSERT(_v.colVar() == _w.colVar(), "colVar of v and w differ.");
1377  }
1378  public:
1379  enum {
1380  ScalarValued = E1::ScalarValued,
1381  ColBlocks = E1::ColBlocks,
1382  Space = E1::Space
1383  }; // == E2::Space
1384 
1385 // const Scalar eval(const index_t k) const { return (_u.eval(k) > 0 ? _v.eval
1386 // (k) : _w.eval(k)); }
1387 
1388  const Temporary_t eval(const index_t k) const {
1389  return (_u.eval(k) > 0 ? _v.eval(k) : _w.eval(k));
1390  }
1391 
1392  // { res = eval_impl(_u,_v,_w,k); return res;}
1393 
1394  index_t rows() const { return _v.rows(); }
1395  index_t cols() const { return _v.cols(); }
1396  void parse(gsExprHelper<Scalar> &evList) const {
1397  _u.parse(evList);
1398  _v.parse(evList);
1399  _w.parse(evList);
1400  }
1401 
1402  const gsFeSpace<Scalar> &rowVar() const { return _v.rowVar(); }
1403  const gsFeSpace<Scalar> &colVar() const { return _v.colVar(); }
1404 
1405  void print(std::ostream &os) const { os << "ternary"; }
1406 
1408 // private:
1409 // typename util::enable_if<E1::ScalarValued,Scalar>::type
1410 // eval_impl(const E0 & u, const E1 & v, const E2 & w, const index_t k) const
1411 // {
1412 // return (u.eval(k) > 0 ? v.eval(k) : w.eval(k));
1413 // }
1414 
1415 // typename util::enable_if<!ScalarValued,gsMatrix<Scalar>>::type
1416 // eval_impl(const E0 & u, const E1 & v, const E2 & w, const index_t k)
1417 // {
1418 // return (u.eval(k).array() > 0 ? v.eval(k) : w.eval(k));
1419 // }
1420 
1421 // private:
1422 // template<class U, class V, class W> static inline
1423 // typename util::enable_if<U::ScalarValued && V::ScalarValued,AutoReturn_t>::type
1424 // eval_impl(const U &u, const V & v, const W & w, const index_t k)
1425 // {
1426 // gsMatrix<Scalar> res(1,1);
1427 // // bool test = u.eval(k) > 0;
1428 // res<<(u.eval(k) > 0 ? v.eval(k) : w.eval(k));
1429 // return res;
1430 // }
1431 
1432 // template<class U, class V, class W> static inline
1433 // typename util::enable_if<U::ScalarValued && !V::ScalarValued,AutoReturn_t>::type
1434 // eval_impl(const U &u, const V & v, const W & w, const index_t k)
1435 // {
1436 // return u.eval(k) > 0 ? v.eval(k) : w.eval(k);
1437 // }
1438 
1439 // template<class U, class V, class W> static inline
1440 // typename util::enable_if<!U::ScalarValued,gsMatrix<Scalar>>::type
1441 // eval_impl(const U &u, const V & v, const W & w, const index_t k)
1442 // {
1443 // GISMO_ERROR("Something went wrong");
1444 // }
1445 };
1446 
1447 /*
1448  Expression for Jacobian matrix for PDE-based parameterization construction
1449 */
1450 template<class E>
1451 class jacScaledLx_expr : public _expr<jacScaledLx_expr<E> > {
1452  public:
1453  typedef typename E::Scalar Scalar;
1454 
1455  private:
1456  typename E::Nested_t _u;
1457  typename gsGeometryMap<Scalar>::Nested_t _G;
1458 
1459  public:
1460  enum { Space = 3, ScalarValued = 0, ColBlocks = 0 };
1461 
1462  jacScaledLx_expr(const E &u, const gsGeometryMap<Scalar> &G) : _u(u), _G(G) {}
1463 
1464  mutable gsMatrix<Scalar> res, derivGeom, deriv2Geom, derivBasis, deriv2Basis;
1465  mutable gsMatrix<Scalar> dg11dx, dg11dy, dg22dx, dg22dy, dg12dx, dg12dy;
1466  mutable gsMatrix<Scalar> commonTerm, dLxdx, dLxdy, dLydx, dLydy;
1467 
1468 // EIGEN_MAKE_ALIGNED_OPERATOR_NEW
1469 
1470  const gsMatrix<Scalar> &eval(const index_t k) const {
1471  gsMatrix<Scalar> basis = _u.data().values[0].col(k);
1472 
1473  derivBasis = _u.data().values[1].col(k).transpose();
1474  deriv2Basis = _u.data().values[2].col(k).transpose();
1475 
1476  derivBasis.blockTransposeInPlace(_u.dim());
1477  deriv2Basis.blockTransposeInPlace(1 + _u.dim());
1478 
1479  derivGeom = _G.data().values[1].col(k);
1480  deriv2Geom = _G.data().values[2].col(k);
1481 
1482  Scalar g11 = derivGeom(0) * derivGeom(0) + derivGeom(2) * derivGeom(2);
1483  Scalar g12 = derivGeom(0) * derivGeom(1) + derivGeom(2) * derivGeom(3);
1484  Scalar g22 = derivGeom(1) * derivGeom(1) + derivGeom(3) * derivGeom(3);
1485 
1486  Scalar scaleFactor = g11 + g22;
1487 
1488  Scalar Lx =
1489  (g22 * deriv2Geom(0) + g11 * deriv2Geom(1) - 2.0 * g12 * deriv2Geom(2))
1490  / scaleFactor;
1491  Scalar Ly =
1492  (g22 * deriv2Geom(3) + g11 * deriv2Geom(4) - 2.0 * g12 * deriv2Geom(5))
1493  / scaleFactor;
1494 
1495  dg11dx.noalias() = 2.0 * derivGeom(0) * derivBasis.row(0);
1496  dg11dy.noalias() = 2.0 * derivGeom(2) * derivBasis.row(0);
1497  dg22dx.noalias() = 2.0 * derivGeom(1) * derivBasis.row(1);
1498  dg22dy.noalias() = 2.0 * derivGeom(3) * derivBasis.row(1);
1499  dg12dx.noalias() =
1500  derivGeom(1) * derivBasis.row(0) + derivGeom(0) * derivBasis.row(1);
1501  dg12dy.noalias() =
1502  derivGeom(3) * derivBasis.row(0) + derivGeom(2) * derivBasis.row(1);
1503 
1504  commonTerm.noalias() = g22 * deriv2Basis.row(0) + g11 * deriv2Basis.row(1)
1505  - 2.0 * g12 * deriv2Basis.row(2);
1506  dLxdx.noalias() = dg22dx * deriv2Geom(0) + dg11dx * deriv2Geom(1)
1507  - 2.0 * dg12dx * deriv2Geom(2) + commonTerm;
1508  dLxdy.noalias() = dg22dy * deriv2Geom(0) + dg11dy * deriv2Geom(1)
1509  - 2.0 * dg12dy * deriv2Geom(2);
1510  dLydx.noalias() = dg22dx * deriv2Geom(3) + dg11dx * deriv2Geom(4)
1511  - 2.0 * dg12dx * deriv2Geom(5);
1512  dLydy.noalias() = dg22dy * deriv2Geom(3) + dg11dy * deriv2Geom(4)
1513  - 2.0 * dg12dy * deriv2Geom(5) + commonTerm;
1514 
1515  dLxdx = (dLxdx - Lx * (dg11dx + dg22dx)) / scaleFactor;
1516  dLxdy = (dLxdy - Lx * (dg11dy + dg22dy)) / scaleFactor;
1517  dLydx = (dLydx - Ly * (dg11dx + dg22dx)) / scaleFactor;
1518  dLydy = (dLydy - Ly * (dg11dy + dg22dy)) / scaleFactor;
1519 
1520  const index_t A = _u.cardinality() / _u.dim(); // _u.data().actives.rows()
1521  res.resize(_u.cardinality(), _u.cardinality());
1522 // res.setZero();
1523  res.topLeftCorner(A, A).noalias() = basis * dLxdx;
1524  res.topRightCorner(A, A).noalias() = basis * dLxdy;
1525  res.bottomLeftCorner(A, A).noalias() = basis * dLydx;
1526  res.bottomRightCorner(A, A).noalias() = basis * dLydy;
1527 
1528  return res;
1529  }
1530 
1531  index_t rows() const { return 1; }
1532 
1533  index_t cols() const { return 1; }
1534 
1535  void parse(gsExprHelper<Scalar> &evList) const {
1536  evList.add(_u);
1537  _u.data().flags |= NEED_VALUE | NEED_GRAD | NEED_DERIV2;
1538 
1539  evList.add(_G);
1540  _G.data().flags |= NEED_DERIV | NEED_DERIV2;
1541  }
1542 
1543  const gsFeSpace<Scalar> &rowVar() const { return _u.rowVar(); }
1544  const gsFeSpace<Scalar> &colVar() const { return _u.rowVar(); }
1545  // TODO: question, what do these parameters mean?
1546  index_t cardinality_impl() const { return _u.cardinality_impl(); }
1547 
1548  void print(std::ostream &os) const {
1549  os << "jacScaledLx(";
1550  _u.print(os);
1551  os << ")";
1552  }
1553 };
1554 
1555 /*
1556  Expression for Jacobian matrix (diagonal part) for PDE-based parameterization construction
1557 */
1558 template<class E>
1559 class jacScaledLxDiag_expr : public _expr<jacScaledLxDiag_expr<E> > {
1560  public:
1561  typedef typename E::Scalar Scalar;
1562 
1563  private:
1564  typename E::Nested_t _u;
1565  typename gsGeometryMap<Scalar>::Nested_t _G;
1566 
1567  public:
1568  enum { Space = 3, ScalarValued = 0, ColBlocks = 0 };
1569 
1570  jacScaledLxDiag_expr(const E &u, const gsGeometryMap<Scalar> &G)
1571  : _u(u), _G(G) {}
1572 
1573  mutable gsMatrix<Scalar> res, derivGeom, deriv2Geom, derivBasis, deriv2Basis;
1574  mutable gsMatrix<Scalar> dg11dx, dg11dy, dg22dx, dg22dy, dg12dx, dg12dy;
1575  mutable gsMatrix<Scalar> commonTerm, dLxdx, dLydy;
1576 
1577 // EIGEN_MAKE_ALIGNED_OPERATOR_NEW
1578 
1579  const gsMatrix<Scalar> &eval(const index_t k) const {
1580  derivBasis = _u.data().values[1].col(k).transpose();
1581  deriv2Basis = _u.data().values[2].col(k).transpose();
1582 
1583  derivBasis.blockTransposeInPlace(_u.dim());
1584  deriv2Basis.blockTransposeInPlace(1 + _u.dim());
1585 
1586  derivGeom = _G.data().values[1].col(k);
1587  deriv2Geom = _G.data().values[2].col(k);
1588 
1589  Scalar g11 = derivGeom(0) * derivGeom(0) + derivGeom(2) * derivGeom(2);
1590  Scalar g12 = derivGeom(0) * derivGeom(1) + derivGeom(2) * derivGeom(3);
1591  Scalar g22 = derivGeom(1) * derivGeom(1) + derivGeom(3) * derivGeom(3);
1592 
1593  Scalar scaleFactor = g11 + g22;
1594 
1595  Scalar Lx =
1596  (g22 * deriv2Geom(0) + g11 * deriv2Geom(1) - 2.0 * g12 * deriv2Geom(2))
1597  / scaleFactor;
1598  Scalar Ly =
1599  (g22 * deriv2Geom(3) + g11 * deriv2Geom(4) - 2.0 * g12 * deriv2Geom(5))
1600  / scaleFactor;
1601 
1602  dg11dx.noalias() = 2.0 * derivGeom(0) * derivBasis.row(0);
1603  dg11dy.noalias() = 2.0 * derivGeom(2) * derivBasis.row(0);
1604  dg22dx.noalias() = 2.0 * derivGeom(1) * derivBasis.row(1);
1605  dg22dy.noalias() = 2.0 * derivGeom(3) * derivBasis.row(1);
1606  dg12dx.noalias() =
1607  derivGeom(1) * derivBasis.row(0) + derivGeom(0) * derivBasis.row(1);
1608  dg12dy.noalias() =
1609  derivGeom(3) * derivBasis.row(0) + derivGeom(2) * derivBasis.row(1);
1610 
1611  commonTerm.noalias() = g22 * deriv2Basis.row(0) + g11 * deriv2Basis.row(1)
1612  - 2.0 * g12 * deriv2Basis.row(2);
1613  dLxdx.noalias() = dg22dx * deriv2Geom(0) + dg11dx * deriv2Geom(1)
1614  - 2.0 * dg12dx * deriv2Geom(2) + commonTerm;
1615  dLydy.noalias() = dg22dy * deriv2Geom(3) + dg11dy * deriv2Geom(4)
1616  - 2.0 * dg12dy * deriv2Geom(5) + commonTerm;
1617 
1618  dLxdx = (dLxdx - Lx * (dg11dx + dg22dx)) / scaleFactor;
1619  dLydy = (dLydy - Ly * (dg11dy + dg22dy)) / scaleFactor;
1620 
1621  const index_t A = _u.cardinality() / _u.dim(); // _u.data().actives.rows()
1622  res.resize(_u.cardinality(), _u.cardinality());
1623  res.setZero();
1624  res.topLeftCorner(A, A) = (_u.data().values[0].col(k).array()
1625  * dLxdx.array()).matrix().asDiagonal();
1626  res.bottomRightCorner(A, A) = (_u.data().values[0].col(k).array()
1627  * dLydy.array()).matrix().asDiagonal();
1628 
1629  return res;
1630  }
1631 
1632  index_t rows() const { return 1; }
1633 
1634  index_t cols() const { return 1; }
1635 
1636  void parse(gsExprHelper<Scalar> &evList) const {
1637  evList.add(_u);
1638  _u.data().flags |= NEED_VALUE | NEED_GRAD | NEED_DERIV2;
1639 
1640  evList.add(_G);
1641  _G.data().flags |= NEED_DERIV | NEED_DERIV2;
1642  }
1643 
1644  const gsFeSpace<Scalar> &rowVar() const { return _u.rowVar(); }
1645  const gsFeSpace<Scalar> &colVar() const { return _u.rowVar(); }
1646  // TODO: question, what do these parameters mean?
1647  index_t cardinality_impl() const { return _u.cardinality_impl(); }
1648 
1649  void print(std::ostream &os) const {
1650  os << "jacScaledLxDiag(";
1651  _u.print(os);
1652  os << ")";
1653  }
1654 };
1655 
1656 /*
1657  Expression for Jacobian matrix (diagonal block) for PDE-based parameterization construction
1658 */
1659 template<class E>
1660 class jacScaledLxDiagBlock_expr : public _expr<jacScaledLxDiagBlock_expr<E> > {
1661  public:
1662  typedef typename E::Scalar Scalar;
1663 
1664  private:
1665  typename E::Nested_t _u;
1666  typename gsGeometryMap<Scalar>::Nested_t _G;
1667 
1668  public:
1669  enum { Space = 3, ScalarValued = 0, ColBlocks = 0 };
1670 
1671  jacScaledLxDiagBlock_expr(const E &u, const gsGeometryMap<Scalar> &G)
1672  : _u(u), _G(G) {}
1673 
1674  mutable gsMatrix<Scalar> res, basis, derivGeom, deriv2Geom, derivBasis,
1675  deriv2Basis;
1676  mutable gsMatrix<Scalar> dg11dx, dg11dy, dg22dx, dg22dy, dg12dx, dg12dy;
1677  mutable gsMatrix<Scalar> commonTerm, dLxdx, dLydy;
1678 
1679 // EIGEN_MAKE_ALIGNED_OPERATOR_NEW
1680 
1681  const gsMatrix<Scalar> &eval(const index_t k) const {
1682 // basis = _u.data().values[0].col(k);
1683 
1684  derivBasis = _u.data().values[1].col(k).transpose();
1685  deriv2Basis = _u.data().values[2].col(k).transpose();
1686 
1687  derivBasis.blockTransposeInPlace(_u.dim());
1688  deriv2Basis.blockTransposeInPlace(1 + _u.dim());
1689 
1690  derivGeom = _G.data().values[1].col(k);
1691  deriv2Geom = _G.data().values[2].col(k);
1692 
1693  Scalar g11 = derivGeom(0) * derivGeom(0) + derivGeom(2) * derivGeom(2);
1694  Scalar g12 = derivGeom(0) * derivGeom(1) + derivGeom(2) * derivGeom(3);
1695  Scalar g22 = derivGeom(1) * derivGeom(1) + derivGeom(3) * derivGeom(3);
1696  Scalar scaleFactor = g11 + g22;
1697 
1698  Scalar Lx =
1699  (g22 * deriv2Geom(0) + g11 * deriv2Geom(1) - 2.0 * g12 * deriv2Geom(2))
1700  / scaleFactor;
1701  Scalar Ly =
1702  (g22 * deriv2Geom(3) + g11 * deriv2Geom(4) - 2.0 * g12 * deriv2Geom(5))
1703  / scaleFactor;
1704 
1705  dg11dx.noalias() = 2.0 * derivGeom(0) * derivBasis.row(0);
1706  dg11dy.noalias() = 2.0 * derivGeom(2) * derivBasis.row(0);
1707  dg22dx.noalias() = 2.0 * derivGeom(1) * derivBasis.row(1);
1708  dg22dy.noalias() = 2.0 * derivGeom(3) * derivBasis.row(1);
1709  dg12dx.noalias() =
1710  derivGeom(1) * derivBasis.row(0) + derivGeom(0) * derivBasis.row(1);
1711  dg12dy.noalias() =
1712  derivGeom(3) * derivBasis.row(0) + derivGeom(2) * derivBasis.row(1);
1713 
1714  commonTerm.noalias() = g22 * deriv2Basis.row(0) + g11 * deriv2Basis.row(1)
1715  - 2.0 * g12 * deriv2Basis.row(2);
1716  dLxdx.noalias() = dg22dx * deriv2Geom(0) + dg11dx * deriv2Geom(1)
1717  - 2.0 * dg12dx * deriv2Geom(2) + commonTerm;
1718  dLydy.noalias() = dg22dy * deriv2Geom(3) + dg11dy * deriv2Geom(4)
1719  - 2.0 * dg12dy * deriv2Geom(5) + commonTerm;
1720 
1721  dLxdx = (dLxdx - Lx * (dg11dx + dg22dx)) / scaleFactor;
1722  dLydy = (dLydy - Ly * (dg11dy + dg22dy)) / scaleFactor;
1723 
1724  const index_t A = _u.cardinality() / _u.dim(); // _u.data().actives.rows()
1725  res.resize(_u.cardinality(), _u.cardinality());
1726  res.setZero();
1727  res.topLeftCorner(A, A).noalias() =
1728  _u.data().values[0].col(k) * dLxdx;
1729  res.bottomRightCorner(A, A).noalias() =
1730  _u.data().values[0].col(k) * dLydy;
1731 
1732  return res;
1733  }
1734 
1735  index_t rows() const { return 1; }
1736 
1737  index_t cols() const { return 1; }
1738 
1739  void parse(gsExprHelper<Scalar> &evList) const {
1740  evList.add(_u);
1741  _u.data().flags |= NEED_VALUE | NEED_GRAD | NEED_DERIV2;
1742 
1743  evList.add(_G);
1744  _G.data().flags |= NEED_DERIV | NEED_DERIV2;
1745  }
1746 
1747  const gsFeSpace<Scalar> &rowVar() const { return _u.rowVar(); }
1748  const gsFeSpace<Scalar> &colVar() const { return _u.rowVar(); }
1749  // TODO: question, what do these parameters mean?
1750  index_t cardinality_impl() const { return _u.cardinality_impl(); }
1751 
1752  void print(std::ostream &os) const {
1753  os << "jacScaledLxDiagBlock(";
1754  _u.print(os);
1755  os << ")";
1756  }
1757 };
1758 
1759 /*
1760  Expression for Jacobian matrix (in H1 space) for PDE-based parameterization construction
1761 */
1762 template<class E>
1763 class jacScaledLxH1_expr : public _expr<jacScaledLxH1_expr<E> > {
1764  public:
1765  typedef typename E::Scalar Scalar;
1766 
1767  private:
1768  typename E::Nested_t _u;
1769  typename gsGeometryMap<Scalar>::Nested_t _G;
1770 
1771  public:
1772  enum { Space = 3, ScalarValued = 0, ColBlocks = 0 };
1773 
1774  jacScaledLxH1_expr(const E &u, const gsGeometryMap<Scalar> &G)
1775  : _u(u), _G(G) {}
1776 
1777  mutable gsMatrix<Scalar> res, derivBasis;
1778 
1779 // EIGEN_MAKE_ALIGNED_OPERATOR_NEW
1780 
1781  const gsMatrix<Scalar> &eval(const index_t k) const {
1782  gsMatrix<Scalar> jacMat = _G.data().values[1].reshapeCol(k,
1783  _G.data().dim.first,
1784  _G.data().dim.second);
1785  gsMatrix<Scalar> invJacMat = jacMat.inverse();
1786 
1787  derivBasis = _u.data().values[1].col(k).transpose();
1788  derivBasis.blockTransposeInPlace(_u.dim());
1789  gsMatrix<Scalar> invJacHatG;
1790  invJacHatG.noalias() = invJacMat * derivBasis;
1791 
1792  const index_t N = _u.cardinality() / _u.dim(); // _u.data().actives.rows()
1793  gsMatrix<Scalar> jacdLxdx(N, N);
1794  gsMatrix<Scalar> jacdLxdy(N, N);
1795  gsMatrix<Scalar> jacdLydx(N, N);
1796  gsMatrix<Scalar> jacdLydy(N, N);
1797 
1798  gsMatrix<> temp(2, N);
1799  for (auto i = 0; i < N; ++i) {
1800  // for x-direction
1801  temp.row(0).noalias() = invJacHatG(0, i) * invJacHatG.row(0);
1802  temp.row(1).noalias() = invJacHatG(1, i) * invJacHatG.row(0);
1803  gsMatrix<Scalar> dinvJacdx(2, 2);
1804  dinvJacdx.row(0).noalias() = invJacHatG(0, i) * invJacMat.row(0);
1805  dinvJacdx.row(1).noalias() = invJacHatG(1, i) * invJacMat.row(0);
1806  jacdLxdx.col(i).noalias() = -temp.transpose() * invJacMat.col(0)
1807  - invJacHatG.transpose() * dinvJacdx.col(0);
1808  jacdLydx.col(i).noalias() = -temp.transpose() * invJacMat.col(1)
1809  - invJacHatG.transpose() * dinvJacdx.col(1);
1810 
1811  // for y-direction
1812  temp.row(0).noalias() = invJacHatG(0, i) * invJacHatG.row(1);
1813  temp.row(1).noalias() = invJacHatG(1, i) * invJacHatG.row(1);
1814  gsMatrix<Scalar> dinvJacdy(2, 2);
1815  dinvJacdy.row(0).noalias() = invJacHatG(0, i) * invJacMat.row(1);
1816  dinvJacdy.row(1).noalias() = invJacHatG(1, i) * invJacMat.row(1);
1817  jacdLxdy.col(i).noalias() = -temp.transpose() * invJacMat.col(0)
1818  - invJacHatG.transpose() * dinvJacdy.col(0);
1819  jacdLydy.col(i).noalias() = -temp.transpose() * invJacMat.col(1)
1820  - invJacHatG.transpose() * dinvJacdy.col(1);
1821  }
1822 
1823  res.resize(_u.cardinality(), _u.cardinality());
1824  res.setZero();
1825  res.topLeftCorner(N, N).noalias() = jacdLxdx;
1826  res.topRightCorner(N, N).noalias() = jacdLxdy;
1827  res.bottomLeftCorner(N, N).noalias() = jacdLydx;
1828  res.bottomRightCorner(N, N).noalias() = jacdLydy;
1829 
1830  return res;
1831  }
1832 
1833  index_t rows() const { return 1; }
1834 
1835  index_t cols() const { return 1; }
1836 
1837  void parse(gsExprHelper<Scalar> &evList) const {
1838  evList.add(_u);
1839  _u.data().flags |= NEED_GRAD;
1840 
1841  evList.add(_G);
1842  _G.data().flags |= NEED_DERIV;
1843  }
1844 
1845  const gsFeSpace<Scalar> &rowVar() const { return _u.rowVar(); }
1846  const gsFeSpace<Scalar> &colVar() const { return _u.rowVar(); }
1847  // TODO: question, what do these parameters mean?
1848  index_t cardinality_impl() const { return _u.cardinality_impl(); }
1849 
1850  void print(std::ostream &os) const {
1851  os << "jacScaledLx(";
1852  _u.print(os);
1853  os << ")";
1854  }
1855 };
1856 
1857 /*
1858  Expression for Jacobian matrix (in H1 space) for PDE-based parameterization construction
1859 */
1860 template<class E>
1861 class jacScaledLxH1DiagBlock_expr
1862  : public _expr<jacScaledLxH1DiagBlock_expr<E> > {
1863  public:
1864  typedef typename E::Scalar Scalar;
1865 
1866  private:
1867  typename E::Nested_t _u;
1868  typename gsGeometryMap<Scalar>::Nested_t _G;
1869 
1870  public:
1871  enum { Space = 3, ScalarValued = 0, ColBlocks = 0 };
1872 
1873  jacScaledLxH1DiagBlock_expr(const E &u, const gsGeometryMap<Scalar> &G)
1874  : _u(u), _G(G) {}
1875 
1876  mutable gsMatrix<Scalar> res, derivBasis;
1877 
1878 // EIGEN_MAKE_ALIGNED_OPERATOR_NEW
1879 
1880  const gsMatrix<Scalar> &eval(const index_t k) const {
1881  gsMatrix<Scalar> jacMat = _G.data().values[1].reshapeCol(k,
1882  _G.data().dim.first,
1883  _G.data().dim.second);
1884  gsMatrix<Scalar> invJacMat = jacMat.inverse();
1885 
1886  derivBasis = _u.data().values[1].col(k).transpose();
1887  derivBasis.blockTransposeInPlace(_u.dim());
1888  gsMatrix<Scalar> invJacHatG = invJacMat * derivBasis;
1889 
1890  const index_t N = _u.cardinality() / _u.dim(); // _u.data().actives.rows()
1891  gsMatrix<Scalar> jacdLxdx(N, N);
1892 // gsMatrix<Scalar> jacdLxdy(N,N);
1893 // gsMatrix<Scalar> jacdLydx(N,N);
1894  gsMatrix<Scalar> jacdLydy(N, N);
1895 
1896  gsMatrix<> temp(2, N);
1897  for (auto i = 0; i < N; ++i) {
1898  // for x-direction
1899  temp.row(0) = invJacHatG(0, i) * invJacHatG.row(0);
1900  temp.row(1) = invJacHatG(1, i) * invJacHatG.row(0);
1901  gsMatrix<Scalar> dinvJacdx(2, 2);
1902  dinvJacdx.row(0) = invJacHatG(0, i) * invJacMat.row(0);
1903  dinvJacdx.row(1) = invJacHatG(1, i) * invJacMat.row(0);
1904  jacdLxdx.col(i) = -temp.transpose() * invJacMat.col(0)
1905  - invJacHatG.transpose() * dinvJacdx.col(0);
1906 // jacdLydx.col(i) = -temp.transpose()*invJacMat.col(1) - invJacHatG.transpose()*dinvJacdx.col(1);
1907 
1908  // for y-direction
1909  temp.row(0) = invJacHatG(0, i) * invJacHatG.row(1);
1910  temp.row(1) = invJacHatG(1, i) * invJacHatG.row(1);
1911  gsMatrix<Scalar> dinvJacdy(2, 2);
1912  dinvJacdy.row(0) = invJacHatG(0, i) * invJacMat.row(1);
1913  dinvJacdy.row(1) = invJacHatG(1, i) * invJacMat.row(1);
1914 // jacdLxdy.col(i) = -temp.transpose()*invJacMat.col(0) - invJacHatG.transpose()*dinvJacdy.col(0);
1915  jacdLydy.col(i) = -temp.transpose() * invJacMat.col(1)
1916  - invJacHatG.transpose() * dinvJacdy.col(1);
1917  }
1918 
1919  res.resize(_u.cardinality(), _u.cardinality());
1920  res.setZero();
1921  res.topLeftCorner(N, N) = jacdLxdx;
1922 // res.topRightCorner(N,N) = jacdLxdy;
1923 // res.bottomLeftCorner(N,N) = jacdLydx;
1924  res.bottomRightCorner(N, N) = jacdLydy;
1925 
1926  return res;
1927  }
1928 
1929  index_t rows() const { return 1; }
1930 
1931  index_t cols() const { return 1; }
1932 
1933  void parse(gsExprHelper<Scalar> &evList) const {
1934  evList.add(_u);
1935  _u.data().flags |= NEED_GRAD;
1936 
1937  evList.add(_G);
1938  _G.data().flags |= NEED_DERIV;
1939  }
1940 
1941  const gsFeSpace<Scalar> &rowVar() const { return _u.rowVar(); }
1942  const gsFeSpace<Scalar> &colVar() const { return _u.rowVar(); }
1943  // TODO: question, what do these parameters mean?
1944  index_t cardinality_impl() const { return _u.cardinality_impl(); }
1945 
1946  void print(std::ostream &os) const {
1947  os << "jacScaledLx(";
1948  _u.print(os);
1949  os << ")";
1950  }
1951 };
1952 
1953 /*
1954 
1955  Expression for the product of the jacob of space and matrix A. Return term
1956  is a column vector (d*n x 1). The result should be the same as the result
1957  from jac(space) % A, whereas avoids redundant multiply with zero-value components.
1958 
1959  For 2D case:
1960  frprod2(space,A) = [a_{11}*dN_1/dxi_1, a_{12}*dN_1/dxi_2, a_{11}*dN_2/dxi_1,
1961  a_{12}*dN_2/dxi_2,...,a_{11}*dN_n/dxi_1, a_{12}*dN_n/dxi_2, a_{21}*dN_1/dxi_1,
1962  a_{22}*dN_1/dxi_2, a_{21}*dN_2/dxi_1, a_{22}*dN_2/dxi_2,...,a_{21}*dN_n/dxi_1,
1963  a_{22}*dN_n/dxi_2]
1964  For 3D case:
1965  frprod2(space,A) = [a_{11}*dN_1/dxi_1, a_{12}*dN_1/dxi_2, a_{13}*dN_1/dxi_3,...
1966  a_{11}*dN_2/dxi_1, a_{12}*dN_2/dxi_1, a_{13}*dN_2/dxi_3,...,
1967  a_{11}*dN_n/dxi_1, a_{12}*dN_n/dxi_1, a_{13}*dN_n/dxi_3,
1968  a_{21}*dN_1/dxi_1, a_{22}*dN_1/dxi_2, a_{23}*dN_1/dxi_3,...
1969  a_{21}*dN_2/dxi_1, a_{22}*dN_2/dxi_1, a_{23}*dN_2/dxi_3,...,
1970  a_{21}*dN_n/dxi_1, a_{22}*dN_n/dxi_1, a_{23}*dN_n/dxi_3,
1971  a_{31}*dN_1/dxi_1, a_{32}*dN_1/dxi_2, a_{33}*dN_1/dxi_3,...
1972  a_{31}*dN_2/dxi_1, a_{32}*dN_2/dxi_1, a_{33}*dN_2/dxi_3,...,
1973  a_{31}*dN_n/dxi_1, a_{32}*dN_n/dxi_1, a_{33}*dN_n/dxi_3,]
1974 
1975 NOTE: _u should be a space, _v should NOT be a space (fix with assert)
1976 
1977  */
1978 template<typename E1, typename E2>
1979 class frprod2_expr : public _expr<frprod2_expr<E1, E2> > {
1980  public:
1981  typedef typename E2::Scalar Scalar;
1982  enum { ScalarValued = 0, Space = E1::Space, ColBlocks = 0 };
1983 
1984  private:
1985  typename E1::Nested_t _u;
1986  typename E2::Nested_t _v;
1987 
1988  mutable gsMatrix<Scalar> res, bGrads;
1989 
1990  public:
1991 
1992  frprod2_expr(_expr<E1> const &u, _expr<E2> const &v)
1993  : _u(u), _v(v) {
1994  // gsInfo << "expression is space ? "<<E1::Space <<"\n"; _u.print(gsInfo);
1995  // GISMO_ASSERT(_u.rows() == _v.rows(),
1996  // "Wrong dimensions "<<_u.rows()<<"!="<<_v.rows()<<" in % operation");
1997  // GISMO_ASSERT(_u.cols() == _v.cols(),
1998  // "Wrong dimensions "<<_u.cols()<<"!="<<_v.cols()<<" in % operation");
1999  }
2000 
2001  const gsMatrix<Scalar> &eval(const index_t k) const //todo: specialize for nb==1
2002  {
2003  auto A = _v.eval(k);
2004  bGrads = _u.data().values[1].col(k).transpose();
2005  bGrads.blockTransposeInPlace(_u.dim());
2006  res.noalias() = A * bGrads;
2007  res.transposeInPlace();
2008  res.resize(1, _u.cardinality());
2009 // res.reshaped(1,_u.cardinality());
2010  res.transposeInPlace();
2011  return res;
2012  }
2013 
2014  index_t rows() const { return 1; }
2015  index_t cols() const { return 1; }
2016 
2017  void parse(gsExprHelper<Scalar> &evList) const {
2018  _v.parse(evList);
2019  evList.add(_u);
2020  _u.data().flags |= NEED_GRAD;
2021  }
2022 
2023  const gsFeSpace<Scalar> &rowVar() const { return _u.rowVar(); }
2024  const gsFeSpace<Scalar> &colVar() const { return _v.rowVar(); }
2025 
2026  void print(std::ostream &os) const {
2027  os << "(";
2028  _u.print(os);
2029  os << " % ";
2030  _v.print(os);
2031  os << ")";
2032  }
2033 };
2034 
2036 template<class E>
2037 EIGEN_STRONG_INLINE
2038 jacScaledLx_expr<E> jacScaledLx(const E &u,
2039  const gsGeometryMap<typename E::Scalar> &G) {
2040  return jacScaledLx_expr<E>(u, G);
2041 }
2042 
2044 template<class E>
2045 EIGEN_STRONG_INLINE
2046 jacScaledLxDiag_expr<E> jacScaledLxDiag(const E &u,
2047  const gsGeometryMap<typename E::Scalar> &G) {
2048  return jacScaledLxDiag_expr<E>(u, G);
2049 }
2050 
2052 template<class E>
2053 EIGEN_STRONG_INLINE
2054 jacScaledLxDiagBlock_expr<E> jacScaledLxDiagBlock(const E &u,
2055  const gsGeometryMap<typename E::Scalar> &G) {
2056  return jacScaledLxDiagBlock_expr<E>(u, G);
2057 }
2058 
2060 template<class E>
2061 EIGEN_STRONG_INLINE
2062 jacScaledLxH1_expr<E> jacScaledLxH1(const E &u,
2063  const gsGeometryMap<typename E::Scalar> &G) {
2064  return jacScaledLxH1_expr<E>(u, G);
2065 }
2066 
2068 template<class E>
2069 EIGEN_STRONG_INLINE
2070 jacScaledLxH1DiagBlock_expr<E> jacScaledLxH1DiagBlock(const E &u,
2071  const gsGeometryMap<
2072  typename E::Scalar> &G) {
2073  return jacScaledLxH1DiagBlock_expr<E>(u, G);
2074 }
2075 
2076 // Frobenious product (also known as double dot product) operator for expressions
2077 template<typename E1, typename E2>
2078 EIGEN_STRONG_INLINE
2079 frprod2_expr<E1, E2> const frprod2(E1 const &u,
2080  E2 const &M) {
2081  return frprod2_expr<E1, E2>(u, M);
2082 }
2083 
2085 template<class E0, class E1, class E2>
2086 EIGEN_STRONG_INLINE
2087 ternary_expr<E0, E1, E2> ternary(const E0 &u,
2088  const E1 &v,
2089  const E2 &w) {
2090  return ternary_expr<E0, E1, E2>(u, v, w);
2091 }
2092 } // namespace expr
2093 
2094 }// namespace gismo
static T computeAreaBoundary(const gsMultiPatch< T > &mp)
Computes the area of a multipatch representing boundary curves.
Definition: gsBarrierCore.hpp:114
EIGEN_STRONG_INLINE ternary_expr< E0, E1, E2 > ternary(const E0 &u, const E1 &v, const E2 &w)
Ternary ternary_expr.
Definition: gsBarrierCore.hpp:2087
T integral(const expr::_expr< E > &expr)
Calculates the integral of the expression expr on the whole integration domain.
Definition: gsExprEvaluator.h:154
Definition: gsExprAssembler.h:30
Scalar epsilon
Definition: gsBarrierCore.h:521
#define GISMO_NO_IMPLEMENTATION
Definition: gsDebug.h:129
Dirichlet type.
Definition: gsBoundaryConditions.h:31
Definition: gsBarrierCore.h:504
Gradient of the object.
Definition: gsForwardDeclarations.h:73
void setIntegrationElements(const gsMultiBasis< T > &mesh)
Sets the domain of integration.
Definition: gsExprEvaluator.h:110
EIGEN_STRONG_INLINE frprod3_expr< E1, E2 > const frprod3(E1 const &u, E2 const &M)
Frobenious product (also known as double dot product) operator for expressions.
Definition: gsBarrierCore.hpp:1339
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
EIGEN_STRONG_INLINE jacScaledLxDiagBlock_expr< E > jacScaledLxDiagBlock(const E &u, const gsGeometryMap< typename E::Scalar > &G)
diagonal block part of jacobian matrix of scaled Lx for PDE-based parameterization construction ...
Definition: gsBarrierCore.hpp:2054
const_biterator bBegin() const
Definition: gsBoxTopology.h:139
const gsVector< Scalar > & compute(const gsVector< Scalar > &u0, const Residual_t &F, const Jacobian_t &Jacobian)
perform Anderson acceleration iteration
Definition: gsBarrierCore.h:618
Anderson acceleration solver and its (preconditioned) variants.
Definition: gsBarrierCore.h:607
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
Second derivatives.
Definition: gsForwardDeclarations.h:80
#define index_t
Definition: gsConfig.h:32
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
Maintains a mapping from patch-local dofs to global dof indices and allows the elimination of individ...
Definition: gsDofMapper.h:68
void verboseLog(const std::string &message, const index_t &verbose)
helper function to set optimizer options
Definition: gsBarrierCore.h:101
static gsMultiPatch< T > computeVHPatch(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
variational harmonic method
const gsMatrix< T > & rhs() const
Returns the right-hand side vector(s)
Definition: gsExprAssembler.h:129
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
const gsSparseMatrix< T > & matrix() const
Returns the left-hand global matrix.
Definition: gsExprAssembler.h:116
std::string to_string(const unsigned &i)
Helper to convert small unsigned to string.
Definition: gsXml.cpp:74
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
static T computeAreaInterior(const gsMultiPatch< T > &mp)
Computes the area of a multipatch representing a full domain.
Definition: gsBarrierCore.hpp:94
#define gsWarn
Definition: gsDebug.h:50
Holds a set of patch-wise bases and their topology information.
Definition: gsMultiBasis.h:36
static gsMultiPatch< T > computeBarrierPatch(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
Barrier function-based method.
static gsMultiPatch< T > compute(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
construct analysis-suitable parameterization
Definition: gsBarrierCore.hpp:129
int updatePreconditionerStep
Definition: gsBarrierCore.h:542
Definition: gsDirichletValues.h:23
geometryMap getMap(const gsFunctionSet< T > &g)
Registers g as an isogeometric geometry map and return a handle to it.
Definition: gsExprAssembler.h:161
Container class for a set of geometry patches and their topology, that is, the interface connections ...
Definition: gsMultiPatch.h:33
Generic evaluator of isogeometric expressions.
Definition: gsExprEvaluator.h:38
space getSpace(const gsFunctionSet< T > &mp, index_t dim=1, index_t id=0)
Definition: gsExprAssembler.h:166
static gsOptionList defaultOptions()
Default options.
Definition: gsBarrierCore.hpp:24
EIGEN_STRONG_INLINE jacScaledLx_expr< E > jacScaledLx(const E &u, const gsGeometryMap< typename E::Scalar > &G)
jacobian matrix of scaled Lx for PDE-based parameterization construction
Definition: gsBarrierCore.hpp:2038
Class containing a set of boundary conditions.
Definition: gsBoundaryConditions.h:341
Value of the object.
Definition: gsForwardDeclarations.h:72
index_t askInt(const std::string &label, const index_t &value=0) const
Reads value for option label from options.
Definition: gsOptionList.cpp:117
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
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
EIGEN_STRONG_INLINE jacScaledLxH1_expr< E > jacScaledLxH1(const E &u, const gsGeometryMap< typename E::Scalar > &G)
jacobian matrix of scaled Lx (in H1 space) for PDE-based parameterization construction ...
Definition: gsBarrierCore.hpp:2062
#define GISMO_ERROR(message)
Definition: gsDebug.h:118
static gsMultiPatch< T > computePenaltyPatch2(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
Penalty function-based method (2)
bool askSwitch(const std::string &label, const bool &value=false) const
Reads value for option label from options.
Definition: gsOptionList.cpp:128
void setGeoMap(const gsFunctionSet< T > &gm)
Set the geometry map to evaluate boundary conditions.
Definition: gsBoundaryConditions.h:916
Class which holds a list of parameters/options, and provides easy access to them. ...
Definition: gsOptionList.h:32
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
void setIntegrationElements(const gsMultiBasis< T > &mesh)
Sets the domain of integration.
Definition: gsExprAssembler.h:136
void derivBasis()
Definition: gsBSplineAlgorithms.h:199
const_biterator bEnd() const
Definition: gsBoxTopology.h:144
void addSwitch(const std::string &label, const std::string &desc, const bool &value)
Adds a option named label, with description desc and value value.
Definition: gsOptionList.cpp:235
EIGEN_STRONG_INLINE jacScaledLxDiag_expr< E > jacScaledLxDiag(const E &u, const gsGeometryMap< typename E::Scalar > &G)
diagonal part of jacobian matrix of scaled Lx for PDE-based parameterization construction ...
Definition: gsBarrierCore.hpp:2046
bool usePreconditioning
Definition: gsBarrierCore.h:552
static T computeArea(const gsMultiPatch< T > &mp)
Compute the area of the computational domain.
Definition: gsBarrierCore.hpp:88
void initSystem(const index_t numRhs=1)
Initializes the sparse system (sparse matrix and rhs)
Definition: gsExprAssembler.h:290
geometryMap getMap(const gsMultiPatch< T > &mp)
Registers mp as an isogeometric geometry map and return a handle to it.
Definition: gsExprEvaluator.h:116
int m
Definition: gsBarrierCore.h:513
static gsMultiPatch< T > computePenaltyPatch(const gsMultiPatch< T > &mp, const gsDofMapper &mapper, const gsOptionList &options)
Penalty function-based method (1)
EIGEN_STRONG_INLINE jacScaledLxH1DiagBlock_expr< E > jacScaledLxH1DiagBlock(const E &u, const gsGeometryMap< typename E::Scalar > &G)
jacobian matrix of scaled Lx (in H1 space) for PDE-based parameterization construction ...
Definition: gsBarrierCore.hpp:2070
void copy(T begin, T end, U *result)
Small wrapper for std::copy mimicking std::copy for a raw pointer destination, copies n positions sta...
Definition: gsMemory.h:391
EIGEN_STRONG_INLINE jac_expr< E > jac(const symbol_expr< E > &u)
The Jacobian matrix of a FE variable.
Definition: gsExpressions.h:4528