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