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: /*@
7: DMPlexMetricSetFromOptions - Configure the Riemannian metric context on a `DMPLEX` from the options database
9: Collective
11: Input Parameter:
12: . dm - The `DM`
14: Options Database Keys:
15: + -dm_plex_metric_isotropic (true|false) - Is the metric isotropic?
16: . -dm_plex_metric_uniform (true|false) - Is the metric uniform?
17: . -dm_plex_metric_restrict_anisotropy_first (true|false) - Restrict anisotropy before normalization
18: . -dm_plex_metric_no_insert (true|false) - Turn off node insertion and deletion during adaptation
19: . -dm_plex_metric_no_swap (true|false) - Turn off facet swapping
20: . -dm_plex_metric_no_move (true|false) - Turn off facet node movement
21: . -dm_plex_metric_no_surf (true|false) - Turn off surface modification
22: . -dm_plex_metric_num_iterations nits - Number of `ParMmg` adaptation iterations
23: . -dm_plex_metric_verbosity verbosity - Verbosity of the metric-based adaptation package (-1 silent, 10 maximum)
24: . -dm_plex_metric_h_min h_min - Minimum tolerated metric magnitude
25: . -dm_plex_metric_h_max h_max - Maximum tolerated metric magnitude
26: . -dm_plex_metric_a_max a_max - Maximum tolerated anisotropy
27: . -dm_plex_metric_p order - L-p normalization order
28: . -dm_plex_metric_target_complexity comp - Target metric complexity
29: . -dm_plex_metric_gradation_factor fact - Metric gradation factor
30: - -dm_plex_metric_hausdorff_number h - Metric Hausdorff number
32: Level: beginner
34: Note:
35: A metric context is created on the `DM` if none exists.
37: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexMetricSetIsotropic()`, `DMPlexMetricSetUniform()`, `DMPlexMetricSetMinimumMagnitude()`, `DMPlexMetricSetMaximumMagnitude()`, `DMPlexMetricSetTargetComplexity()`
38: @*/
39: PetscErrorCode DMPlexMetricSetFromOptions(DM dm)
40: {
41: DM_Plex *plex = (DM_Plex *)dm->data;
42: MPI_Comm comm;
43: PetscBool isotropic = PETSC_FALSE, uniform = PETSC_FALSE, restrictAnisotropyFirst = PETSC_FALSE;
44: PetscBool noInsert = PETSC_FALSE, noSwap = PETSC_FALSE, noMove = PETSC_FALSE, noSurf = PETSC_FALSE;
45: PetscInt verbosity = -1, numIter = 3;
46: 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;
48: PetscFunctionBegin;
49: if (!plex->metricCtx) PetscCall(PetscNew(&plex->metricCtx));
50: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
51: PetscOptionsBegin(comm, "", "Riemannian metric options", "DMPlexMetric");
52: PetscCall(PetscOptionsBool("-dm_plex_metric_isotropic", "Is the metric isotropic?", "DMPlexMetricCreateIsotropic", isotropic, &isotropic, NULL));
53: PetscCall(DMPlexMetricSetIsotropic(dm, isotropic));
54: PetscCall(PetscOptionsBool("-dm_plex_metric_uniform", "Is the metric uniform?", "DMPlexMetricCreateUniform", uniform, &uniform, NULL));
55: PetscCall(DMPlexMetricSetUniform(dm, uniform));
56: PetscCall(PetscOptionsBool("-dm_plex_metric_restrict_anisotropy_first", "Should anisotropy be restricted before normalization?", "DMPlexNormalize", restrictAnisotropyFirst, &restrictAnisotropyFirst, NULL));
57: PetscCall(DMPlexMetricSetRestrictAnisotropyFirst(dm, restrictAnisotropyFirst));
58: PetscCall(PetscOptionsBool("-dm_plex_metric_no_insert", "Turn off node insertion and deletion", "DMAdaptMetric", noInsert, &noInsert, NULL));
59: PetscCall(DMPlexMetricSetNoInsertion(dm, noInsert));
60: PetscCall(PetscOptionsBool("-dm_plex_metric_no_swap", "Turn off facet swapping", "DMAdaptMetric", noSwap, &noSwap, NULL));
61: PetscCall(DMPlexMetricSetNoSwapping(dm, noSwap));
62: PetscCall(PetscOptionsBool("-dm_plex_metric_no_move", "Turn off facet node movement", "DMAdaptMetric", noMove, &noMove, NULL));
63: PetscCall(DMPlexMetricSetNoMovement(dm, noMove));
64: PetscCall(PetscOptionsBool("-dm_plex_metric_no_surf", "Turn off surface modification", "DMAdaptMetric", noSurf, &noSurf, NULL));
65: PetscCall(DMPlexMetricSetNoSurf(dm, noSurf));
66: PetscCall(PetscOptionsBoundedInt("-dm_plex_metric_num_iterations", "Number of ParMmg adaptation iterations", "DMAdaptMetric", numIter, &numIter, NULL, 0));
67: PetscCall(DMPlexMetricSetNumIterations(dm, numIter));
68: PetscCall(PetscOptionsRangeInt("-dm_plex_metric_verbosity", "Verbosity of metric-based mesh adaptation package (-1 = silent, 10 = maximum)", "DMAdaptMetric", verbosity, &verbosity, NULL, -1, 10));
69: PetscCall(DMPlexMetricSetVerbosity(dm, verbosity));
70: PetscCall(PetscOptionsReal("-dm_plex_metric_h_min", "Minimum tolerated metric magnitude", "DMPlexMetricEnforceSPD", h_min, &h_min, NULL));
71: PetscCall(DMPlexMetricSetMinimumMagnitude(dm, h_min));
72: PetscCall(PetscOptionsReal("-dm_plex_metric_h_max", "Maximum tolerated metric magnitude", "DMPlexMetricEnforceSPD", h_max, &h_max, NULL));
73: PetscCall(DMPlexMetricSetMaximumMagnitude(dm, h_max));
74: PetscCall(PetscOptionsReal("-dm_plex_metric_a_max", "Maximum tolerated anisotropy", "DMPlexMetricEnforceSPD", a_max, &a_max, NULL));
75: PetscCall(DMPlexMetricSetMaximumAnisotropy(dm, a_max));
76: PetscCall(PetscOptionsReal("-dm_plex_metric_p", "L-p normalization order", "DMPlexMetricNormalize", p, &p, NULL));
77: PetscCall(DMPlexMetricSetNormalizationOrder(dm, p));
78: PetscCall(PetscOptionsReal("-dm_plex_metric_target_complexity", "Target metric complexity", "DMPlexMetricNormalize", target, &target, NULL));
79: PetscCall(DMPlexMetricSetTargetComplexity(dm, target));
80: PetscCall(PetscOptionsReal("-dm_plex_metric_gradation_factor", "Metric gradation factor", "DMAdaptMetric", beta, &beta, NULL));
81: PetscCall(DMPlexMetricSetGradationFactor(dm, beta));
82: PetscCall(PetscOptionsReal("-dm_plex_metric_hausdorff_number", "Metric Hausdorff number", "DMAdaptMetric", hausd, &hausd, NULL));
83: PetscCall(DMPlexMetricSetHausdorffNumber(dm, hausd));
84: PetscOptionsEnd();
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
88: /*@
89: DMPlexMetricSetIsotropic - Record whether a metric is isotropic
91: Input Parameters:
92: + dm - The `DM`
93: - isotropic - Is the metric isotropic?
95: Level: beginner
97: .seealso: `DMPLEX`, `DMPlexMetricIsIsotropic()`, `DMPlexMetricSetUniform()`, `DMPlexMetricSetRestrictAnisotropyFirst()`
98: @*/
99: PetscErrorCode DMPlexMetricSetIsotropic(DM dm, PetscBool isotropic)
100: {
101: DM_Plex *plex = (DM_Plex *)dm->data;
103: PetscFunctionBegin;
104: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
105: plex->metricCtx->isotropic = isotropic;
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: /*@
110: DMPlexMetricIsIsotropic - Is a metric isotropic?
112: Input Parameters:
113: . dm - The `DM`
115: Output Parameters:
116: . isotropic - Is the metric isotropic?
118: Level: beginner
120: .seealso: `DMPLEX`, `DMPlexMetricSetIsotropic()`, `DMPlexMetricIsUniform()`, `DMPlexMetricRestrictAnisotropyFirst()`
121: @*/
122: PetscErrorCode DMPlexMetricIsIsotropic(DM dm, PetscBool *isotropic)
123: {
124: DM_Plex *plex = (DM_Plex *)dm->data;
126: PetscFunctionBegin;
127: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
128: *isotropic = plex->metricCtx->isotropic;
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: /*@
133: DMPlexMetricSetUniform - Record whether a metric is uniform
135: Input Parameters:
136: + dm - The `DM`
137: - uniform - Is the metric uniform?
139: Level: beginner
141: Note:
142: If the metric is specified as uniform then it is assumed to be isotropic, too.
144: .seealso: `DMPLEX`, `DMPlexMetricIsUniform()`, `DMPlexMetricSetIsotropic()`, `DMPlexMetricSetRestrictAnisotropyFirst()`
145: @*/
146: PetscErrorCode DMPlexMetricSetUniform(DM dm, PetscBool uniform)
147: {
148: DM_Plex *plex = (DM_Plex *)dm->data;
150: PetscFunctionBegin;
151: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
152: plex->metricCtx->uniform = uniform;
153: if (uniform) plex->metricCtx->isotropic = uniform;
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: /*@
158: DMPlexMetricIsUniform - Is a metric uniform?
160: Input Parameters:
161: . dm - The `DM`
163: Output Parameters:
164: . uniform - Is the metric uniform?
166: Level: beginner
168: .seealso: `DMPLEX`, `DMPlexMetricSetUniform()`, `DMPlexMetricIsIsotropic()`, `DMPlexMetricRestrictAnisotropyFirst()`
169: @*/
170: PetscErrorCode DMPlexMetricIsUniform(DM dm, PetscBool *uniform)
171: {
172: DM_Plex *plex = (DM_Plex *)dm->data;
174: PetscFunctionBegin;
175: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
176: *uniform = plex->metricCtx->uniform;
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: /*@
181: DMPlexMetricSetRestrictAnisotropyFirst - Record whether anisotropy should be restricted before normalization
183: Input Parameters:
184: + dm - The `DM`
185: - restrictAnisotropyFirst - Should anisotropy be normalized first?
187: Level: beginner
189: .seealso: `DMPLEX`, `DMPlexMetricSetIsotropic()`, `DMPlexMetricRestrictAnisotropyFirst()`
190: @*/
191: PetscErrorCode DMPlexMetricSetRestrictAnisotropyFirst(DM dm, PetscBool restrictAnisotropyFirst)
192: {
193: DM_Plex *plex = (DM_Plex *)dm->data;
195: PetscFunctionBegin;
196: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
197: plex->metricCtx->restrictAnisotropyFirst = restrictAnisotropyFirst;
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: /*@
202: DMPlexMetricRestrictAnisotropyFirst - Is anisotropy restricted before normalization or after?
204: Input Parameters:
205: . dm - The `DM`
207: Output Parameters:
208: . restrictAnisotropyFirst - Is anisotropy be normalized first?
210: Level: beginner
212: .seealso: `DMPLEX`, `DMPlexMetricIsIsotropic()`, `DMPlexMetricSetRestrictAnisotropyFirst()`
213: @*/
214: PetscErrorCode DMPlexMetricRestrictAnisotropyFirst(DM dm, PetscBool *restrictAnisotropyFirst)
215: {
216: DM_Plex *plex = (DM_Plex *)dm->data;
218: PetscFunctionBegin;
219: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
220: *restrictAnisotropyFirst = plex->metricCtx->restrictAnisotropyFirst;
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: /*@
225: DMPlexMetricSetNoInsertion - Should node insertion and deletion be turned off?
227: Input Parameters:
228: + dm - The `DM`
229: - noInsert - Should node insertion and deletion be turned off?
231: Level: beginner
233: Note:
234: This is only used by Mmg and ParMmg (not Pragmatic).
236: .seealso: `DMPLEX`, `DMPlexMetricNoInsertion()`, `DMPlexMetricSetNoSwapping()`, `DMPlexMetricSetNoMovement()`, `DMPlexMetricSetNoSurf()`
237: @*/
238: PetscErrorCode DMPlexMetricSetNoInsertion(DM dm, PetscBool noInsert)
239: {
240: DM_Plex *plex = (DM_Plex *)dm->data;
242: PetscFunctionBegin;
243: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
244: plex->metricCtx->noInsert = noInsert;
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: /*@
249: DMPlexMetricNoInsertion - Are node insertion and deletion turned off?
251: Input Parameters:
252: . dm - The `DM`
254: Output Parameters:
255: . noInsert - Are node insertion and deletion turned off?
257: Level: beginner
259: Note:
260: This is only used by Mmg and ParMmg (not Pragmatic).
262: .seealso: `DMPLEX`, `DMPlexMetricSetNoInsertion()`, `DMPlexMetricNoSwapping()`, `DMPlexMetricNoMovement()`, `DMPlexMetricNoSurf()`
263: @*/
264: PetscErrorCode DMPlexMetricNoInsertion(DM dm, PetscBool *noInsert)
265: {
266: DM_Plex *plex = (DM_Plex *)dm->data;
268: PetscFunctionBegin;
269: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
270: *noInsert = plex->metricCtx->noInsert;
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: /*@
275: DMPlexMetricSetNoSwapping - Should facet swapping be turned off?
277: Input Parameters:
278: + dm - The `DM`
279: - noSwap - Should facet swapping be turned off?
281: Level: beginner
283: Note:
284: This is only used by Mmg and ParMmg (not Pragmatic).
286: .seealso: `DMPLEX`, `DMPlexMetricNoSwapping()`, `DMPlexMetricSetNoInsertion()`, `DMPlexMetricSetNoMovement()`, `DMPlexMetricSetNoSurf()`
287: @*/
288: PetscErrorCode DMPlexMetricSetNoSwapping(DM dm, PetscBool noSwap)
289: {
290: DM_Plex *plex = (DM_Plex *)dm->data;
292: PetscFunctionBegin;
293: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
294: plex->metricCtx->noSwap = noSwap;
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: /*@
299: DMPlexMetricNoSwapping - Is facet swapping turned off?
301: Input Parameters:
302: . dm - The `DM`
304: Output Parameters:
305: . noSwap - Is facet swapping turned off?
307: Level: beginner
309: Note:
310: This is only used by Mmg and ParMmg (not Pragmatic).
312: .seealso: `DMPLEX`, `DMPlexMetricSetNoSwapping()`, `DMPlexMetricNoInsertion()`, `DMPlexMetricNoMovement()`, `DMPlexMetricNoSurf()`
313: @*/
314: PetscErrorCode DMPlexMetricNoSwapping(DM dm, PetscBool *noSwap)
315: {
316: DM_Plex *plex = (DM_Plex *)dm->data;
318: PetscFunctionBegin;
319: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
320: *noSwap = plex->metricCtx->noSwap;
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: /*@
325: DMPlexMetricSetNoMovement - Should node movement be turned off?
327: Input Parameters:
328: + dm - The `DM`
329: - noMove - Should node movement be turned off?
331: Level: beginner
333: Note:
334: This is only used by Mmg and ParMmg (not Pragmatic).
336: .seealso: `DMPLEX`, `DMPlexMetricNoMovement()`, `DMPlexMetricSetNoInsertion()`, `DMPlexMetricSetNoSwapping()`, `DMPlexMetricSetNoSurf()`
337: @*/
338: PetscErrorCode DMPlexMetricSetNoMovement(DM dm, PetscBool noMove)
339: {
340: DM_Plex *plex = (DM_Plex *)dm->data;
342: PetscFunctionBegin;
343: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
344: plex->metricCtx->noMove = noMove;
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: /*@
349: DMPlexMetricNoMovement - Is node movement turned off?
351: Input Parameters:
352: . dm - The `DM`
354: Output Parameters:
355: . noMove - Is node movement turned off?
357: Level: beginner
359: Note:
360: This is only used by Mmg and ParMmg (not Pragmatic).
362: .seealso: `DMPLEX`, `DMPlexMetricSetNoMovement()`, `DMPlexMetricNoInsertion()`, `DMPlexMetricNoSwapping()`, `DMPlexMetricNoSurf()`
363: @*/
364: PetscErrorCode DMPlexMetricNoMovement(DM dm, PetscBool *noMove)
365: {
366: DM_Plex *plex = (DM_Plex *)dm->data;
368: PetscFunctionBegin;
369: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
370: *noMove = plex->metricCtx->noMove;
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: /*@
375: DMPlexMetricSetNoSurf - Should surface modification be turned off?
377: Input Parameters:
378: + dm - The `DM`
379: - noSurf - Should surface modification be turned off?
381: Level: beginner
383: Note:
384: This is only used by Mmg and ParMmg (not Pragmatic).
386: .seealso: `DMPLEX`, `DMPlexMetricNoSurf()`, `DMPlexMetricSetNoMovement()`, `DMPlexMetricSetNoInsertion()`, `DMPlexMetricSetNoSwapping()`
387: @*/
388: PetscErrorCode DMPlexMetricSetNoSurf(DM dm, PetscBool noSurf)
389: {
390: DM_Plex *plex = (DM_Plex *)dm->data;
392: PetscFunctionBegin;
393: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
394: plex->metricCtx->noSurf = noSurf;
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: /*@
399: DMPlexMetricNoSurf - Is surface modification turned off?
401: Input Parameters:
402: . dm - The `DM`
404: Output Parameters:
405: . noSurf - Is surface modification turned off?
407: Level: beginner
409: Note:
410: This is only used by Mmg and ParMmg (not Pragmatic).
412: .seealso: `DMPLEX`, `DMPlexMetricSetNoSurf()`, `DMPlexMetricNoMovement()`, `DMPlexMetricNoInsertion()`, `DMPlexMetricNoSwapping()`
413: @*/
414: PetscErrorCode DMPlexMetricNoSurf(DM dm, PetscBool *noSurf)
415: {
416: DM_Plex *plex = (DM_Plex *)dm->data;
418: PetscFunctionBegin;
419: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
420: *noSurf = plex->metricCtx->noSurf;
421: PetscFunctionReturn(PETSC_SUCCESS);
422: }
424: /*@
425: DMPlexMetricSetMinimumMagnitude - Set the minimum tolerated metric magnitude
427: Input Parameters:
428: + dm - The `DM`
429: - h_min - The minimum tolerated metric magnitude
431: Level: beginner
433: .seealso: `DMPLEX`, `DMPlexMetricGetMinimumMagnitude()`, `DMPlexMetricSetMaximumMagnitude()`
434: @*/
435: PetscErrorCode DMPlexMetricSetMinimumMagnitude(DM dm, PetscReal h_min)
436: {
437: DM_Plex *plex = (DM_Plex *)dm->data;
439: PetscFunctionBegin;
440: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
441: PetscCheck(h_min > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Metric magnitudes must be in (0, inf)");
442: plex->metricCtx->h_min = h_min;
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: /*@
447: DMPlexMetricGetMinimumMagnitude - Get the minimum tolerated metric magnitude
449: Input Parameters:
450: . dm - The `DM`
452: Output Parameters:
453: . h_min - The minimum tolerated metric magnitude
455: Level: beginner
457: .seealso: `DMPLEX`, `DMPlexMetricSetMinimumMagnitude()`, `DMPlexMetricGetMaximumMagnitude()`
458: @*/
459: PetscErrorCode DMPlexMetricGetMinimumMagnitude(DM dm, PetscReal *h_min)
460: {
461: DM_Plex *plex = (DM_Plex *)dm->data;
463: PetscFunctionBegin;
464: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
465: *h_min = plex->metricCtx->h_min;
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: /*@
470: DMPlexMetricSetMaximumMagnitude - Set the maximum tolerated metric magnitude
472: Input Parameters:
473: + dm - The `DM`
474: - h_max - The maximum tolerated metric magnitude
476: Level: beginner
478: .seealso: `DMPLEX`, `DMPlexMetricGetMaximumMagnitude()`, `DMPlexMetricSetMinimumMagnitude()`
479: @*/
480: PetscErrorCode DMPlexMetricSetMaximumMagnitude(DM dm, PetscReal h_max)
481: {
482: DM_Plex *plex = (DM_Plex *)dm->data;
484: PetscFunctionBegin;
485: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
486: PetscCheck(h_max > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Metric magnitudes must be in (0, inf)");
487: plex->metricCtx->h_max = h_max;
488: PetscFunctionReturn(PETSC_SUCCESS);
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: `DMPLEX`, `DMPlexMetricSetMaximumMagnitude()`, `DMPlexMetricGetMinimumMagnitude()`
503: @*/
504: PetscErrorCode DMPlexMetricGetMaximumMagnitude(DM dm, PetscReal *h_max)
505: {
506: DM_Plex *plex = (DM_Plex *)dm->data;
508: PetscFunctionBegin;
509: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
510: *h_max = plex->metricCtx->h_max;
511: PetscFunctionReturn(PETSC_SUCCESS);
512: }
514: /*@
515: DMPlexMetricSetMaximumAnisotropy - Set the maximum tolerated metric anisotropy
517: Input Parameters:
518: + dm - The `DM`
519: - a_max - The maximum tolerated metric anisotropy
521: Level: beginner
523: Note:
524: If the value zero is given then anisotropy will not be restricted. Otherwise, it should be at least one.
526: .seealso: `DMPLEX`, `DMPlexMetricGetMaximumAnisotropy()`, `DMPlexMetricSetMaximumMagnitude()`
527: @*/
528: PetscErrorCode DMPlexMetricSetMaximumAnisotropy(DM dm, PetscReal a_max)
529: {
530: DM_Plex *plex = (DM_Plex *)dm->data;
532: PetscFunctionBegin;
533: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
534: PetscCheck(a_max >= 1.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Anisotropy must be in [1, inf)");
535: plex->metricCtx->a_max = a_max;
536: PetscFunctionReturn(PETSC_SUCCESS);
537: }
539: /*@
540: DMPlexMetricGetMaximumAnisotropy - Get the maximum tolerated metric anisotropy
542: Input Parameters:
543: . dm - The `DM`
545: Output Parameters:
546: . a_max - The maximum tolerated metric anisotropy
548: Level: beginner
550: .seealso: `DMPLEX`, `DMPlexMetricSetMaximumAnisotropy()`, `DMPlexMetricGetMaximumMagnitude()`
551: @*/
552: PetscErrorCode DMPlexMetricGetMaximumAnisotropy(DM dm, PetscReal *a_max)
553: {
554: DM_Plex *plex = (DM_Plex *)dm->data;
556: PetscFunctionBegin;
557: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
558: *a_max = plex->metricCtx->a_max;
559: PetscFunctionReturn(PETSC_SUCCESS);
560: }
562: /*@
563: DMPlexMetricSetTargetComplexity - Set the target metric complexity
565: Input Parameters:
566: + dm - The `DM`
567: - targetComplexity - The target metric complexity
569: Level: beginner
571: .seealso: `DMPLEX`, `DMPlexMetricGetTargetComplexity()`, `DMPlexMetricSetNormalizationOrder()`
572: @*/
573: PetscErrorCode DMPlexMetricSetTargetComplexity(DM dm, PetscReal targetComplexity)
574: {
575: DM_Plex *plex = (DM_Plex *)dm->data;
577: PetscFunctionBegin;
578: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
579: PetscCheck(targetComplexity > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Metric complexity must be in (0, inf)");
580: plex->metricCtx->targetComplexity = targetComplexity;
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: /*@
585: DMPlexMetricGetTargetComplexity - Get the target metric complexity
587: Input Parameters:
588: . dm - The `DM`
590: Output Parameters:
591: . targetComplexity - The target metric complexity
593: Level: beginner
595: .seealso: `DMPLEX`, `DMPlexMetricSetTargetComplexity()`, `DMPlexMetricGetNormalizationOrder()`
596: @*/
597: PetscErrorCode DMPlexMetricGetTargetComplexity(DM dm, PetscReal *targetComplexity)
598: {
599: DM_Plex *plex = (DM_Plex *)dm->data;
601: PetscFunctionBegin;
602: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
603: *targetComplexity = plex->metricCtx->targetComplexity;
604: PetscFunctionReturn(PETSC_SUCCESS);
605: }
607: /*@
608: DMPlexMetricSetNormalizationOrder - Set the order p for L-p normalization
610: Input Parameters:
611: + dm - The `DM`
612: - p - The normalization order
614: Level: beginner
616: .seealso: `DMPLEX`, `DMPlexMetricGetNormalizationOrder()`, `DMPlexMetricSetTargetComplexity()`
617: @*/
618: PetscErrorCode DMPlexMetricSetNormalizationOrder(DM dm, PetscReal p)
619: {
620: DM_Plex *plex = (DM_Plex *)dm->data;
622: PetscFunctionBegin;
623: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
624: PetscCheck(p >= 1.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Normalization order must be in [1, inf)");
625: plex->metricCtx->p = p;
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }
629: /*@
630: DMPlexMetricGetNormalizationOrder - Get the order p for L-p normalization
632: Input Parameters:
633: . dm - The `DM`
635: Output Parameters:
636: . p - The normalization order
638: Level: beginner
640: .seealso: `DMPLEX`, `DMPlexMetricSetNormalizationOrder()`, `DMPlexMetricGetTargetComplexity()`
641: @*/
642: PetscErrorCode DMPlexMetricGetNormalizationOrder(DM dm, PetscReal *p)
643: {
644: DM_Plex *plex = (DM_Plex *)dm->data;
646: PetscFunctionBegin;
647: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
648: *p = plex->metricCtx->p;
649: PetscFunctionReturn(PETSC_SUCCESS);
650: }
652: /*@
653: DMPlexMetricSetGradationFactor - Set the metric gradation factor
655: Input Parameters:
656: + dm - The `DM`
657: - beta - The metric gradation factor
659: Level: beginner
661: Notes:
662: The gradation factor is the maximum tolerated length ratio between adjacent edges.
664: Turn off gradation by passing the value -1. Otherwise, pass a positive value.
666: This is only used by Mmg and ParMmg (not Pragmatic).
668: .seealso: `DMPLEX`, `DMPlexMetricGetGradationFactor()`, `DMPlexMetricSetHausdorffNumber()`
669: @*/
670: PetscErrorCode DMPlexMetricSetGradationFactor(DM dm, PetscReal beta)
671: {
672: DM_Plex *plex = (DM_Plex *)dm->data;
674: PetscFunctionBegin;
675: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
676: PetscCheck(beta > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Metric gradation factor must be in (0, inf)");
677: plex->metricCtx->gradationFactor = beta;
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: /*@
682: DMPlexMetricGetGradationFactor - Get the metric gradation factor
684: Input Parameters:
685: . dm - The `DM`
687: Output Parameters:
688: . beta - The metric gradation factor
690: Level: beginner
692: Notes:
693: The gradation factor is the maximum tolerated length ratio between adjacent edges.
695: The value -1 implies that gradation is turned off.
697: This is only used by Mmg and ParMmg (not Pragmatic).
699: .seealso: `DMPLEX`, `DMPlexMetricSetGradationFactor()`, `DMPlexMetricGetHausdorffNumber()`
700: @*/
701: PetscErrorCode DMPlexMetricGetGradationFactor(DM dm, PetscReal *beta)
702: {
703: DM_Plex *plex = (DM_Plex *)dm->data;
705: PetscFunctionBegin;
706: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
707: *beta = plex->metricCtx->gradationFactor;
708: PetscFunctionReturn(PETSC_SUCCESS);
709: }
711: /*@
712: DMPlexMetricSetHausdorffNumber - Set the metric Hausdorff number
714: Input Parameters:
715: + dm - The `DM`
716: - hausd - The metric Hausdorff number
718: Level: beginner
720: Notes:
721: The Hausdorff number imposes the maximal distance between the piecewise linear approximation of the
722: boundary and the reconstructed ideal boundary. Thus, a low Hausdorff parameter leads to refine the
723: high curvature areas. By default, the Hausdorff value is set to 0.01, which is a suitable value for
724: an object of size 1 in each direction. For smaller (resp. larger) objects, you may need to decrease
725: (resp. increase) the Hausdorff parameter. (Taken from
726: https://www.mmgtools.org/mmg-remesher-try-mmg/mmg-remesher-options/mmg-remesher-option-hausd).
728: This is only used by Mmg and ParMmg (not Pragmatic).
730: .seealso: `DMPLEX`, `DMPlexMetricSetGradationFactor()`, `DMPlexMetricGetHausdorffNumber()`
731: @*/
732: PetscErrorCode DMPlexMetricSetHausdorffNumber(DM dm, PetscReal hausd)
733: {
734: DM_Plex *plex = (DM_Plex *)dm->data;
736: PetscFunctionBegin;
737: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
738: PetscCheck(hausd > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Metric Hausdorff number must be in (0, inf)");
739: plex->metricCtx->hausdorffNumber = hausd;
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
743: /*@
744: DMPlexMetricGetHausdorffNumber - Get the metric Hausdorff number
746: Input Parameters:
747: . dm - The `DM`
749: Output Parameters:
750: . hausd - The metric Hausdorff number
752: Level: beginner
754: Notes:
755: The Hausdorff number imposes the maximal distance between the piecewise linear approximation of the
756: boundary and the reconstructed ideal boundary. Thus, a low Hausdorff parameter leads to refine the
757: high curvature areas. By default, the Hausdorff value is set to 0.01, which is a suitable value for
758: an object of size 1 in each direction. For smaller (resp. larger) objects, you may need to decrease
759: (resp. increase) the Hausdorff parameter. (Taken from
760: https://www.mmgtools.org/mmg-remesher-try-mmg/mmg-remesher-options/mmg-remesher-option-hausd).
762: This is only used by Mmg and ParMmg (not Pragmatic).
764: .seealso: `DMPLEX`, `DMPlexMetricGetGradationFactor()`, `DMPlexMetricSetHausdorffNumber()`
765: @*/
766: PetscErrorCode DMPlexMetricGetHausdorffNumber(DM dm, PetscReal *hausd)
767: {
768: DM_Plex *plex = (DM_Plex *)dm->data;
770: PetscFunctionBegin;
771: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
772: *hausd = plex->metricCtx->hausdorffNumber;
773: PetscFunctionReturn(PETSC_SUCCESS);
774: }
776: /*@
777: DMPlexMetricSetVerbosity - Set the verbosity of the mesh adaptation package
779: Input Parameters:
780: + dm - The `DM`
781: - verbosity - The verbosity, where -1 is silent and 10 is maximum
783: Level: beginner
785: Note:
786: This is only used by Mmg and ParMmg (not Pragmatic).
788: .seealso: `DMPLEX`, `DMPlexMetricGetVerbosity()`, `DMPlexMetricSetNumIterations()`
789: @*/
790: PetscErrorCode DMPlexMetricSetVerbosity(DM dm, PetscInt verbosity)
791: {
792: DM_Plex *plex = (DM_Plex *)dm->data;
794: PetscFunctionBegin;
795: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
796: plex->metricCtx->verbosity = verbosity;
797: PetscFunctionReturn(PETSC_SUCCESS);
798: }
800: /*@
801: DMPlexMetricGetVerbosity - Get the verbosity of the mesh adaptation package
803: Input Parameters:
804: . dm - The `DM`
806: Output Parameters:
807: . verbosity - The verbosity, where -1 is silent and 10 is maximum
809: Level: beginner
811: Note:
812: This is only used by Mmg and ParMmg (not Pragmatic).
814: .seealso: `DMPLEX`, `DMPlexMetricSetVerbosity()`, `DMPlexMetricGetNumIterations()`
815: @*/
816: PetscErrorCode DMPlexMetricGetVerbosity(DM dm, PetscInt *verbosity)
817: {
818: DM_Plex *plex = (DM_Plex *)dm->data;
820: PetscFunctionBegin;
821: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
822: *verbosity = plex->metricCtx->verbosity;
823: PetscFunctionReturn(PETSC_SUCCESS);
824: }
826: /*@
827: DMPlexMetricSetNumIterations - Set the number of parallel adaptation iterations
829: Input Parameters:
830: + dm - The `DM`
831: - numIter - the number of parallel adaptation iterations
833: Level: beginner
835: Note:
836: This is only used by ParMmg (not Pragmatic or Mmg).
838: .seealso: `DMPLEX`, `DMPlexMetricSetVerbosity()`, `DMPlexMetricGetNumIterations()`
839: @*/
840: PetscErrorCode DMPlexMetricSetNumIterations(DM dm, PetscInt numIter)
841: {
842: DM_Plex *plex = (DM_Plex *)dm->data;
844: PetscFunctionBegin;
845: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
846: plex->metricCtx->numIter = numIter;
847: PetscFunctionReturn(PETSC_SUCCESS);
848: }
850: /*@
851: DMPlexMetricGetNumIterations - Get the number of parallel adaptation iterations
853: Input Parameters:
854: . dm - The `DM`
856: Output Parameters:
857: . numIter - the number of parallel adaptation iterations
859: Level: beginner
861: Note:
862: This is only used by Mmg and ParMmg (not Pragmatic or Mmg).
864: .seealso: `DMPLEX`, `DMPlexMetricSetNumIterations()`, `DMPlexMetricGetVerbosity()`
865: @*/
866: PetscErrorCode DMPlexMetricGetNumIterations(DM dm, PetscInt *numIter)
867: {
868: DM_Plex *plex = (DM_Plex *)dm->data;
870: PetscFunctionBegin;
871: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
872: *numIter = plex->metricCtx->numIter;
873: PetscFunctionReturn(PETSC_SUCCESS);
874: }
876: static PetscErrorCode DMPlexP1FieldCreate_Private(DM dm, PetscInt f, PetscInt size, Vec *metric)
877: {
878: MPI_Comm comm;
879: PetscFE fe;
880: PetscInt dim;
882: PetscFunctionBegin;
883: /* Extract metadata from dm */
884: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
885: PetscCall(DMGetDimension(dm, &dim));
887: /* Create a P1 field of the requested size */
888: PetscCall(PetscFECreateLagrange(comm, dim, size, PETSC_TRUE, 1, PETSC_DETERMINE, &fe));
889: PetscCall(DMSetField(dm, f, NULL, (PetscObject)fe));
890: PetscCall(DMCreateDS(dm));
891: PetscCall(PetscFEDestroy(&fe));
892: PetscCall(DMCreateLocalVector(dm, metric));
893: PetscFunctionReturn(PETSC_SUCCESS);
894: }
896: /*@
897: DMPlexMetricCreate - Create a Riemannian metric field
899: Input Parameters:
900: + dm - The `DM`
901: - f - The field number to use
903: Output Parameter:
904: . metric - The metric
906: Options Database Key:
907: . -dm_adaptor (pragmatic|mmg|parmmg) - specify `DMAdapterType` to use
909: Options Database Keys for Mmg and ParMmg:
910: + -dm_plex_metric_gradation_factor - Maximum ratio by which edge lengths may grow during gradation
911: . -dm_plex_metric_num_iterations - Number of parallel mesh adaptation iterations for ParMmg
912: . -dm_plex_metric_no_insert - Should node insertion/deletion be turned off?
913: . -dm_plex_metric_no_swap - Should facet swapping be turned off?
914: . -dm_plex_metric_no_move - Should node movement be turned off?
915: - -dm_plex_metric_verbosity - Choose a verbosity level from -1 (silent) to 10 (maximum).
917: Options Database Keys for Riemannian metrics:
918: + -dm_plex_metric_isotropic - Is the metric isotropic?
919: . -dm_plex_metric_uniform - Is the metric uniform?
920: . -dm_plex_metric_restrict_anisotropy_first - Should anisotropy be restricted before normalization?
921: . -dm_plex_metric_h_min - Minimum tolerated metric magnitude
922: . -dm_plex_metric_h_max - Maximum tolerated metric magnitude
923: . -dm_plex_metric_a_max - Maximum tolerated anisotropy
924: . -dm_plex_metric_p - L-p normalization order
925: - -dm_plex_metric_target_complexity - Target metric complexity
927: Level: beginner
929: Note:
930: It is assumed that the `DM` is comprised of simplices.
932: .seealso: `DMPLEX`, `DMPlexMetricCreateUniform()`, `DMPlexMetricCreateIsotropic()`
933: @*/
934: PetscErrorCode DMPlexMetricCreate(DM dm, PetscInt f, Vec *metric)
935: {
936: PetscBool isotropic, uniform;
937: PetscInt coordDim, Nd;
939: PetscFunctionBegin;
940: PetscCall(DMGetCoordinateDim(dm, &coordDim));
941: Nd = coordDim * coordDim;
942: PetscCall(DMPlexMetricIsUniform(dm, &uniform));
943: PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic));
944: if (uniform) {
945: MPI_Comm comm;
947: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
948: PetscCheck(isotropic, comm, PETSC_ERR_SUP, "Uniform anisotropic metrics not supported");
949: PetscCall(VecCreate(comm, metric));
950: PetscCall(VecSetSizes(*metric, 1, PETSC_DECIDE));
951: PetscCall(VecSetFromOptions(*metric));
952: } else if (isotropic) PetscCall(DMPlexP1FieldCreate_Private(dm, f, 1, metric));
953: else PetscCall(DMPlexP1FieldCreate_Private(dm, f, Nd, metric));
954: PetscFunctionReturn(PETSC_SUCCESS);
955: }
957: /*@
958: DMPlexMetricCreateUniform - Construct a uniform isotropic metric
960: Input Parameters:
961: + dm - The `DM`
962: . f - The field number to use
963: - alpha - Scaling parameter for the diagonal
965: Output Parameter:
966: . metric - The uniform metric
968: Level: beginner
970: Note:
971: In this case, the metric is constant in space and so is only specified for a single vertex.
973: .seealso: `DMPLEX`, `DMPlexMetricCreate()`, `DMPlexMetricCreateIsotropic()`
974: @*/
975: PetscErrorCode DMPlexMetricCreateUniform(DM dm, PetscInt f, PetscReal alpha, Vec *metric)
976: {
977: PetscFunctionBegin;
978: PetscCall(DMPlexMetricSetUniform(dm, PETSC_TRUE));
979: PetscCall(DMPlexMetricCreate(dm, f, metric));
980: PetscCheck(alpha, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Uniform metric scaling is undefined");
981: PetscCheck(alpha > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Uniform metric scaling should be in (0, inf)");
982: PetscCall(VecSet(*metric, alpha));
983: PetscCall(VecAssemblyBegin(*metric));
984: PetscCall(VecAssemblyEnd(*metric));
985: PetscFunctionReturn(PETSC_SUCCESS);
986: }
988: 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[])
989: {
990: f0[0] = u[0];
991: }
993: /*@
994: DMPlexMetricCreateIsotropic - Construct an isotropic metric from an error indicator
996: Input Parameters:
997: + dm - The `DM`
998: . f - The field number to use
999: - indicator - The error indicator
1001: Output Parameter:
1002: . metric - The isotropic metric
1004: Level: beginner
1006: Notes:
1007: It is assumed that the `DM` is comprised of simplices.
1009: The indicator needs to be a scalar field. If it is not defined vertex-wise, then it is projected appropriately.
1011: .seealso: `DMPLEX`, `DMPlexMetricCreate()`, `DMPlexMetricCreateUniform()`
1012: @*/
1013: PetscErrorCode DMPlexMetricCreateIsotropic(DM dm, PetscInt f, Vec indicator, Vec *metric)
1014: {
1015: PetscInt m, n;
1017: PetscFunctionBegin;
1018: PetscCall(DMPlexMetricSetIsotropic(dm, PETSC_TRUE));
1019: PetscCall(DMPlexMetricCreate(dm, f, metric));
1020: PetscCall(VecGetSize(indicator, &m));
1021: PetscCall(VecGetSize(*metric, &n));
1022: if (m == n) PetscCall(VecCopy(indicator, *metric));
1023: else {
1024: DM dmIndi;
1025: 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[]);
1027: PetscCall(VecGetDM(indicator, &dmIndi));
1028: funcs[0] = identity;
1029: PetscCall(DMProjectFieldLocal(dmIndi, 0.0, indicator, funcs, INSERT_VALUES, *metric));
1030: }
1031: PetscFunctionReturn(PETSC_SUCCESS);
1032: }
1034: /*@
1035: DMPlexMetricDeterminantCreate - Create the determinant field for a Riemannian metric
1037: Input Parameters:
1038: + dm - The `DM` of the metric field
1039: - f - The field number to use
1041: Output Parameters:
1042: + determinant - The determinant field
1043: - dmDet - The corresponding `DM`
1045: Level: beginner
1047: .seealso: `DMPLEX`, `DMPlexMetricCreateUniform()`, `DMPlexMetricCreateIsotropic()`, `DMPlexMetricCreate()`
1048: @*/
1049: PetscErrorCode DMPlexMetricDeterminantCreate(DM dm, PetscInt f, Vec *determinant, DM *dmDet)
1050: {
1051: PetscBool uniform;
1053: PetscFunctionBegin;
1054: PetscCall(DMPlexMetricIsUniform(dm, &uniform));
1055: PetscCall(DMClone(dm, dmDet));
1056: if (uniform) {
1057: MPI_Comm comm;
1059: PetscCall(PetscObjectGetComm((PetscObject)*dmDet, &comm));
1060: PetscCall(VecCreate(comm, determinant));
1061: PetscCall(VecSetSizes(*determinant, 1, PETSC_DECIDE));
1062: PetscCall(VecSetFromOptions(*determinant));
1063: } else PetscCall(DMPlexP1FieldCreate_Private(*dmDet, f, 1, determinant));
1064: PetscFunctionReturn(PETSC_SUCCESS);
1065: }
1067: static PetscErrorCode LAPACKsyevFail(PetscInt dim, PetscScalar Mpos[])
1068: {
1069: PetscInt i, j;
1071: PetscFunctionBegin;
1072: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Failed to apply LAPACKsyev to the matrix\n"));
1073: for (i = 0; i < dim; ++i) {
1074: if (i == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " [["));
1075: else PetscCall(PetscPrintf(PETSC_COMM_SELF, " ["));
1076: for (j = 0; j < dim; ++j) {
1077: if (j < dim - 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%15.8e, ", (double)PetscAbsScalar(Mpos[i * dim + j])));
1078: else PetscCall(PetscPrintf(PETSC_COMM_SELF, "%15.8e", (double)PetscAbsScalar(Mpos[i * dim + j])));
1079: }
1080: if (i < dim - 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, "]\n"));
1081: else PetscCall(PetscPrintf(PETSC_COMM_SELF, "]]\n"));
1082: }
1083: PetscFunctionReturn(PETSC_SUCCESS);
1084: }
1086: static PetscErrorCode DMPlexMetricModify_Private(PetscInt dim, PetscReal h_min, PetscReal h_max, PetscReal a_max, PetscScalar Mp[], PetscScalar *detMp)
1087: {
1088: PetscInt i, j, k;
1089: 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);
1090: PetscScalar *Mpos;
1092: PetscFunctionBegin;
1093: PetscCall(PetscMalloc2(dim * dim, &Mpos, dim, &eigs));
1095: /* Symmetrize */
1096: for (i = 0; i < dim; ++i) {
1097: Mpos[i * dim + i] = Mp[i * dim + i];
1098: for (j = i + 1; j < dim; ++j) {
1099: Mpos[i * dim + j] = 0.5 * (Mp[i * dim + j] + Mp[j * dim + i]);
1100: Mpos[j * dim + i] = Mpos[i * dim + j];
1101: }
1102: }
1104: /* Compute eigendecomposition */
1105: if (dim == 1) {
1106: /* Isotropic case */
1107: eigs[0] = PetscRealPart(Mpos[0]);
1108: Mpos[0] = 1.0;
1109: } else {
1110: /* Anisotropic case */
1111: PetscScalar *work;
1112: PetscBLASInt lwork;
1114: PetscCall(PetscBLASIntCast(5 * dim, &lwork));
1115: PetscCall(PetscMalloc1(5 * dim, &work));
1116: {
1117: PetscBLASInt lierr;
1118: PetscBLASInt nb;
1120: PetscCall(PetscBLASIntCast(dim, &nb));
1121: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1122: #if defined(PETSC_USE_COMPLEX)
1123: {
1124: PetscReal *rwork;
1125: PetscCall(PetscMalloc1(3 * dim, &rwork));
1126: PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, Mpos, &nb, eigs, work, &lwork, rwork, &lierr));
1127: PetscCall(PetscFree(rwork));
1128: }
1129: #else
1130: PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, Mpos, &nb, eigs, work, &lwork, &lierr));
1131: #endif
1132: if (lierr) {
1133: for (i = 0; i < dim; ++i) {
1134: Mpos[i * dim + i] = Mp[i * dim + i];
1135: for (j = i + 1; j < dim; ++j) {
1136: Mpos[i * dim + j] = 0.5 * (Mp[i * dim + j] + Mp[j * dim + i]);
1137: Mpos[j * dim + i] = Mpos[i * dim + j];
1138: }
1139: }
1140: PetscCall(LAPACKsyevFail(dim, Mpos));
1141: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %" PetscBLASInt_FMT, lierr);
1142: }
1143: PetscCall(PetscFPTrapPop());
1144: }
1145: PetscCall(PetscFree(work));
1146: }
1148: /* Reflect to positive orthant and enforce maximum and minimum size */
1149: max_eig = 0.0;
1150: for (i = 0; i < dim; ++i) {
1151: eigs[i] = PetscMin(l_max, PetscMax(l_min, PetscAbsReal(eigs[i])));
1152: max_eig = PetscMax(eigs[i], max_eig);
1153: }
1155: /* Enforce maximum anisotropy and compute determinant */
1156: *detMp = 1.0;
1157: for (i = 0; i < dim; ++i) {
1158: if (a_max >= 1.0) eigs[i] = PetscMax(eigs[i], max_eig * la_min);
1159: *detMp *= eigs[i];
1160: }
1162: /* Reconstruct Hessian */
1163: for (i = 0; i < dim; ++i) {
1164: for (j = 0; j < dim; ++j) {
1165: Mp[i * dim + j] = 0.0;
1166: for (k = 0; k < dim; ++k) Mp[i * dim + j] += Mpos[k * dim + i] * eigs[k] * Mpos[k * dim + j];
1167: }
1168: }
1169: PetscCall(PetscFree2(Mpos, eigs));
1170: PetscFunctionReturn(PETSC_SUCCESS);
1171: }
1173: /*@
1174: DMPlexMetricEnforceSPD - Enforce symmetric positive-definiteness of a metric
1176: Input Parameters:
1177: + dm - The `DM`
1178: . metricIn - The metric
1179: . restrictSizes - Should maximum/minimum metric magnitudes be enforced?
1180: - restrictAnisotropy - Should maximum anisotropy be enforced?
1182: Output Parameters:
1183: + metricOut - The metric
1184: - determinant - Its determinant
1186: Options Database Keys:
1187: + -dm_plex_metric_isotropic - Is the metric isotropic?
1188: . -dm_plex_metric_uniform - Is the metric uniform?
1189: . -dm_plex_metric_h_min - Minimum tolerated metric magnitude
1190: . -dm_plex_metric_h_max - Maximum tolerated metric magnitude
1191: - -dm_plex_metric_a_max - Maximum tolerated anisotropy
1193: Level: beginner
1195: .seealso: `DMPLEX`, `DMPlexMetricNormalize()`, `DMPlexMetricIntersection()`
1196: @*/
1197: PetscErrorCode DMPlexMetricEnforceSPD(DM dm, Vec metricIn, PetscBool restrictSizes, PetscBool restrictAnisotropy, Vec metricOut, Vec determinant)
1198: {
1199: DM dmDet;
1200: PetscBool isotropic, uniform;
1201: PetscInt dim, vStart, vEnd, v;
1202: PetscScalar *met, *det;
1203: PetscReal h_min = 1.0e-30, h_max = 1.0e+30, a_max = 1.0e+15;
1205: PetscFunctionBegin;
1206: PetscCall(PetscLogEventBegin(DMPLEX_MetricEnforceSPD, 0, 0, 0, 0));
1208: /* Extract metadata from dm */
1209: PetscCall(DMGetDimension(dm, &dim));
1210: if (restrictSizes) {
1211: PetscCall(DMPlexMetricGetMinimumMagnitude(dm, &h_min));
1212: PetscCall(DMPlexMetricGetMaximumMagnitude(dm, &h_max));
1213: h_min = PetscMax(h_min, 1.0e-30);
1214: h_max = PetscMin(h_max, 1.0e+30);
1215: PetscCheck(h_min < h_max, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Minimum metric magnitude should be smaller than maximum metric magnitude");
1216: }
1217: if (restrictAnisotropy) {
1218: PetscCall(DMPlexMetricGetMaximumAnisotropy(dm, &a_max));
1219: a_max = PetscMin(a_max, 1.0e+30);
1220: }
1222: /* Setup output metric */
1223: PetscCall(VecCopy(metricIn, metricOut));
1225: /* Enforce SPD and extract determinant */
1226: PetscCall(VecGetArray(metricOut, &met));
1227: PetscCall(DMPlexMetricIsUniform(dm, &uniform));
1228: PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic));
1229: if (uniform) {
1230: PetscCheck(isotropic, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Uniform anisotropic metrics cannot exist");
1232: /* Uniform case */
1233: PetscCall(VecGetArray(determinant, &det));
1234: PetscCall(DMPlexMetricModify_Private(1, h_min, h_max, a_max, met, det));
1235: PetscCall(VecRestoreArray(determinant, &det));
1236: } else {
1237: /* Spatially varying case */
1238: PetscInt nrow;
1240: if (isotropic) nrow = 1;
1241: else nrow = dim;
1242: PetscCall(VecGetDM(determinant, &dmDet));
1243: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1244: PetscCall(VecGetArray(determinant, &det));
1245: for (v = vStart; v < vEnd; ++v) {
1246: PetscScalar *vmet, *vdet;
1247: PetscCall(DMPlexPointLocalRef(dm, v, met, &vmet));
1248: PetscCall(DMPlexPointLocalRef(dmDet, v, det, &vdet));
1249: PetscCall(DMPlexMetricModify_Private(nrow, h_min, h_max, a_max, vmet, vdet));
1250: }
1251: PetscCall(VecRestoreArray(determinant, &det));
1252: }
1253: PetscCall(VecRestoreArray(metricOut, &met));
1255: PetscCall(PetscLogEventEnd(DMPLEX_MetricEnforceSPD, 0, 0, 0, 0));
1256: PetscFunctionReturn(PETSC_SUCCESS);
1257: }
1259: 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[])
1260: {
1261: const PetscScalar p = constants[0];
1263: f0[0] = PetscPowScalar(u[0], p / (2.0 * p + dim));
1264: }
1266: /*@
1267: DMPlexMetricNormalize - Apply L-p normalization to a metric
1269: Input Parameters:
1270: + dm - The `DM`
1271: . metricIn - The unnormalized metric
1272: . restrictSizes - Should maximum/minimum metric magnitudes be enforced?
1273: - restrictAnisotropy - Should maximum metric anisotropy be enforced?
1275: Output Parameters:
1276: + metricOut - The normalized metric
1277: - determinant - computed determinant
1279: Options Database Keys:
1280: + -dm_plex_metric_isotropic - Is the metric isotropic?
1281: . -dm_plex_metric_uniform - Is the metric uniform?
1282: . -dm_plex_metric_restrict_anisotropy_first - Should anisotropy be restricted before normalization?
1283: . -dm_plex_metric_h_min - Minimum tolerated metric magnitude
1284: . -dm_plex_metric_h_max - Maximum tolerated metric magnitude
1285: . -dm_plex_metric_a_max - Maximum tolerated anisotropy
1286: . -dm_plex_metric_p - L-p normalization order
1287: - -dm_plex_metric_target_complexity - Target metric complexity
1289: Level: beginner
1291: .seealso: `DMPLEX`, `DMPlexMetricEnforceSPD()`, `DMPlexMetricIntersection()`
1292: @*/
1293: PetscErrorCode DMPlexMetricNormalize(DM dm, Vec metricIn, PetscBool restrictSizes, PetscBool restrictAnisotropy, Vec metricOut, Vec determinant)
1294: {
1295: DM dmDet;
1296: MPI_Comm comm;
1297: PetscBool restrictAnisotropyFirst, isotropic, uniform;
1298: PetscDS ds;
1299: PetscInt dim, Nd, vStart, vEnd, v, i;
1300: PetscScalar *met, *det, integral, constants[1];
1301: PetscReal p, h_min = 1.0e-30, h_max = 1.0e+30, a_max = 0.0, factGlob, fact, target, realIntegral;
1303: PetscFunctionBegin;
1304: PetscCall(PetscLogEventBegin(DMPLEX_MetricNormalize, 0, 0, 0, 0));
1306: /* Extract metadata from dm */
1307: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1308: PetscCall(DMGetDimension(dm, &dim));
1309: PetscCall(DMPlexMetricIsUniform(dm, &uniform));
1310: PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic));
1311: if (isotropic) Nd = 1;
1312: else Nd = dim * dim;
1314: /* Set up metric and ensure it is SPD */
1315: PetscCall(DMPlexMetricRestrictAnisotropyFirst(dm, &restrictAnisotropyFirst));
1316: PetscCall(DMPlexMetricEnforceSPD(dm, metricIn, PETSC_FALSE, (PetscBool)(restrictAnisotropy && restrictAnisotropyFirst), metricOut, determinant));
1318: /* Compute global normalization factor */
1319: PetscCall(DMPlexMetricGetTargetComplexity(dm, &target));
1320: PetscCall(DMPlexMetricGetNormalizationOrder(dm, &p));
1321: constants[0] = p;
1322: if (uniform) {
1323: PetscCheck(isotropic, comm, PETSC_ERR_SUP, "Uniform anisotropic metrics not supported");
1324: DM dmTmp;
1325: Vec tmp;
1327: PetscCall(DMClone(dm, &dmTmp));
1328: PetscCall(DMPlexP1FieldCreate_Private(dmTmp, 0, 1, &tmp));
1329: PetscCall(VecGetArray(determinant, &det));
1330: PetscCall(VecSet(tmp, det[0]));
1331: PetscCall(VecRestoreArray(determinant, &det));
1332: PetscCall(DMGetDS(dmTmp, &ds));
1333: PetscCall(PetscDSSetConstants(ds, 1, constants));
1334: PetscCall(PetscDSSetObjective(ds, 0, detMFunc));
1335: PetscCall(DMPlexComputeIntegralFEM(dmTmp, tmp, &integral, NULL));
1336: PetscCall(VecDestroy(&tmp));
1337: PetscCall(DMDestroy(&dmTmp));
1338: } else {
1339: PetscCall(VecGetDM(determinant, &dmDet));
1340: PetscCall(DMGetDS(dmDet, &ds));
1341: PetscCall(PetscDSSetConstants(ds, 1, constants));
1342: PetscCall(PetscDSSetObjective(ds, 0, detMFunc));
1343: PetscCall(DMPlexComputeIntegralFEM(dmDet, determinant, &integral, NULL));
1344: }
1345: realIntegral = PetscRealPart(integral);
1346: PetscCheck(realIntegral > 1.0e-30, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Global metric normalization factor must be in (0, inf). Is the input metric positive-definite?");
1347: factGlob = PetscPowReal(target / realIntegral, 2.0 / dim);
1349: /* Apply local scaling */
1350: if (restrictSizes) {
1351: PetscCall(DMPlexMetricGetMinimumMagnitude(dm, &h_min));
1352: PetscCall(DMPlexMetricGetMaximumMagnitude(dm, &h_max));
1353: h_min = PetscMax(h_min, 1.0e-30);
1354: h_max = PetscMin(h_max, 1.0e+30);
1355: PetscCheck(h_min < h_max, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Minimum metric magnitude should be smaller than maximum metric magnitude");
1356: }
1357: if (restrictAnisotropy && !restrictAnisotropyFirst) {
1358: PetscCall(DMPlexMetricGetMaximumAnisotropy(dm, &a_max));
1359: a_max = PetscMin(a_max, 1.0e+30);
1360: }
1361: PetscCall(VecGetArray(metricOut, &met));
1362: PetscCall(VecGetArray(determinant, &det));
1363: if (uniform) {
1364: /* Uniform case */
1365: met[0] *= factGlob * PetscPowReal(PetscAbsScalar(det[0]), -1.0 / (2 * p + dim));
1366: if (restrictSizes) PetscCall(DMPlexMetricModify_Private(1, h_min, h_max, a_max, met, det));
1367: } else {
1368: /* Spatially varying case */
1369: PetscInt nrow;
1371: if (isotropic) nrow = 1;
1372: else nrow = dim;
1373: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1374: PetscCall(VecGetDM(determinant, &dmDet));
1375: for (v = vStart; v < vEnd; ++v) {
1376: PetscScalar *Mp, *detM;
1378: PetscCall(DMPlexPointLocalRef(dm, v, met, &Mp));
1379: PetscCall(DMPlexPointLocalRef(dmDet, v, det, &detM));
1380: fact = factGlob * PetscPowReal(PetscAbsScalar(detM[0]), -1.0 / (2 * p + dim));
1381: for (i = 0; i < Nd; ++i) Mp[i] *= fact;
1382: if (restrictSizes) PetscCall(DMPlexMetricModify_Private(nrow, h_min, h_max, a_max, Mp, detM));
1383: }
1384: }
1385: PetscCall(VecRestoreArray(determinant, &det));
1386: PetscCall(VecRestoreArray(metricOut, &met));
1388: PetscCall(PetscLogEventEnd(DMPLEX_MetricNormalize, 0, 0, 0, 0));
1389: PetscFunctionReturn(PETSC_SUCCESS);
1390: }
1392: /*@
1393: DMPlexMetricAverage - Compute the average of a list of metrics
1395: Input Parameters:
1396: + dm - The `DM`
1397: . numMetrics - The number of metrics to be averaged
1398: . weights - Weights for the average
1399: - metrics - The metrics to be averaged
1401: Output Parameter:
1402: . metricAvg - The averaged metric
1404: Level: beginner
1406: Notes:
1407: The weights should sum to unity.
1409: If weights are not provided then an unweighted average is used.
1411: .seealso: `DMPLEX`, `DMPlexMetricAverage2()`, `DMPlexMetricAverage3()`, `DMPlexMetricIntersection()`
1412: @*/
1413: PetscErrorCode DMPlexMetricAverage(DM dm, PetscInt numMetrics, PetscReal weights[], Vec metrics[], Vec metricAvg)
1414: {
1415: PetscBool haveWeights = PETSC_TRUE;
1416: PetscInt i, m, n;
1417: PetscReal sum = 0.0, tol = 1.0e-10;
1419: PetscFunctionBegin;
1420: PetscCall(PetscLogEventBegin(DMPLEX_MetricAverage, 0, 0, 0, 0));
1421: PetscCheck(numMetrics >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot average %" PetscInt_FMT " < 1 metrics", numMetrics);
1422: PetscCall(VecSet(metricAvg, 0.0));
1423: PetscCall(VecGetSize(metricAvg, &m));
1424: for (i = 0; i < numMetrics; ++i) {
1425: PetscCall(VecGetSize(metrics[i], &n));
1426: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Averaging different metric types not implemented");
1427: }
1429: /* Default to the unweighted case */
1430: if (!weights) {
1431: PetscCall(PetscMalloc1(numMetrics, &weights));
1432: haveWeights = PETSC_FALSE;
1433: for (i = 0; i < numMetrics; ++i) weights[i] = 1.0 / numMetrics;
1434: }
1436: /* Check weights sum to unity */
1437: for (i = 0; i < numMetrics; ++i) sum += weights[i];
1438: PetscCheck(PetscAbsReal(sum - 1) <= tol, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Weights do not sum to unity");
1440: /* Compute metric average */
1441: for (i = 0; i < numMetrics; ++i) PetscCall(VecAXPY(metricAvg, weights[i], metrics[i]));
1442: if (!haveWeights) PetscCall(PetscFree(weights));
1444: PetscCall(PetscLogEventEnd(DMPLEX_MetricAverage, 0, 0, 0, 0));
1445: PetscFunctionReturn(PETSC_SUCCESS);
1446: }
1448: /*@
1449: DMPlexMetricAverage2 - Compute the unweighted average of two metrics
1451: Input Parameters:
1452: + dm - The `DM`
1453: . metric1 - The first metric to be averaged
1454: - metric2 - The second metric to be averaged
1456: Output Parameter:
1457: . metricAvg - The averaged metric
1459: Level: beginner
1461: .seealso: `DMPLEX`, `DMPlexMetricAverage()`, `DMPlexMetricAverage3()`
1462: @*/
1463: PetscErrorCode DMPlexMetricAverage2(DM dm, Vec metric1, Vec metric2, Vec metricAvg)
1464: {
1465: PetscReal weights[2] = {0.5, 0.5};
1466: Vec metrics[2] = {metric1, metric2};
1468: PetscFunctionBegin;
1469: PetscCall(DMPlexMetricAverage(dm, 2, weights, metrics, metricAvg));
1470: PetscFunctionReturn(PETSC_SUCCESS);
1471: }
1473: /*@
1474: DMPlexMetricAverage3 - Compute the unweighted average of three metrics
1476: Input Parameters:
1477: + dm - The `DM`
1478: . metric1 - The first metric to be averaged
1479: . metric2 - The second metric to be averaged
1480: - metric3 - The third metric to be averaged
1482: Output Parameter:
1483: . metricAvg - The averaged metric
1485: Level: beginner
1487: .seealso: `DMPLEX`, `DMPlexMetricAverage()`, `DMPlexMetricAverage2()`
1488: @*/
1489: PetscErrorCode DMPlexMetricAverage3(DM dm, Vec metric1, Vec metric2, Vec metric3, Vec metricAvg)
1490: {
1491: PetscReal weights[3] = {1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0};
1492: Vec metrics[3] = {metric1, metric2, metric3};
1494: PetscFunctionBegin;
1495: PetscCall(DMPlexMetricAverage(dm, 3, weights, metrics, metricAvg));
1496: PetscFunctionReturn(PETSC_SUCCESS);
1497: }
1499: static PetscErrorCode DMPlexMetricIntersection_Private(PetscInt dim, PetscScalar M1[], PetscScalar M2[])
1500: {
1501: PetscInt i, j, k, l, m;
1502: PetscReal *evals;
1503: PetscScalar *evecs, *sqrtM1, *isqrtM1;
1505: PetscFunctionBegin;
1506: /* Isotropic case */
1507: if (dim == 1) {
1508: M2[0] = PetscMax(PetscRealPart(M1[0]), PetscRealPart(M2[0]));
1509: PetscFunctionReturn(PETSC_SUCCESS);
1510: }
1512: /* Anisotropic case */
1513: PetscCall(PetscMalloc4(dim * dim, &evecs, dim * dim, &sqrtM1, dim * dim, &isqrtM1, dim, &evals));
1514: for (i = 0; i < dim; ++i) {
1515: for (j = 0; j < dim; ++j) evecs[i * dim + j] = M1[i * dim + j];
1516: }
1517: {
1518: PetscScalar *work;
1519: PetscBLASInt lwork;
1521: PetscCall(PetscBLASIntCast(5 * dim, &lwork));
1522: PetscCall(PetscMalloc1(5 * dim, &work));
1523: {
1524: PetscBLASInt lierr, nb;
1525: PetscReal sqrtj;
1527: /* Compute eigendecomposition of M1 */
1528: PetscCall(PetscBLASIntCast(dim, &nb));
1529: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1530: #if defined(PETSC_USE_COMPLEX)
1531: {
1532: PetscReal *rwork;
1533: PetscCall(PetscMalloc1(3 * dim, &rwork));
1534: PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals, work, &lwork, rwork, &lierr));
1535: PetscCall(PetscFree(rwork));
1536: }
1537: #else
1538: PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals, work, &lwork, &lierr));
1539: #endif
1540: if (lierr) {
1541: PetscCall(LAPACKsyevFail(dim, M1));
1542: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %" PetscBLASInt_FMT, lierr);
1543: }
1544: PetscCall(PetscFPTrapPop());
1546: /* Compute square root and the reciprocal thereof */
1547: for (i = 0; i < dim; ++i) {
1548: for (k = 0; k < dim; ++k) {
1549: sqrtM1[i * dim + k] = 0.0;
1550: isqrtM1[i * dim + k] = 0.0;
1551: for (j = 0; j < dim; ++j) {
1552: sqrtj = PetscSqrtReal(evals[j]);
1553: sqrtM1[i * dim + k] += evecs[j * dim + i] * sqrtj * evecs[j * dim + k];
1554: isqrtM1[i * dim + k] += evecs[j * dim + i] * (1.0 / sqrtj) * evecs[j * dim + k];
1555: }
1556: }
1557: }
1559: /* Map M2 into the space spanned by the eigenvectors of M1 */
1560: for (i = 0; i < dim; ++i) {
1561: for (l = 0; l < dim; ++l) {
1562: evecs[i * dim + l] = 0.0;
1563: for (j = 0; j < dim; ++j) {
1564: for (k = 0; k < dim; ++k) evecs[i * dim + l] += isqrtM1[j * dim + i] * M2[j * dim + k] * isqrtM1[k * dim + l];
1565: }
1566: }
1567: }
1569: /* Compute eigendecomposition */
1570: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1571: #if defined(PETSC_USE_COMPLEX)
1572: {
1573: PetscReal *rwork;
1574: PetscCall(PetscMalloc1(3 * dim, &rwork));
1575: PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals, work, &lwork, rwork, &lierr));
1576: PetscCall(PetscFree(rwork));
1577: }
1578: #else
1579: PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals, work, &lwork, &lierr));
1580: #endif
1581: if (lierr) {
1582: for (i = 0; i < dim; ++i) {
1583: for (l = 0; l < dim; ++l) {
1584: evecs[i * dim + l] = 0.0;
1585: for (j = 0; j < dim; ++j) {
1586: for (k = 0; k < dim; ++k) evecs[i * dim + l] += isqrtM1[j * dim + i] * M2[j * dim + k] * isqrtM1[k * dim + l];
1587: }
1588: }
1589: }
1590: PetscCall(LAPACKsyevFail(dim, evecs));
1591: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %" PetscBLASInt_FMT, lierr);
1592: }
1593: PetscCall(PetscFPTrapPop());
1595: /* Modify eigenvalues */
1596: for (i = 0; i < dim; ++i) evals[i] = PetscMax(evals[i], 1.0);
1598: /* Map back to get the intersection */
1599: for (i = 0; i < dim; ++i) {
1600: for (m = 0; m < dim; ++m) {
1601: M2[i * dim + m] = 0.0;
1602: for (j = 0; j < dim; ++j) {
1603: for (k = 0; k < dim; ++k) {
1604: for (l = 0; l < dim; ++l) M2[i * dim + m] += sqrtM1[j * dim + i] * evecs[k * dim + j] * evals[k] * evecs[k * dim + l] * sqrtM1[l * dim + m];
1605: }
1606: }
1607: }
1608: }
1609: }
1610: PetscCall(PetscFree(work));
1611: }
1612: PetscCall(PetscFree4(evecs, sqrtM1, isqrtM1, evals));
1613: PetscFunctionReturn(PETSC_SUCCESS);
1614: }
1616: /*@
1617: DMPlexMetricIntersection - Compute the intersection of a list of metrics
1619: Input Parameters:
1620: + dm - The `DM`
1621: . numMetrics - The number of metrics to be intersected
1622: - metrics - The metrics to be intersected
1624: Output Parameter:
1625: . metricInt - The intersected metric
1627: Level: beginner
1629: Notes:
1630: The intersection of a list of metrics has the minimal ellipsoid which fits within the ellipsoids of the component metrics.
1632: The implementation used here is only consistent with the minimal ellipsoid definition in the case numMetrics = 2.
1634: .seealso: `DMPLEX`, `DMPlexMetricIntersection2()`, `DMPlexMetricIntersection3()`, `DMPlexMetricAverage()`
1635: @*/
1636: PetscErrorCode DMPlexMetricIntersection(DM dm, PetscInt numMetrics, Vec metrics[], Vec metricInt)
1637: {
1638: PetscBool isotropic, uniform;
1639: PetscInt v, i, m, n;
1640: PetscScalar *met, *meti;
1642: PetscFunctionBegin;
1643: PetscCall(PetscLogEventBegin(DMPLEX_MetricIntersection, 0, 0, 0, 0));
1644: PetscCheck(numMetrics >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot intersect %" PetscInt_FMT " < 1 metrics", numMetrics);
1646: /* Copy over the first metric */
1647: PetscCall(VecCopy(metrics[0], metricInt));
1648: if (numMetrics == 1) PetscFunctionReturn(PETSC_SUCCESS);
1649: PetscCall(VecGetSize(metricInt, &m));
1650: for (i = 0; i < numMetrics; ++i) {
1651: PetscCall(VecGetSize(metrics[i], &n));
1652: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Intersecting different metric types not implemented");
1653: }
1655: /* Intersect subsequent metrics in turn */
1656: PetscCall(DMPlexMetricIsUniform(dm, &uniform));
1657: PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic));
1658: if (uniform) {
1659: /* Uniform case */
1660: PetscCall(VecGetArray(metricInt, &met));
1661: for (i = 1; i < numMetrics; ++i) {
1662: PetscCall(VecGetArray(metrics[i], &meti));
1663: PetscCall(DMPlexMetricIntersection_Private(1, meti, met));
1664: PetscCall(VecRestoreArray(metrics[i], &meti));
1665: }
1666: PetscCall(VecRestoreArray(metricInt, &met));
1667: } else {
1668: /* Spatially varying case */
1669: PetscInt dim, vStart, vEnd, nrow;
1670: PetscScalar *M, *Mi;
1672: PetscCall(DMGetDimension(dm, &dim));
1673: if (isotropic) nrow = 1;
1674: else nrow = dim;
1675: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1676: PetscCall(VecGetArray(metricInt, &met));
1677: for (i = 1; i < numMetrics; ++i) {
1678: PetscCall(VecGetArray(metrics[i], &meti));
1679: for (v = vStart; v < vEnd; ++v) {
1680: PetscCall(DMPlexPointLocalRef(dm, v, met, &M));
1681: PetscCall(DMPlexPointLocalRef(dm, v, meti, &Mi));
1682: PetscCall(DMPlexMetricIntersection_Private(nrow, Mi, M));
1683: }
1684: PetscCall(VecRestoreArray(metrics[i], &meti));
1685: }
1686: PetscCall(VecRestoreArray(metricInt, &met));
1687: }
1689: PetscCall(PetscLogEventEnd(DMPLEX_MetricIntersection, 0, 0, 0, 0));
1690: PetscFunctionReturn(PETSC_SUCCESS);
1691: }
1693: /*@
1694: DMPlexMetricIntersection2 - Compute the intersection of two metrics
1696: Input Parameters:
1697: + dm - The `DM`
1698: . metric1 - The first metric to be intersected
1699: - metric2 - The second metric to be intersected
1701: Output Parameter:
1702: . metricInt - The intersected metric
1704: Level: beginner
1706: .seealso: `DMPLEX`, `DMPlexMetricIntersection()`, `DMPlexMetricIntersection3()`
1707: @*/
1708: PetscErrorCode DMPlexMetricIntersection2(DM dm, Vec metric1, Vec metric2, Vec metricInt)
1709: {
1710: Vec metrics[2] = {metric1, metric2};
1712: PetscFunctionBegin;
1713: PetscCall(DMPlexMetricIntersection(dm, 2, metrics, metricInt));
1714: PetscFunctionReturn(PETSC_SUCCESS);
1715: }
1717: /*@
1718: DMPlexMetricIntersection3 - Compute the intersection of three metrics
1720: Input Parameters:
1721: + dm - The `DM`
1722: . metric1 - The first metric to be intersected
1723: . metric2 - The second metric to be intersected
1724: - metric3 - The third metric to be intersected
1726: Output Parameter:
1727: . metricInt - The intersected metric
1729: Level: beginner
1731: .seealso: `DMPLEX`, `DMPlexMetricIntersection()`, `DMPlexMetricIntersection2()`
1732: @*/
1733: PetscErrorCode DMPlexMetricIntersection3(DM dm, Vec metric1, Vec metric2, Vec metric3, Vec metricInt)
1734: {
1735: Vec metrics[3] = {metric1, metric2, metric3};
1737: PetscFunctionBegin;
1738: PetscCall(DMPlexMetricIntersection(dm, 3, metrics, metricInt));
1739: PetscFunctionReturn(PETSC_SUCCESS);
1740: }