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