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

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

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

 61:   Level: beginner

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

 69:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
 70:   plex->metricCtx->isotropic = isotropic;
 71:   return 0;
 72: }

 74: /*@
 75:   DMPlexMetricIsIsotropic - Is a metric isotropic?

 77:   Input parameters:
 78: . dm        - The DM

 80:   Output parameters:
 81: . isotropic - Is the metric isotropic?

 83:   Level: beginner

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

 91:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
 92:   *isotropic = plex->metricCtx->isotropic;
 93:   return 0;
 94: }

 96: /*@
 97:   DMPlexMetricSetUniform - Record whether a metric is uniform

 99:   Input parameters:
100: + dm      - The DM
101: - uniform - Is the metric uniform?

103:   Level: beginner

105:   Notes:

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

109: .seealso: `DMPlexMetricIsUniform()`, `DMPlexMetricSetIsotropic()`, `DMPlexMetricSetRestrictAnisotropyFirst()`
110: @*/
111: PetscErrorCode DMPlexMetricSetUniform(DM dm, PetscBool uniform)
112: {
113:   DM_Plex *plex = (DM_Plex *)dm->data;

115:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
116:   plex->metricCtx->uniform = uniform;
117:   if (uniform) plex->metricCtx->isotropic = uniform;
118:   return 0;
119: }

121: /*@
122:   DMPlexMetricIsUniform - Is a metric uniform?

124:   Input parameters:
125: . dm      - The DM

127:   Output parameters:
128: . uniform - Is the metric uniform?

130:   Level: beginner

132: .seealso: `DMPlexMetricSetUniform()`, `DMPlexMetricIsIsotropic()`, `DMPlexMetricRestrictAnisotropyFirst()`
133: @*/
134: PetscErrorCode DMPlexMetricIsUniform(DM dm, PetscBool *uniform)
135: {
136:   DM_Plex *plex = (DM_Plex *)dm->data;

138:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
139:   *uniform = plex->metricCtx->uniform;
140:   return 0;
141: }

143: /*@
144:   DMPlexMetricSetRestrictAnisotropyFirst - Record whether anisotropy should be restricted before normalization

146:   Input parameters:
147: + dm                      - The DM
148: - restrictAnisotropyFirst - Should anisotropy be normalized first?

150:   Level: beginner

152: .seealso: `DMPlexMetricSetIsotropic()`, `DMPlexMetricRestrictAnisotropyFirst()`
153: @*/
154: PetscErrorCode DMPlexMetricSetRestrictAnisotropyFirst(DM dm, PetscBool restrictAnisotropyFirst)
155: {
156:   DM_Plex *plex = (DM_Plex *)dm->data;

158:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
159:   plex->metricCtx->restrictAnisotropyFirst = restrictAnisotropyFirst;
160:   return 0;
161: }

163: /*@
164:   DMPlexMetricRestrictAnisotropyFirst - Is anisotropy restricted before normalization or after?

166:   Input parameters:
167: . dm                      - The DM

169:   Output parameters:
170: . restrictAnisotropyFirst - Is anisotropy be normalized first?

172:   Level: beginner

174: .seealso: `DMPlexMetricIsIsotropic()`, `DMPlexMetricSetRestrictAnisotropyFirst()`
175: @*/
176: PetscErrorCode DMPlexMetricRestrictAnisotropyFirst(DM dm, PetscBool *restrictAnisotropyFirst)
177: {
178:   DM_Plex *plex = (DM_Plex *)dm->data;

180:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
181:   *restrictAnisotropyFirst = plex->metricCtx->restrictAnisotropyFirst;
182:   return 0;
183: }

185: /*@
186:   DMPlexMetricSetNoInsertion - Should node insertion and deletion be turned off?

188:   Input parameters:
189: + dm       - The DM
190: - noInsert - Should node insertion and deletion be turned off?

192:   Level: beginner

194:   Notes:
195:   This is only used by Mmg and ParMmg (not Pragmatic).

197: .seealso: `DMPlexMetricNoInsertion()`, `DMPlexMetricSetNoSwapping()`, `DMPlexMetricSetNoMovement()`, `DMPlexMetricSetNoSurf()`
198: @*/
199: PetscErrorCode DMPlexMetricSetNoInsertion(DM dm, PetscBool noInsert)
200: {
201:   DM_Plex *plex = (DM_Plex *)dm->data;

203:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
204:   plex->metricCtx->noInsert = noInsert;
205:   return 0;
206: }

208: /*@
209:   DMPlexMetricNoInsertion - Are node insertion and deletion turned off?

211:   Input parameters:
212: . dm       - The DM

214:   Output parameters:
215: . noInsert - Are node insertion and deletion turned off?

217:   Level: beginner

219:   Notes:
220:   This is only used by Mmg and ParMmg (not Pragmatic).

222: .seealso: `DMPlexMetricSetNoInsertion()`, `DMPlexMetricNoSwapping()`, `DMPlexMetricNoMovement()`, `DMPlexMetricNoSurf()`
223: @*/
224: PetscErrorCode DMPlexMetricNoInsertion(DM dm, PetscBool *noInsert)
225: {
226:   DM_Plex *plex = (DM_Plex *)dm->data;

228:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
229:   *noInsert = plex->metricCtx->noInsert;
230:   return 0;
231: }

233: /*@
234:   DMPlexMetricSetNoSwapping - Should facet swapping be turned off?

236:   Input parameters:
237: + dm     - The DM
238: - noSwap - Should facet swapping be turned off?

240:   Level: beginner

242:   Notes:
243:   This is only used by Mmg and ParMmg (not Pragmatic).

245: .seealso: `DMPlexMetricNoSwapping()`, `DMPlexMetricSetNoInsertion()`, `DMPlexMetricSetNoMovement()`, `DMPlexMetricSetNoSurf()`
246: @*/
247: PetscErrorCode DMPlexMetricSetNoSwapping(DM dm, PetscBool noSwap)
248: {
249:   DM_Plex *plex = (DM_Plex *)dm->data;

251:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
252:   plex->metricCtx->noSwap = noSwap;
253:   return 0;
254: }

256: /*@
257:   DMPlexMetricNoSwapping - Is facet swapping turned off?

259:   Input parameters:
260: . dm     - The DM

262:   Output parameters:
263: . noSwap - Is facet swapping turned off?

265:   Level: beginner

267:   Notes:
268:   This is only used by Mmg and ParMmg (not Pragmatic).

270: .seealso: `DMPlexMetricSetNoSwapping()`, `DMPlexMetricNoInsertion()`, `DMPlexMetricNoMovement()`, `DMPlexMetricNoSurf()`
271: @*/
272: PetscErrorCode DMPlexMetricNoSwapping(DM dm, PetscBool *noSwap)
273: {
274:   DM_Plex *plex = (DM_Plex *)dm->data;

276:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
277:   *noSwap = plex->metricCtx->noSwap;
278:   return 0;
279: }

281: /*@
282:   DMPlexMetricSetNoMovement - Should node movement be turned off?

284:   Input parameters:
285: + dm     - The DM
286: - noMove - Should node movement be turned off?

288:   Level: beginner

290:   Notes:
291:   This is only used by Mmg and ParMmg (not Pragmatic).

293: .seealso: `DMPlexMetricNoMovement()`, `DMPlexMetricSetNoInsertion()`, `DMPlexMetricSetNoSwapping()`, `DMPlexMetricSetNoSurf()`
294: @*/
295: PetscErrorCode DMPlexMetricSetNoMovement(DM dm, PetscBool noMove)
296: {
297:   DM_Plex *plex = (DM_Plex *)dm->data;

299:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
300:   plex->metricCtx->noMove = noMove;
301:   return 0;
302: }

304: /*@
305:   DMPlexMetricNoMovement - Is node movement turned off?

307:   Input parameters:
308: . dm     - The DM

310:   Output parameters:
311: . noMove - Is node movement turned off?

313:   Level: beginner

315:   Notes:
316:   This is only used by Mmg and ParMmg (not Pragmatic).

318: .seealso: `DMPlexMetricSetNoMovement()`, `DMPlexMetricNoInsertion()`, `DMPlexMetricNoSwapping()`, `DMPlexMetricNoSurf()`
319: @*/
320: PetscErrorCode DMPlexMetricNoMovement(DM dm, PetscBool *noMove)
321: {
322:   DM_Plex *plex = (DM_Plex *)dm->data;

324:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
325:   *noMove = plex->metricCtx->noMove;
326:   return 0;
327: }

329: /*@
330:   DMPlexMetricSetNoSurf - Should surface modification be turned off?

332:   Input parameters:
333: + dm     - The DM
334: - noSurf - Should surface modification be turned off?

336:   Level: beginner

338:   Notes:
339:   This is only used by Mmg and ParMmg (not Pragmatic).

341: .seealso: `DMPlexMetricNoSurf()`, `DMPlexMetricSetNoMovement()`, `DMPlexMetricSetNoInsertion()`, `DMPlexMetricSetNoSwapping()`
342: @*/
343: PetscErrorCode DMPlexMetricSetNoSurf(DM dm, PetscBool noSurf)
344: {
345:   DM_Plex *plex = (DM_Plex *)dm->data;

347:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
348:   plex->metricCtx->noSurf = noSurf;
349:   return 0;
350: }

352: /*@
353:   DMPlexMetricNoSurf - Is surface modification turned off?

355:   Input parameters:
356: . dm     - The DM

358:   Output parameters:
359: . noSurf - Is surface modification turned off?

361:   Level: beginner

363:   Notes:
364:   This is only used by Mmg and ParMmg (not Pragmatic).

366: .seealso: `DMPlexMetricSetNoSurf()`, `DMPlexMetricNoMovement()`, `DMPlexMetricNoInsertion()`, `DMPlexMetricNoSwapping()`
367: @*/
368: PetscErrorCode DMPlexMetricNoSurf(DM dm, PetscBool *noSurf)
369: {
370:   DM_Plex *plex = (DM_Plex *)dm->data;

372:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
373:   *noSurf = plex->metricCtx->noSurf;
374:   return 0;
375: }

377: /*@
378:   DMPlexMetricSetMinimumMagnitude - Set the minimum tolerated metric magnitude

380:   Input parameters:
381: + dm    - The DM
382: - h_min - The minimum tolerated metric magnitude

384:   Level: beginner

386: .seealso: `DMPlexMetricGetMinimumMagnitude()`, `DMPlexMetricSetMaximumMagnitude()`
387: @*/
388: PetscErrorCode DMPlexMetricSetMinimumMagnitude(DM dm, PetscReal h_min)
389: {
390:   DM_Plex *plex = (DM_Plex *)dm->data;

392:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
394:   plex->metricCtx->h_min = h_min;
395:   return 0;
396: }

398: /*@
399:   DMPlexMetricGetMinimumMagnitude - Get the minimum tolerated metric magnitude

401:   Input parameters:
402: . dm    - The DM

404:   Output parameters:
405: . h_min - The minimum tolerated metric magnitude

407:   Level: beginner

409: .seealso: `DMPlexMetricSetMinimumMagnitude()`, `DMPlexMetricGetMaximumMagnitude()`
410: @*/
411: PetscErrorCode DMPlexMetricGetMinimumMagnitude(DM dm, PetscReal *h_min)
412: {
413:   DM_Plex *plex = (DM_Plex *)dm->data;

415:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
416:   *h_min = plex->metricCtx->h_min;
417:   return 0;
418: }

420: /*@
421:   DMPlexMetricSetMaximumMagnitude - Set the maximum tolerated metric magnitude

423:   Input parameters:
424: + dm    - The DM
425: - h_max - The maximum tolerated metric magnitude

427:   Level: beginner

429: .seealso: `DMPlexMetricGetMaximumMagnitude()`, `DMPlexMetricSetMinimumMagnitude()`
430: @*/
431: PetscErrorCode DMPlexMetricSetMaximumMagnitude(DM dm, PetscReal h_max)
432: {
433:   DM_Plex *plex = (DM_Plex *)dm->data;

435:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
437:   plex->metricCtx->h_max = h_max;
438:   return 0;
439: }

441: /*@
442:   DMPlexMetricGetMaximumMagnitude - Get the maximum tolerated metric magnitude

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

447:   Output parameters:
448: . h_max - The maximum tolerated metric magnitude

450:   Level: beginner

452: .seealso: `DMPlexMetricSetMaximumMagnitude()`, `DMPlexMetricGetMinimumMagnitude()`
453: @*/
454: PetscErrorCode DMPlexMetricGetMaximumMagnitude(DM dm, PetscReal *h_max)
455: {
456:   DM_Plex *plex = (DM_Plex *)dm->data;

458:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
459:   *h_max = plex->metricCtx->h_max;
460:   return 0;
461: }

463: /*@
464:   DMPlexMetricSetMaximumAnisotropy - Set the maximum tolerated metric anisotropy

466:   Input parameters:
467: + dm    - The DM
468: - a_max - The maximum tolerated metric anisotropy

470:   Level: beginner

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

474: .seealso: `DMPlexMetricGetMaximumAnisotropy()`, `DMPlexMetricSetMaximumMagnitude()`
475: @*/
476: PetscErrorCode DMPlexMetricSetMaximumAnisotropy(DM dm, PetscReal a_max)
477: {
478:   DM_Plex *plex = (DM_Plex *)dm->data;

480:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
482:   plex->metricCtx->a_max = a_max;
483:   return 0;
484: }

486: /*@
487:   DMPlexMetricGetMaximumAnisotropy - Get the maximum tolerated metric anisotropy

489:   Input parameters:
490: . dm    - The DM

492:   Output parameters:
493: . a_max - The maximum tolerated metric anisotropy

495:   Level: beginner

497: .seealso: `DMPlexMetricSetMaximumAnisotropy()`, `DMPlexMetricGetMaximumMagnitude()`
498: @*/
499: PetscErrorCode DMPlexMetricGetMaximumAnisotropy(DM dm, PetscReal *a_max)
500: {
501:   DM_Plex *plex = (DM_Plex *)dm->data;

503:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
504:   *a_max = plex->metricCtx->a_max;
505:   return 0;
506: }

508: /*@
509:   DMPlexMetricSetTargetComplexity - Set the target metric complexity

511:   Input parameters:
512: + dm               - The DM
513: - targetComplexity - The target metric complexity

515:   Level: beginner

517: .seealso: `DMPlexMetricGetTargetComplexity()`, `DMPlexMetricSetNormalizationOrder()`
518: @*/
519: PetscErrorCode DMPlexMetricSetTargetComplexity(DM dm, PetscReal targetComplexity)
520: {
521:   DM_Plex *plex = (DM_Plex *)dm->data;

523:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
525:   plex->metricCtx->targetComplexity = targetComplexity;
526:   return 0;
527: }

529: /*@
530:   DMPlexMetricGetTargetComplexity - Get the target metric complexity

532:   Input parameters:
533: . dm               - The DM

535:   Output parameters:
536: . targetComplexity - The target metric complexity

538:   Level: beginner

540: .seealso: `DMPlexMetricSetTargetComplexity()`, `DMPlexMetricGetNormalizationOrder()`
541: @*/
542: PetscErrorCode DMPlexMetricGetTargetComplexity(DM dm, PetscReal *targetComplexity)
543: {
544:   DM_Plex *plex = (DM_Plex *)dm->data;

546:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
547:   *targetComplexity = plex->metricCtx->targetComplexity;
548:   return 0;
549: }

551: /*@
552:   DMPlexMetricSetNormalizationOrder - Set the order p for L-p normalization

554:   Input parameters:
555: + dm - The DM
556: - p  - The normalization order

558:   Level: beginner

560: .seealso: `DMPlexMetricGetNormalizationOrder()`, `DMPlexMetricSetTargetComplexity()`
561: @*/
562: PetscErrorCode DMPlexMetricSetNormalizationOrder(DM dm, PetscReal p)
563: {
564:   DM_Plex *plex = (DM_Plex *)dm->data;

566:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
568:   plex->metricCtx->p = p;
569:   return 0;
570: }

572: /*@
573:   DMPlexMetricGetNormalizationOrder - Get the order p for L-p normalization

575:   Input parameters:
576: . dm - The DM

578:   Output parameters:
579: . p - The normalization order

581:   Level: beginner

583: .seealso: `DMPlexMetricSetNormalizationOrder()`, `DMPlexMetricGetTargetComplexity()`
584: @*/
585: PetscErrorCode DMPlexMetricGetNormalizationOrder(DM dm, PetscReal *p)
586: {
587:   DM_Plex *plex = (DM_Plex *)dm->data;

589:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
590:   *p = plex->metricCtx->p;
591:   return 0;
592: }

594: /*@
595:   DMPlexMetricSetGradationFactor - Set the metric gradation factor

597:   Input parameters:
598: + dm   - The DM
599: - beta - The metric gradation factor

601:   Level: beginner

603:   Notes:

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

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

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

611: .seealso: `DMPlexMetricGetGradationFactor()`, `DMPlexMetricSetHausdorffNumber()`
612: @*/
613: PetscErrorCode DMPlexMetricSetGradationFactor(DM dm, PetscReal beta)
614: {
615:   DM_Plex *plex = (DM_Plex *)dm->data;

617:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
619:   plex->metricCtx->gradationFactor = beta;
620:   return 0;
621: }

623: /*@
624:   DMPlexMetricGetGradationFactor - Get the metric gradation factor

626:   Input parameters:
627: . dm   - The DM

629:   Output parameters:
630: . beta - The metric gradation factor

632:   Level: beginner

634:   Notes:

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

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

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

642: .seealso: `DMPlexMetricSetGradationFactor()`, `DMPlexMetricGetHausdorffNumber()`
643: @*/
644: PetscErrorCode DMPlexMetricGetGradationFactor(DM dm, PetscReal *beta)
645: {
646:   DM_Plex *plex = (DM_Plex *)dm->data;

648:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
649:   *beta = plex->metricCtx->gradationFactor;
650:   return 0;
651: }

653: /*@
654:   DMPlexMetricSetHausdorffNumber - Set the metric Hausdorff number

656:   Input parameters:
657: + dm    - The DM
658: - hausd - The metric Hausdorff number

660:   Level: beginner

662:   Notes:

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

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

673: .seealso: `DMPlexMetricSetGradationFactor()`, `DMPlexMetricGetHausdorffNumber()`
674: @*/
675: PetscErrorCode DMPlexMetricSetHausdorffNumber(DM dm, PetscReal hausd)
676: {
677:   DM_Plex *plex = (DM_Plex *)dm->data;

679:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
681:   plex->metricCtx->hausdorffNumber = hausd;
682:   return 0;
683: }

685: /*@
686:   DMPlexMetricGetHausdorffNumber - Get the metric Hausdorff number

688:   Input parameters:
689: . dm    - The DM

691:   Output parameters:
692: . hausd - The metric Hausdorff number

694:   Level: beginner

696:   Notes:

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

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

707: .seealso: `DMPlexMetricGetGradationFactor()`, `DMPlexMetricSetHausdorffNumber()`
708: @*/
709: PetscErrorCode DMPlexMetricGetHausdorffNumber(DM dm, PetscReal *hausd)
710: {
711:   DM_Plex *plex = (DM_Plex *)dm->data;

713:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
714:   *hausd = plex->metricCtx->hausdorffNumber;
715:   return 0;
716: }

718: /*@
719:   DMPlexMetricSetVerbosity - Set the verbosity of the mesh adaptation package

721:   Input parameters:
722: + dm        - The DM
723: - verbosity - The verbosity, where -1 is silent and 10 is maximum

725:   Level: beginner

727:   Notes:
728:   This is only used by Mmg and ParMmg (not Pragmatic).

730: .seealso: `DMPlexMetricGetVerbosity()`, `DMPlexMetricSetNumIterations()`
731: @*/
732: PetscErrorCode DMPlexMetricSetVerbosity(DM dm, PetscInt verbosity)
733: {
734:   DM_Plex *plex = (DM_Plex *)dm->data;

736:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
737:   plex->metricCtx->verbosity = verbosity;
738:   return 0;
739: }

741: /*@
742:   DMPlexMetricGetVerbosity - Get the verbosity of the mesh adaptation package

744:   Input parameters:
745: . dm        - The DM

747:   Output parameters:
748: . verbosity - The verbosity, where -1 is silent and 10 is maximum

750:   Level: beginner

752:   Notes:
753:   This is only used by Mmg and ParMmg (not Pragmatic).

755: .seealso: `DMPlexMetricSetVerbosity()`, `DMPlexMetricGetNumIterations()`
756: @*/
757: PetscErrorCode DMPlexMetricGetVerbosity(DM dm, PetscInt *verbosity)
758: {
759:   DM_Plex *plex = (DM_Plex *)dm->data;

761:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
762:   *verbosity = plex->metricCtx->verbosity;
763:   return 0;
764: }

766: /*@
767:   DMPlexMetricSetNumIterations - Set the number of parallel adaptation iterations

769:   Input parameters:
770: + dm      - The DM
771: - numIter - the number of parallel adaptation iterations

773:   Level: beginner

775:   Notes:
776:   This is only used by ParMmg (not Pragmatic or Mmg).

778: .seealso: `DMPlexMetricSetVerbosity()`, `DMPlexMetricGetNumIterations()`
779: @*/
780: PetscErrorCode DMPlexMetricSetNumIterations(DM dm, PetscInt numIter)
781: {
782:   DM_Plex *plex = (DM_Plex *)dm->data;

784:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
785:   plex->metricCtx->numIter = numIter;
786:   return 0;
787: }

789: /*@
790:   DMPlexMetricGetNumIterations - Get the number of parallel adaptation iterations

792:   Input parameters:
793: . dm      - The DM

795:   Output parameters:
796: . numIter - the number of parallel adaptation iterations

798:   Level: beginner

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

803: .seealso: `DMPlexMetricSetNumIterations()`, `DMPlexMetricGetVerbosity()`
804: @*/
805: PetscErrorCode DMPlexMetricGetNumIterations(DM dm, PetscInt *numIter)
806: {
807:   DM_Plex *plex = (DM_Plex *)dm->data;

809:   if (!plex->metricCtx) DMPlexMetricSetFromOptions(dm);
810:   *numIter = plex->metricCtx->numIter;
811:   return 0;
812: }

814: PetscErrorCode DMPlexP1FieldCreate_Private(DM dm, PetscInt f, PetscInt size, Vec *metric)
815: {
816:   MPI_Comm comm;
817:   PetscFE  fe;
818:   PetscInt dim;


821:   /* Extract metadata from dm */
822:   PetscObjectGetComm((PetscObject)dm, &comm);
823:   DMGetDimension(dm, &dim);

825:   /* Create a P1 field of the requested size */
826:   PetscFECreateLagrange(comm, dim, size, PETSC_TRUE, 1, PETSC_DETERMINE, &fe);
827:   DMSetField(dm, f, NULL, (PetscObject)fe);
828:   DMCreateDS(dm);
829:   PetscFEDestroy(&fe);
830:   DMCreateLocalVector(dm, metric);

832:   return 0;
833: }

835: /*@
836:   DMPlexMetricCreate - Create a Riemannian metric field

838:   Input parameters:
839: + dm     - The DM
840: - f      - The field number to use

842:   Output parameter:
843: . metric - The metric

845:   Level: beginner

847:   Notes:

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

851:   Command line options for Riemannian metrics:

853: + -dm_plex_metric_isotropic                 - Is the metric isotropic?
854: . -dm_plex_metric_uniform                   - Is the metric uniform?
855: . -dm_plex_metric_restrict_anisotropy_first - Should anisotropy be restricted before normalization?
856: . -dm_plex_metric_h_min                     - Minimum tolerated metric magnitude
857: . -dm_plex_metric_h_max                     - Maximum tolerated metric magnitude
858: . -dm_plex_metric_a_max                     - Maximum tolerated anisotropy
859: . -dm_plex_metric_p                         - L-p normalization order
860: - -dm_plex_metric_target_complexity         - Target metric complexity

862:   Switching between remeshers can be achieved using

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

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

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

875: .seealso: `DMPlexMetricCreateUniform()`, `DMPlexMetricCreateIsotropic()`
876: @*/
877: PetscErrorCode DMPlexMetricCreate(DM dm, PetscInt f, Vec *metric)
878: {
879:   PetscBool isotropic, uniform;
880:   PetscInt  coordDim, Nd;

882:   DMGetCoordinateDim(dm, &coordDim);
883:   Nd = coordDim * coordDim;
884:   DMPlexMetricIsUniform(dm, &uniform);
885:   DMPlexMetricIsIsotropic(dm, &isotropic);
886:   if (uniform) {
887:     MPI_Comm comm;

889:     PetscObjectGetComm((PetscObject)dm, &comm);
891:     VecCreate(comm, metric);
892:     VecSetSizes(*metric, 1, PETSC_DECIDE);
893:     VecSetFromOptions(*metric);
894:   } else if (isotropic) DMPlexP1FieldCreate_Private(dm, f, 1, metric);
895:   else DMPlexP1FieldCreate_Private(dm, f, Nd, metric);
896:   return 0;
897: }

899: /*@
900:   DMPlexMetricCreateUniform - Construct a uniform isotropic metric

902:   Input parameters:
903: + dm     - The DM
904: . f      - The field number to use
905: - alpha  - Scaling parameter for the diagonal

907:   Output parameter:
908: . metric - The uniform metric

910:   Level: beginner

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

914: .seealso: `DMPlexMetricCreate()`, `DMPlexMetricCreateIsotropic()`
915: @*/
916: PetscErrorCode DMPlexMetricCreateUniform(DM dm, PetscInt f, PetscReal alpha, Vec *metric)
917: {
918:   DMPlexMetricSetUniform(dm, PETSC_TRUE);
919:   DMPlexMetricCreate(dm, f, metric);
922:   VecSet(*metric, alpha);
923:   VecAssemblyBegin(*metric);
924:   VecAssemblyEnd(*metric);
925:   return 0;
926: }

928: static void identity(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
929: {
930:   f0[0] = u[0];
931: }

933: /*@
934:   DMPlexMetricCreateIsotropic - Construct an isotropic metric from an error indicator

936:   Input parameters:
937: + dm        - The DM
938: . f         - The field number to use
939: - indicator - The error indicator

941:   Output parameter:
942: . metric    - The isotropic metric

944:   Level: beginner

946:   Notes:

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

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

952: .seealso: `DMPlexMetricCreate()`, `DMPlexMetricCreateUniform()`
953: @*/
954: PetscErrorCode DMPlexMetricCreateIsotropic(DM dm, PetscInt f, Vec indicator, Vec *metric)
955: {
956:   PetscInt m, n;

958:   DMPlexMetricSetIsotropic(dm, PETSC_TRUE);
959:   DMPlexMetricCreate(dm, f, metric);
960:   VecGetSize(indicator, &m);
961:   VecGetSize(*metric, &n);
962:   if (m == n) VecCopy(indicator, *metric);
963:   else {
964:     DM dmIndi;
965:     void (*funcs[1])(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]);

967:     VecGetDM(indicator, &dmIndi);
968:     funcs[0] = identity;
969:     DMProjectFieldLocal(dmIndi, 0.0, indicator, funcs, INSERT_VALUES, *metric);
970:   }
971:   return 0;
972: }

974: /*@
975:   DMPlexMetricDeterminantCreate - Create the determinant field for a Riemannian metric

977:   Input parameters:
978: + dm          - The DM of the metric field
979: - f           - The field number to use

981:   Output parameter:
982: + determinant - The determinant field
983: - dmDet       - The corresponding DM

985:   Level: beginner

987: .seealso: DMPlexMetricCreateUniform(), DMPlexMetricCreateIsotropic(), DMPlexMetricCreate()
988: @*/
989: PetscErrorCode DMPlexMetricDeterminantCreate(DM dm, PetscInt f, Vec *determinant, DM *dmDet)
990: {
991:   PetscBool uniform;

993:   DMPlexMetricIsUniform(dm, &uniform);
994:   DMClone(dm, dmDet);
995:   if (uniform) {
996:     MPI_Comm comm;

998:     PetscObjectGetComm((PetscObject)*dmDet, &comm);
999:     VecCreate(comm, determinant);
1000:     VecSetSizes(*determinant, 1, PETSC_DECIDE);
1001:     VecSetFromOptions(*determinant);
1002:   } else DMPlexP1FieldCreate_Private(*dmDet, f, 1, determinant);
1003:   return 0;
1004: }

1006: static PetscErrorCode LAPACKsyevFail(PetscInt dim, PetscScalar Mpos[])
1007: {
1008:   PetscInt i, j;

1010:   PetscPrintf(PETSC_COMM_SELF, "Failed to apply LAPACKsyev to the matrix\n");
1011:   for (i = 0; i < dim; ++i) {
1012:     if (i == 0) PetscPrintf(PETSC_COMM_SELF, "    [[");
1013:     else PetscPrintf(PETSC_COMM_SELF, "     [");
1014:     for (j = 0; j < dim; ++j) {
1015:       if (j < dim - 1) PetscPrintf(PETSC_COMM_SELF, "%15.8e, ", (double)PetscAbsScalar(Mpos[i * dim + j]));
1016:       else PetscPrintf(PETSC_COMM_SELF, "%15.8e", (double)PetscAbsScalar(Mpos[i * dim + j]));
1017:     }
1018:     if (i < dim - 1) PetscPrintf(PETSC_COMM_SELF, "]\n");
1019:     else PetscPrintf(PETSC_COMM_SELF, "]]\n");
1020:   }
1021:   return 0;
1022: }

1024: static PetscErrorCode DMPlexMetricModify_Private(PetscInt dim, PetscReal h_min, PetscReal h_max, PetscReal a_max, PetscScalar Mp[], PetscScalar *detMp)
1025: {
1026:   PetscInt     i, j, k;
1027:   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);
1028:   PetscScalar *Mpos;

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

1032:   /* Symmetrize */
1033:   for (i = 0; i < dim; ++i) {
1034:     Mpos[i * dim + i] = Mp[i * dim + i];
1035:     for (j = i + 1; j < dim; ++j) {
1036:       Mpos[i * dim + j] = 0.5 * (Mp[i * dim + j] + Mp[j * dim + i]);
1037:       Mpos[j * dim + i] = Mpos[i * dim + j];
1038:     }
1039:   }

1041:   /* Compute eigendecomposition */
1042:   if (dim == 1) {
1043:     /* Isotropic case */
1044:     eigs[0] = PetscRealPart(Mpos[0]);
1045:     Mpos[0] = 1.0;
1046:   } else {
1047:     /* Anisotropic case */
1048:     PetscScalar *work;
1049:     PetscBLASInt lwork;

1051:     lwork = 5 * dim;
1052:     PetscMalloc1(5 * dim, &work);
1053:     {
1054:       PetscBLASInt lierr;
1055:       PetscBLASInt nb;

1057:       PetscBLASIntCast(dim, &nb);
1058:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1059: #if defined(PETSC_USE_COMPLEX)
1060:       {
1061:         PetscReal *rwork;
1062:         PetscMalloc1(3 * dim, &rwork);
1063:         PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, Mpos, &nb, eigs, work, &lwork, rwork, &lierr));
1064:         PetscFree(rwork);
1065:       }
1066: #else
1067:       PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, Mpos, &nb, eigs, work, &lwork, &lierr));
1068: #endif
1069:       if (lierr) {
1070:         for (i = 0; i < dim; ++i) {
1071:           Mpos[i * dim + i] = Mp[i * dim + i];
1072:           for (j = i + 1; j < dim; ++j) {
1073:             Mpos[i * dim + j] = 0.5 * (Mp[i * dim + j] + Mp[j * dim + i]);
1074:             Mpos[j * dim + i] = Mpos[i * dim + j];
1075:           }
1076:         }
1077:         LAPACKsyevFail(dim, Mpos);
1078:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int)lierr);
1079:       }
1080:       PetscFPTrapPop();
1081:     }
1082:     PetscFree(work);
1083:   }

1085:   /* Reflect to positive orthant and enforce maximum and minimum size */
1086:   max_eig = 0.0;
1087:   for (i = 0; i < dim; ++i) {
1088:     eigs[i] = PetscMin(l_max, PetscMax(l_min, PetscAbsReal(eigs[i])));
1089:     max_eig = PetscMax(eigs[i], max_eig);
1090:   }

1092:   /* Enforce maximum anisotropy and compute determinant */
1093:   *detMp = 1.0;
1094:   for (i = 0; i < dim; ++i) {
1095:     if (a_max > 1.0) eigs[i] = PetscMax(eigs[i], max_eig * la_min);
1096:     *detMp *= eigs[i];
1097:   }

1099:   /* Reconstruct Hessian */
1100:   for (i = 0; i < dim; ++i) {
1101:     for (j = 0; j < dim; ++j) {
1102:       Mp[i * dim + j] = 0.0;
1103:       for (k = 0; k < dim; ++k) Mp[i * dim + j] += Mpos[k * dim + i] * eigs[k] * Mpos[k * dim + j];
1104:     }
1105:   }
1106:   PetscFree2(Mpos, eigs);

1108:   return 0;
1109: }

1111: /*@
1112:   DMPlexMetricEnforceSPD - Enforce symmetric positive-definiteness of a metric

1114:   Input parameters:
1115: + dm                 - The DM
1116: . metricIn           - The metric
1117: . restrictSizes      - Should maximum/minimum metric magnitudes be enforced?
1118: - restrictAnisotropy - Should maximum anisotropy be enforced?

1120:   Output parameter:
1121: + metricOut          - The metric
1122: - determinant        - Its determinant

1124:   Level: beginner

1126:   Notes:

1128:   Relevant command line options:

1130: + -dm_plex_metric_isotropic - Is the metric isotropic?
1131: . -dm_plex_metric_uniform   - Is the metric uniform?
1132: . -dm_plex_metric_h_min     - Minimum tolerated metric magnitude
1133: . -dm_plex_metric_h_max     - Maximum tolerated metric magnitude
1134: - -dm_plex_metric_a_max     - Maximum tolerated anisotropy

1136: .seealso: `DMPlexMetricNormalize()`, `DMPlexMetricIntersection()`
1137: @*/
1138: PetscErrorCode DMPlexMetricEnforceSPD(DM dm, Vec metricIn, PetscBool restrictSizes, PetscBool restrictAnisotropy, Vec metricOut, Vec determinant)
1139: {
1140:   DM           dmDet;
1141:   PetscBool    isotropic, uniform;
1142:   PetscInt     dim, vStart, vEnd, v;
1143:   PetscScalar *met, *det;
1144:   PetscReal    h_min = 1.0e-30, h_max = 1.0e+30, a_max = 0.0;

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

1148:   /* Extract metadata from dm */
1149:   DMGetDimension(dm, &dim);
1150:   if (restrictSizes) {
1151:     DMPlexMetricGetMinimumMagnitude(dm, &h_min);
1152:     DMPlexMetricGetMaximumMagnitude(dm, &h_max);
1153:     h_min = PetscMax(h_min, 1.0e-30);
1154:     h_max = PetscMin(h_max, 1.0e+30);
1156:   }
1157:   if (restrictAnisotropy) {
1158:     DMPlexMetricGetMaximumAnisotropy(dm, &a_max);
1159:     a_max = PetscMin(a_max, 1.0e+30);
1160:   }

1162:   /* Setup output metric */
1163:   VecCopy(metricIn, metricOut);

1165:   /* Enforce SPD and extract determinant */
1166:   VecGetArray(metricOut, &met);
1167:   DMPlexMetricIsUniform(dm, &uniform);
1168:   DMPlexMetricIsIsotropic(dm, &isotropic);
1169:   if (uniform) {

1172:     /* Uniform case */
1173:     VecGetArray(determinant, &det);
1174:     DMPlexMetricModify_Private(1, h_min, h_max, a_max, met, det);
1175:     VecRestoreArray(determinant, &det);
1176:   } else {
1177:     /* Spatially varying case */
1178:     PetscInt nrow;

1180:     if (isotropic) nrow = 1;
1181:     else nrow = dim;
1182:     VecGetDM(determinant, &dmDet);
1183:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1184:     VecGetArray(determinant, &det);
1185:     for (v = vStart; v < vEnd; ++v) {
1186:       PetscScalar *vmet, *vdet;
1187:       DMPlexPointLocalRef(dm, v, met, &vmet);
1188:       DMPlexPointLocalRef(dmDet, v, det, &vdet);
1189:       DMPlexMetricModify_Private(nrow, h_min, h_max, a_max, vmet, vdet);
1190:     }
1191:     VecRestoreArray(determinant, &det);
1192:   }
1193:   VecRestoreArray(metricOut, &met);

1195:   PetscLogEventEnd(DMPLEX_MetricEnforceSPD, 0, 0, 0, 0);
1196:   return 0;
1197: }

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

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

1206: /*@
1207:   DMPlexMetricNormalize - Apply L-p normalization to a metric

1209:   Input parameters:
1210: + dm                 - The DM
1211: . metricIn           - The unnormalized metric
1212: . restrictSizes      - Should maximum/minimum metric magnitudes be enforced?
1213: - restrictAnisotropy - Should maximum metric anisotropy be enforced?

1215:   Output parameter:
1216: . metricOut          - The normalized metric

1218:   Level: beginner

1220:   Notes:

1222:   Relevant command line options:

1224: + -dm_plex_metric_isotropic                 - Is the metric isotropic?
1225: . -dm_plex_metric_uniform                   - Is the metric uniform?
1226: . -dm_plex_metric_restrict_anisotropy_first - Should anisotropy be restricted before normalization?
1227: . -dm_plex_metric_h_min                     - Minimum tolerated metric magnitude
1228: . -dm_plex_metric_h_max                     - Maximum tolerated metric magnitude
1229: . -dm_plex_metric_a_max                     - Maximum tolerated anisotropy
1230: . -dm_plex_metric_p                         - L-p normalization order
1231: - -dm_plex_metric_target_complexity         - Target metric complexity

1233: .seealso: `DMPlexMetricEnforceSPD()`, `DMPlexMetricIntersection()`
1234: @*/
1235: PetscErrorCode DMPlexMetricNormalize(DM dm, Vec metricIn, PetscBool restrictSizes, PetscBool restrictAnisotropy, Vec metricOut, Vec determinant)
1236: {
1237:   DM           dmDet;
1238:   MPI_Comm     comm;
1239:   PetscBool    restrictAnisotropyFirst, isotropic, uniform;
1240:   PetscDS      ds;
1241:   PetscInt     dim, Nd, vStart, vEnd, v, i;
1242:   PetscScalar *met, *det, integral, constants[1];
1243:   PetscReal    p, h_min = 1.0e-30, h_max = 1.0e+30, a_max = 0.0, factGlob, fact, target, realIntegral;

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

1247:   /* Extract metadata from dm */
1248:   PetscObjectGetComm((PetscObject)dm, &comm);
1249:   DMGetDimension(dm, &dim);
1250:   DMPlexMetricIsUniform(dm, &uniform);
1251:   DMPlexMetricIsIsotropic(dm, &isotropic);
1252:   if (isotropic) Nd = 1;
1253:   else Nd = dim * dim;

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

1259:   /* Compute global normalization factor */
1260:   DMPlexMetricGetTargetComplexity(dm, &target);
1261:   DMPlexMetricGetNormalizationOrder(dm, &p);
1262:   constants[0] = p;
1263:   if (uniform) {
1265:     DM  dmTmp;
1266:     Vec tmp;

1268:     DMClone(dm, &dmTmp);
1269:     DMPlexP1FieldCreate_Private(dmTmp, 0, 1, &tmp);
1270:     VecGetArray(determinant, &det);
1271:     VecSet(tmp, det[0]);
1272:     VecRestoreArray(determinant, &det);
1273:     DMGetDS(dmTmp, &ds);
1274:     PetscDSSetConstants(ds, 1, constants);
1275:     PetscDSSetObjective(ds, 0, detMFunc);
1276:     DMPlexComputeIntegralFEM(dmTmp, tmp, &integral, NULL);
1277:     VecDestroy(&tmp);
1278:     DMDestroy(&dmTmp);
1279:   } else {
1280:     VecGetDM(determinant, &dmDet);
1281:     DMGetDS(dmDet, &ds);
1282:     PetscDSSetConstants(ds, 1, constants);
1283:     PetscDSSetObjective(ds, 0, detMFunc);
1284:     DMPlexComputeIntegralFEM(dmDet, determinant, &integral, NULL);
1285:   }
1286:   realIntegral = PetscRealPart(integral);
1288:   factGlob = PetscPowReal(target / realIntegral, 2.0 / dim);

1290:   /* Apply local scaling */
1291:   if (restrictSizes) {
1292:     DMPlexMetricGetMinimumMagnitude(dm, &h_min);
1293:     DMPlexMetricGetMaximumMagnitude(dm, &h_max);
1294:     h_min = PetscMax(h_min, 1.0e-30);
1295:     h_max = PetscMin(h_max, 1.0e+30);
1297:   }
1298:   if (restrictAnisotropy && !restrictAnisotropyFirst) {
1299:     DMPlexMetricGetMaximumAnisotropy(dm, &a_max);
1300:     a_max = PetscMin(a_max, 1.0e+30);
1301:   }
1302:   VecGetArray(metricOut, &met);
1303:   VecGetArray(determinant, &det);
1304:   if (uniform) {
1305:     /* Uniform case */
1306:     met[0] *= factGlob * PetscPowReal(PetscAbsScalar(det[0]), -1.0 / (2 * p + dim));
1307:     if (restrictSizes) DMPlexMetricModify_Private(1, h_min, h_max, a_max, met, det);
1308:   } else {
1309:     /* Spatially varying case */
1310:     PetscInt nrow;

1312:     if (isotropic) nrow = 1;
1313:     else nrow = dim;
1314:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1315:     VecGetDM(determinant, &dmDet);
1316:     for (v = vStart; v < vEnd; ++v) {
1317:       PetscScalar *Mp, *detM;

1319:       DMPlexPointLocalRef(dm, v, met, &Mp);
1320:       DMPlexPointLocalRef(dmDet, v, det, &detM);
1321:       fact = factGlob * PetscPowReal(PetscAbsScalar(detM[0]), -1.0 / (2 * p + dim));
1322:       for (i = 0; i < Nd; ++i) Mp[i] *= fact;
1323:       if (restrictSizes) DMPlexMetricModify_Private(nrow, h_min, h_max, a_max, Mp, detM);
1324:     }
1325:   }
1326:   VecRestoreArray(determinant, &det);
1327:   VecRestoreArray(metricOut, &met);

1329:   PetscLogEventEnd(DMPLEX_MetricNormalize, 0, 0, 0, 0);
1330:   return 0;
1331: }

1333: /*@
1334:   DMPlexMetricAverage - Compute the average of a list of metrics

1336:   Input Parameters:
1337: + dm         - The DM
1338: . numMetrics - The number of metrics to be averaged
1339: . weights    - Weights for the average
1340: - metrics    - The metrics to be averaged

1342:   Output Parameter:
1343: . metricAvg  - The averaged metric

1345:   Level: beginner

1347:   Notes:
1348:   The weights should sum to unity.

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

1352: .seealso: `DMPlexMetricAverage2()`, `DMPlexMetricAverage3()`, `DMPlexMetricIntersection()`
1353: @*/
1354: PetscErrorCode DMPlexMetricAverage(DM dm, PetscInt numMetrics, PetscReal weights[], Vec metrics[], Vec metricAvg)
1355: {
1356:   PetscBool haveWeights = PETSC_TRUE;
1357:   PetscInt  i, m, n;
1358:   PetscReal sum = 0.0, tol = 1.0e-10;

1360:   PetscLogEventBegin(DMPLEX_MetricAverage, 0, 0, 0, 0);
1362:   VecSet(metricAvg, 0.0);
1363:   VecGetSize(metricAvg, &m);
1364:   for (i = 0; i < numMetrics; ++i) {
1365:     VecGetSize(metrics[i], &n);
1367:   }

1369:   /* Default to the unweighted case */
1370:   if (!weights) {
1371:     PetscMalloc1(numMetrics, &weights);
1372:     haveWeights = PETSC_FALSE;
1373:     for (i = 0; i < numMetrics; ++i) weights[i] = 1.0 / numMetrics;
1374:   }

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

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

1384:   PetscLogEventEnd(DMPLEX_MetricAverage, 0, 0, 0, 0);
1385:   return 0;
1386: }

1388: /*@
1389:   DMPlexMetricAverage2 - Compute the unweighted average of two metrics

1391:   Input Parameters:
1392: + dm         - The DM
1393: . metric1    - The first metric to be averaged
1394: - metric2    - The second metric to be averaged

1396:   Output Parameter:
1397: . metricAvg  - The averaged metric

1399:   Level: beginner

1401: .seealso: `DMPlexMetricAverage()`, `DMPlexMetricAverage3()`
1402: @*/
1403: PetscErrorCode DMPlexMetricAverage2(DM dm, Vec metric1, Vec metric2, Vec metricAvg)
1404: {
1405:   PetscReal weights[2] = {0.5, 0.5};
1406:   Vec       metrics[2] = {metric1, metric2};

1408:   DMPlexMetricAverage(dm, 2, weights, metrics, metricAvg);
1409:   return 0;
1410: }

1412: /*@
1413:   DMPlexMetricAverage3 - Compute the unweighted average of three metrics

1415:   Input Parameters:
1416: + dm         - The DM
1417: . metric1    - The first metric to be averaged
1418: . metric2    - The second metric to be averaged
1419: - metric3    - The third metric to be averaged

1421:   Output Parameter:
1422: . metricAvg  - The averaged metric

1424:   Level: beginner

1426: .seealso: `DMPlexMetricAverage()`, `DMPlexMetricAverage2()`
1427: @*/
1428: PetscErrorCode DMPlexMetricAverage3(DM dm, Vec metric1, Vec metric2, Vec metric3, Vec metricAvg)
1429: {
1430:   PetscReal weights[3] = {1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0};
1431:   Vec       metrics[3] = {metric1, metric2, metric3};

1433:   DMPlexMetricAverage(dm, 3, weights, metrics, metricAvg);
1434:   return 0;
1435: }

1437: static PetscErrorCode DMPlexMetricIntersection_Private(PetscInt dim, PetscScalar M1[], PetscScalar M2[])
1438: {
1439:   PetscInt     i, j, k, l, m;
1440:   PetscReal   *evals;
1441:   PetscScalar *evecs, *sqrtM1, *isqrtM1;


1444:   /* Isotropic case */
1445:   if (dim == 1) {
1446:     M2[0] = (PetscScalar)PetscMax(PetscRealPart(M1[0]), PetscRealPart(M2[0]));
1447:     return 0;
1448:   }

1450:   /* Anisotropic case */
1451:   PetscMalloc4(dim * dim, &evecs, dim * dim, &sqrtM1, dim * dim, &isqrtM1, dim, &evals);
1452:   for (i = 0; i < dim; ++i) {
1453:     for (j = 0; j < dim; ++j) evecs[i * dim + j] = M1[i * dim + j];
1454:   }
1455:   {
1456:     PetscScalar *work;
1457:     PetscBLASInt lwork;

1459:     lwork = 5 * dim;
1460:     PetscMalloc1(5 * dim, &work);
1461:     {
1462:       PetscBLASInt lierr, nb;
1463:       PetscReal    sqrtj;

1465:       /* Compute eigendecomposition of M1 */
1466:       PetscBLASIntCast(dim, &nb);
1467:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1468: #if defined(PETSC_USE_COMPLEX)
1469:       {
1470:         PetscReal *rwork;
1471:         PetscMalloc1(3 * dim, &rwork);
1472:         PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals, work, &lwork, rwork, &lierr));
1473:         PetscFree(rwork);
1474:       }
1475: #else
1476:       PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals, work, &lwork, &lierr));
1477: #endif
1478:       if (lierr) {
1479:         LAPACKsyevFail(dim, M1);
1480:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int)lierr);
1481:       }
1482:       PetscFPTrapPop();

1484:       /* Compute square root and the reciprocal thereof */
1485:       for (i = 0; i < dim; ++i) {
1486:         for (k = 0; k < dim; ++k) {
1487:           sqrtM1[i * dim + k]  = 0.0;
1488:           isqrtM1[i * dim + k] = 0.0;
1489:           for (j = 0; j < dim; ++j) {
1490:             sqrtj = PetscSqrtReal(evals[j]);
1491:             sqrtM1[i * dim + k] += evecs[j * dim + i] * sqrtj * evecs[j * dim + k];
1492:             isqrtM1[i * dim + k] += evecs[j * dim + i] * (1.0 / sqrtj) * evecs[j * dim + k];
1493:           }
1494:         }
1495:       }

1497:       /* Map M2 into the space spanned by the eigenvectors of M1 */
1498:       for (i = 0; i < dim; ++i) {
1499:         for (l = 0; l < dim; ++l) {
1500:           evecs[i * dim + l] = 0.0;
1501:           for (j = 0; j < dim; ++j) {
1502:             for (k = 0; k < dim; ++k) evecs[i * dim + l] += isqrtM1[j * dim + i] * M2[j * dim + k] * isqrtM1[k * dim + l];
1503:           }
1504:         }
1505:       }

1507:       /* Compute eigendecomposition */
1508:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
1509: #if defined(PETSC_USE_COMPLEX)
1510:       {
1511:         PetscReal *rwork;
1512:         PetscMalloc1(3 * dim, &rwork);
1513:         PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals, work, &lwork, rwork, &lierr));
1514:         PetscFree(rwork);
1515:       }
1516: #else
1517:       PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals, work, &lwork, &lierr));
1518: #endif
1519:       if (lierr) {
1520:         for (i = 0; i < dim; ++i) {
1521:           for (l = 0; l < dim; ++l) {
1522:             evecs[i * dim + l] = 0.0;
1523:             for (j = 0; j < dim; ++j) {
1524:               for (k = 0; k < dim; ++k) evecs[i * dim + l] += isqrtM1[j * dim + i] * M2[j * dim + k] * isqrtM1[k * dim + l];
1525:             }
1526:           }
1527:         }
1528:         LAPACKsyevFail(dim, evecs);
1529:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int)lierr);
1530:       }
1531:       PetscFPTrapPop();

1533:       /* Modify eigenvalues */
1534:       for (i = 0; i < dim; ++i) evals[i] = PetscMax(evals[i], 1.0);

1536:       /* Map back to get the intersection */
1537:       for (i = 0; i < dim; ++i) {
1538:         for (m = 0; m < dim; ++m) {
1539:           M2[i * dim + m] = 0.0;
1540:           for (j = 0; j < dim; ++j) {
1541:             for (k = 0; k < dim; ++k) {
1542:               for (l = 0; l < dim; ++l) M2[i * dim + m] += sqrtM1[j * dim + i] * evecs[j * dim + k] * evals[k] * evecs[l * dim + k] * sqrtM1[l * dim + m];
1543:             }
1544:           }
1545:         }
1546:       }
1547:     }
1548:     PetscFree(work);
1549:   }
1550:   PetscFree4(evecs, sqrtM1, isqrtM1, evals);
1551:   return 0;
1552: }

1554: /*@
1555:   DMPlexMetricIntersection - Compute the intersection of a list of metrics

1557:   Input Parameters:
1558: + dm         - The DM
1559: . numMetrics - The number of metrics to be intersected
1560: - metrics    - The metrics to be intersected

1562:   Output Parameter:
1563: . metricInt  - The intersected metric

1565:   Level: beginner

1567:   Notes:

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

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

1573: .seealso: `DMPlexMetricIntersection2()`, `DMPlexMetricIntersection3()`, `DMPlexMetricAverage()`
1574: @*/
1575: PetscErrorCode DMPlexMetricIntersection(DM dm, PetscInt numMetrics, Vec metrics[], Vec metricInt)
1576: {
1577:   PetscBool    isotropic, uniform;
1578:   PetscInt     v, i, m, n;
1579:   PetscScalar *met, *meti;

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

1584:   /* Copy over the first metric */
1585:   VecCopy(metrics[0], metricInt);
1586:   if (numMetrics == 1) return 0;
1587:   VecGetSize(metricInt, &m);
1588:   for (i = 0; i < numMetrics; ++i) {
1589:     VecGetSize(metrics[i], &n);
1591:   }

1593:   /* Intersect subsequent metrics in turn */
1594:   DMPlexMetricIsUniform(dm, &uniform);
1595:   DMPlexMetricIsIsotropic(dm, &isotropic);
1596:   if (uniform) {
1597:     /* Uniform case */
1598:     VecGetArray(metricInt, &met);
1599:     for (i = 1; i < numMetrics; ++i) {
1600:       VecGetArray(metrics[i], &meti);
1601:       DMPlexMetricIntersection_Private(1, meti, met);
1602:       VecRestoreArray(metrics[i], &meti);
1603:     }
1604:     VecRestoreArray(metricInt, &met);
1605:   } else {
1606:     /* Spatially varying case */
1607:     PetscInt     dim, vStart, vEnd, nrow;
1608:     PetscScalar *M, *Mi;

1610:     DMGetDimension(dm, &dim);
1611:     if (isotropic) nrow = 1;
1612:     else nrow = dim;
1613:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1614:     VecGetArray(metricInt, &met);
1615:     for (i = 1; i < numMetrics; ++i) {
1616:       VecGetArray(metrics[i], &meti);
1617:       for (v = vStart; v < vEnd; ++v) {
1618:         DMPlexPointLocalRef(dm, v, met, &M);
1619:         DMPlexPointLocalRef(dm, v, meti, &Mi);
1620:         DMPlexMetricIntersection_Private(nrow, Mi, M);
1621:       }
1622:       VecRestoreArray(metrics[i], &meti);
1623:     }
1624:     VecRestoreArray(metricInt, &met);
1625:   }

1627:   PetscLogEventEnd(DMPLEX_MetricIntersection, 0, 0, 0, 0);
1628:   return 0;
1629: }

1631: /*@
1632:   DMPlexMetricIntersection2 - Compute the intersection of two metrics

1634:   Input Parameters:
1635: + dm        - The DM
1636: . metric1   - The first metric to be intersected
1637: - metric2   - The second metric to be intersected

1639:   Output Parameter:
1640: . metricInt - The intersected metric

1642:   Level: beginner

1644: .seealso: `DMPlexMetricIntersection()`, `DMPlexMetricIntersection3()`
1645: @*/
1646: PetscErrorCode DMPlexMetricIntersection2(DM dm, Vec metric1, Vec metric2, Vec metricInt)
1647: {
1648:   Vec metrics[2] = {metric1, metric2};

1650:   DMPlexMetricIntersection(dm, 2, metrics, metricInt);
1651:   return 0;
1652: }

1654: /*@
1655:   DMPlexMetricIntersection3 - Compute the intersection of three metrics

1657:   Input Parameters:
1658: + dm        - The DM
1659: . metric1   - The first metric to be intersected
1660: . metric2   - The second metric to be intersected
1661: - metric3   - The third metric to be intersected

1663:   Output Parameter:
1664: . metricInt - The intersected metric

1666:   Level: beginner

1668: .seealso: `DMPlexMetricIntersection()`, `DMPlexMetricIntersection2()`
1669: @*/
1670: PetscErrorCode DMPlexMetricIntersection3(DM dm, Vec metric1, Vec metric2, Vec metric3, Vec metricInt)
1671: {
1672:   Vec metrics[3] = {metric1, metric2, metric3};

1674:   DMPlexMetricIntersection(dm, 3, metrics, metricInt);
1675:   return 0;
1676: }