Actual source code: ex22.c
2: const char help[] = "Test DMPlexCoordinatesToReference().\n\n";
4: #include <petscdm.h>
5: #include <petscdmplex.h>
7: static PetscErrorCode testIdentity(DM dm, PetscBool dmIsSimplicial, PetscInt cell, PetscRandom randCtx, PetscInt numPoints, PetscReal tol)
8: {
9: PetscInt i, j, dimC, dimR;
10: PetscReal *preimage, *mapped, *inverted;
12: PetscFunctionBegin;
13: PetscCall(DMGetDimension(dm, &dimR));
14: PetscCall(DMGetCoordinateDim(dm, &dimC));
16: PetscCall(DMGetWorkArray(dm, dimR * numPoints, MPIU_REAL, &preimage));
17: PetscCall(DMGetWorkArray(dm, dimC * numPoints, MPIU_REAL, &mapped));
18: PetscCall(DMGetWorkArray(dm, dimR * numPoints, MPIU_REAL, &inverted));
20: for (i = 0; i < dimR * numPoints; i++) PetscCall(PetscRandomGetValueReal(randCtx, &preimage[i]));
21: if (dmIsSimplicial && dimR > 1) {
22: for (i = 0; i < numPoints; i++) {
23: for (j = 0; j < dimR; j++) {
24: PetscReal x = preimage[i * dimR + j];
25: PetscReal y = preimage[i * dimR + ((j + 1) % dimR)];
27: preimage[i * dimR + ((j + 1) % dimR)] = -1. + (y + 1.) * 0.5 * (1. - x);
28: }
29: }
30: }
32: PetscCall(DMPlexReferenceToCoordinates(dm, cell, numPoints, preimage, mapped));
33: PetscCall(DMPlexCoordinatesToReference(dm, cell, numPoints, mapped, inverted));
35: for (i = 0; i < numPoints; i++) {
36: PetscReal max = 0.;
38: for (j = 0; j < dimR; j++) max = PetscMax(max, PetscAbsReal(preimage[i * dimR + j] - inverted[i * dimR + j]));
39: if (max > tol) {
40: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Bad inversion for cell %" PetscInt_FMT " with error %g (tol %g): (", cell, (double)max, (double)tol));
41: for (j = 0; j < dimR; j++) {
42: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%+f", (double)preimage[i * dimR + j]));
43: if (j < dimR - 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, ","));
44: }
45: PetscCall(PetscPrintf(PETSC_COMM_SELF, ") --> ("));
46: for (j = 0; j < dimC; j++) {
47: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%+f", (double)mapped[i * dimC + j]));
48: if (j < dimC - 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, ","));
49: }
50: PetscCall(PetscPrintf(PETSC_COMM_SELF, ") --> ("));
51: for (j = 0; j < dimR; j++) {
52: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%+f", (double)inverted[i * dimR + j]));
53: if (j < dimR - 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, ","));
54: }
55: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
56: } else {
57: char strBuf[BUFSIZ] = {'\0'};
58: size_t offset = 0, count;
60: PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, "Good inversion for cell %" PetscInt_FMT ": (", &count, cell));
61: offset += count;
62: for (j = 0; j < dimR; j++) {
63: PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, "%+f", &count, (double)preimage[i * dimR + j]));
64: offset += count;
65: if (j < dimR - 1) {
66: PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, ",", &count));
67: offset += count;
68: }
69: }
70: PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, ") --> (", &count));
71: offset += count;
72: for (j = 0; j < dimC; j++) {
73: PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, "%+f", &count, (double)mapped[i * dimC + j]));
74: offset += count;
75: if (j < dimC - 1) {
76: PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, ",", &count));
77: offset += count;
78: }
79: }
80: PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, ") --> (", &count));
81: offset += count;
82: for (j = 0; j < dimR; j++) {
83: PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, "%+f", &count, (double)inverted[i * dimR + j]));
84: offset += count;
85: if (j < dimR - 1) {
86: PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, ",", &count));
87: offset += count;
88: }
89: }
90: PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, ")\n", &count));
91: PetscCall(PetscInfo(dm, "%s", strBuf));
92: }
93: }
95: PetscCall(DMRestoreWorkArray(dm, dimR * numPoints, MPIU_REAL, &inverted));
96: PetscCall(DMRestoreWorkArray(dm, dimC * numPoints, MPIU_REAL, &mapped));
97: PetscCall(DMRestoreWorkArray(dm, dimR * numPoints, MPIU_REAL, &preimage));
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: static PetscErrorCode identityEmbedding(PetscInt dim, PetscReal time, const PetscReal *x, PetscInt Nf, PetscScalar *u, void *ctx)
102: {
103: PetscInt i;
105: for (i = 0; i < dim; i++) u[i] = x[i];
106: return PETSC_SUCCESS;
107: }
109: int main(int argc, char **argv)
110: {
111: PetscRandom randCtx;
112: PetscInt dim, dimC, isSimplex, isFE, numTests = 10;
113: PetscReal perturb = 0.1, tol = 10. * PETSC_SMALL;
115: PetscFunctionBeginUser;
116: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
117: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &randCtx));
118: PetscCall(PetscRandomSetInterval(randCtx, -1., 1.));
119: PetscCall(PetscRandomSetFromOptions(randCtx));
120: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex21", NULL);
121: PetscCall(PetscOptionsReal("-vertex_perturbation", "scale of random vertex distortion", NULL, perturb, &perturb, NULL));
122: PetscCall(PetscOptionsBoundedInt("-num_test_points", "number of points to test", NULL, numTests, &numTests, NULL, 0));
123: PetscOptionsEnd();
124: for (dim = 1; dim <= 3; dim++) {
125: for (dimC = dim; dimC <= PetscMin(3, dim + 1); dimC++) {
126: for (isSimplex = 0; isSimplex < 2; isSimplex++) {
127: for (isFE = 0; isFE < 2; isFE++) {
128: DM dm;
129: Vec coords;
130: PetscScalar *coordArray;
131: PetscReal noise;
132: PetscInt i, n, order = 1;
134: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex ? PETSC_TRUE : PETSC_FALSE), &dm));
135: if (isFE) {
136: DM dmCoord;
137: PetscSpace sp;
138: PetscFE fe;
139: Vec localCoords;
140: PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *) = {identityEmbedding};
141: void *ctxs[1] = {NULL};
143: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, dim, isSimplex ? PETSC_TRUE : PETSC_FALSE, isSimplex ? NULL : "tensor_", PETSC_DEFAULT, &fe));
144: PetscCall(PetscFEGetBasisSpace(fe, &sp));
145: PetscCall(PetscSpaceGetDegree(sp, &order, NULL));
146: PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
147: PetscCall(DMCreateDS(dm));
148: PetscCall(DMCreateLocalVector(dm, &localCoords));
149: PetscCall(DMProjectFunctionLocal(dm, 0, funcs, ctxs, INSERT_VALUES, localCoords));
150: PetscCall(VecSetDM(localCoords, NULL)); /* This is necessary to prevent a reference loop */
151: PetscCall(DMClone(dm, &dmCoord));
152: PetscCall(DMSetField(dmCoord, 0, NULL, (PetscObject)fe));
153: PetscCall(PetscFEDestroy(&fe));
154: PetscCall(DMCreateDS(dmCoord));
155: PetscCall(DMSetCoordinateDM(dm, dmCoord));
156: PetscCall(DMDestroy(&dmCoord));
157: PetscCall(DMSetCoordinatesLocal(dm, localCoords));
158: PetscCall(VecDestroy(&localCoords));
159: }
160: PetscCall(PetscInfo(dm, "Testing %s%" PetscInt_FMT " %" PetscInt_FMT "D mesh embedded in %" PetscInt_FMT "D\n", isSimplex ? "P" : "Q", order, dim, dimC));
161: PetscCall(DMGetCoordinatesLocal(dm, &coords));
162: PetscCall(VecGetLocalSize(coords, &n));
163: if (dimC > dim) { /* reembed in higher dimension */
164: PetscSection sec, newSec;
165: PetscInt pStart, pEnd, p, i, newN;
166: Vec newVec;
167: DM coordDM;
168: PetscScalar *newCoordArray;
170: PetscCall(DMGetCoordinateSection(dm, &sec));
171: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)sec), &newSec));
172: PetscCall(PetscSectionSetNumFields(newSec, 1));
173: PetscCall(PetscSectionGetChart(sec, &pStart, &pEnd));
174: PetscCall(PetscSectionSetChart(newSec, pStart, pEnd));
175: for (p = pStart; p < pEnd; p++) {
176: PetscInt nDof;
178: PetscCall(PetscSectionGetDof(sec, p, &nDof));
179: PetscCheck(nDof % dim == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Coordinate section point %" PetscInt_FMT " has %" PetscInt_FMT " dofs != 0 mod %" PetscInt_FMT, p, nDof, dim);
180: PetscCall(PetscSectionSetDof(newSec, p, (nDof / dim) * dimC));
181: PetscCall(PetscSectionSetFieldDof(newSec, p, 0, (nDof / dim) * dimC));
182: }
183: PetscCall(PetscSectionSetUp(newSec));
184: PetscCall(PetscSectionGetStorageSize(newSec, &newN));
185: PetscCall(VecCreateSeq(PETSC_COMM_SELF, newN, &newVec));
186: PetscCall(VecGetArray(newVec, &newCoordArray));
187: PetscCall(VecGetArray(coords, &coordArray));
188: for (i = 0; i < n / dim; i++) {
189: PetscInt j;
191: for (j = 0; j < dim; j++) newCoordArray[i * dimC + j] = coordArray[i * dim + j];
192: for (; j < dimC; j++) newCoordArray[i * dimC + j] = 0.;
193: }
194: PetscCall(VecRestoreArray(coords, &coordArray));
195: PetscCall(VecRestoreArray(newVec, &newCoordArray));
196: PetscCall(DMSetCoordinateDim(dm, dimC));
197: PetscCall(DMSetCoordinateSection(dm, dimC, newSec));
198: PetscCall(DMSetCoordinatesLocal(dm, newVec));
199: PetscCall(VecDestroy(&newVec));
200: PetscCall(PetscSectionDestroy(&newSec));
201: PetscCall(DMGetCoordinatesLocal(dm, &coords));
202: PetscCall(DMGetCoordinateDM(dm, &coordDM));
203: if (isFE) {
204: PetscFE fe;
206: PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, dimC, isSimplex ? PETSC_TRUE : PETSC_FALSE, isSimplex ? NULL : "tensor_", PETSC_DEFAULT, &fe));
207: PetscCall(DMSetField(coordDM, 0, NULL, (PetscObject)fe));
208: PetscCall(PetscFEDestroy(&fe));
209: PetscCall(DMCreateDS(coordDM));
210: }
211: PetscCall(DMSetCoordinateDim(coordDM, dimC));
212: PetscCall(VecGetLocalSize(coords, &n));
213: }
214: PetscCall(VecGetArray(coords, &coordArray));
215: for (i = 0; i < n; i++) {
216: PetscCall(PetscRandomGetValueReal(randCtx, &noise));
217: coordArray[i] += noise * perturb;
218: }
219: PetscCall(VecRestoreArray(coords, &coordArray));
220: PetscCall(DMSetCoordinatesLocal(dm, coords));
221: PetscCall(testIdentity(dm, isSimplex ? PETSC_TRUE : PETSC_FALSE, 0, randCtx, numTests, tol));
222: PetscCall(DMDestroy(&dm));
223: }
224: }
225: }
226: }
227: PetscCall(PetscRandomDestroy(&randCtx));
228: PetscCall(PetscFinalize());
229: return 0;
230: }
232: /*TEST
234: test:
235: suffix: 0
236: args: -petscspace_degree 2 -tensor_petscspace_degree 2
238: TEST*/