Actual source code: ex76.c
1: static char help[] = "Test for DMPlexMetricIntersection using two constant 3x3 metrics.\n\n";
3: #include <petscdmplex.h>
4: #include <petscblaslapack.h>
6: /* Euler z-x-z (extrinsic) rotation same as used in DMPlexCreateBasisRotation
7: R = Rz(alpha) * Rx(beta) * Rz(gamma) */
9: static void RotationZ(PetscScalar a, PetscScalar R[3][3])
10: {
11: PetscScalar c = PetscCosScalar(a), s = PetscSinScalar(a);
12: R[0][0] = c;
13: R[0][1] = -s;
14: R[0][2] = 0.;
15: R[1][0] = s;
16: R[1][1] = c;
17: R[1][2] = 0.;
18: R[2][0] = 0.;
19: R[2][1] = 0.;
20: R[2][2] = 1.;
21: }
23: static void RotationX(PetscScalar b, PetscScalar R[3][3])
24: {
25: PetscScalar c = PetscCosScalar(b), s = PetscSinScalar(b);
26: R[0][0] = 1.;
27: R[0][1] = 0.;
28: R[0][2] = 0.;
29: R[1][0] = 0.;
30: R[1][1] = c;
31: R[1][2] = -s;
32: R[2][0] = 0.;
33: R[2][1] = s;
34: R[2][2] = c;
35: }
37: static void MatMult3(const PetscScalar A[3][3], const PetscScalar B[3][3], PetscScalar C[3][3])
38: {
39: for (PetscInt i = 0; i < 3; i++)
40: for (PetscInt j = 0; j < 3; j++) {
41: C[i][j] = 0.0;
42: for (PetscInt k = 0; k < 3; k++) C[i][j] += A[i][k] * B[k][j];
43: }
44: }
46: static void EulerZXZ(PetscScalar a, PetscScalar b, PetscScalar g, PetscScalar Q[3][3])
47: {
48: PetscScalar Rz1[3][3], Rx[3][3], Rz2[3][3], T[3][3];
49: RotationZ(a, Rz1);
50: RotationX(b, Rx);
51: RotationZ(g, Rz2);
52: MatMult3(Rz1, Rx, T);
53: MatMult3(T, Rz2, Q);
54: }
56: /* Build metric M = Q^T D Q given Euler and eigenvalues */
57: static void BuildMetric(PetscScalar a, PetscScalar b, PetscScalar g, const PetscScalar lam[3], PetscScalar M[3][3])
58: {
59: PetscScalar Q[3][3], QT[3][3], DQ[3][3];
60: EulerZXZ(a, b, g, Q);
61: for (PetscInt i = 0; i < 3; i++)
62: for (PetscInt j = 0; j < 3; j++) QT[i][j] = Q[j][i];
63: for (PetscInt i = 0; i < 3; i++)
64: for (PetscInt j = 0; j < 3; j++) DQ[i][j] = lam[i] * Q[i][j];
65: MatMult3(QT, DQ, M);
66: }
68: /* Write a constant 3x3 metric into the metric Vec */
69: static PetscErrorCode SetConstantMetricVertices(DM dm, Vec metric, const PetscScalar M[3][3])
70: {
71: PetscFunctionBeginUser;
72: PetscInt vStart, vEnd;
73: PetscSection sec;
74: Vec loc;
75: PetscScalar *lptr;
77: PetscCall(DMGetLocalSection(dm, &sec));
78: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
80: PetscCall(DMGetLocalVector(dm, &loc));
81: PetscCall(VecZeroEntries(loc));
82: PetscCall(VecGetArray(loc, &lptr));
84: for (PetscInt p = vStart; p < vEnd; ++p) {
85: PetscInt dof, off;
86: PetscCall(PetscSectionGetDof(sec, p, &dof));
87: PetscCall(PetscSectionGetOffset(sec, p, &off));
88: if (dof == 9) {
89: PetscInt t = 0;
90: for (PetscInt i = 0; i < 3; i++)
91: for (PetscInt j = 0; j < 3; j++) lptr[off + t++] = M[i][j];
92: } else if (dof > 0) {
93: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected dof count %" PetscInt_FMT, dof);
94: }
95: }
97: PetscCall(VecRestoreArray(loc, &lptr));
98: PetscCall(DMLocalToGlobal(dm, loc, INSERT_VALUES, metric));
99: PetscCall(DMRestoreLocalVector(dm, &loc));
100: PetscFunctionReturn(PETSC_SUCCESS);
101: }
103: /* Loewner (semidefinite) inequality check via eigenvalues of A-B, using LAPACK syev
104: Loewner inequality: A >= B iff A-B is positive semi-definite */
105: static PetscBool LoewnerGE(const PetscScalar A[3][3], const PetscScalar B[3][3], PetscScalar tol)
106: {
107: /* Build C = A - B in flat buffer */
108: PetscScalar C[9];
109: {
110: PetscInt t = 0;
111: for (PetscInt i = 0; i < 3; i++)
112: for (PetscInt j = 0; j < 3; j++) C[t++] = A[i][j] - B[i][j];
113: }
115: PetscScalar w[3]; /* eigenvalues */
116: PetscScalar work[9]; /* minimal workspace for 3x3 */
117: PetscBLASInt n = 3, lda = 3, lwork = 9;
119: /* jobz='N' (eigenvalues only), uplo='L' (lower-triangular part).
120: Column-major vs row-major does not matter for symmetric matrices. */
121: PetscCallLAPACKInfo("LAPACKsyev", LAPACKsyev_("N", "L", &n, C, &lda, w, work, &lwork, &info));
123: /* Smallest eigenvalue >= -tol? */
124: return (w[0] >= -tol) ? PETSC_TRUE : PETSC_FALSE;
125: }
127: int main(int argc, char **argv)
128: {
129: DM dm;
130: Vec m1, m2, mcap;
131: PetscMPIInt rank;
133: PetscFunctionBeginUser;
134: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
135: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
137: /* 1) Tiny 3D simplex mesh */
138: PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, 3, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, 0, PETSC_FALSE, &dm));
139: PetscCall(DMSetUp(dm));
141: /* Allocate metric vectors with the correct layout for this DM */
142: PetscCall(DMPlexMetricCreate(dm, 0, &m1));
143: PetscCall(DMPlexMetricCreate(dm, 0, &m2));
144: PetscCall(DMPlexMetricCreate(dm, 0, &mcap));
146: /* 2) Two constant metrics M1, M2 via Q^T D Q (aligned eigenbasis) */
147: PetscScalar a = 0.3, b = 0.7, g = 1.1; /* Euler angles (z-x-z extrinsic) */
148: PetscScalar lam1[3] = {1.0, 4.0, 9.0};
149: PetscScalar lam2[3] = {2.0, 3.0, 16.0};
151: PetscScalar M1mat[3][3], M2mat[3][3], Mcap_expected[3][3];
152: BuildMetric(a, b, g, lam1, M1mat);
153: BuildMetric(a, b, g, lam2, M2mat);
155: /* Expected intersection for aligned eigenbasis: diag(max(lam1, lam2)) */
156: PetscScalar lamcap[3] = {PetscMax(lam1[0], lam2[0]), PetscMax(lam1[1], lam2[1]), PetscMax(lam1[2], lam2[2])};
157: BuildMetric(a, b, g, lamcap, Mcap_expected);
159: PetscCall(VecZeroEntries(m1));
160: PetscCall(VecZeroEntries(m2));
161: PetscCall(SetConstantMetricVertices(dm, m1, M1mat));
162: PetscCall(SetConstantMetricVertices(dm, m2, M2mat));
164: /* 3) Intersect */
165: Vec metrics[2] = {m1, m2};
166: PetscCall(DMPlexMetricIntersection(dm, 2, metrics, mcap));
168: /* 4) Verify Loewner order and aligned-eigenbasis equality */
169: PetscInt vStart, vEnd;
170: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
171: const PetscScalar *arr;
172: PetscCall(VecGetArrayRead(mcap, &arr));
173: PetscSection sec;
174: PetscCall(DMGetLocalSection(dm, &sec));
176: for (PetscInt p = vStart; p < vEnd; ++p) {
177: PetscInt dof, off;
178: PetscCall(PetscSectionGetDof(sec, p, &dof));
179: PetscCall(PetscSectionGetOffset(sec, p, &off));
180: if (dof == 9) {
181: /* Inline read: directly from arr into Mread */
182: PetscScalar Mread[3][3];
183: for (PetscInt i = 0, t = 0; i < 3; i++)
184: for (PetscInt j = 0; j < 3; j++) Mread[i][j] = arr[off + t++];
186: /* Loewner checks */
187: PetscBool ge1 = LoewnerGE(Mread, M1mat, 1e-10);
188: PetscBool ge2 = LoewnerGE(Mread, M2mat, 1e-10);
189: PetscCheck(ge1 && ge2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Intersection not >= inputs in Loewner order");
191: /* Exact match for aligned eigenbasis */
192: PetscScalar maxdiff = 0.0;
193: for (PetscInt i = 0; i < 3; i++)
194: for (PetscInt j = 0; j < 3; j++) {
195: PetscScalar d = fabs(Mread[i][j] - Mcap_expected[i][j]);
196: if (d > maxdiff) maxdiff = d;
197: }
198: PetscCheck(maxdiff < 1e-12, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Aligned-eigenbasis formula failed at vertex (maxdiff=%.3e)", maxdiff);
199: }
200: }
201: PetscCall(VecRestoreArrayRead(mcap, &arr));
203: if (!rank) PetscCall(PetscPrintf(PETSC_COMM_SELF, "OK: Intersection passed (3x3 @ vertices, aligned case).\n"));
205: /* misaligned subtest: change Q of M2 */
206: {
207: PetscScalar a2 = a + 0.2, b2 = b - 0.1, g2 = g + 0.3;
208: BuildMetric(a2, b2, g2, lam2, M2mat);
209: PetscCall(VecZeroEntries(m2));
210: PetscCall(SetConstantMetricVertices(dm, m2, M2mat));
212: PetscCall(DMPlexMetricIntersection(dm, 2, metrics, mcap));
214: const PetscScalar *ar2;
215: PetscCall(VecGetArrayRead(mcap, &ar2));
216: for (PetscInt p = vStart; p < vEnd; ++p) {
217: PetscInt dof, off;
218: PetscCall(PetscSectionGetDof(sec, p, &dof));
219: PetscCall(PetscSectionGetOffset(sec, p, &off));
220: if (dof == 9) {
221: PetscScalar Mread[3][3];
222: for (PetscInt i = 0, t = 0; i < 3; i++)
223: for (PetscInt j = 0; j < 3; j++) Mread[i][j] = ar2[off + t++];
225: PetscBool ge1 = LoewnerGE(Mread, M1mat, 1e-10);
226: PetscBool ge2 = LoewnerGE(Mread, M2mat, 1e-10);
227: PetscCheck(ge1 && ge2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Intersection not >= inputs (misaligned)");
228: }
229: }
230: PetscCall(VecRestoreArrayRead(mcap, &ar2));
231: if (!rank) PetscCall(PetscPrintf(PETSC_COMM_SELF, "OK: Intersection passed Loewner checks (misaligned Q).\n"));
232: }
234: PetscCall(VecDestroy(&m1));
235: PetscCall(VecDestroy(&m2));
236: PetscCall(VecDestroy(&mcap));
237: PetscCall(DMDestroy(&dm));
238: PetscCall(PetscFinalize());
239: return 0;
240: }
242: /*TEST
243: build:
244: requires: !complex
245: test:
246: requires: ctetgen
248: TEST*/