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: }