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