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;
120: /* jobz='N' (eigenvalues only), uplo='L' (lower-triangular part).
121: Column-major vs row-major does not matter for symmetric matrices. */
122: PetscCallLAPACKInfo("LAPACKsyev", LAPACKsyev_("N", "L", &n, C, &lda, w, work, &lwork, &info));
124: /* Smallest eigenvalue >= -tol? */
125: return (w[0] >= -tol) ? PETSC_TRUE : PETSC_FALSE;
126: }
128: int main(int argc, char **argv)
129: {
130: DM dm;
131: Vec m1, m2, mcap;
132: PetscMPIInt rank;
134: PetscFunctionBeginUser;
135: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
136: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
138: /* 1) Tiny 3D simplex mesh */
139: PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, 3, PETSC_TRUE, NULL, NULL, NULL, NULL, PETSC_TRUE, 0, PETSC_FALSE, &dm));
140: PetscCall(DMSetUp(dm));
142: /* Allocate metric vectors with the correct layout for this DM */
143: PetscCall(DMPlexMetricCreate(dm, 0, &m1));
144: PetscCall(DMPlexMetricCreate(dm, 0, &m2));
145: PetscCall(DMPlexMetricCreate(dm, 0, &mcap));
147: /* 2) Two constant metrics M1, M2 via Q^T D Q (aligned eigenbasis) */
148: PetscScalar a = 0.3, b = 0.7, g = 1.1; /* Euler angles (z-x-z extrinsic) */
149: PetscScalar lam1[3] = {1.0, 4.0, 9.0};
150: PetscScalar lam2[3] = {2.0, 3.0, 16.0};
152: PetscScalar M1mat[3][3], M2mat[3][3], Mcap_expected[3][3];
153: BuildMetric(a, b, g, lam1, M1mat);
154: BuildMetric(a, b, g, lam2, M2mat);
156: /* Expected intersection for aligned eigenbasis: diag(max(lam1, lam2)) */
157: PetscScalar lamcap[3] = {PetscMax(lam1[0], lam2[0]), PetscMax(lam1[1], lam2[1]), PetscMax(lam1[2], lam2[2])};
158: BuildMetric(a, b, g, lamcap, Mcap_expected);
160: PetscCall(VecZeroEntries(m1));
161: PetscCall(VecZeroEntries(m2));
162: PetscCall(SetConstantMetricVertices(dm, m1, M1mat));
163: PetscCall(SetConstantMetricVertices(dm, m2, M2mat));
165: /* 3) Intersect */
166: Vec metrics[2] = {m1, m2};
167: PetscCall(DMPlexMetricIntersection(dm, 2, metrics, mcap));
169: /* 4) Verify Loewner order and aligned-eigenbasis equality */
170: PetscInt vStart, vEnd;
171: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
172: const PetscScalar *arr;
173: PetscCall(VecGetArrayRead(mcap, &arr));
174: PetscSection sec;
175: PetscCall(DMGetLocalSection(dm, &sec));
177: for (PetscInt p = vStart; p < vEnd; ++p) {
178: PetscInt dof, off;
179: PetscCall(PetscSectionGetDof(sec, p, &dof));
180: PetscCall(PetscSectionGetOffset(sec, p, &off));
181: if (dof == 9) {
182: /* Inline read: directly from arr into Mread */
183: PetscScalar Mread[3][3];
184: for (PetscInt i = 0, t = 0; i < 3; i++)
185: for (PetscInt j = 0; j < 3; j++) Mread[i][j] = arr[off + t++];
187: /* Loewner checks */
188: PetscBool ge1 = LoewnerGE(Mread, M1mat, 1e-10);
189: PetscBool ge2 = LoewnerGE(Mread, M2mat, 1e-10);
190: PetscCheck(ge1 && ge2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Intersection not >= inputs in Loewner order");
192: /* Exact match for aligned eigenbasis */
193: PetscScalar maxdiff = 0.0;
194: for (PetscInt i = 0; i < 3; i++)
195: for (PetscInt j = 0; j < 3; j++) {
196: PetscScalar d = fabs(Mread[i][j] - Mcap_expected[i][j]);
197: if (d > maxdiff) maxdiff = d;
198: }
199: PetscCheck(maxdiff < 1e-12, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Aligned-eigenbasis formula failed at vertex (maxdiff=%.3e)", maxdiff);
200: }
201: }
202: PetscCall(VecRestoreArrayRead(mcap, &arr));
204: if (!rank) PetscCall(PetscPrintf(PETSC_COMM_SELF, "OK: Intersection passed (3x3 @ vertices, aligned case).\n"));
206: /* misaligned subtest: change Q of M2 */
207: {
208: PetscScalar a2 = a + 0.2, b2 = b - 0.1, g2 = g + 0.3;
209: BuildMetric(a2, b2, g2, lam2, M2mat);
210: PetscCall(VecZeroEntries(m2));
211: PetscCall(SetConstantMetricVertices(dm, m2, M2mat));
213: PetscCall(DMPlexMetricIntersection(dm, 2, metrics, mcap));
215: const PetscScalar *ar2;
216: PetscCall(VecGetArrayRead(mcap, &ar2));
217: for (PetscInt p = vStart; p < vEnd; ++p) {
218: PetscInt dof, off;
219: PetscCall(PetscSectionGetDof(sec, p, &dof));
220: PetscCall(PetscSectionGetOffset(sec, p, &off));
221: if (dof == 9) {
222: PetscScalar Mread[3][3];
223: for (PetscInt i = 0, t = 0; i < 3; i++)
224: for (PetscInt j = 0; j < 3; j++) Mread[i][j] = ar2[off + t++];
226: PetscBool ge1 = LoewnerGE(Mread, M1mat, 1e-10);
227: PetscBool ge2 = LoewnerGE(Mread, M2mat, 1e-10);
228: PetscCheck(ge1 && ge2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Intersection not >= inputs (misaligned)");
229: }
230: }
231: PetscCall(VecRestoreArrayRead(mcap, &ar2));
232: if (!rank) PetscCall(PetscPrintf(PETSC_COMM_SELF, "OK: Intersection passed Loewner checks (misaligned Q).\n"));
233: }
235: PetscCall(VecDestroy(&m1));
236: PetscCall(VecDestroy(&m2));
237: PetscCall(VecDestroy(&mcap));
238: PetscCall(DMDestroy(&dm));
239: PetscCall(PetscFinalize());
240: return 0;
241: }
243: /*TEST
244: build:
245: requires: !complex
246: test:
247: requires: ctetgen
249: TEST*/