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: PetscFunctionBegin;
16: if (!plex->metricCtx) PetscCall(PetscNew(&plex->metricCtx));
17: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
18: PetscOptionsBegin(comm, "", "Riemannian metric options", "DMPlexMetric");
19: PetscCall(PetscOptionsBool("-dm_plex_metric_isotropic", "Is the metric isotropic?", "DMPlexMetricCreateIsotropic", isotropic, &isotropic, NULL));
20: PetscCall(DMPlexMetricSetIsotropic(dm, isotropic));
21: PetscCall(PetscOptionsBool("-dm_plex_metric_uniform", "Is the metric uniform?", "DMPlexMetricCreateUniform", uniform, &uniform, NULL));
22: PetscCall(DMPlexMetricSetUniform(dm, uniform));
23: PetscCall(PetscOptionsBool("-dm_plex_metric_restrict_anisotropy_first", "Should anisotropy be restricted before normalization?", "DMPlexNormalize", restrictAnisotropyFirst, &restrictAnisotropyFirst, NULL));
24: PetscCall(DMPlexMetricSetRestrictAnisotropyFirst(dm, restrictAnisotropyFirst));
25: PetscCall(PetscOptionsBool("-dm_plex_metric_no_insert", "Turn off node insertion and deletion", "DMAdaptMetric", noInsert, &noInsert, NULL));
26: PetscCall(DMPlexMetricSetNoInsertion(dm, noInsert));
27: PetscCall(PetscOptionsBool("-dm_plex_metric_no_swap", "Turn off facet swapping", "DMAdaptMetric", noSwap, &noSwap, NULL));
28: PetscCall(DMPlexMetricSetNoSwapping(dm, noSwap));
29: PetscCall(PetscOptionsBool("-dm_plex_metric_no_move", "Turn off facet node movement", "DMAdaptMetric", noMove, &noMove, NULL));
30: PetscCall(DMPlexMetricSetNoMovement(dm, noMove));
31: PetscCall(PetscOptionsBool("-dm_plex_metric_no_surf", "Turn off surface modification", "DMAdaptMetric", noSurf, &noSurf, NULL));
32: PetscCall(DMPlexMetricSetNoSurf(dm, noSurf));
33: PetscCall(PetscOptionsBoundedInt("-dm_plex_metric_num_iterations", "Number of ParMmg adaptation iterations", "DMAdaptMetric", numIter, &numIter, NULL, 0));
34: PetscCall(DMPlexMetricSetNumIterations(dm, numIter));
35: PetscCall(PetscOptionsRangeInt("-dm_plex_metric_verbosity", "Verbosity of metric-based mesh adaptation package (-1 = silent, 10 = maximum)", "DMAdaptMetric", verbosity, &verbosity, NULL, -1, 10));
36: PetscCall(DMPlexMetricSetVerbosity(dm, verbosity));
37: PetscCall(PetscOptionsReal("-dm_plex_metric_h_min", "Minimum tolerated metric magnitude", "DMPlexMetricEnforceSPD", h_min, &h_min, NULL));
38: PetscCall(DMPlexMetricSetMinimumMagnitude(dm, h_min));
39: PetscCall(PetscOptionsReal("-dm_plex_metric_h_max", "Maximum tolerated metric magnitude", "DMPlexMetricEnforceSPD", h_max, &h_max, NULL));
40: PetscCall(DMPlexMetricSetMaximumMagnitude(dm, h_max));
41: PetscCall(PetscOptionsReal("-dm_plex_metric_a_max", "Maximum tolerated anisotropy", "DMPlexMetricEnforceSPD", a_max, &a_max, NULL));
42: PetscCall(DMPlexMetricSetMaximumAnisotropy(dm, a_max));
43: PetscCall(PetscOptionsReal("-dm_plex_metric_p", "L-p normalization order", "DMPlexMetricNormalize", p, &p, NULL));
44: PetscCall(DMPlexMetricSetNormalizationOrder(dm, p));
45: PetscCall(PetscOptionsReal("-dm_plex_metric_target_complexity", "Target metric complexity", "DMPlexMetricNormalize", target, &target, NULL));
46: PetscCall(DMPlexMetricSetTargetComplexity(dm, target));
47: PetscCall(PetscOptionsReal("-dm_plex_metric_gradation_factor", "Metric gradation factor", "DMAdaptMetric", beta, &beta, NULL));
48: PetscCall(DMPlexMetricSetGradationFactor(dm, beta));
49: PetscCall(PetscOptionsReal("-dm_plex_metric_hausdorff_number", "Metric Hausdorff number", "DMAdaptMetric", hausd, &hausd, NULL));
50: PetscCall(DMPlexMetricSetHausdorffNumber(dm, hausd));
51: PetscOptionsEnd();
52: PetscFunctionReturn(PETSC_SUCCESS);
53: }
55: /*@
56: DMPlexMetricSetIsotropic - Record whether a metric is isotropic
58: Input Parameters:
59: + dm - The `DM`
60: - isotropic - Is the metric isotropic?
62: Level: beginner
64: .seealso: `DMPLEX`, `DMPlexMetricIsIsotropic()`, `DMPlexMetricSetUniform()`, `DMPlexMetricSetRestrictAnisotropyFirst()`
65: @*/
66: PetscErrorCode DMPlexMetricSetIsotropic(DM dm, PetscBool isotropic)
67: {
68: DM_Plex *plex = (DM_Plex *)dm->data;
70: PetscFunctionBegin;
71: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
72: plex->metricCtx->isotropic = isotropic;
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: /*@
77: DMPlexMetricIsIsotropic - Is a metric isotropic?
79: Input Parameters:
80: . dm - The `DM`
82: Output Parameters:
83: . isotropic - Is the metric isotropic?
85: Level: beginner
87: .seealso: `DMPLEX`, `DMPlexMetricSetIsotropic()`, `DMPlexMetricIsUniform()`, `DMPlexMetricRestrictAnisotropyFirst()`
88: @*/
89: PetscErrorCode DMPlexMetricIsIsotropic(DM dm, PetscBool *isotropic)
90: {
91: DM_Plex *plex = (DM_Plex *)dm->data;
93: PetscFunctionBegin;
94: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
95: *isotropic = plex->metricCtx->isotropic;
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: /*@
100: DMPlexMetricSetUniform - Record whether a metric is uniform
102: Input Parameters:
103: + dm - The `DM`
104: - uniform - Is the metric uniform?
106: Level: beginner
108: Note:
109: If the metric is specified as uniform then it is assumed to be isotropic, too.
111: .seealso: `DMPLEX`, `DMPlexMetricIsUniform()`, `DMPlexMetricSetIsotropic()`, `DMPlexMetricSetRestrictAnisotropyFirst()`
112: @*/
113: PetscErrorCode DMPlexMetricSetUniform(DM dm, PetscBool uniform)
114: {
115: DM_Plex *plex = (DM_Plex *)dm->data;
117: PetscFunctionBegin;
118: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
119: plex->metricCtx->uniform = uniform;
120: if (uniform) plex->metricCtx->isotropic = uniform;
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
124: /*@
125: DMPlexMetricIsUniform - Is a metric uniform?
127: Input Parameters:
128: . dm - The `DM`
130: Output Parameters:
131: . uniform - Is the metric uniform?
133: Level: beginner
135: .seealso: `DMPLEX`, `DMPlexMetricSetUniform()`, `DMPlexMetricIsIsotropic()`, `DMPlexMetricRestrictAnisotropyFirst()`
136: @*/
137: PetscErrorCode DMPlexMetricIsUniform(DM dm, PetscBool *uniform)
138: {
139: DM_Plex *plex = (DM_Plex *)dm->data;
141: PetscFunctionBegin;
142: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
143: *uniform = plex->metricCtx->uniform;
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: /*@
148: DMPlexMetricSetRestrictAnisotropyFirst - Record whether anisotropy should be restricted before normalization
150: Input Parameters:
151: + dm - The `DM`
152: - restrictAnisotropyFirst - Should anisotropy be normalized first?
154: Level: beginner
156: .seealso: `DMPLEX`, `DMPlexMetricSetIsotropic()`, `DMPlexMetricRestrictAnisotropyFirst()`
157: @*/
158: PetscErrorCode DMPlexMetricSetRestrictAnisotropyFirst(DM dm, PetscBool restrictAnisotropyFirst)
159: {
160: DM_Plex *plex = (DM_Plex *)dm->data;
162: PetscFunctionBegin;
163: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
164: plex->metricCtx->restrictAnisotropyFirst = restrictAnisotropyFirst;
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: /*@
169: DMPlexMetricRestrictAnisotropyFirst - Is anisotropy restricted before normalization or after?
171: Input Parameters:
172: . dm - The `DM`
174: Output Parameters:
175: . restrictAnisotropyFirst - Is anisotropy be normalized first?
177: Level: beginner
179: .seealso: `DMPLEX`, `DMPlexMetricIsIsotropic()`, `DMPlexMetricSetRestrictAnisotropyFirst()`
180: @*/
181: PetscErrorCode DMPlexMetricRestrictAnisotropyFirst(DM dm, PetscBool *restrictAnisotropyFirst)
182: {
183: DM_Plex *plex = (DM_Plex *)dm->data;
185: PetscFunctionBegin;
186: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
187: *restrictAnisotropyFirst = plex->metricCtx->restrictAnisotropyFirst;
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: /*@
192: DMPlexMetricSetNoInsertion - Should node insertion and deletion be turned off?
194: Input Parameters:
195: + dm - The `DM`
196: - noInsert - Should node insertion and deletion be turned off?
198: Level: beginner
200: Note:
201: This is only used by Mmg and ParMmg (not Pragmatic).
203: .seealso: `DMPLEX`, `DMPlexMetricNoInsertion()`, `DMPlexMetricSetNoSwapping()`, `DMPlexMetricSetNoMovement()`, `DMPlexMetricSetNoSurf()`
204: @*/
205: PetscErrorCode DMPlexMetricSetNoInsertion(DM dm, PetscBool noInsert)
206: {
207: DM_Plex *plex = (DM_Plex *)dm->data;
209: PetscFunctionBegin;
210: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
211: plex->metricCtx->noInsert = noInsert;
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: /*@
216: DMPlexMetricNoInsertion - Are node insertion and deletion turned off?
218: Input Parameters:
219: . dm - The `DM`
221: Output Parameters:
222: . noInsert - Are node insertion and deletion turned off?
224: Level: beginner
226: Note:
227: This is only used by Mmg and ParMmg (not Pragmatic).
229: .seealso: `DMPLEX`, `DMPlexMetricSetNoInsertion()`, `DMPlexMetricNoSwapping()`, `DMPlexMetricNoMovement()`, `DMPlexMetricNoSurf()`
230: @*/
231: PetscErrorCode DMPlexMetricNoInsertion(DM dm, PetscBool *noInsert)
232: {
233: DM_Plex *plex = (DM_Plex *)dm->data;
235: PetscFunctionBegin;
236: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
237: *noInsert = plex->metricCtx->noInsert;
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: /*@
242: DMPlexMetricSetNoSwapping - Should facet swapping be turned off?
244: Input Parameters:
245: + dm - The `DM`
246: - noSwap - Should facet swapping be turned off?
248: Level: beginner
250: Note:
251: This is only used by Mmg and ParMmg (not Pragmatic).
253: .seealso: `DMPLEX`, `DMPlexMetricNoSwapping()`, `DMPlexMetricSetNoInsertion()`, `DMPlexMetricSetNoMovement()`, `DMPlexMetricSetNoSurf()`
254: @*/
255: PetscErrorCode DMPlexMetricSetNoSwapping(DM dm, PetscBool noSwap)
256: {
257: DM_Plex *plex = (DM_Plex *)dm->data;
259: PetscFunctionBegin;
260: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
261: plex->metricCtx->noSwap = noSwap;
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: /*@
266: DMPlexMetricNoSwapping - Is facet swapping turned off?
268: Input Parameters:
269: . dm - The `DM`
271: Output Parameters:
272: . noSwap - Is facet swapping turned off?
274: Level: beginner
276: Note:
277: This is only used by Mmg and ParMmg (not Pragmatic).
279: .seealso: `DMPLEX`, `DMPlexMetricSetNoSwapping()`, `DMPlexMetricNoInsertion()`, `DMPlexMetricNoMovement()`, `DMPlexMetricNoSurf()`
280: @*/
281: PetscErrorCode DMPlexMetricNoSwapping(DM dm, PetscBool *noSwap)
282: {
283: DM_Plex *plex = (DM_Plex *)dm->data;
285: PetscFunctionBegin;
286: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
287: *noSwap = plex->metricCtx->noSwap;
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: /*@
292: DMPlexMetricSetNoMovement - Should node movement be turned off?
294: Input Parameters:
295: + dm - The `DM`
296: - noMove - Should node movement be turned off?
298: Level: beginner
300: Note:
301: This is only used by Mmg and ParMmg (not Pragmatic).
303: .seealso: `DMPLEX`, `DMPlexMetricNoMovement()`, `DMPlexMetricSetNoInsertion()`, `DMPlexMetricSetNoSwapping()`, `DMPlexMetricSetNoSurf()`
304: @*/
305: PetscErrorCode DMPlexMetricSetNoMovement(DM dm, PetscBool noMove)
306: {
307: DM_Plex *plex = (DM_Plex *)dm->data;
309: PetscFunctionBegin;
310: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
311: plex->metricCtx->noMove = noMove;
312: PetscFunctionReturn(PETSC_SUCCESS);
313: }
315: /*@
316: DMPlexMetricNoMovement - Is node movement turned off?
318: Input Parameters:
319: . dm - The `DM`
321: Output Parameters:
322: . noMove - Is node movement turned off?
324: Level: beginner
326: Note:
327: This is only used by Mmg and ParMmg (not Pragmatic).
329: .seealso: `DMPLEX`, `DMPlexMetricSetNoMovement()`, `DMPlexMetricNoInsertion()`, `DMPlexMetricNoSwapping()`, `DMPlexMetricNoSurf()`
330: @*/
331: PetscErrorCode DMPlexMetricNoMovement(DM dm, PetscBool *noMove)
332: {
333: DM_Plex *plex = (DM_Plex *)dm->data;
335: PetscFunctionBegin;
336: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
337: *noMove = plex->metricCtx->noMove;
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: /*@
342: DMPlexMetricSetNoSurf - Should surface modification be turned off?
344: Input Parameters:
345: + dm - The `DM`
346: - noSurf - Should surface modification be turned off?
348: Level: beginner
350: Note:
351: This is only used by Mmg and ParMmg (not Pragmatic).
353: .seealso: `DMPLEX`, `DMPlexMetricNoSurf()`, `DMPlexMetricSetNoMovement()`, `DMPlexMetricSetNoInsertion()`, `DMPlexMetricSetNoSwapping()`
354: @*/
355: PetscErrorCode DMPlexMetricSetNoSurf(DM dm, PetscBool noSurf)
356: {
357: DM_Plex *plex = (DM_Plex *)dm->data;
359: PetscFunctionBegin;
360: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
361: plex->metricCtx->noSurf = noSurf;
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: /*@
366: DMPlexMetricNoSurf - Is surface modification turned off?
368: Input Parameters:
369: . dm - The `DM`
371: Output Parameters:
372: . noSurf - Is surface modification turned off?
374: Level: beginner
376: Note:
377: This is only used by Mmg and ParMmg (not Pragmatic).
379: .seealso: `DMPLEX`, `DMPlexMetricSetNoSurf()`, `DMPlexMetricNoMovement()`, `DMPlexMetricNoInsertion()`, `DMPlexMetricNoSwapping()`
380: @*/
381: PetscErrorCode DMPlexMetricNoSurf(DM dm, PetscBool *noSurf)
382: {
383: DM_Plex *plex = (DM_Plex *)dm->data;
385: PetscFunctionBegin;
386: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
387: *noSurf = plex->metricCtx->noSurf;
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: /*@
392: DMPlexMetricSetMinimumMagnitude - Set the minimum tolerated metric magnitude
394: Input Parameters:
395: + dm - The `DM`
396: - h_min - The minimum tolerated metric magnitude
398: Level: beginner
400: .seealso: `DMPLEX`, `DMPlexMetricGetMinimumMagnitude()`, `DMPlexMetricSetMaximumMagnitude()`
401: @*/
402: PetscErrorCode DMPlexMetricSetMinimumMagnitude(DM dm, PetscReal h_min)
403: {
404: DM_Plex *plex = (DM_Plex *)dm->data;
406: PetscFunctionBegin;
407: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
408: PetscCheck(h_min > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Metric magnitudes must be in (0, inf)");
409: plex->metricCtx->h_min = h_min;
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: /*@
414: DMPlexMetricGetMinimumMagnitude - Get the minimum tolerated metric magnitude
416: Input Parameters:
417: . dm - The `DM`
419: Output Parameters:
420: . h_min - The minimum tolerated metric magnitude
422: Level: beginner
424: .seealso: `DMPLEX`, `DMPlexMetricSetMinimumMagnitude()`, `DMPlexMetricGetMaximumMagnitude()`
425: @*/
426: PetscErrorCode DMPlexMetricGetMinimumMagnitude(DM dm, PetscReal *h_min)
427: {
428: DM_Plex *plex = (DM_Plex *)dm->data;
430: PetscFunctionBegin;
431: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
432: *h_min = plex->metricCtx->h_min;
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: /*@
437: DMPlexMetricSetMaximumMagnitude - Set the maximum tolerated metric magnitude
439: Input Parameters:
440: + dm - The `DM`
441: - h_max - The maximum tolerated metric magnitude
443: Level: beginner
445: .seealso: `DMPLEX`, `DMPlexMetricGetMaximumMagnitude()`, `DMPlexMetricSetMinimumMagnitude()`
446: @*/
447: PetscErrorCode DMPlexMetricSetMaximumMagnitude(DM dm, PetscReal h_max)
448: {
449: DM_Plex *plex = (DM_Plex *)dm->data;
451: PetscFunctionBegin;
452: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
453: PetscCheck(h_max > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Metric magnitudes must be in (0, inf)");
454: plex->metricCtx->h_max = h_max;
455: PetscFunctionReturn(PETSC_SUCCESS);
456: }
458: /*@
459: DMPlexMetricGetMaximumMagnitude - Get the maximum tolerated metric magnitude
461: Input Parameters:
462: . dm - The `DM`
464: Output Parameters:
465: . h_max - The maximum tolerated metric magnitude
467: Level: beginner
469: .seealso: `DMPLEX`, `DMPlexMetricSetMaximumMagnitude()`, `DMPlexMetricGetMinimumMagnitude()`
470: @*/
471: PetscErrorCode DMPlexMetricGetMaximumMagnitude(DM dm, PetscReal *h_max)
472: {
473: DM_Plex *plex = (DM_Plex *)dm->data;
475: PetscFunctionBegin;
476: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
477: *h_max = plex->metricCtx->h_max;
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: /*@
482: DMPlexMetricSetMaximumAnisotropy - Set the maximum tolerated metric anisotropy
484: Input Parameters:
485: + dm - The `DM`
486: - a_max - The maximum tolerated metric anisotropy
488: Level: beginner
490: Note:
491: If the value zero is given then anisotropy will not be restricted. Otherwise, it should be at least one.
493: .seealso: `DMPLEX`, `DMPlexMetricGetMaximumAnisotropy()`, `DMPlexMetricSetMaximumMagnitude()`
494: @*/
495: PetscErrorCode DMPlexMetricSetMaximumAnisotropy(DM dm, PetscReal a_max)
496: {
497: DM_Plex *plex = (DM_Plex *)dm->data;
499: PetscFunctionBegin;
500: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
501: PetscCheck(a_max >= 1.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Anisotropy must be in [1, inf)");
502: plex->metricCtx->a_max = a_max;
503: PetscFunctionReturn(PETSC_SUCCESS);
504: }
506: /*@
507: DMPlexMetricGetMaximumAnisotropy - Get the maximum tolerated metric anisotropy
509: Input Parameters:
510: . dm - The `DM`
512: Output Parameters:
513: . a_max - The maximum tolerated metric anisotropy
515: Level: beginner
517: .seealso: `DMPLEX`, `DMPlexMetricSetMaximumAnisotropy()`, `DMPlexMetricGetMaximumMagnitude()`
518: @*/
519: PetscErrorCode DMPlexMetricGetMaximumAnisotropy(DM dm, PetscReal *a_max)
520: {
521: DM_Plex *plex = (DM_Plex *)dm->data;
523: PetscFunctionBegin;
524: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
525: *a_max = plex->metricCtx->a_max;
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: /*@
530: DMPlexMetricSetTargetComplexity - Set the target metric complexity
532: Input Parameters:
533: + dm - The `DM`
534: - targetComplexity - The target metric complexity
536: Level: beginner
538: .seealso: `DMPLEX`, `DMPlexMetricGetTargetComplexity()`, `DMPlexMetricSetNormalizationOrder()`
539: @*/
540: PetscErrorCode DMPlexMetricSetTargetComplexity(DM dm, PetscReal targetComplexity)
541: {
542: DM_Plex *plex = (DM_Plex *)dm->data;
544: PetscFunctionBegin;
545: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
546: PetscCheck(targetComplexity > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Metric complexity must be in (0, inf)");
547: plex->metricCtx->targetComplexity = targetComplexity;
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }
551: /*@
552: DMPlexMetricGetTargetComplexity - Get the target metric complexity
554: Input Parameters:
555: . dm - The `DM`
557: Output Parameters:
558: . targetComplexity - The target metric complexity
560: Level: beginner
562: .seealso: `DMPLEX`, `DMPlexMetricSetTargetComplexity()`, `DMPlexMetricGetNormalizationOrder()`
563: @*/
564: PetscErrorCode DMPlexMetricGetTargetComplexity(DM dm, PetscReal *targetComplexity)
565: {
566: DM_Plex *plex = (DM_Plex *)dm->data;
568: PetscFunctionBegin;
569: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
570: *targetComplexity = plex->metricCtx->targetComplexity;
571: PetscFunctionReturn(PETSC_SUCCESS);
572: }
574: /*@
575: DMPlexMetricSetNormalizationOrder - Set the order p for L-p normalization
577: Input Parameters:
578: + dm - The `DM`
579: - p - The normalization order
581: Level: beginner
583: .seealso: `DMPLEX`, `DMPlexMetricGetNormalizationOrder()`, `DMPlexMetricSetTargetComplexity()`
584: @*/
585: PetscErrorCode DMPlexMetricSetNormalizationOrder(DM dm, PetscReal p)
586: {
587: DM_Plex *plex = (DM_Plex *)dm->data;
589: PetscFunctionBegin;
590: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
591: PetscCheck(p >= 1.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Normalization order must be in [1, inf)");
592: plex->metricCtx->p = p;
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
596: /*@
597: DMPlexMetricGetNormalizationOrder - Get the order p for L-p normalization
599: Input Parameters:
600: . dm - The `DM`
602: Output Parameters:
603: . p - The normalization order
605: Level: beginner
607: .seealso: `DMPLEX`, `DMPlexMetricSetNormalizationOrder()`, `DMPlexMetricGetTargetComplexity()`
608: @*/
609: PetscErrorCode DMPlexMetricGetNormalizationOrder(DM dm, PetscReal *p)
610: {
611: DM_Plex *plex = (DM_Plex *)dm->data;
613: PetscFunctionBegin;
614: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
615: *p = plex->metricCtx->p;
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: /*@
620: DMPlexMetricSetGradationFactor - Set the metric gradation factor
622: Input Parameters:
623: + dm - The `DM`
624: - beta - The metric gradation factor
626: Level: beginner
628: Notes:
629: The gradation factor is the maximum tolerated length ratio between adjacent edges.
631: Turn off gradation by passing the value -1. Otherwise, pass a positive value.
633: This is only used by Mmg and ParMmg (not Pragmatic).
635: .seealso: `DMPLEX`, `DMPlexMetricGetGradationFactor()`, `DMPlexMetricSetHausdorffNumber()`
636: @*/
637: PetscErrorCode DMPlexMetricSetGradationFactor(DM dm, PetscReal beta)
638: {
639: DM_Plex *plex = (DM_Plex *)dm->data;
641: PetscFunctionBegin;
642: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
643: PetscCheck(beta > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Metric gradation factor must be in (0, inf)");
644: plex->metricCtx->gradationFactor = beta;
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
648: /*@
649: DMPlexMetricGetGradationFactor - Get the metric gradation factor
651: Input Parameters:
652: . dm - The `DM`
654: Output Parameters:
655: . beta - The metric gradation factor
657: Level: beginner
659: Notes:
660: The gradation factor is the maximum tolerated length ratio between adjacent edges.
662: The value -1 implies that gradation is turned off.
664: This is only used by Mmg and ParMmg (not Pragmatic).
666: .seealso: `DMPLEX`, `DMPlexMetricSetGradationFactor()`, `DMPlexMetricGetHausdorffNumber()`
667: @*/
668: PetscErrorCode DMPlexMetricGetGradationFactor(DM dm, PetscReal *beta)
669: {
670: DM_Plex *plex = (DM_Plex *)dm->data;
672: PetscFunctionBegin;
673: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
674: *beta = plex->metricCtx->gradationFactor;
675: PetscFunctionReturn(PETSC_SUCCESS);
676: }
678: /*@
679: DMPlexMetricSetHausdorffNumber - Set the metric Hausdorff number
681: Input Parameters:
682: + dm - The `DM`
683: - hausd - The metric Hausdorff number
685: Level: beginner
687: Notes:
688: The Hausdorff number imposes the maximal distance between the piecewise linear approximation of the
689: boundary and the reconstructed ideal boundary. Thus, a low Hausdorff parameter leads to refine the
690: high curvature areas. By default, the Hausdorff value is set to 0.01, which is a suitable value for
691: an object of size 1 in each direction. For smaller (resp. larger) objects, you may need to decrease
692: (resp. increase) the Hausdorff parameter. (Taken from
693: https://www.mmgtools.org/mmg-remesher-try-mmg/mmg-remesher-options/mmg-remesher-option-hausd).
695: This is only used by Mmg and ParMmg (not Pragmatic).
697: .seealso: `DMPLEX`, `DMPlexMetricSetGradationFactor()`, `DMPlexMetricGetHausdorffNumber()`
698: @*/
699: PetscErrorCode DMPlexMetricSetHausdorffNumber(DM dm, PetscReal hausd)
700: {
701: DM_Plex *plex = (DM_Plex *)dm->data;
703: PetscFunctionBegin;
704: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
705: PetscCheck(hausd > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Metric Hausdorff number must be in (0, inf)");
706: plex->metricCtx->hausdorffNumber = hausd;
707: PetscFunctionReturn(PETSC_SUCCESS);
708: }
710: /*@
711: DMPlexMetricGetHausdorffNumber - Get the metric Hausdorff number
713: Input Parameters:
714: . dm - The `DM`
716: Output Parameters:
717: . hausd - The metric Hausdorff number
719: Level: beginner
721: Notes:
722: The Hausdorff number imposes the maximal distance between the piecewise linear approximation of the
723: boundary and the reconstructed ideal boundary. Thus, a low Hausdorff parameter leads to refine the
724: high curvature areas. By default, the Hausdorff value is set to 0.01, which is a suitable value for
725: an object of size 1 in each direction. For smaller (resp. larger) objects, you may need to decrease
726: (resp. increase) the Hausdorff parameter. (Taken from
727: https://www.mmgtools.org/mmg-remesher-try-mmg/mmg-remesher-options/mmg-remesher-option-hausd).
729: This is only used by Mmg and ParMmg (not Pragmatic).
731: .seealso: `DMPLEX`, `DMPlexMetricGetGradationFactor()`, `DMPlexMetricSetHausdorffNumber()`
732: @*/
733: PetscErrorCode DMPlexMetricGetHausdorffNumber(DM dm, PetscReal *hausd)
734: {
735: DM_Plex *plex = (DM_Plex *)dm->data;
737: PetscFunctionBegin;
738: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
739: *hausd = plex->metricCtx->hausdorffNumber;
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
743: /*@
744: DMPlexMetricSetVerbosity - Set the verbosity of the mesh adaptation package
746: Input Parameters:
747: + dm - The `DM`
748: - verbosity - The verbosity, where -1 is silent and 10 is maximum
750: Level: beginner
752: Note:
753: This is only used by Mmg and ParMmg (not Pragmatic).
755: .seealso: `DMPLEX`, `DMPlexMetricGetVerbosity()`, `DMPlexMetricSetNumIterations()`
756: @*/
757: PetscErrorCode DMPlexMetricSetVerbosity(DM dm, PetscInt verbosity)
758: {
759: DM_Plex *plex = (DM_Plex *)dm->data;
761: PetscFunctionBegin;
762: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
763: plex->metricCtx->verbosity = verbosity;
764: PetscFunctionReturn(PETSC_SUCCESS);
765: }
767: /*@
768: DMPlexMetricGetVerbosity - Get the verbosity of the mesh adaptation package
770: Input Parameters:
771: . dm - The `DM`
773: Output Parameters:
774: . verbosity - The verbosity, where -1 is silent and 10 is maximum
776: Level: beginner
778: Note:
779: This is only used by Mmg and ParMmg (not Pragmatic).
781: .seealso: `DMPLEX`, `DMPlexMetricSetVerbosity()`, `DMPlexMetricGetNumIterations()`
782: @*/
783: PetscErrorCode DMPlexMetricGetVerbosity(DM dm, PetscInt *verbosity)
784: {
785: DM_Plex *plex = (DM_Plex *)dm->data;
787: PetscFunctionBegin;
788: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
789: *verbosity = plex->metricCtx->verbosity;
790: PetscFunctionReturn(PETSC_SUCCESS);
791: }
793: /*@
794: DMPlexMetricSetNumIterations - Set the number of parallel adaptation iterations
796: Input Parameters:
797: + dm - The `DM`
798: - numIter - the number of parallel adaptation iterations
800: Level: beginner
802: Note:
803: This is only used by ParMmg (not Pragmatic or Mmg).
805: .seealso: `DMPLEX`, `DMPlexMetricSetVerbosity()`, `DMPlexMetricGetNumIterations()`
806: @*/
807: PetscErrorCode DMPlexMetricSetNumIterations(DM dm, PetscInt numIter)
808: {
809: DM_Plex *plex = (DM_Plex *)dm->data;
811: PetscFunctionBegin;
812: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
813: plex->metricCtx->numIter = numIter;
814: PetscFunctionReturn(PETSC_SUCCESS);
815: }
817: /*@
818: DMPlexMetricGetNumIterations - Get the number of parallel adaptation iterations
820: Input Parameters:
821: . dm - The `DM`
823: Output Parameters:
824: . numIter - the number of parallel adaptation iterations
826: Level: beginner
828: Note:
829: This is only used by Mmg and ParMmg (not Pragmatic or Mmg).
831: .seealso: `DMPLEX`, `DMPlexMetricSetNumIterations()`, `DMPlexMetricGetVerbosity()`
832: @*/
833: PetscErrorCode DMPlexMetricGetNumIterations(DM dm, PetscInt *numIter)
834: {
835: DM_Plex *plex = (DM_Plex *)dm->data;
837: PetscFunctionBegin;
838: if (!plex->metricCtx) PetscCall(DMPlexMetricSetFromOptions(dm));
839: *numIter = plex->metricCtx->numIter;
840: PetscFunctionReturn(PETSC_SUCCESS);
841: }
843: static PetscErrorCode DMPlexP1FieldCreate_Private(DM dm, PetscInt f, PetscInt size, Vec *metric)
844: {
845: MPI_Comm comm;
846: PetscFE fe;
847: PetscInt dim;
849: PetscFunctionBegin;
851: /* Extract metadata from dm */
852: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
853: PetscCall(DMGetDimension(dm, &dim));
855: /* Create a P1 field of the requested size */
856: PetscCall(PetscFECreateLagrange(comm, dim, size, PETSC_TRUE, 1, PETSC_DETERMINE, &fe));
857: PetscCall(DMSetField(dm, f, NULL, (PetscObject)fe));
858: PetscCall(DMCreateDS(dm));
859: PetscCall(PetscFEDestroy(&fe));
860: PetscCall(DMCreateLocalVector(dm, metric));
862: PetscFunctionReturn(PETSC_SUCCESS);
863: }
865: /*@
866: DMPlexMetricCreate - Create a Riemannian metric field
868: Input Parameters:
869: + dm - The `DM`
870: - f - The field number to use
872: Output Parameter:
873: . metric - The metric
875: Options Database Key:
876: . -dm_adaptor <pragmatic/mmg/parmmg> - specify dm adaptor to use
878: Options Database Keys for Mmg and ParMmg:
879: + -dm_plex_metric_gradation_factor - Maximum ratio by which edge lengths may grow during gradation
880: . -dm_plex_metric_num_iterations - Number of parallel mesh adaptation iterations for ParMmg
881: . -dm_plex_metric_no_insert - Should node insertion/deletion be turned off?
882: . -dm_plex_metric_no_swap - Should facet swapping be turned off?
883: . -dm_plex_metric_no_move - Should node movement be turned off?
884: - -dm_plex_metric_verbosity - Choose a verbosity level from -1 (silent) to 10 (maximum).
886: Options Database Keys for Riemannian metrics:
887: + -dm_plex_metric_isotropic - Is the metric isotropic?
888: . -dm_plex_metric_uniform - Is the metric uniform?
889: . -dm_plex_metric_restrict_anisotropy_first - Should anisotropy be restricted before normalization?
890: . -dm_plex_metric_h_min - Minimum tolerated metric magnitude
891: . -dm_plex_metric_h_max - Maximum tolerated metric magnitude
892: . -dm_plex_metric_a_max - Maximum tolerated anisotropy
893: . -dm_plex_metric_p - L-p normalization order
894: - -dm_plex_metric_target_complexity - Target metric complexity
896: Level: beginner
898: Note:
899: It is assumed that the `DM` is comprised of simplices.
901: .seealso: `DMPLEX`, `DMPlexMetricCreateUniform()`, `DMPlexMetricCreateIsotropic()`
902: @*/
903: PetscErrorCode DMPlexMetricCreate(DM dm, PetscInt f, Vec *metric)
904: {
905: PetscBool isotropic, uniform;
906: PetscInt coordDim, Nd;
908: PetscFunctionBegin;
909: PetscCall(DMGetCoordinateDim(dm, &coordDim));
910: Nd = coordDim * coordDim;
911: PetscCall(DMPlexMetricIsUniform(dm, &uniform));
912: PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic));
913: if (uniform) {
914: MPI_Comm comm;
916: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
917: PetscCheck(isotropic, comm, PETSC_ERR_SUP, "Uniform anisotropic metrics not supported");
918: PetscCall(VecCreate(comm, metric));
919: PetscCall(VecSetSizes(*metric, 1, PETSC_DECIDE));
920: PetscCall(VecSetFromOptions(*metric));
921: } else if (isotropic) PetscCall(DMPlexP1FieldCreate_Private(dm, f, 1, metric));
922: else PetscCall(DMPlexP1FieldCreate_Private(dm, f, Nd, metric));
923: PetscFunctionReturn(PETSC_SUCCESS);
924: }
926: /*@
927: DMPlexMetricCreateUniform - Construct a uniform isotropic metric
929: Input Parameters:
930: + dm - The `DM`
931: . f - The field number to use
932: - alpha - Scaling parameter for the diagonal
934: Output Parameter:
935: . metric - The uniform metric
937: Level: beginner
939: Note:
940: In this case, the metric is constant in space and so is only specified for a single vertex.
942: .seealso: `DMPLEX`, `DMPlexMetricCreate()`, `DMPlexMetricCreateIsotropic()`
943: @*/
944: PetscErrorCode DMPlexMetricCreateUniform(DM dm, PetscInt f, PetscReal alpha, Vec *metric)
945: {
946: PetscFunctionBegin;
947: PetscCall(DMPlexMetricSetUniform(dm, PETSC_TRUE));
948: PetscCall(DMPlexMetricCreate(dm, f, metric));
949: PetscCheck(alpha, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Uniform metric scaling is undefined");
950: PetscCheck(alpha > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Uniform metric scaling should be in (0, inf)");
951: PetscCall(VecSet(*metric, alpha));
952: PetscCall(VecAssemblyBegin(*metric));
953: PetscCall(VecAssemblyEnd(*metric));
954: PetscFunctionReturn(PETSC_SUCCESS);
955: }
957: 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[])
958: {
959: f0[0] = u[0];
960: }
962: /*@
963: DMPlexMetricCreateIsotropic - Construct an isotropic metric from an error indicator
965: Input Parameters:
966: + dm - The `DM`
967: . f - The field number to use
968: - indicator - The error indicator
970: Output Parameter:
971: . metric - The isotropic metric
973: Level: beginner
975: Notes:
976: It is assumed that the `DM` is comprised of simplices.
978: The indicator needs to be a scalar field. If it is not defined vertex-wise, then it is projected appropriately.
980: .seealso: `DMPLEX`, `DMPlexMetricCreate()`, `DMPlexMetricCreateUniform()`
981: @*/
982: PetscErrorCode DMPlexMetricCreateIsotropic(DM dm, PetscInt f, Vec indicator, Vec *metric)
983: {
984: PetscInt m, n;
986: PetscFunctionBegin;
987: PetscCall(DMPlexMetricSetIsotropic(dm, PETSC_TRUE));
988: PetscCall(DMPlexMetricCreate(dm, f, metric));
989: PetscCall(VecGetSize(indicator, &m));
990: PetscCall(VecGetSize(*metric, &n));
991: if (m == n) PetscCall(VecCopy(indicator, *metric));
992: else {
993: DM dmIndi;
994: 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[]);
996: PetscCall(VecGetDM(indicator, &dmIndi));
997: funcs[0] = identity;
998: PetscCall(DMProjectFieldLocal(dmIndi, 0.0, indicator, funcs, INSERT_VALUES, *metric));
999: }
1000: PetscFunctionReturn(PETSC_SUCCESS);
1001: }
1003: /*@
1004: DMPlexMetricDeterminantCreate - Create the determinant field for a Riemannian metric
1006: Input Parameters:
1007: + dm - The `DM` of the metric field
1008: - f - The field number to use
1010: Output Parameters:
1011: + determinant - The determinant field
1012: - dmDet - The corresponding `DM`
1014: Level: beginner
1016: .seealso: `DMPLEX`, `DMPlexMetricCreateUniform()`, `DMPlexMetricCreateIsotropic()`, `DMPlexMetricCreate()`
1017: @*/
1018: PetscErrorCode DMPlexMetricDeterminantCreate(DM dm, PetscInt f, Vec *determinant, DM *dmDet)
1019: {
1020: PetscBool uniform;
1022: PetscFunctionBegin;
1023: PetscCall(DMPlexMetricIsUniform(dm, &uniform));
1024: PetscCall(DMClone(dm, dmDet));
1025: if (uniform) {
1026: MPI_Comm comm;
1028: PetscCall(PetscObjectGetComm((PetscObject)*dmDet, &comm));
1029: PetscCall(VecCreate(comm, determinant));
1030: PetscCall(VecSetSizes(*determinant, 1, PETSC_DECIDE));
1031: PetscCall(VecSetFromOptions(*determinant));
1032: } else PetscCall(DMPlexP1FieldCreate_Private(*dmDet, f, 1, determinant));
1033: PetscFunctionReturn(PETSC_SUCCESS);
1034: }
1036: static PetscErrorCode LAPACKsyevFail(PetscInt dim, PetscScalar Mpos[])
1037: {
1038: PetscInt i, j;
1040: PetscFunctionBegin;
1041: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Failed to apply LAPACKsyev to the matrix\n"));
1042: for (i = 0; i < dim; ++i) {
1043: if (i == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " [["));
1044: else PetscCall(PetscPrintf(PETSC_COMM_SELF, " ["));
1045: for (j = 0; j < dim; ++j) {
1046: if (j < dim - 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%15.8e, ", (double)PetscAbsScalar(Mpos[i * dim + j])));
1047: else PetscCall(PetscPrintf(PETSC_COMM_SELF, "%15.8e", (double)PetscAbsScalar(Mpos[i * dim + j])));
1048: }
1049: if (i < dim - 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, "]\n"));
1050: else PetscCall(PetscPrintf(PETSC_COMM_SELF, "]]\n"));
1051: }
1052: PetscFunctionReturn(PETSC_SUCCESS);
1053: }
1055: static PetscErrorCode DMPlexMetricModify_Private(PetscInt dim, PetscReal h_min, PetscReal h_max, PetscReal a_max, PetscScalar Mp[], PetscScalar *detMp)
1056: {
1057: PetscInt i, j, k;
1058: 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);
1059: PetscScalar *Mpos;
1061: PetscFunctionBegin;
1062: PetscCall(PetscMalloc2(dim * dim, &Mpos, dim, &eigs));
1064: /* Symmetrize */
1065: for (i = 0; i < dim; ++i) {
1066: Mpos[i * dim + i] = Mp[i * dim + i];
1067: for (j = i + 1; j < dim; ++j) {
1068: Mpos[i * dim + j] = 0.5 * (Mp[i * dim + j] + Mp[j * dim + i]);
1069: Mpos[j * dim + i] = Mpos[i * dim + j];
1070: }
1071: }
1073: /* Compute eigendecomposition */
1074: if (dim == 1) {
1075: /* Isotropic case */
1076: eigs[0] = PetscRealPart(Mpos[0]);
1077: Mpos[0] = 1.0;
1078: } else {
1079: /* Anisotropic case */
1080: PetscScalar *work;
1081: PetscBLASInt lwork;
1083: lwork = 5 * dim;
1084: PetscCall(PetscMalloc1(5 * dim, &work));
1085: {
1086: PetscBLASInt lierr;
1087: PetscBLASInt nb;
1089: PetscCall(PetscBLASIntCast(dim, &nb));
1090: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1091: #if defined(PETSC_USE_COMPLEX)
1092: {
1093: PetscReal *rwork;
1094: PetscCall(PetscMalloc1(3 * dim, &rwork));
1095: PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, Mpos, &nb, eigs, work, &lwork, rwork, &lierr));
1096: PetscCall(PetscFree(rwork));
1097: }
1098: #else
1099: PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, Mpos, &nb, eigs, work, &lwork, &lierr));
1100: #endif
1101: if (lierr) {
1102: for (i = 0; i < dim; ++i) {
1103: Mpos[i * dim + i] = Mp[i * dim + i];
1104: for (j = i + 1; j < dim; ++j) {
1105: Mpos[i * dim + j] = 0.5 * (Mp[i * dim + j] + Mp[j * dim + i]);
1106: Mpos[j * dim + i] = Mpos[i * dim + j];
1107: }
1108: }
1109: PetscCall(LAPACKsyevFail(dim, Mpos));
1110: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int)lierr);
1111: }
1112: PetscCall(PetscFPTrapPop());
1113: }
1114: PetscCall(PetscFree(work));
1115: }
1117: /* Reflect to positive orthant and enforce maximum and minimum size */
1118: max_eig = 0.0;
1119: for (i = 0; i < dim; ++i) {
1120: eigs[i] = PetscMin(l_max, PetscMax(l_min, PetscAbsReal(eigs[i])));
1121: max_eig = PetscMax(eigs[i], max_eig);
1122: }
1124: /* Enforce maximum anisotropy and compute determinant */
1125: *detMp = 1.0;
1126: for (i = 0; i < dim; ++i) {
1127: if (a_max > 1.0) eigs[i] = PetscMax(eigs[i], max_eig * la_min);
1128: *detMp *= eigs[i];
1129: }
1131: /* Reconstruct Hessian */
1132: for (i = 0; i < dim; ++i) {
1133: for (j = 0; j < dim; ++j) {
1134: Mp[i * dim + j] = 0.0;
1135: for (k = 0; k < dim; ++k) Mp[i * dim + j] += Mpos[k * dim + i] * eigs[k] * Mpos[k * dim + j];
1136: }
1137: }
1138: PetscCall(PetscFree2(Mpos, eigs));
1140: PetscFunctionReturn(PETSC_SUCCESS);
1141: }
1143: /*@
1144: DMPlexMetricEnforceSPD - Enforce symmetric positive-definiteness of a metric
1146: Input Parameters:
1147: + dm - The `DM`
1148: . metricIn - The metric
1149: . restrictSizes - Should maximum/minimum metric magnitudes be enforced?
1150: - restrictAnisotropy - Should maximum anisotropy be enforced?
1152: Output Parameters:
1153: + metricOut - The metric
1154: - determinant - Its determinant
1156: Options Database Keys:
1157: + -dm_plex_metric_isotropic - Is the metric isotropic?
1158: . -dm_plex_metric_uniform - Is the metric uniform?
1159: . -dm_plex_metric_h_min - Minimum tolerated metric magnitude
1160: . -dm_plex_metric_h_max - Maximum tolerated metric magnitude
1161: - -dm_plex_metric_a_max - Maximum tolerated anisotropy
1163: Level: beginner
1165: .seealso: `DMPLEX`, `DMPlexMetricNormalize()`, `DMPlexMetricIntersection()`
1166: @*/
1167: PetscErrorCode DMPlexMetricEnforceSPD(DM dm, Vec metricIn, PetscBool restrictSizes, PetscBool restrictAnisotropy, Vec metricOut, Vec determinant)
1168: {
1169: DM dmDet;
1170: PetscBool isotropic, uniform;
1171: PetscInt dim, vStart, vEnd, v;
1172: PetscScalar *met, *det;
1173: PetscReal h_min = 1.0e-30, h_max = 1.0e+30, a_max = 1.0e+15;
1175: PetscFunctionBegin;
1176: PetscCall(PetscLogEventBegin(DMPLEX_MetricEnforceSPD, 0, 0, 0, 0));
1178: /* Extract metadata from dm */
1179: PetscCall(DMGetDimension(dm, &dim));
1180: if (restrictSizes) {
1181: PetscCall(DMPlexMetricGetMinimumMagnitude(dm, &h_min));
1182: PetscCall(DMPlexMetricGetMaximumMagnitude(dm, &h_max));
1183: h_min = PetscMax(h_min, 1.0e-30);
1184: h_max = PetscMin(h_max, 1.0e+30);
1185: PetscCheck(h_min < h_max, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Minimum metric magnitude should be smaller than maximum metric magnitude");
1186: }
1187: if (restrictAnisotropy) {
1188: PetscCall(DMPlexMetricGetMaximumAnisotropy(dm, &a_max));
1189: a_max = PetscMin(a_max, 1.0e+30);
1190: }
1192: /* Setup output metric */
1193: PetscCall(VecCopy(metricIn, metricOut));
1195: /* Enforce SPD and extract determinant */
1196: PetscCall(VecGetArray(metricOut, &met));
1197: PetscCall(DMPlexMetricIsUniform(dm, &uniform));
1198: PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic));
1199: if (uniform) {
1200: PetscCheck(isotropic, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Uniform anisotropic metrics cannot exist");
1202: /* Uniform case */
1203: PetscCall(VecGetArray(determinant, &det));
1204: PetscCall(DMPlexMetricModify_Private(1, h_min, h_max, a_max, met, det));
1205: PetscCall(VecRestoreArray(determinant, &det));
1206: } else {
1207: /* Spatially varying case */
1208: PetscInt nrow;
1210: if (isotropic) nrow = 1;
1211: else nrow = dim;
1212: PetscCall(VecGetDM(determinant, &dmDet));
1213: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1214: PetscCall(VecGetArray(determinant, &det));
1215: for (v = vStart; v < vEnd; ++v) {
1216: PetscScalar *vmet, *vdet;
1217: PetscCall(DMPlexPointLocalRef(dm, v, met, &vmet));
1218: PetscCall(DMPlexPointLocalRef(dmDet, v, det, &vdet));
1219: PetscCall(DMPlexMetricModify_Private(nrow, h_min, h_max, a_max, vmet, vdet));
1220: }
1221: PetscCall(VecRestoreArray(determinant, &det));
1222: }
1223: PetscCall(VecRestoreArray(metricOut, &met));
1225: PetscCall(PetscLogEventEnd(DMPLEX_MetricEnforceSPD, 0, 0, 0, 0));
1226: PetscFunctionReturn(PETSC_SUCCESS);
1227: }
1229: 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[])
1230: {
1231: const PetscScalar p = constants[0];
1233: f0[0] = PetscPowScalar(u[0], p / (2.0 * p + dim));
1234: }
1236: /*@
1237: DMPlexMetricNormalize - Apply L-p normalization to a metric
1239: Input Parameters:
1240: + dm - The `DM`
1241: . metricIn - The unnormalized metric
1242: . restrictSizes - Should maximum/minimum metric magnitudes be enforced?
1243: - restrictAnisotropy - Should maximum metric anisotropy be enforced?
1245: Output Parameters:
1246: + metricOut - The normalized metric
1247: - determinant - computed determinant
1249: Options Database Keys:
1250: + -dm_plex_metric_isotropic - Is the metric isotropic?
1251: . -dm_plex_metric_uniform - Is the metric uniform?
1252: . -dm_plex_metric_restrict_anisotropy_first - Should anisotropy be restricted before normalization?
1253: . -dm_plex_metric_h_min - Minimum tolerated metric magnitude
1254: . -dm_plex_metric_h_max - Maximum tolerated metric magnitude
1255: . -dm_plex_metric_a_max - Maximum tolerated anisotropy
1256: . -dm_plex_metric_p - L-p normalization order
1257: - -dm_plex_metric_target_complexity - Target metric complexity
1259: Level: beginner
1261: .seealso: `DMPLEX`, `DMPlexMetricEnforceSPD()`, `DMPlexMetricIntersection()`
1262: @*/
1263: PetscErrorCode DMPlexMetricNormalize(DM dm, Vec metricIn, PetscBool restrictSizes, PetscBool restrictAnisotropy, Vec metricOut, Vec determinant)
1264: {
1265: DM dmDet;
1266: MPI_Comm comm;
1267: PetscBool restrictAnisotropyFirst, isotropic, uniform;
1268: PetscDS ds;
1269: PetscInt dim, Nd, vStart, vEnd, v, i;
1270: PetscScalar *met, *det, integral, constants[1];
1271: PetscReal p, h_min = 1.0e-30, h_max = 1.0e+30, a_max = 0.0, factGlob, fact, target, realIntegral;
1273: PetscFunctionBegin;
1274: PetscCall(PetscLogEventBegin(DMPLEX_MetricNormalize, 0, 0, 0, 0));
1276: /* Extract metadata from dm */
1277: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1278: PetscCall(DMGetDimension(dm, &dim));
1279: PetscCall(DMPlexMetricIsUniform(dm, &uniform));
1280: PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic));
1281: if (isotropic) Nd = 1;
1282: else Nd = dim * dim;
1284: /* Set up metric and ensure it is SPD */
1285: PetscCall(DMPlexMetricRestrictAnisotropyFirst(dm, &restrictAnisotropyFirst));
1286: PetscCall(DMPlexMetricEnforceSPD(dm, metricIn, PETSC_FALSE, (PetscBool)(restrictAnisotropy && restrictAnisotropyFirst), metricOut, determinant));
1288: /* Compute global normalization factor */
1289: PetscCall(DMPlexMetricGetTargetComplexity(dm, &target));
1290: PetscCall(DMPlexMetricGetNormalizationOrder(dm, &p));
1291: constants[0] = p;
1292: if (uniform) {
1293: PetscCheck(isotropic, comm, PETSC_ERR_SUP, "Uniform anisotropic metrics not supported");
1294: DM dmTmp;
1295: Vec tmp;
1297: PetscCall(DMClone(dm, &dmTmp));
1298: PetscCall(DMPlexP1FieldCreate_Private(dmTmp, 0, 1, &tmp));
1299: PetscCall(VecGetArray(determinant, &det));
1300: PetscCall(VecSet(tmp, det[0]));
1301: PetscCall(VecRestoreArray(determinant, &det));
1302: PetscCall(DMGetDS(dmTmp, &ds));
1303: PetscCall(PetscDSSetConstants(ds, 1, constants));
1304: PetscCall(PetscDSSetObjective(ds, 0, detMFunc));
1305: PetscCall(DMPlexComputeIntegralFEM(dmTmp, tmp, &integral, NULL));
1306: PetscCall(VecDestroy(&tmp));
1307: PetscCall(DMDestroy(&dmTmp));
1308: } else {
1309: PetscCall(VecGetDM(determinant, &dmDet));
1310: PetscCall(DMGetDS(dmDet, &ds));
1311: PetscCall(PetscDSSetConstants(ds, 1, constants));
1312: PetscCall(PetscDSSetObjective(ds, 0, detMFunc));
1313: PetscCall(DMPlexComputeIntegralFEM(dmDet, determinant, &integral, NULL));
1314: }
1315: realIntegral = PetscRealPart(integral);
1316: 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?");
1317: factGlob = PetscPowReal(target / realIntegral, 2.0 / dim);
1319: /* Apply local scaling */
1320: if (restrictSizes) {
1321: PetscCall(DMPlexMetricGetMinimumMagnitude(dm, &h_min));
1322: PetscCall(DMPlexMetricGetMaximumMagnitude(dm, &h_max));
1323: h_min = PetscMax(h_min, 1.0e-30);
1324: h_max = PetscMin(h_max, 1.0e+30);
1325: PetscCheck(h_min < h_max, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Minimum metric magnitude should be smaller than maximum metric magnitude");
1326: }
1327: if (restrictAnisotropy && !restrictAnisotropyFirst) {
1328: PetscCall(DMPlexMetricGetMaximumAnisotropy(dm, &a_max));
1329: a_max = PetscMin(a_max, 1.0e+30);
1330: }
1331: PetscCall(VecGetArray(metricOut, &met));
1332: PetscCall(VecGetArray(determinant, &det));
1333: if (uniform) {
1334: /* Uniform case */
1335: met[0] *= factGlob * PetscPowReal(PetscAbsScalar(det[0]), -1.0 / (2 * p + dim));
1336: if (restrictSizes) PetscCall(DMPlexMetricModify_Private(1, h_min, h_max, a_max, met, det));
1337: } else {
1338: /* Spatially varying case */
1339: PetscInt nrow;
1341: if (isotropic) nrow = 1;
1342: else nrow = dim;
1343: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1344: PetscCall(VecGetDM(determinant, &dmDet));
1345: for (v = vStart; v < vEnd; ++v) {
1346: PetscScalar *Mp, *detM;
1348: PetscCall(DMPlexPointLocalRef(dm, v, met, &Mp));
1349: PetscCall(DMPlexPointLocalRef(dmDet, v, det, &detM));
1350: fact = factGlob * PetscPowReal(PetscAbsScalar(detM[0]), -1.0 / (2 * p + dim));
1351: for (i = 0; i < Nd; ++i) Mp[i] *= fact;
1352: if (restrictSizes) PetscCall(DMPlexMetricModify_Private(nrow, h_min, h_max, a_max, Mp, detM));
1353: }
1354: }
1355: PetscCall(VecRestoreArray(determinant, &det));
1356: PetscCall(VecRestoreArray(metricOut, &met));
1358: PetscCall(PetscLogEventEnd(DMPLEX_MetricNormalize, 0, 0, 0, 0));
1359: PetscFunctionReturn(PETSC_SUCCESS);
1360: }
1362: /*@
1363: DMPlexMetricAverage - Compute the average of a list of metrics
1365: Input Parameters:
1366: + dm - The `DM`
1367: . numMetrics - The number of metrics to be averaged
1368: . weights - Weights for the average
1369: - metrics - The metrics to be averaged
1371: Output Parameter:
1372: . metricAvg - The averaged metric
1374: Level: beginner
1376: Notes:
1377: The weights should sum to unity.
1379: If weights are not provided then an unweighted average is used.
1381: .seealso: `DMPLEX`, `DMPlexMetricAverage2()`, `DMPlexMetricAverage3()`, `DMPlexMetricIntersection()`
1382: @*/
1383: PetscErrorCode DMPlexMetricAverage(DM dm, PetscInt numMetrics, PetscReal weights[], Vec metrics[], Vec metricAvg)
1384: {
1385: PetscBool haveWeights = PETSC_TRUE;
1386: PetscInt i, m, n;
1387: PetscReal sum = 0.0, tol = 1.0e-10;
1389: PetscFunctionBegin;
1390: PetscCall(PetscLogEventBegin(DMPLEX_MetricAverage, 0, 0, 0, 0));
1391: PetscCheck(numMetrics >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot average %" PetscInt_FMT " < 1 metrics", numMetrics);
1392: PetscCall(VecSet(metricAvg, 0.0));
1393: PetscCall(VecGetSize(metricAvg, &m));
1394: for (i = 0; i < numMetrics; ++i) {
1395: PetscCall(VecGetSize(metrics[i], &n));
1396: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Averaging different metric types not implemented");
1397: }
1399: /* Default to the unweighted case */
1400: if (!weights) {
1401: PetscCall(PetscMalloc1(numMetrics, &weights));
1402: haveWeights = PETSC_FALSE;
1403: for (i = 0; i < numMetrics; ++i) weights[i] = 1.0 / numMetrics;
1404: }
1406: /* Check weights sum to unity */
1407: for (i = 0; i < numMetrics; ++i) sum += weights[i];
1408: PetscCheck(PetscAbsReal(sum - 1) <= tol, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Weights do not sum to unity");
1410: /* Compute metric average */
1411: for (i = 0; i < numMetrics; ++i) PetscCall(VecAXPY(metricAvg, weights[i], metrics[i]));
1412: if (!haveWeights) PetscCall(PetscFree(weights));
1414: PetscCall(PetscLogEventEnd(DMPLEX_MetricAverage, 0, 0, 0, 0));
1415: PetscFunctionReturn(PETSC_SUCCESS);
1416: }
1418: /*@
1419: DMPlexMetricAverage2 - Compute the unweighted average of two metrics
1421: Input Parameters:
1422: + dm - The `DM`
1423: . metric1 - The first metric to be averaged
1424: - metric2 - The second metric to be averaged
1426: Output Parameter:
1427: . metricAvg - The averaged metric
1429: Level: beginner
1431: .seealso: `DMPLEX`, `DMPlexMetricAverage()`, `DMPlexMetricAverage3()`
1432: @*/
1433: PetscErrorCode DMPlexMetricAverage2(DM dm, Vec metric1, Vec metric2, Vec metricAvg)
1434: {
1435: PetscReal weights[2] = {0.5, 0.5};
1436: Vec metrics[2] = {metric1, metric2};
1438: PetscFunctionBegin;
1439: PetscCall(DMPlexMetricAverage(dm, 2, weights, metrics, metricAvg));
1440: PetscFunctionReturn(PETSC_SUCCESS);
1441: }
1443: /*@
1444: DMPlexMetricAverage3 - Compute the unweighted average of three metrics
1446: Input Parameters:
1447: + dm - The `DM`
1448: . metric1 - The first metric to be averaged
1449: . metric2 - The second metric to be averaged
1450: - metric3 - The third metric to be averaged
1452: Output Parameter:
1453: . metricAvg - The averaged metric
1455: Level: beginner
1457: .seealso: `DMPLEX`, `DMPlexMetricAverage()`, `DMPlexMetricAverage2()`
1458: @*/
1459: PetscErrorCode DMPlexMetricAverage3(DM dm, Vec metric1, Vec metric2, Vec metric3, Vec metricAvg)
1460: {
1461: PetscReal weights[3] = {1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0};
1462: Vec metrics[3] = {metric1, metric2, metric3};
1464: PetscFunctionBegin;
1465: PetscCall(DMPlexMetricAverage(dm, 3, weights, metrics, metricAvg));
1466: PetscFunctionReturn(PETSC_SUCCESS);
1467: }
1469: static PetscErrorCode DMPlexMetricIntersection_Private(PetscInt dim, PetscScalar M1[], PetscScalar M2[])
1470: {
1471: PetscInt i, j, k, l, m;
1472: PetscReal *evals;
1473: PetscScalar *evecs, *sqrtM1, *isqrtM1;
1475: PetscFunctionBegin;
1477: /* Isotropic case */
1478: if (dim == 1) {
1479: M2[0] = (PetscScalar)PetscMax(PetscRealPart(M1[0]), PetscRealPart(M2[0]));
1480: PetscFunctionReturn(PETSC_SUCCESS);
1481: }
1483: /* Anisotropic case */
1484: PetscCall(PetscMalloc4(dim * dim, &evecs, dim * dim, &sqrtM1, dim * dim, &isqrtM1, dim, &evals));
1485: for (i = 0; i < dim; ++i) {
1486: for (j = 0; j < dim; ++j) evecs[i * dim + j] = M1[i * dim + j];
1487: }
1488: {
1489: PetscScalar *work;
1490: PetscBLASInt lwork;
1492: lwork = 5 * dim;
1493: PetscCall(PetscMalloc1(5 * dim, &work));
1494: {
1495: PetscBLASInt lierr, nb;
1496: PetscReal sqrtj;
1498: /* Compute eigendecomposition of M1 */
1499: PetscCall(PetscBLASIntCast(dim, &nb));
1500: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1501: #if defined(PETSC_USE_COMPLEX)
1502: {
1503: PetscReal *rwork;
1504: PetscCall(PetscMalloc1(3 * dim, &rwork));
1505: PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals, work, &lwork, rwork, &lierr));
1506: PetscCall(PetscFree(rwork));
1507: }
1508: #else
1509: PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals, work, &lwork, &lierr));
1510: #endif
1511: if (lierr) {
1512: PetscCall(LAPACKsyevFail(dim, M1));
1513: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int)lierr);
1514: }
1515: PetscCall(PetscFPTrapPop());
1517: /* Compute square root and the reciprocal thereof */
1518: for (i = 0; i < dim; ++i) {
1519: for (k = 0; k < dim; ++k) {
1520: sqrtM1[i * dim + k] = 0.0;
1521: isqrtM1[i * dim + k] = 0.0;
1522: for (j = 0; j < dim; ++j) {
1523: sqrtj = PetscSqrtReal(evals[j]);
1524: sqrtM1[i * dim + k] += evecs[j * dim + i] * sqrtj * evecs[j * dim + k];
1525: isqrtM1[i * dim + k] += evecs[j * dim + i] * (1.0 / sqrtj) * evecs[j * dim + k];
1526: }
1527: }
1528: }
1530: /* Map M2 into the space spanned by the eigenvectors of M1 */
1531: for (i = 0; i < dim; ++i) {
1532: for (l = 0; l < dim; ++l) {
1533: evecs[i * dim + l] = 0.0;
1534: for (j = 0; j < dim; ++j) {
1535: for (k = 0; k < dim; ++k) evecs[i * dim + l] += isqrtM1[j * dim + i] * M2[j * dim + k] * isqrtM1[k * dim + l];
1536: }
1537: }
1538: }
1540: /* Compute eigendecomposition */
1541: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1542: #if defined(PETSC_USE_COMPLEX)
1543: {
1544: PetscReal *rwork;
1545: PetscCall(PetscMalloc1(3 * dim, &rwork));
1546: PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals, work, &lwork, rwork, &lierr));
1547: PetscCall(PetscFree(rwork));
1548: }
1549: #else
1550: PetscCallBLAS("LAPACKsyev", LAPACKsyev_("V", "U", &nb, evecs, &nb, evals, work, &lwork, &lierr));
1551: #endif
1552: if (lierr) {
1553: for (i = 0; i < dim; ++i) {
1554: for (l = 0; l < dim; ++l) {
1555: evecs[i * dim + l] = 0.0;
1556: for (j = 0; j < dim; ++j) {
1557: for (k = 0; k < dim; ++k) evecs[i * dim + l] += isqrtM1[j * dim + i] * M2[j * dim + k] * isqrtM1[k * dim + l];
1558: }
1559: }
1560: }
1561: PetscCall(LAPACKsyevFail(dim, evecs));
1562: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int)lierr);
1563: }
1564: PetscCall(PetscFPTrapPop());
1566: /* Modify eigenvalues */
1567: for (i = 0; i < dim; ++i) evals[i] = PetscMax(evals[i], 1.0);
1569: /* Map back to get the intersection */
1570: for (i = 0; i < dim; ++i) {
1571: for (m = 0; m < dim; ++m) {
1572: M2[i * dim + m] = 0.0;
1573: for (j = 0; j < dim; ++j) {
1574: for (k = 0; k < dim; ++k) {
1575: 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];
1576: }
1577: }
1578: }
1579: }
1580: }
1581: PetscCall(PetscFree(work));
1582: }
1583: PetscCall(PetscFree4(evecs, sqrtM1, isqrtM1, evals));
1584: PetscFunctionReturn(PETSC_SUCCESS);
1585: }
1587: /*@
1588: DMPlexMetricIntersection - Compute the intersection of a list of metrics
1590: Input Parameters:
1591: + dm - The `DM`
1592: . numMetrics - The number of metrics to be intersected
1593: - metrics - The metrics to be intersected
1595: Output Parameter:
1596: . metricInt - The intersected metric
1598: Level: beginner
1600: Notes:
1601: The intersection of a list of metrics has the minimal ellipsoid which fits within the ellipsoids of the component metrics.
1603: The implementation used here is only consistent with the minimal ellipsoid definition in the case numMetrics = 2.
1605: .seealso: `DMPLEX`, `DMPlexMetricIntersection2()`, `DMPlexMetricIntersection3()`, `DMPlexMetricAverage()`
1606: @*/
1607: PetscErrorCode DMPlexMetricIntersection(DM dm, PetscInt numMetrics, Vec metrics[], Vec metricInt)
1608: {
1609: PetscBool isotropic, uniform;
1610: PetscInt v, i, m, n;
1611: PetscScalar *met, *meti;
1613: PetscFunctionBegin;
1614: PetscCall(PetscLogEventBegin(DMPLEX_MetricIntersection, 0, 0, 0, 0));
1615: PetscCheck(numMetrics >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot intersect %" PetscInt_FMT " < 1 metrics", numMetrics);
1617: /* Copy over the first metric */
1618: PetscCall(VecCopy(metrics[0], metricInt));
1619: if (numMetrics == 1) PetscFunctionReturn(PETSC_SUCCESS);
1620: PetscCall(VecGetSize(metricInt, &m));
1621: for (i = 0; i < numMetrics; ++i) {
1622: PetscCall(VecGetSize(metrics[i], &n));
1623: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Intersecting different metric types not implemented");
1624: }
1626: /* Intersect subsequent metrics in turn */
1627: PetscCall(DMPlexMetricIsUniform(dm, &uniform));
1628: PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic));
1629: if (uniform) {
1630: /* Uniform case */
1631: PetscCall(VecGetArray(metricInt, &met));
1632: for (i = 1; i < numMetrics; ++i) {
1633: PetscCall(VecGetArray(metrics[i], &meti));
1634: PetscCall(DMPlexMetricIntersection_Private(1, meti, met));
1635: PetscCall(VecRestoreArray(metrics[i], &meti));
1636: }
1637: PetscCall(VecRestoreArray(metricInt, &met));
1638: } else {
1639: /* Spatially varying case */
1640: PetscInt dim, vStart, vEnd, nrow;
1641: PetscScalar *M, *Mi;
1643: PetscCall(DMGetDimension(dm, &dim));
1644: if (isotropic) nrow = 1;
1645: else nrow = dim;
1646: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1647: PetscCall(VecGetArray(metricInt, &met));
1648: for (i = 1; i < numMetrics; ++i) {
1649: PetscCall(VecGetArray(metrics[i], &meti));
1650: for (v = vStart; v < vEnd; ++v) {
1651: PetscCall(DMPlexPointLocalRef(dm, v, met, &M));
1652: PetscCall(DMPlexPointLocalRef(dm, v, meti, &Mi));
1653: PetscCall(DMPlexMetricIntersection_Private(nrow, Mi, M));
1654: }
1655: PetscCall(VecRestoreArray(metrics[i], &meti));
1656: }
1657: PetscCall(VecRestoreArray(metricInt, &met));
1658: }
1660: PetscCall(PetscLogEventEnd(DMPLEX_MetricIntersection, 0, 0, 0, 0));
1661: PetscFunctionReturn(PETSC_SUCCESS);
1662: }
1664: /*@
1665: DMPlexMetricIntersection2 - Compute the intersection of two metrics
1667: Input Parameters:
1668: + dm - The `DM`
1669: . metric1 - The first metric to be intersected
1670: - metric2 - The second metric to be intersected
1672: Output Parameter:
1673: . metricInt - The intersected metric
1675: Level: beginner
1677: .seealso: `DMPLEX`, `DMPlexMetricIntersection()`, `DMPlexMetricIntersection3()`
1678: @*/
1679: PetscErrorCode DMPlexMetricIntersection2(DM dm, Vec metric1, Vec metric2, Vec metricInt)
1680: {
1681: Vec metrics[2] = {metric1, metric2};
1683: PetscFunctionBegin;
1684: PetscCall(DMPlexMetricIntersection(dm, 2, metrics, metricInt));
1685: PetscFunctionReturn(PETSC_SUCCESS);
1686: }
1688: /*@
1689: DMPlexMetricIntersection3 - Compute the intersection of three metrics
1691: Input Parameters:
1692: + dm - The `DM`
1693: . metric1 - The first metric to be intersected
1694: . metric2 - The second metric to be intersected
1695: - metric3 - The third metric to be intersected
1697: Output Parameter:
1698: . metricInt - The intersected metric
1700: Level: beginner
1702: .seealso: `DMPLEX`, `DMPlexMetricIntersection()`, `DMPlexMetricIntersection2()`
1703: @*/
1704: PetscErrorCode DMPlexMetricIntersection3(DM dm, Vec metric1, Vec metric2, Vec metric3, Vec metricInt)
1705: {
1706: Vec metrics[3] = {metric1, metric2, metric3};
1708: PetscFunctionBegin;
1709: PetscCall(DMPlexMetricIntersection(dm, 3, metrics, metricInt));
1710: PetscFunctionReturn(PETSC_SUCCESS);
1711: }