Actual source code: plexmetric.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <petscblaslapack.h>

  4: PetscLogEvent DMPLEX_MetricEnforceSPD, DMPLEX_MetricNormalize, DMPLEX_MetricAverage, DMPLEX_MetricIntersection;

  6: PetscErrorCode DMPlexMetricSetFromOptions(DM dm)
  7: {
  8:   MPI_Comm       comm;
  9:   PetscBool      isotropic = PETSC_FALSE, uniform = PETSC_FALSE, restrictAnisotropyFirst = PETSC_FALSE;
 10:   PetscBool      noInsert = PETSC_FALSE, noSwap = PETSC_FALSE, noMove = PETSC_FALSE, noSurf = PETSC_FALSE;
 12:   PetscInt       verbosity = -1, numIter = 3;
 13:   PetscReal      h_min = 1.0e-30, h_max = 1.0e+30, a_max = 1.0e+05, p = 1.0, target = 1000.0, beta = 1.3, hausd = 0.01;

 15:   PetscObjectGetComm((PetscObject) dm, &comm);
 16:   PetscOptionsBegin(comm, "", "Riemannian metric options", "DMPlexMetric");
 17:   PetscOptionsBool("-dm_plex_metric_isotropic", "Is the metric isotropic?", "DMPlexMetricCreateIsotropic", isotropic, &isotropic, NULL);
 18:   DMPlexMetricSetIsotropic(dm, isotropic);
 19:   PetscOptionsBool("-dm_plex_metric_uniform", "Is the metric uniform?", "DMPlexMetricCreateUniform", uniform, &uniform, NULL);
 20:   DMPlexMetricSetUniform(dm, uniform);
 21:   PetscOptionsBool("-dm_plex_metric_restrict_anisotropy_first", "Should anisotropy be restricted before normalization?", "DMPlexNormalize", restrictAnisotropyFirst, &restrictAnisotropyFirst, NULL);
 22:   DMPlexMetricSetRestrictAnisotropyFirst(dm, restrictAnisotropyFirst);
 23:   PetscOptionsBool("-dm_plex_metric_no_insert", "Turn off node insertion and deletion", "DMAdaptMetric", noInsert, &noInsert, NULL);
 24:   DMPlexMetricSetNoInsertion(dm, noInsert);
 25:   PetscOptionsBool("-dm_plex_metric_no_swap", "Turn off facet swapping", "DMAdaptMetric", noSwap, &noSwap, NULL);
 26:   DMPlexMetricSetNoSwapping(dm, noSwap);
 27:   PetscOptionsBool("-dm_plex_metric_no_move", "Turn off facet node movement", "DMAdaptMetric", noMove, &noMove, NULL);
 28:   DMPlexMetricSetNoMovement(dm, noMove);
 29:   PetscOptionsBool("-dm_plex_metric_no_surf", "Turn off surface modification", "DMAdaptMetric", noSurf, &noSurf, NULL);
 30:   DMPlexMetricSetNoSurf(dm, noSurf);
 31:   PetscOptionsBoundedInt("-dm_plex_metric_num_iterations", "Number of ParMmg adaptation iterations", "DMAdaptMetric", numIter, &numIter, NULL, 0);
 32:   DMPlexMetricSetNumIterations(dm, numIter);
 33:   PetscOptionsRangeInt("-dm_plex_metric_verbosity", "Verbosity of metric-based mesh adaptation package (-1 = silent, 10 = maximum)", "DMAdaptMetric", verbosity, &verbosity, NULL, -1, 10);
 34:   DMPlexMetricSetVerbosity(dm, verbosity);
 35:   PetscOptionsReal("-dm_plex_metric_h_min", "Minimum tolerated metric magnitude", "DMPlexMetricEnforceSPD", h_min, &h_min, NULL);
 36:   DMPlexMetricSetMinimumMagnitude(dm, h_min);
 37:   PetscOptionsReal("-dm_plex_metric_h_max", "Maximum tolerated metric magnitude", "DMPlexMetricEnforceSPD", h_max, &h_max, NULL);
 38:   DMPlexMetricSetMaximumMagnitude(dm, h_max);
 39:   PetscOptionsReal("-dm_plex_metric_a_max", "Maximum tolerated anisotropy", "DMPlexMetricEnforceSPD", a_max, &a_max, NULL);
 40:   DMPlexMetricSetMaximumAnisotropy(dm, a_max);
 41:   PetscOptionsReal("-dm_plex_metric_p", "L-p normalization order", "DMPlexMetricNormalize", p, &p, NULL);
 42:   DMPlexMetricSetNormalizationOrder(dm, p);
 43:   PetscOptionsReal("-dm_plex_metric_target_complexity", "Target metric complexity", "DMPlexMetricNormalize", target, &target, NULL);
 44:   DMPlexMetricSetTargetComplexity(dm, target);
 45:   PetscOptionsReal("-dm_plex_metric_gradation_factor", "Metric gradation factor", "DMAdaptMetric", beta, &beta, NULL);
 46:   DMPlexMetricSetGradationFactor(dm, beta);
 47:   PetscOptionsReal("-dm_plex_metric_hausdorff_number", "Metric Hausdorff number", "DMAdaptMetric", hausd, &hausd, NULL);
 48:   DMPlexMetricSetHausdorffNumber(dm, hausd);
 49:   PetscOptionsEnd();
 50:   return 0;
 51: }

 53: /*@
 54:   DMPlexMetricSetIsotropic - Record whether a metric is isotropic

 56:   Input parameters:
 57: + dm        - The DM
 58: - isotropic - Is the metric isotropic?

 60:   Level: beginner

 62: .seealso: DMPlexMetricIsIsotropic(), DMPlexMetricSetUniform(), DMPlexMetricSetRestrictAnisotropyFirst()
 63: @*/
 64: PetscErrorCode DMPlexMetricSetIsotropic(DM dm, PetscBool isotropic)
 65: {
 66:   DM_Plex       *plex = (DM_Plex *) dm->data;

 68:   if (!plex->metricCtx) {
 69:     PetscNew(&plex->metricCtx);
 70:     DMPlexMetricSetFromOptions(dm);
 71:   }
 72:   plex->metricCtx->isotropic = isotropic;
 73:   return 0;
 74: }

 76: /*@
 77:   DMPlexMetricIsIsotropic - Is a metric isotropic?

 79:   Input parameters:
 80: . dm        - The DM

 82:   Output parameters:
 83: . isotropic - Is the metric isotropic?

 85:   Level: beginner

 87: .seealso: DMPlexMetricSetIsotropic(), DMPlexMetricIsUniform(), DMPlexMetricRestrictAnisotropyFirst()
 88: @*/
 89: PetscErrorCode DMPlexMetricIsIsotropic(DM dm, PetscBool *isotropic)
 90: {
 91:   DM_Plex       *plex = (DM_Plex *) dm->data;

 93:   if (!plex->metricCtx) {
 94:     PetscNew(&plex->metricCtx);
 95:     DMPlexMetricSetFromOptions(dm);
 96:   }
 97:   *isotropic = plex->metricCtx->isotropic;
 98:   return 0;
 99: }

101: /*@
102:   DMPlexMetricSetUniform - Record whether a metric is uniform

104:   Input parameters:
105: + dm      - The DM
106: - uniform - Is the metric uniform?

108:   Level: beginner

110:   Notes:

112:   If the metric is specified as uniform then it is assumed to be isotropic, too.

114: .seealso: DMPlexMetricIsUniform(), DMPlexMetricSetIsotropic(), DMPlexMetricSetRestrictAnisotropyFirst()
115: @*/
116: PetscErrorCode DMPlexMetricSetUniform(DM dm, PetscBool uniform)
117: {
118:   DM_Plex       *plex = (DM_Plex *) dm->data;

120:   if (!plex->metricCtx) {
121:     PetscNew(&plex->metricCtx);
122:     DMPlexMetricSetFromOptions(dm);
123:   }
124:   plex->metricCtx->uniform = uniform;
125:   if (uniform) plex->metricCtx->isotropic = uniform;
126:   return 0;
127: }

129: /*@
130:   DMPlexMetricIsUniform - Is a metric uniform?

132:   Input parameters:
133: . dm      - The DM

135:   Output parameters:
136: . uniform - Is the metric uniform?

138:   Level: beginner

140: .seealso: DMPlexMetricSetUniform(), DMPlexMetricIsIsotropic(), DMPlexMetricRestrictAnisotropyFirst()
141: @*/
142: PetscErrorCode DMPlexMetricIsUniform(DM dm, PetscBool *uniform)
143: {
144:   DM_Plex       *plex = (DM_Plex *) dm->data;

146:   if (!plex->metricCtx) {
147:     PetscNew(&plex->metricCtx);
148:     DMPlexMetricSetFromOptions(dm);
149:   }
150:   *uniform = plex->metricCtx->uniform;
151:   return 0;
152: }

154: /*@
155:   DMPlexMetricSetRestrictAnisotropyFirst - Record whether anisotropy should be restricted before normalization

157:   Input parameters:
158: + dm                      - The DM
159: - restrictAnisotropyFirst - Should anisotropy be normalized first?

161:   Level: beginner

163: .seealso: DMPlexMetricSetIsotropic(), DMPlexMetricRestrictAnisotropyFirst()
164: @*/
165: PetscErrorCode DMPlexMetricSetRestrictAnisotropyFirst(DM dm, PetscBool restrictAnisotropyFirst)
166: {
167:   DM_Plex       *plex = (DM_Plex *) dm->data;

169:   if (!plex->metricCtx) {
170:     PetscNew(&plex->metricCtx);
171:     DMPlexMetricSetFromOptions(dm);
172:   }
173:   plex->metricCtx->restrictAnisotropyFirst = restrictAnisotropyFirst;
174:   return 0;
175: }

177: /*@
178:   DMPlexMetricRestrictAnisotropyFirst - Is anisotropy restricted before normalization or after?

180:   Input parameters:
181: . dm                      - The DM

183:   Output parameters:
184: . restrictAnisotropyFirst - Is anisotropy be normalized first?

186:   Level: beginner

188: .seealso: DMPlexMetricIsIsotropic(), DMPlexMetricSetRestrictAnisotropyFirst()
189: @*/
190: PetscErrorCode DMPlexMetricRestrictAnisotropyFirst(DM dm, PetscBool *restrictAnisotropyFirst)
191: {
192:   DM_Plex       *plex = (DM_Plex *) dm->data;

194:   if (!plex->metricCtx) {
195:     PetscNew(&plex->metricCtx);
196:     DMPlexMetricSetFromOptions(dm);
197:   }
198:   *restrictAnisotropyFirst = plex->metricCtx->restrictAnisotropyFirst;
199:   return 0;
200: }

202: /*@
203:   DMPlexMetricSetNoInsertion - Should node insertion and deletion be turned off?

205:   Input parameters:
206: + dm       - The DM
207: - noInsert - Should node insertion and deletion be turned off?

209:   Level: beginner

211:   Notes:
212:   This is only used by Mmg and ParMmg (not Pragmatic).

214: .seealso: DMPlexMetricNoInsertion(), DMPlexMetricSetNoSwapping(), DMPlexMetricSetNoMovement(), DMPlexMetricSetNoSurf()
215: @*/
216: PetscErrorCode DMPlexMetricSetNoInsertion(DM dm, PetscBool noInsert)
217: {
218:   DM_Plex       *plex = (DM_Plex *) dm->data;

220:   if (!plex->metricCtx) {
221:     PetscNew(&plex->metricCtx);
222:     DMPlexMetricSetFromOptions(dm);
223:   }
224:   plex->metricCtx->noInsert = noInsert;
225:   return 0;
226: }

228: /*@
229:   DMPlexMetricNoInsertion - Are node insertion and deletion turned off?

231:   Input parameters:
232: . dm       - The DM

234:   Output parameters:
235: . noInsert - Are node insertion and deletion turned off?

237:   Level: beginner

239:   Notes:
240:   This is only used by Mmg and ParMmg (not Pragmatic).

242: .seealso: DMPlexMetricSetNoInsertion(), DMPlexMetricNoSwapping(), DMPlexMetricNoMovement(), DMPlexMetricNoSurf()
243: @*/
244: PetscErrorCode DMPlexMetricNoInsertion(DM dm, PetscBool *noInsert)
245: {
246:   DM_Plex       *plex = (DM_Plex *) dm->data;

248:   if (!plex->metricCtx) {
249:     PetscNew(&plex->metricCtx);
250:     DMPlexMetricSetFromOptions(dm);
251:   }
252:   *noInsert = plex->metricCtx->noInsert;
253:   return 0;
254: }

256: /*@
257:   DMPlexMetricSetNoSwapping - Should facet swapping be turned off?

259:   Input parameters:
260: + dm     - The DM
261: - noSwap - Should facet swapping be turned off?

263:   Level: beginner

265:   Notes:
266:   This is only used by Mmg and ParMmg (not Pragmatic).

268: .seealso: DMPlexMetricNoSwapping(), DMPlexMetricSetNoInsertion(), DMPlexMetricSetNoMovement(), DMPlexMetricSetNoSurf()
269: @*/
270: PetscErrorCode DMPlexMetricSetNoSwapping(DM dm, PetscBool noSwap)
271: {
272:   DM_Plex       *plex = (DM_Plex *) dm->data;

274:   if (!plex->metricCtx) {
275:     PetscNew(&plex->metricCtx);
276:     DMPlexMetricSetFromOptions(dm);
277:   }
278:   plex->metricCtx->noSwap = noSwap;
279:   return 0;
280: }

282: /*@
283:   DMPlexMetricNoSwapping - Is facet swapping turned off?

285:   Input parameters:
286: . dm     - The DM

288:   Output parameters:
289: . noSwap - Is facet swapping turned off?

291:   Level: beginner

293:   Notes:
294:   This is only used by Mmg and ParMmg (not Pragmatic).

296: .seealso: DMPlexMetricSetNoSwapping(), DMPlexMetricNoInsertion(), DMPlexMetricNoMovement(), DMPlexMetricNoSurf()
297: @*/
298: PetscErrorCode DMPlexMetricNoSwapping(DM dm, PetscBool *noSwap)
299: {
300:   DM_Plex       *plex = (DM_Plex *) dm->data;

302:   if (!plex->metricCtx) {
303:     PetscNew(&plex->metricCtx);
304:     DMPlexMetricSetFromOptions(dm);
305:   }
306:   *noSwap = plex->metricCtx->noSwap;
307:   return 0;
308: }

310: /*@
311:   DMPlexMetricSetNoMovement - Should node movement be turned off?

313:   Input parameters:
314: + dm     - The DM
315: - noMove - Should node movement be turned off?

317:   Level: beginner

319:   Notes:
320:   This is only used by Mmg and ParMmg (not Pragmatic).

322: .seealso: DMPlexMetricNoMovement(), DMPlexMetricSetNoInsertion(), DMPlexMetricSetNoSwapping(), DMPlexMetricSetNoSurf()
323: @*/
324: PetscErrorCode DMPlexMetricSetNoMovement(DM dm, PetscBool noMove)
325: {
326:   DM_Plex       *plex = (DM_Plex *) dm->data;

328:   if (!plex->metricCtx) {
329:     PetscNew(&plex->metricCtx);
330:     DMPlexMetricSetFromOptions(dm);
331:   }
332:   plex->metricCtx->noMove = noMove;
333:   return 0;
334: }

336: /*@
337:   DMPlexMetricNoMovement - Is node movement turned off?

339:   Input parameters:
340: . dm     - The DM

342:   Output parameters:
343: . noMove - Is node movement turned off?

345:   Level: beginner

347:   Notes:
348:   This is only used by Mmg and ParMmg (not Pragmatic).

350: .seealso: DMPlexMetricSetNoMovement(), DMPlexMetricNoInsertion(), DMPlexMetricNoSwapping(), DMPlexMetricNoSurf()
351: @*/
352: PetscErrorCode DMPlexMetricNoMovement(DM dm, PetscBool *noMove)
353: {
354:   DM_Plex       *plex = (DM_Plex *) dm->data;

356:   if (!plex->metricCtx) {
357:     PetscNew(&plex->metricCtx);
358:     DMPlexMetricSetFromOptions(dm);
359:   }
360:   *noMove = plex->metricCtx->noMove;
361:   return 0;
362: }

364: /*@
365:   DMPlexMetricSetNoSurf - Should surface modification be turned off?

367:   Input parameters:
368: + dm     - The DM
369: - noSurf - Should surface modification be turned off?

371:   Level: beginner

373:   Notes:
374:   This is only used by Mmg and ParMmg (not Pragmatic).

376: .seealso: DMPlexMetricNoSurf(), DMPlexMetricSetNoMovement(), DMPlexMetricSetNoInsertion(), DMPlexMetricSetNoSwapping()
377: @*/
378: PetscErrorCode DMPlexMetricSetNoSurf(DM dm, PetscBool noSurf)
379: {
380:   DM_Plex       *plex = (DM_Plex *) dm->data;

382:   if (!plex->metricCtx) {
383:     PetscNew(&plex->metricCtx);
384:     DMPlexMetricSetFromOptions(dm);
385:   }
386:   plex->metricCtx->noSurf = noSurf;
387:   return 0;
388: }

390: /*@
391:   DMPlexMetricNoSurf - Is surface modification turned off?

393:   Input parameters:
394: . dm     - The DM

396:   Output parameters:
397: . noSurf - Is surface modification turned off?

399:   Level: beginner

401:   Notes:
402:   This is only used by Mmg and ParMmg (not Pragmatic).

404: .seealso: DMPlexMetricSetNoSurf(), DMPlexMetricNoMovement(), DMPlexMetricNoInsertion(), DMPlexMetricNoSwapping()
405: @*/
406: PetscErrorCode DMPlexMetricNoSurf(DM dm, PetscBool *noSurf)
407: {
408:   DM_Plex       *plex = (DM_Plex *) dm->data;

410:   if (!plex->metricCtx) {
411:     PetscNew(&plex->metricCtx);
412:     DMPlexMetricSetFromOptions(dm);
413:   }
414:   *noSurf = plex->metricCtx->noSurf;
415:   return 0;
416: }

418: /*@
419:   DMPlexMetricSetMinimumMagnitude - Set the minimum tolerated metric magnitude

421:   Input parameters:
422: + dm    - The DM
423: - h_min - The minimum tolerated metric magnitude

425:   Level: beginner

427: .seealso: DMPlexMetricGetMinimumMagnitude(), DMPlexMetricSetMaximumMagnitude()
428: @*/
429: PetscErrorCode DMPlexMetricSetMinimumMagnitude(DM dm, PetscReal h_min)
430: {
431:   DM_Plex       *plex = (DM_Plex *) dm->data;

433:   if (!plex->metricCtx) {
434:     PetscNew(&plex->metricCtx);
435:     DMPlexMetricSetFromOptions(dm);
436:   }
438:   plex->metricCtx->h_min = h_min;
439:   return 0;
440: }

442: /*@
443:   DMPlexMetricGetMinimumMagnitude - Get the minimum tolerated metric magnitude

445:   Input parameters:
446: . dm    - The DM

448:   Output parameters:
449: . h_min - The minimum tolerated metric magnitude

451:   Level: beginner

453: .seealso: DMPlexMetricSetMinimumMagnitude(), DMPlexMetricGetMaximumMagnitude()
454: @*/
455: PetscErrorCode DMPlexMetricGetMinimumMagnitude(DM dm, PetscReal *h_min)
456: {
457:   DM_Plex       *plex = (DM_Plex *) dm->data;

459:   if (!plex->metricCtx) {
460:     PetscNew(&plex->metricCtx);
461:     DMPlexMetricSetFromOptions(dm);
462:   }
463:   *h_min = plex->metricCtx->h_min;
464:   return 0;
465: }

467: /*@
468:   DMPlexMetricSetMaximumMagnitude - Set the maximum tolerated metric magnitude

470:   Input parameters:
471: + dm    - The DM
472: - h_max - The maximum tolerated metric magnitude

474:   Level: beginner

476: .seealso: DMPlexMetricGetMaximumMagnitude(), DMPlexMetricSetMinimumMagnitude()
477: @*/
478: PetscErrorCode DMPlexMetricSetMaximumMagnitude(DM dm, PetscReal h_max)
479: {
480:   DM_Plex       *plex = (DM_Plex *) dm->data;

482:   if (!plex->metricCtx) {
483:     PetscNew(&plex->metricCtx);
484:     DMPlexMetricSetFromOptions(dm);
485:   }
487:   plex->metricCtx->h_max = h_max;
488:   return 0;
489: }

491: /*@
492:   DMPlexMetricGetMaximumMagnitude - Get the maximum tolerated metric magnitude

494:   Input parameters:
495: . dm    - The DM

497:   Output parameters:
498: . h_max - The maximum tolerated metric magnitude

500:   Level: beginner

502: .seealso: DMPlexMetricSetMaximumMagnitude(), DMPlexMetricGetMinimumMagnitude()
503: @*/
504: PetscErrorCode DMPlexMetricGetMaximumMagnitude(DM dm, PetscReal *h_max)
505: {
506:   DM_Plex       *plex = (DM_Plex *) dm->data;

508:   if (!plex->metricCtx) {
509:     PetscNew(&plex->metricCtx);
510:     DMPlexMetricSetFromOptions(dm);
511:   }
512:   *h_max = plex->metricCtx->h_max;
513:   return 0;
514: }

516: /*@
517:   DMPlexMetricSetMaximumAnisotropy - Set the maximum tolerated metric anisotropy

519:   Input parameters:
520: + dm    - The DM
521: - a_max - The maximum tolerated metric anisotropy

523:   Level: beginner

525:   Note: If the value zero is given then anisotropy will not be restricted. Otherwise, it should be at least one.

527: .seealso: DMPlexMetricGetMaximumAnisotropy(), DMPlexMetricSetMaximumMagnitude()
528: @*/
529: PetscErrorCode DMPlexMetricSetMaximumAnisotropy(DM dm, PetscReal a_max)
530: {
531:   DM_Plex       *plex = (DM_Plex *) dm->data;

533:   if (!plex->metricCtx) {
534:     PetscNew(&plex->metricCtx);
535:     DMPlexMetricSetFromOptions(dm);
536:   }
538:   plex->metricCtx->a_max = a_max;
539:   return 0;
540: }

542: /*@
543:   DMPlexMetricGetMaximumAnisotropy - Get the maximum tolerated metric anisotropy

545:   Input parameters:
546: . dm    - The DM

548:   Output parameters:
549: . a_max - The maximum tolerated metric anisotropy

551:   Level: beginner

553: .seealso: DMPlexMetricSetMaximumAnisotropy(), DMPlexMetricGetMaximumMagnitude()
554: @*/
555: PetscErrorCode DMPlexMetricGetMaximumAnisotropy(DM dm, PetscReal *a_max)
556: {
557:   DM_Plex       *plex = (DM_Plex *) dm->data;

559:   if (!plex->metricCtx) {
560:     PetscNew(&plex->metricCtx);
561:     DMPlexMetricSetFromOptions(dm);
562:   }
563:   *a_max = plex->metricCtx->a_max;
564:   return 0;
565: }

567: /*@
568:   DMPlexMetricSetTargetComplexity - Set the target metric complexity

570:   Input parameters:
571: + dm               - The DM
572: - targetComplexity - The target metric complexity

574:   Level: beginner

576: .seealso: DMPlexMetricGetTargetComplexity(), DMPlexMetricSetNormalizationOrder()
577: @*/
578: PetscErrorCode DMPlexMetricSetTargetComplexity(DM dm, PetscReal targetComplexity)
579: {
580:   DM_Plex       *plex = (DM_Plex *) dm->data;

582:   if (!plex->metricCtx) {
583:     PetscNew(&plex->metricCtx);
584:     DMPlexMetricSetFromOptions(dm);
585:   }
587:   plex->metricCtx->targetComplexity = targetComplexity;
588:   return 0;
589: }

591: /*@
592:   DMPlexMetricGetTargetComplexity - Get the target metric complexity

594:   Input parameters:
595: . dm               - The DM

597:   Output parameters:
598: . targetComplexity - The target metric complexity

600:   Level: beginner

602: .seealso: DMPlexMetricSetTargetComplexity(), DMPlexMetricGetNormalizationOrder()
603: @*/
604: PetscErrorCode DMPlexMetricGetTargetComplexity(DM dm, PetscReal *targetComplexity)
605: {
606:   DM_Plex       *plex = (DM_Plex *) dm->data;

608:   if (!plex->metricCtx) {
609:     PetscNew(&plex->metricCtx);
610:     DMPlexMetricSetFromOptions(dm);
611:   }
612:   *targetComplexity = plex->metricCtx->targetComplexity;
613:   return 0;
614: }

616: /*@
617:   DMPlexMetricSetNormalizationOrder - Set the order p for L-p normalization

619:   Input parameters:
620: + dm - The DM
621: - p  - The normalization order

623:   Level: beginner

625: .seealso: DMPlexMetricGetNormalizationOrder(), DMPlexMetricSetTargetComplexity()
626: @*/
627: PetscErrorCode DMPlexMetricSetNormalizationOrder(DM dm, PetscReal p)
628: {
629:   DM_Plex       *plex = (DM_Plex *) dm->data;

631:   if (!plex->metricCtx) {
632:     PetscNew(&plex->metricCtx);
633:     DMPlexMetricSetFromOptions(dm);
634:   }
636:   plex->metricCtx->p = p;
637:   return 0;
638: }

640: /*@
641:   DMPlexMetricGetNormalizationOrder - Get the order p for L-p normalization

643:   Input parameters:
644: . dm - The DM

646:   Output parameters:
647: . p - The normalization order

649:   Level: beginner

651: .seealso: DMPlexMetricSetNormalizationOrder(), DMPlexMetricGetTargetComplexity()
652: @*/
653: PetscErrorCode DMPlexMetricGetNormalizationOrder(DM dm, PetscReal *p)
654: {
655:   DM_Plex       *plex = (DM_Plex *) dm->data;

657:   if (!plex->metricCtx) {
658:     PetscNew(&plex->metricCtx);
659:     DMPlexMetricSetFromOptions(dm);
660:   }
661:   *p = plex->metricCtx->p;
662:   return 0;
663: }

665: /*@
666:   DMPlexMetricSetGradationFactor - Set the metric gradation factor

668:   Input parameters:
669: + dm   - The DM
670: - beta - The metric gradation factor

672:   Level: beginner

674:   Notes:

676:   The gradation factor is the maximum tolerated length ratio between adjacent edges.

678:   Turn off gradation by passing the value -1. Otherwise, pass a positive value.

680:   This is only used by Mmg and ParMmg (not Pragmatic).

682: .seealso: DMPlexMetricGetGradationFactor(), DMPlexMetricSetHausdorffNumber()
683: @*/
684: PetscErrorCode DMPlexMetricSetGradationFactor(DM dm, PetscReal beta)
685: {
686:   DM_Plex       *plex = (DM_Plex *) dm->data;

688:   if (!plex->metricCtx) {
689:     PetscNew(&plex->metricCtx);
690:     DMPlexMetricSetFromOptions(dm);
691:   }
692:   plex->metricCtx->gradationFactor = beta;
693:   return 0;
694: }

696: /*@
697:   DMPlexMetricGetGradationFactor - Get the metric gradation factor

699:   Input parameters:
700: . dm   - The DM

702:   Output parameters:
703: . beta - The metric gradation factor

705:   Level: beginner

707:   Notes:

709:   The gradation factor is the maximum tolerated length ratio between adjacent edges.

711:   The value -1 implies that gradation is turned off.

713:   This is only used by Mmg and ParMmg (not Pragmatic).

715: .seealso: DMPlexMetricSetGradationFactor(), DMPlexMetricGetHausdorffNumber()
716: @*/
717: PetscErrorCode DMPlexMetricGetGradationFactor(DM dm, PetscReal *beta)
718: {
719:   DM_Plex       *plex = (DM_Plex *) dm->data;

721:   if (!plex->metricCtx) {
722:     PetscNew(&plex->metricCtx);
723:     DMPlexMetricSetFromOptions(dm);
724:   }
725:   *beta = plex->metricCtx->gradationFactor;
726:   return 0;
727: }

729: /*@
730:   DMPlexMetricSetHausdorffNumber - Set the metric Hausdorff number

732:   Input parameters:
733: + dm    - The DM
734: - hausd - The metric Hausdorff number

736:   Level: beginner

738:   Notes:

740:   The Hausdorff number imposes the maximal distance between the piecewise linear approximation of the
741:   boundary and the reconstructed ideal boundary. Thus, a low Hausdorff parameter leads to refine the
742:   high curvature areas. By default, the Hausdorff value is set to 0.01, which is a suitable value for
743:   an object of size 1 in each direction. For smaller (resp. larger) objects, you may need to decrease
744:   (resp. increase) the Hausdorff parameter. (Taken from
745:   https://www.mmgtools.org/mmg-remesher-try-mmg/mmg-remesher-options/mmg-remesher-option-hausd).

747:   This is only used by Mmg and ParMmg (not Pragmatic).

749: .seealso: DMPlexMetricSetGradationFactor(), DMPlexMetricGetHausdorffNumber()
750: @*/
751: PetscErrorCode DMPlexMetricSetHausdorffNumber(DM dm, PetscReal hausd)
752: {
753:   DM_Plex       *plex = (DM_Plex *) dm->data;

755:   if (!plex->metricCtx) {
756:     PetscNew(&plex->metricCtx);
757:     DMPlexMetricSetFromOptions(dm);
758:   }
759:   plex->metricCtx->hausdorffNumber = hausd;
760:   return 0;
761: }

763: /*@
764:   DMPlexMetricGetHausdorffNumber - Get the metric Hausdorff number

766:   Input parameters:
767: . dm    - The DM

769:   Output parameters:
770: . hausd - The metric Hausdorff number

772:   Level: beginner

774:   Notes:

776:   The Hausdorff number imposes the maximal distance between the piecewise linear approximation of the
777:   boundary and the reconstructed ideal boundary. Thus, a low Hausdorff parameter leads to refine the
778:   high curvature areas. By default, the Hausdorff value is set to 0.01, which is a suitable value for
779:   an object of size 1 in each direction. For smaller (resp. larger) objects, you may need to decrease
780:   (resp. increase) the Hausdorff parameter. (Taken from
781:   https://www.mmgtools.org/mmg-remesher-try-mmg/mmg-remesher-options/mmg-remesher-option-hausd).

783:   This is only used by Mmg and ParMmg (not Pragmatic).

785: .seealso: DMPlexMetricGetGradationFactor(), DMPlexMetricSetHausdorffNumber()
786: @*/
787: PetscErrorCode DMPlexMetricGetHausdorffNumber(DM dm, PetscReal *hausd)
788: {
789:   DM_Plex       *plex = (DM_Plex *) dm->data;

791:   if (!plex->metricCtx) {
792:     PetscNew(&plex->metricCtx);
793:     DMPlexMetricSetFromOptions(dm);
794:   }
795:   *hausd = plex->metricCtx->hausdorffNumber;
796:   return 0;
797: }

799: /*@
800:   DMPlexMetricSetVerbosity - Set the verbosity of the mesh adaptation package

802:   Input parameters:
803: + dm        - The DM
804: - verbosity - The verbosity, where -1 is silent and 10 is maximum

806:   Level: beginner

808:   Notes:
809:   This is only used by Mmg and ParMmg (not Pragmatic).

811: .seealso: DMPlexMetricGetVerbosity(), DMPlexMetricSetNumIterations()
812: @*/
813: PetscErrorCode DMPlexMetricSetVerbosity(DM dm, PetscInt verbosity)
814: {
815:   DM_Plex       *plex = (DM_Plex *) dm->data;

817:   if (!plex->metricCtx) {
818:     PetscNew(&plex->metricCtx);
819:     DMPlexMetricSetFromOptions(dm);
820:   }
821:   plex->metricCtx->verbosity = verbosity;
822:   return 0;
823: }

825: /*@
826:   DMPlexMetricGetVerbosity - Get the verbosity of the mesh adaptation package

828:   Input parameters:
829: . dm        - The DM

831:   Output parameters:
832: . verbosity - The verbosity, where -1 is silent and 10 is maximum

834:   Level: beginner

836:   Notes:
837:   This is only used by Mmg and ParMmg (not Pragmatic).

839: .seealso: DMPlexMetricSetVerbosity(), DMPlexMetricGetNumIterations()
840: @*/
841: PetscErrorCode DMPlexMetricGetVerbosity(DM dm, PetscInt *verbosity)
842: {
843:   DM_Plex       *plex = (DM_Plex *) dm->data;

845:   if (!plex->metricCtx) {
846:     PetscNew(&plex->metricCtx);
847:     DMPlexMetricSetFromOptions(dm);
848:   }
849:   *verbosity = plex->metricCtx->verbosity;
850:   return 0;
851: }

853: /*@
854:   DMPlexMetricSetNumIterations - Set the number of parallel adaptation iterations

856:   Input parameters:
857: + dm      - The DM
858: - numIter - the number of parallel adaptation iterations

860:   Level: beginner

862:   Notes:
863:   This is only used by ParMmg (not Pragmatic or Mmg).

865: .seealso: DMPlexMetricSetVerbosity(), DMPlexMetricGetNumIterations()
866: @*/
867: PetscErrorCode DMPlexMetricSetNumIterations(DM dm, PetscInt numIter)
868: {
869:   DM_Plex       *plex = (DM_Plex *) dm->data;

871:   if (!plex->metricCtx) {
872:     PetscNew(&plex->metricCtx);
873:     DMPlexMetricSetFromOptions(dm);
874:   }
875:   plex->metricCtx->numIter = numIter;
876:   return 0;
877: }

879: /*@
880:   DMPlexMetricGetNumIterations - Get the number of parallel adaptation iterations

882:   Input parameters:
883: . dm      - The DM

885:   Output parameters:
886: . numIter - the number of parallel adaptation iterations

888:   Level: beginner

890:   Notes:
891:   This is only used by Mmg and ParMmg (not Pragmatic or Mmg).

893: .seealso: DMPlexMetricSetNumIterations(), DMPlexMetricGetVerbosity()
894: @*/
895: PetscErrorCode DMPlexMetricGetNumIterations(DM dm, PetscInt *numIter)
896: {
897:   DM_Plex       *plex = (DM_Plex *) dm->data;

899:   if (!plex->metricCtx) {
900:     PetscNew(&plex->metricCtx);
901:     DMPlexMetricSetFromOptions(dm);
902:   }
903:   *numIter = plex->metricCtx->numIter;
904:   return 0;
905: }

907: PetscErrorCode DMPlexP1FieldCreate_Private(DM dm, PetscInt f, PetscInt size, Vec *metric)
908: {
909:   MPI_Comm       comm;
910:   PetscFE        fe;
911:   PetscInt       dim;


914:   /* Extract metadata from dm */
915:   PetscObjectGetComm((PetscObject) dm, &comm);
916:   DMGetDimension(dm, &dim);

918:   /* Create a P1 field of the requested size */
919:   PetscFECreateLagrange(comm, dim, size, PETSC_TRUE, 1, PETSC_DETERMINE, &fe);
920:   DMSetField(dm, f, NULL, (PetscObject)fe);
921:   DMCreateDS(dm);
922:   PetscFEDestroy(&fe);
923:   DMCreateLocalVector(dm, metric);

925:   return 0;
926: }

928: /*@
929:   DMPlexMetricCreate - Create a Riemannian metric field

931:   Input parameters:
932: + dm     - The DM
933: - f      - The field number to use

935:   Output parameter:
936: . metric - The metric

938:   Level: beginner

940:   Notes:

942:   It is assumed that the DM is comprised of simplices.

944:   Command line options for Riemannian metrics:

946: + -dm_plex_metric_isotropic                 - Is the metric isotropic?
947: . -dm_plex_metric_uniform                   - Is the metric uniform?
948: . -dm_plex_metric_restrict_anisotropy_first - Should anisotropy be restricted before normalization?
949: . -dm_plex_metric_h_min                     - Minimum tolerated metric magnitude
950: . -dm_plex_metric_h_max                     - Maximum tolerated metric magnitude
951: . -dm_plex_metric_a_max                     - Maximum tolerated anisotropy
952: . -dm_plex_metric_p                         - L-p normalization order
953: - -dm_plex_metric_target_complexity         - Target metric complexity

955:   Switching between remeshers can be achieved using

957: . -dm_adaptor <pragmatic/mmg/parmmg>        - specify dm adaptor to use

959:   Further options that are only relevant to Mmg and ParMmg:

961: + -dm_plex_metric_gradation_factor          - Maximum ratio by which edge lengths may grow during gradation
962: . -dm_plex_metric_num_iterations            - Number of parallel mesh adaptation iterations for ParMmg
963: . -dm_plex_metric_no_insert                 - Should node insertion/deletion be turned off?
964: . -dm_plex_metric_no_swap                   - Should facet swapping be turned off?
965: . -dm_plex_metric_no_move                   - Should node movement be turned off?
966: - -dm_plex_metric_verbosity                 - Choose a verbosity level from -1 (silent) to 10 (maximum).

968: .seealso: DMPlexMetricCreateUniform(), DMPlexMetricCreateIsotropic()
969: @*/
970: PetscErrorCode DMPlexMetricCreate(DM dm, PetscInt f, Vec *metric)
971: {
972:   DM_Plex       *plex = (DM_Plex *) dm->data;
973:   PetscBool      isotropic, uniform;
974:   PetscInt       coordDim, Nd;

976:   if (!plex->metricCtx) {
977:     PetscNew(&plex->metricCtx);
978:     DMPlexMetricSetFromOptions(dm);
979:   }
980:   DMGetCoordinateDim(dm, &coordDim);
981:   Nd = coordDim*coordDim;
982:   DMPlexMetricIsUniform(dm, &uniform);
983:   DMPlexMetricIsIsotropic(dm, &isotropic);
984:   if (uniform) {
985:     MPI_Comm comm;

987:     PetscObjectGetComm((PetscObject) dm, &comm);
989:     VecCreate(comm, metric);
990:     VecSetSizes(*metric, 1, PETSC_DECIDE);
991:     VecSetFromOptions(*metric);
992:   } else if (isotropic) {
993:     DMPlexP1FieldCreate_Private(dm, f, 1, metric);
994:   } else {
995:     DMPlexP1FieldCreate_Private(dm, f, Nd, metric);
996:   }
997:   return 0;
998: }

1000: /*@
1001:   DMPlexMetricCreateUniform - Construct a uniform isotropic metric

1003:   Input parameters:
1004: + dm     - The DM
1005: . f      - The field number to use
1006: - alpha  - Scaling parameter for the diagonal

1008:   Output parameter:
1009: . metric - The uniform metric

1011:   Level: beginner

1013:   Note: In this case, the metric is constant in space and so is only specified for a single vertex.

1015: .seealso: DMPlexMetricCreate(), DMPlexMetricCreateIsotropic()
1016: @*/
1017: PetscErrorCode DMPlexMetricCreateUniform(DM dm, PetscInt f, PetscReal alpha, Vec *metric)
1018: {
1019:   DMPlexMetricSetUniform(dm, PETSC_TRUE);
1020:   DMPlexMetricCreate(dm, f, metric);
1023:   VecSet(*metric, alpha);
1024:   VecAssemblyBegin(*metric);
1025:   VecAssemblyEnd(*metric);
1026:   return 0;
1027: }

1029: static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1030:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1031:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1032:                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1033: {
1034:   f0[0] = u[0];
1035: }

1037: /*@
1038:   DMPlexMetricCreateIsotropic - Construct an isotropic metric from an error indicator

1040:   Input parameters:
1041: + dm        - The DM
1042: . f         - The field number to use
1043: - indicator - The error indicator

1045:   Output parameter:
1046: . metric    - The isotropic metric

1048:   Level: beginner

1050:   Notes:

1052:   It is assumed that the DM is comprised of simplices.

1054:   The indicator needs to be a scalar field. If it is not defined vertex-wise, then it is projected appropriately.

1056: .seealso: DMPlexMetricCreate(), DMPlexMetricCreateUniform()
1057: @*/
1058: PetscErrorCode DMPlexMetricCreateIsotropic(DM dm, PetscInt f, Vec indicator, Vec *metric)
1059: {
1060:   PetscInt       m, n;

1062:   DMPlexMetricSetIsotropic(dm, PETSC_TRUE);
1063:   DMPlexMetricCreate(dm, f, metric);
1064:   VecGetSize(indicator, &m);
1065:   VecGetSize(*metric, &n);
1066:   if (m == n) VecCopy(indicator, *metric);
1067:   else {
1068:     DM     dmIndi;
1069:     void (*funcs[1])(PetscInt, PetscInt, PetscInt,
1070:                      const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1071:                      const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
1072:                      PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);

1074:     VecGetDM(indicator, &dmIndi);
1075:     funcs[0] = identity;
1076:     DMProjectFieldLocal(dmIndi, 0.0, indicator, funcs, INSERT_VALUES, *metric);
1077:   }
1078:   return 0;
1079: }

1081: static PetscErrorCode LAPACKsyevFail(PetscInt dim, PetscScalar Mpos[])
1082: {
1083:   PetscInt i, j;

1085:   PetscPrintf(PETSC_COMM_SELF, "Failed to apply LAPACKsyev to the matrix\n");
1086:   for (i = 0; i < dim; ++i) {
1087:     if (i == 0) PetscPrintf(PETSC_COMM_SELF, "    [[");
1088:     else        PetscPrintf(PETSC_COMM_SELF, "     [");
1089:     for (j = 0; j < dim; ++j) {
1090:       if (j < dim-1) PetscPrintf(PETSC_COMM_SELF, "%15.8e, ", Mpos[i*dim+j]);
1091:       else           PetscPrintf(PETSC_COMM_SELF, "%15.8e", Mpos[i*dim+j]);
1092:     }
1093:     if (i < dim-1) PetscPrintf(PETSC_COMM_SELF, "]\n");
1094:     else           PetscPrintf(PETSC_COMM_SELF, "]]\n");
1095:   }
1096:   return 0;
1097: }

1099: static PetscErrorCode DMPlexMetricModify_Private(PetscInt dim, PetscReal h_min, PetscReal h_max, PetscReal a_max, PetscScalar Mp[], PetscScalar *detMp)
1100: {
1101:   PetscInt       i, j, k;
1102:   PetscReal     *eigs, max_eig, l_min = 1.0/(h_max*h_max), l_max = 1.0/(h_min*h_min), la_min = 1.0/(a_max*a_max);
1103:   PetscScalar   *Mpos;

1105:   PetscMalloc2(dim*dim, &Mpos, dim, &eigs);

1107:   /* Symmetrize */
1108:   for (i = 0; i < dim; ++i) {
1109:     Mpos[i*dim+i] = Mp[i*dim+i];
1110:     for (j = i+1; j < dim; ++j) {
1111:       Mpos[i*dim+j] = 0.5*(Mp[i*dim+j] + Mp[j*dim+i]);
1112:       Mpos[j*dim+i] = Mpos[i*dim+j];
1113:     }
1114:   }

1116:   /* Compute eigendecomposition */
1117:   if (dim == 1) {

1119:     /* Isotropic case */
1120:     eigs[0] = PetscRealPart(Mpos[0]);
1121:     Mpos[0] = 1.0;
1122:   } else {

1124:     /* Anisotropic case */
1125:     PetscScalar  *work;
1126:     PetscBLASInt lwork;

1128:     lwork = 5*dim;
1129:     PetscMalloc1(5*dim, &work);
1130:     {
1131:       PetscBLASInt lierr;
1132:       PetscBLASInt nb;

1134:       PetscBLASIntCast(dim, &nb);
1135:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1136: #if defined(PETSC_USE_COMPLEX)
1137:       {
1138:         PetscReal *rwork;
1139:         PetscMalloc1(3*dim, &rwork);
1140:         PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&nb,Mpos,&nb,eigs,work,&lwork,rwork,&lierr));
1141:         PetscFree(rwork);
1142:       }
1143: #else
1144:       PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&nb,Mpos,&nb,eigs,work,&lwork,&lierr));
1145: #endif
1146:       if (lierr) {
1147:         for (i = 0; i < dim; ++i) {
1148:           Mpos[i*dim+i] = Mp[i*dim+i];
1149:           for (j = i+1; j < dim; ++j) {
1150:             Mpos[i*dim+j] = 0.5*(Mp[i*dim+j] + Mp[j*dim+i]);
1151:             Mpos[j*dim+i] = Mpos[i*dim+j];
1152:           }
1153:         }
1154:         LAPACKsyevFail(dim, Mpos);
1155:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int) lierr);
1156:       }
1157:       PetscFPTrapPop();
1158:     }
1159:     PetscFree(work);
1160:   }

1162:   /* Reflect to positive orthant and enforce maximum and minimum size */
1163:   max_eig = 0.0;
1164:   for (i = 0; i < dim; ++i) {
1165:     eigs[i] = PetscMin(l_max, PetscMax(l_min, PetscAbsReal(eigs[i])));
1166:     max_eig = PetscMax(eigs[i], max_eig);
1167:   }

1169:   /* Enforce maximum anisotropy and compute determinant */
1170:   *detMp = 1.0;
1171:   for (i = 0; i < dim; ++i) {
1172:     if (a_max > 1.0) eigs[i] = PetscMax(eigs[i], max_eig*la_min);
1173:     *detMp *= eigs[i];
1174:   }

1176:   /* Reconstruct Hessian */
1177:   for (i = 0; i < dim; ++i) {
1178:     for (j = 0; j < dim; ++j) {
1179:       Mp[i*dim+j] = 0.0;
1180:       for (k = 0; k < dim; ++k) {
1181:         Mp[i*dim+j] += Mpos[k*dim+i] * eigs[k] * Mpos[k*dim+j];
1182:       }
1183:     }
1184:   }
1185:   PetscFree2(Mpos, eigs);

1187:   return 0;
1188: }

1190: /*@
1191:   DMPlexMetricEnforceSPD - Enforce symmetric positive-definiteness of a metric

1193:   Input parameters:
1194: + dm                 - The DM
1195: . metricIn           - The metric
1196: . restrictSizes      - Should maximum/minimum metric magnitudes be enforced?
1197: - restrictAnisotropy - Should maximum anisotropy be enforced?

1199:   Output parameter:
1200: + metricOut          - The metric
1201: - determinant        - Its determinant

1203:   Level: beginner

1205:   Notes:

1207:   Relevant command line options:

1209: + -dm_plex_metric_isotropic - Is the metric isotropic?
1210: . -dm_plex_metric_uniform   - Is the metric uniform?
1211: . -dm_plex_metric_h_min     - Minimum tolerated metric magnitude
1212: . -dm_plex_metric_h_max     - Maximum tolerated metric magnitude
1213: - -dm_plex_metric_a_max     - Maximum tolerated anisotropy

1215: .seealso: DMPlexMetricNormalize(), DMPlexMetricIntersection()
1216: @*/
1217: PetscErrorCode DMPlexMetricEnforceSPD(DM dm, Vec metricIn, PetscBool restrictSizes, PetscBool restrictAnisotropy, Vec *metricOut, Vec *determinant)
1218: {
1219:   DM             dmDet;
1220:   PetscBool      isotropic, uniform;
1221:   PetscInt       dim, vStart, vEnd, v;
1222:   PetscScalar   *met, *det;
1223:   PetscReal      h_min = 1.0e-30, h_max = 1.0e+30, a_max = 0.0;

1225:   PetscLogEventBegin(DMPLEX_MetricEnforceSPD,0,0,0,0);

1227:   /* Extract metadata from dm */
1228:   DMGetDimension(dm, &dim);
1229:   if (restrictSizes) {
1230:     DMPlexMetricGetMinimumMagnitude(dm, &h_min);
1231:     DMPlexMetricGetMaximumMagnitude(dm, &h_max);
1232:     h_min = PetscMax(h_min, 1.0e-30);
1233:     h_max = PetscMin(h_max, 1.0e+30);
1235:   }
1236:   if (restrictAnisotropy) {
1237:     DMPlexMetricGetMaximumAnisotropy(dm, &a_max);
1238:     a_max = PetscMin(a_max, 1.0e+30);
1239:   }

1241:   /* Setup output metric */
1242:   DMPlexMetricCreate(dm, 0, metricOut);
1243:   VecCopy(metricIn, *metricOut);

1245:   /* Enforce SPD and extract determinant */
1246:   VecGetArray(*metricOut, &met);
1247:   DMPlexMetricIsUniform(dm, &uniform);
1248:   DMPlexMetricIsIsotropic(dm, &isotropic);
1249:   if (uniform) {

1252:     /* Uniform case */
1253:     VecDuplicate(metricIn, determinant);
1254:     VecGetArray(*determinant, &det);
1255:     DMPlexMetricModify_Private(1, h_min, h_max, a_max, met, det);
1256:     VecRestoreArray(*determinant, &det);
1257:   } else {

1259:     /* Spatially varying case */
1260:     PetscInt nrow;

1262:     if (isotropic) nrow = 1;
1263:     else nrow = dim;
1264:     DMClone(dm, &dmDet);
1265:     DMPlexP1FieldCreate_Private(dmDet, 0, 1, determinant);
1266:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1267:     VecGetArray(*determinant, &det);
1268:     for (v = vStart; v < vEnd; ++v) {
1269:       PetscScalar *vmet, *vdet;
1270:       DMPlexPointLocalRef(dm, v, met, &vmet);
1271:       DMPlexPointLocalRef(dmDet, v, det, &vdet);
1272:       DMPlexMetricModify_Private(nrow, h_min, h_max, a_max, vmet, vdet);
1273:     }
1274:     VecRestoreArray(*determinant, &det);
1275:   }
1276:   VecRestoreArray(*metricOut, &met);

1278:   PetscLogEventEnd(DMPLEX_MetricEnforceSPD,0,0,0,0);
1279:   return 0;
1280: }

1282: static void detMFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1283:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1284:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1285:                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1286: {
1287:   const PetscScalar p = constants[0];

1289:   f0[0] = PetscPowScalar(u[0], p/(2.0*p + dim));
1290: }

1292: /*@
1293:   DMPlexMetricNormalize - Apply L-p normalization to a metric

1295:   Input parameters:
1296: + dm                 - The DM
1297: . metricIn           - The unnormalized metric
1298: . restrictSizes      - Should maximum/minimum metric magnitudes be enforced?
1299: - restrictAnisotropy - Should maximum metric anisotropy be enforced?

1301:   Output parameter:
1302: . metricOut          - The normalized metric

1304:   Level: beginner

1306:   Notes:

1308:   Relevant command line options:

1310: + -dm_plex_metric_isotropic                 - Is the metric isotropic?
1311: . -dm_plex_metric_uniform                   - Is the metric uniform?
1312: . -dm_plex_metric_restrict_anisotropy_first - Should anisotropy be restricted before normalization?
1313: . -dm_plex_metric_h_min                     - Minimum tolerated metric magnitude
1314: . -dm_plex_metric_h_max                     - Maximum tolerated metric magnitude
1315: . -dm_plex_metric_a_max                     - Maximum tolerated anisotropy
1316: . -dm_plex_metric_p                         - L-p normalization order
1317: - -dm_plex_metric_target_complexity         - Target metric complexity

1319: .seealso: DMPlexMetricEnforceSPD(), DMPlexMetricIntersection()
1320: @*/
1321: PetscErrorCode DMPlexMetricNormalize(DM dm, Vec metricIn, PetscBool restrictSizes, PetscBool restrictAnisotropy, Vec *metricOut)
1322: {
1323:   DM               dmDet;
1324:   MPI_Comm         comm;
1325:   PetscBool        restrictAnisotropyFirst, isotropic, uniform;
1326:   PetscDS          ds;
1327:   PetscInt         dim, Nd, vStart, vEnd, v, i;
1328:   PetscScalar     *met, *det, integral, constants[1];
1329:   PetscReal        p, h_min = 1.0e-30, h_max = 1.0e+30, a_max = 0.0, factGlob, fact, target, realIntegral;
1330:   Vec              determinant;

1332:   PetscLogEventBegin(DMPLEX_MetricNormalize,0,0,0,0);

1334:   /* Extract metadata from dm */
1335:   PetscObjectGetComm((PetscObject) dm, &comm);
1336:   DMGetDimension(dm, &dim);
1337:   DMPlexMetricIsUniform(dm, &uniform);
1338:   DMPlexMetricIsIsotropic(dm, &isotropic);
1339:   if (isotropic) Nd = 1;
1340:   else Nd = dim*dim;

1342:   /* Set up metric and ensure it is SPD */
1343:   DMPlexMetricRestrictAnisotropyFirst(dm, &restrictAnisotropyFirst);
1344:   DMPlexMetricEnforceSPD(dm, metricIn, PETSC_FALSE, (PetscBool)(restrictAnisotropy && restrictAnisotropyFirst), metricOut, &determinant);

1346:   /* Compute global normalization factor */
1347:   DMPlexMetricGetTargetComplexity(dm, &target);
1348:   DMPlexMetricGetNormalizationOrder(dm, &p);
1349:   constants[0] = p;
1350:   if (uniform) {
1352:     DM  dmTmp;
1353:     Vec tmp;

1355:     DMClone(dm, &dmTmp);
1356:     DMPlexP1FieldCreate_Private(dmTmp, 0, 1, &tmp);
1357:     VecGetArray(determinant, &det);
1358:     VecSet(tmp, det[0]);
1359:     VecRestoreArray(determinant, &det);
1360:     DMGetDS(dmTmp, &ds);
1361:     PetscDSSetConstants(ds, 1, constants);
1362:     PetscDSSetObjective(ds, 0, detMFunc);
1363:     DMPlexComputeIntegralFEM(dmTmp, tmp, &integral, NULL);
1364:     VecDestroy(&tmp);
1365:     DMDestroy(&dmTmp);
1366:   } else {
1367:     VecGetDM(determinant, &dmDet);
1368:     DMGetDS(dmDet, &ds);
1369:     PetscDSSetConstants(ds, 1, constants);
1370:     PetscDSSetObjective(ds, 0, detMFunc);
1371:     DMPlexComputeIntegralFEM(dmDet, determinant, &integral, NULL);
1372:   }
1373:   realIntegral = PetscRealPart(integral);
1375:   factGlob = PetscPowReal(target/realIntegral, 2.0/dim);

1377:   /* Apply local scaling */
1378:   if (restrictSizes) {
1379:     DMPlexMetricGetMinimumMagnitude(dm, &h_min);
1380:     DMPlexMetricGetMaximumMagnitude(dm, &h_max);
1381:     h_min = PetscMax(h_min, 1.0e-30);
1382:     h_max = PetscMin(h_max, 1.0e+30);
1384:   }
1385:   if (restrictAnisotropy && !restrictAnisotropyFirst) {
1386:     DMPlexMetricGetMaximumAnisotropy(dm, &a_max);
1387:     a_max = PetscMin(a_max, 1.0e+30);
1388:   }
1389:   VecGetArray(*metricOut, &met);
1390:   VecGetArray(determinant, &det);
1391:   if (uniform) {

1393:     /* Uniform case */
1394:     met[0] *= factGlob * PetscPowReal(PetscAbsScalar(det[0]), -1.0/(2*p+dim));
1395:     if (restrictSizes) DMPlexMetricModify_Private(1, h_min, h_max, a_max, met, det);
1396:   } else {

1398:     /* Spatially varying case */
1399:     PetscInt nrow;

1401:     if (isotropic) nrow = 1;
1402:     else nrow = dim;
1403:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1404:     for (v = vStart; v < vEnd; ++v) {
1405:       PetscScalar *Mp, *detM;

1407:       DMPlexPointLocalRef(dm, v, met, &Mp);
1408:       DMPlexPointLocalRef(dmDet, v, det, &detM);
1409:       fact = factGlob * PetscPowReal(PetscAbsScalar(detM[0]), -1.0/(2*p+dim));
1410:       for (i = 0; i < Nd; ++i) Mp[i] *= fact;
1411:       if (restrictSizes) DMPlexMetricModify_Private(nrow, h_min, h_max, a_max, Mp, detM);
1412:     }
1413:   }
1414:   VecRestoreArray(determinant, &det);
1415:   VecRestoreArray(*metricOut, &met);
1416:   VecDestroy(&determinant);
1417:   if (!uniform) DMDestroy(&dmDet);

1419:   PetscLogEventEnd(DMPLEX_MetricNormalize,0,0,0,0);
1420:   return 0;
1421: }

1423: /*@
1424:   DMPlexMetricAverage - Compute the average of a list of metrics

1426:   Input Parameters:
1427: + dm         - The DM
1428: . numMetrics - The number of metrics to be averaged
1429: . weights    - Weights for the average
1430: - metrics    - The metrics to be averaged

1432:   Output Parameter:
1433: . metricAvg  - The averaged metric

1435:   Level: beginner

1437:   Notes:
1438:   The weights should sum to unity.

1440:   If weights are not provided then an unweighted average is used.

1442: .seealso: DMPlexMetricAverage2(), DMPlexMetricAverage3(), DMPlexMetricIntersection()
1443: @*/
1444: PetscErrorCode DMPlexMetricAverage(DM dm, PetscInt numMetrics, PetscReal weights[], Vec metrics[], Vec *metricAvg)
1445: {
1446:   PetscBool haveWeights = PETSC_TRUE;
1447:   PetscInt  i, m, n;
1448:   PetscReal sum = 0.0, tol = 1.0e-10;

1450:   PetscLogEventBegin(DMPLEX_MetricAverage,0,0,0,0);
1452:   DMPlexMetricCreate(dm, 0, metricAvg);
1453:   VecSet(*metricAvg, 0.0);
1454:   VecGetSize(*metricAvg, &m);
1455:   for (i = 0; i < numMetrics; ++i) {
1456:     VecGetSize(metrics[i], &n);
1458:   }

1460:   /* Default to the unweighted case */
1461:   if (!weights) {
1462:     PetscMalloc1(numMetrics, &weights);
1463:     haveWeights = PETSC_FALSE;
1464:     for (i = 0; i < numMetrics; ++i) {weights[i] = 1.0/numMetrics; }
1465:   }

1467:   /* Check weights sum to unity */
1468:   for (i = 0; i < numMetrics; ++i) sum += weights[i];

1471:   /* Compute metric average */
1472:   for (i = 0; i < numMetrics; ++i) VecAXPY(*metricAvg, weights[i], metrics[i]);
1473:   if (!haveWeights) PetscFree(weights);

1475:   PetscLogEventEnd(DMPLEX_MetricAverage,0,0,0,0);
1476:   return 0;
1477: }

1479: /*@
1480:   DMPlexMetricAverage2 - Compute the unweighted average of two metrics

1482:   Input Parameters:
1483: + dm         - The DM
1484: . metric1    - The first metric to be averaged
1485: - metric2    - The second metric to be averaged

1487:   Output Parameter:
1488: . metricAvg  - The averaged metric

1490:   Level: beginner

1492: .seealso: DMPlexMetricAverage(), DMPlexMetricAverage3()
1493: @*/
1494: PetscErrorCode DMPlexMetricAverage2(DM dm, Vec metric1, Vec metric2, Vec *metricAvg)
1495: {
1496:   PetscReal      weights[2] = {0.5, 0.5};
1497:   Vec            metrics[2] = {metric1, metric2};

1499:   DMPlexMetricAverage(dm, 2, weights, metrics, metricAvg);
1500:   return 0;
1501: }

1503: /*@
1504:   DMPlexMetricAverage3 - Compute the unweighted average of three metrics

1506:   Input Parameters:
1507: + dm         - The DM
1508: . metric1    - The first metric to be averaged
1509: . metric2    - The second metric to be averaged
1510: - metric3    - The third metric to be averaged

1512:   Output Parameter:
1513: . metricAvg  - The averaged metric

1515:   Level: beginner

1517: .seealso: DMPlexMetricAverage(), DMPlexMetricAverage2()
1518: @*/
1519: PetscErrorCode DMPlexMetricAverage3(DM dm, Vec metric1, Vec metric2, Vec metric3, Vec *metricAvg)
1520: {
1521:   PetscReal      weights[3] = {1.0/3.0, 1.0/3.0, 1.0/3.0};
1522:   Vec            metrics[3] = {metric1, metric2, metric3};

1524:   DMPlexMetricAverage(dm, 3, weights, metrics, metricAvg);
1525:   return 0;
1526: }

1528: static PetscErrorCode DMPlexMetricIntersection_Private(PetscInt dim, PetscScalar M1[], PetscScalar M2[])
1529: {
1530:   PetscInt       i, j, k, l, m;
1531:   PetscReal     *evals, *evals1;
1532:   PetscScalar   *evecs, *sqrtM1, *isqrtM1;


1535:   /* Isotropic case */
1536:   if (dim == 1) {
1537:     M2[0] = (PetscScalar)PetscMin(PetscRealPart(M1[0]), PetscRealPart(M2[0]));
1538:     return 0;
1539:   }

1541:   /* Anisotropic case */
1542:   PetscMalloc5(dim*dim, &evecs, dim*dim, &sqrtM1, dim*dim, &isqrtM1, dim, &evals, dim, &evals1);
1543:   for (i = 0; i < dim; ++i) {
1544:     for (j = 0; j < dim; ++j) {
1545:       evecs[i*dim+j] = M1[i*dim+j];
1546:     }
1547:   }
1548:   {
1549:     PetscScalar *work;
1550:     PetscBLASInt lwork;

1552:     lwork = 5*dim;
1553:     PetscMalloc1(5*dim, &work);
1554:     {
1555:       PetscBLASInt lierr, nb;
1556:       PetscReal    sqrtk;

1558:       /* Compute eigendecomposition of M1 */
1559:       PetscBLASIntCast(dim, &nb);
1560:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1561: #if defined(PETSC_USE_COMPLEX)
1562:       {
1563:         PetscReal *rwork;
1564:         PetscMalloc1(3*dim, &rwork);
1565:         PetscStackCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals1, work, &lwork, rwork, &lierr));
1566:         PetscFree(rwork);
1567:       }
1568: #else
1569:       PetscStackCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals1, work, &lwork, &lierr));
1570: #endif
1571:       if (lierr) {
1572:         LAPACKsyevFail(dim, M1);
1573:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int) lierr);
1574:       }
1575:       PetscFPTrapPop();

1577:       /* Compute square root and reciprocal */
1578:       for (i = 0; i < dim; ++i) {
1579:         for (j = 0; j < dim; ++j) {
1580:           sqrtM1[i*dim+j] = 0.0;
1581:           isqrtM1[i*dim+j] = 0.0;
1582:           for (k = 0; k < dim; ++k) {
1583:             sqrtk = PetscSqrtReal(evals1[k]);
1584:             sqrtM1[i*dim+j] += evecs[k*dim+i] * sqrtk * evecs[k*dim+j];
1585:             isqrtM1[i*dim+j] += evecs[k*dim+i] * (1.0/sqrtk) * evecs[k*dim+j];
1586:           }
1587:         }
1588:       }

1590:       /* Map into the space spanned by the eigenvectors of M1 */
1591:       for (i = 0; i < dim; ++i) {
1592:         for (j = 0; j < dim; ++j) {
1593:           evecs[i*dim+j] = 0.0;
1594:           for (k = 0; k < dim; ++k) {
1595:             for (l = 0; l < dim; ++l) {
1596:               evecs[i*dim+j] += isqrtM1[i*dim+k] * M2[l*dim+k] * isqrtM1[j*dim+l];
1597:             }
1598:           }
1599:         }
1600:       }

1602:       /* Compute eigendecomposition */
1603:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1604: #if defined(PETSC_USE_COMPLEX)
1605:       {
1606:         PetscReal *rwork;
1607:         PetscMalloc1(3*dim, &rwork);
1608:         PetscStackCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals, work, &lwork, rwork, &lierr));
1609:         PetscFree(rwork);
1610:       }
1611: #else
1612:       PetscStackCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals, work, &lwork, &lierr));
1613: #endif
1614:       if (lierr) {
1615:         for (i = 0; i < dim; ++i) {
1616:           for (j = 0; j < dim; ++j) {
1617:             evecs[i*dim+j] = 0.0;
1618:             for (k = 0; k < dim; ++k) {
1619:               for (l = 0; l < dim; ++l) {
1620:                 evecs[i*dim+j] += isqrtM1[i*dim+k] * M2[l*dim+k] * isqrtM1[j*dim+l];
1621:               }
1622:             }
1623:           }
1624:         }
1625:         LAPACKsyevFail(dim, evecs);
1626:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int) lierr);
1627:       }
1628:       PetscFPTrapPop();

1630:       /* Modify eigenvalues */
1631:       for (i = 0; i < dim; ++i) evals[i] = PetscMin(evals[i], evals1[i]);

1633:       /* Map back to get the intersection */
1634:       for (i = 0; i < dim; ++i) {
1635:         for (j = 0; j < dim; ++j) {
1636:           M2[i*dim+j] = 0.0;
1637:           for (k = 0; k < dim; ++k) {
1638:             for (l = 0; l < dim; ++l) {
1639:               for (m = 0; m < dim; ++m) {
1640:                 M2[i*dim+j] += sqrtM1[i*dim+k] * evecs[l*dim+k] * evals[l] * evecs[l*dim+m] * sqrtM1[j*dim+m];
1641:               }
1642:             }
1643:           }
1644:         }
1645:       }
1646:     }
1647:     PetscFree(work);
1648:   }
1649:   PetscFree5(evecs, sqrtM1, isqrtM1, evals, evals1);
1650:   return 0;
1651: }

1653: /*@
1654:   DMPlexMetricIntersection - Compute the intersection of a list of metrics

1656:   Input Parameters:
1657: + dm         - The DM
1658: . numMetrics - The number of metrics to be intersected
1659: - metrics    - The metrics to be intersected

1661:   Output Parameter:
1662: . metricInt  - The intersected metric

1664:   Level: beginner

1666:   Notes:

1668:   The intersection of a list of metrics has the maximal ellipsoid which fits within the ellipsoids of the component metrics.

1670:   The implementation used here is only consistent with the maximal ellipsoid definition in the case numMetrics = 2.

1672: .seealso: DMPlexMetricIntersection2(), DMPlexMetricIntersection3(), DMPlexMetricAverage()
1673: @*/
1674: PetscErrorCode DMPlexMetricIntersection(DM dm, PetscInt numMetrics, Vec metrics[], Vec *metricInt)
1675: {
1676:   PetscBool      isotropic, uniform;
1677:   PetscInt       v, i, m, n;
1678:   PetscScalar   *met, *meti;

1680:   PetscLogEventBegin(DMPLEX_MetricIntersection,0,0,0,0);

1683:   /* Copy over the first metric */
1684:   DMPlexMetricCreate(dm, 0, metricInt);
1685:   VecCopy(metrics[0], *metricInt);
1686:   if (numMetrics == 1) return 0;
1687:   VecGetSize(*metricInt, &m);
1688:   for (i = 0; i < numMetrics; ++i) {
1689:     VecGetSize(metrics[i], &n);
1691:   }

1693:   /* Intersect subsequent metrics in turn */
1694:   DMPlexMetricIsUniform(dm, &uniform);
1695:   DMPlexMetricIsIsotropic(dm, &isotropic);
1696:   if (uniform) {

1698:     /* Uniform case */
1699:     VecGetArray(*metricInt, &met);
1700:     for (i = 1; i < numMetrics; ++i) {
1701:       VecGetArray(metrics[i], &meti);
1702:       DMPlexMetricIntersection_Private(1, meti, met);
1703:       VecRestoreArray(metrics[i], &meti);
1704:     }
1705:     VecRestoreArray(*metricInt, &met);
1706:   } else {

1708:     /* Spatially varying case */
1709:     PetscInt     dim, vStart, vEnd, nrow;
1710:     PetscScalar *M, *Mi;

1712:     DMGetDimension(dm, &dim);
1713:     if (isotropic) nrow = 1;
1714:     else nrow = dim;
1715:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1716:     VecGetArray(*metricInt, &met);
1717:     for (i = 1; i < numMetrics; ++i) {
1718:       VecGetArray(metrics[i], &meti);
1719:       for (v = vStart; v < vEnd; ++v) {
1720:         DMPlexPointLocalRef(dm, v, met, &M);
1721:         DMPlexPointLocalRef(dm, v, meti, &Mi);
1722:         DMPlexMetricIntersection_Private(nrow, Mi, M);
1723:       }
1724:       VecRestoreArray(metrics[i], &meti);
1725:     }
1726:     VecRestoreArray(*metricInt, &met);
1727:   }

1729:   PetscLogEventEnd(DMPLEX_MetricIntersection,0,0,0,0);
1730:   return 0;
1731: }

1733: /*@
1734:   DMPlexMetricIntersection2 - Compute the intersection of two metrics

1736:   Input Parameters:
1737: + dm        - The DM
1738: . metric1   - The first metric to be intersected
1739: - metric2   - The second metric to be intersected

1741:   Output Parameter:
1742: . metricInt - The intersected metric

1744:   Level: beginner

1746: .seealso: DMPlexMetricIntersection(), DMPlexMetricIntersection3()
1747: @*/
1748: PetscErrorCode DMPlexMetricIntersection2(DM dm, Vec metric1, Vec metric2, Vec *metricInt)
1749: {
1750:   Vec            metrics[2] = {metric1, metric2};

1752:   DMPlexMetricIntersection(dm, 2, metrics, metricInt);
1753:   return 0;
1754: }

1756: /*@
1757:   DMPlexMetricIntersection3 - Compute the intersection of three metrics

1759:   Input Parameters:
1760: + dm        - The DM
1761: . metric1   - The first metric to be intersected
1762: . metric2   - The second metric to be intersected
1763: - metric3   - The third metric to be intersected

1765:   Output Parameter:
1766: . metricInt - The intersected metric

1768:   Level: beginner

1770: .seealso: DMPlexMetricIntersection(), DMPlexMetricIntersection2()
1771: @*/
1772: PetscErrorCode DMPlexMetricIntersection3(DM dm, Vec metric1, Vec metric2, Vec metric3, Vec *metricInt)
1773: {
1774:   Vec            metrics[3] = {metric1, metric2, metric3};

1776:   DMPlexMetricIntersection(dm, 3, metrics, metricInt);
1777:   return 0;
1778: }