1 const char help[] = "Test DMPlexCoordinatesToReference().\n\n";
2
3 #include <petscdm.h>
4 #include <petscdmplex.h>
5
testIdentity(DM dm,PetscBool dmIsSimplicial,PetscInt cell,PetscRandom randCtx,PetscInt numPoints,PetscReal tol)6 static PetscErrorCode testIdentity(DM dm, PetscBool dmIsSimplicial, PetscInt cell, PetscRandom randCtx, PetscInt numPoints, PetscReal tol)
7 {
8 PetscInt i, j, dimC, dimR;
9 PetscReal *preimage, *mapped, *inverted;
10
11 PetscFunctionBegin;
12 PetscCall(DMGetDimension(dm, &dimR));
13 PetscCall(DMGetCoordinateDim(dm, &dimC));
14
15 PetscCall(DMGetWorkArray(dm, dimR * numPoints, MPIU_REAL, &preimage));
16 PetscCall(DMGetWorkArray(dm, dimC * numPoints, MPIU_REAL, &mapped));
17 PetscCall(DMGetWorkArray(dm, dimR * numPoints, MPIU_REAL, &inverted));
18
19 for (i = 0; i < dimR * numPoints; i++) PetscCall(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)];
25
26 preimage[i * dimR + ((j + 1) % dimR)] = -1. + (y + 1.) * 0.5 * (1. - x);
27 }
28 }
29 }
30
31 PetscCall(DMPlexReferenceToCoordinates(dm, cell, numPoints, preimage, mapped));
32 PetscCall(DMPlexCoordinatesToReference(dm, cell, numPoints, mapped, inverted));
33
34 for (i = 0; i < numPoints; i++) {
35 PetscReal max = 0.;
36
37 for (j = 0; j < dimR; j++) max = PetscMax(max, PetscAbsReal(preimage[i * dimR + j] - inverted[i * dimR + j]));
38 if (max > tol) {
39 PetscCall(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 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%+f", (double)preimage[i * dimR + j]));
42 if (j < dimR - 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, ","));
43 }
44 PetscCall(PetscPrintf(PETSC_COMM_SELF, ") --> ("));
45 for (j = 0; j < dimC; j++) {
46 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%+f", (double)mapped[i * dimC + j]));
47 if (j < dimC - 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, ","));
48 }
49 PetscCall(PetscPrintf(PETSC_COMM_SELF, ") --> ("));
50 for (j = 0; j < dimR; j++) {
51 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%+f", (double)inverted[i * dimR + j]));
52 if (j < dimR - 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, ","));
53 }
54 PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
55 } else {
56 char strBuf[BUFSIZ] = {'\0'};
57 size_t offset = 0, count;
58
59 PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, "Good inversion for cell %" PetscInt_FMT ": (", &count, cell));
60 offset += count;
61 for (j = 0; j < dimR; j++) {
62 PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, "%+f", &count, (double)preimage[i * dimR + j]));
63 offset += count;
64 if (j < dimR - 1) {
65 PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, ",", &count));
66 offset += count;
67 }
68 }
69 PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, ") --> (", &count));
70 offset += count;
71 for (j = 0; j < dimC; j++) {
72 PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, "%+f", &count, (double)mapped[i * dimC + j]));
73 offset += count;
74 if (j < dimC - 1) {
75 PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, ",", &count));
76 offset += count;
77 }
78 }
79 PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, ") --> (", &count));
80 offset += count;
81 for (j = 0; j < dimR; j++) {
82 PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, "%+f", &count, (double)inverted[i * dimR + j]));
83 offset += count;
84 if (j < dimR - 1) {
85 PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, ",", &count));
86 offset += count;
87 }
88 }
89 PetscCall(PetscSNPrintfCount(strBuf + offset, BUFSIZ - offset, ")", &count));
90 PetscCall(PetscInfo(dm, "%s\n", strBuf));
91 }
92 }
93
94 PetscCall(DMRestoreWorkArray(dm, dimR * numPoints, MPIU_REAL, &inverted));
95 PetscCall(DMRestoreWorkArray(dm, dimC * numPoints, MPIU_REAL, &mapped));
96 PetscCall(DMRestoreWorkArray(dm, dimR * numPoints, MPIU_REAL, &preimage));
97 PetscFunctionReturn(PETSC_SUCCESS);
98 }
99
identityEmbedding(PetscInt dim,PetscReal time,const PetscReal * x,PetscInt Nf,PetscScalar * u,PetscCtx ctx)100 static PetscErrorCode identityEmbedding(PetscInt dim, PetscReal time, const PetscReal *x, PetscInt Nf, PetscScalar *u, PetscCtx ctx)
101 {
102 PetscInt i;
103
104 for (i = 0; i < dim; i++) u[i] = x[i];
105 return PETSC_SUCCESS;
106 }
107
main(int argc,char ** argv)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;
113
114 PetscFunctionBeginUser;
115 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
116 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &randCtx));
117 PetscCall(PetscRandomSetInterval(randCtx, -1., 1.));
118 PetscCall(PetscRandomSetFromOptions(randCtx));
119 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex21", NULL);
120 PetscCall(PetscOptionsReal("-vertex_perturbation", "scale of random vertex distortion", NULL, perturb, &perturb, NULL));
121 PetscCall(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;
132
133 PetscCall(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 PetscCtx ctxs[1] = {NULL};
141
142 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, dim, isSimplex ? PETSC_TRUE : PETSC_FALSE, isSimplex ? NULL : "tensor_", PETSC_DEFAULT, &fe));
143 PetscCall(PetscFEGetBasisSpace(fe, &sp));
144 PetscCall(PetscSpaceGetDegree(sp, &order, NULL));
145 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
146 PetscCall(DMCreateDS(dm));
147 PetscCall(DMCreateLocalVector(dm, &localCoords));
148 PetscCall(DMProjectFunctionLocal(dm, 0, funcs, ctxs, INSERT_VALUES, localCoords));
149 PetscCall(VecSetDM(localCoords, NULL)); /* This is necessary to prevent a reference loop */
150 PetscCall(DMClone(dm, &dmCoord));
151 PetscCall(DMSetField(dmCoord, 0, NULL, (PetscObject)fe));
152 PetscCall(PetscFEDestroy(&fe));
153 PetscCall(DMCreateDS(dmCoord));
154 PetscCall(DMSetCoordinateDM(dm, dmCoord));
155 PetscCall(DMDestroy(&dmCoord));
156 PetscCall(DMSetCoordinatesLocal(dm, localCoords));
157 PetscCall(VecDestroy(&localCoords));
158 }
159 PetscCall(PetscInfo(dm, "Testing %s%" PetscInt_FMT " %" PetscInt_FMT "D mesh embedded in %" PetscInt_FMT "D\n", isSimplex ? "P" : "Q", order, dim, dimC));
160 PetscCall(DMGetCoordinatesLocal(dm, &coords));
161 PetscCall(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;
168
169 PetscCall(DMGetCoordinateSection(dm, &sec));
170 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)sec), &newSec));
171 PetscCall(PetscSectionSetNumFields(newSec, 1));
172 PetscCall(PetscSectionGetChart(sec, &pStart, &pEnd));
173 PetscCall(PetscSectionSetChart(newSec, pStart, pEnd));
174 for (p = pStart; p < pEnd; p++) {
175 PetscInt nDof;
176
177 PetscCall(PetscSectionGetDof(sec, p, &nDof));
178 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);
179 PetscCall(PetscSectionSetDof(newSec, p, (nDof / dim) * dimC));
180 PetscCall(PetscSectionSetFieldDof(newSec, p, 0, (nDof / dim) * dimC));
181 }
182 PetscCall(PetscSectionSetUp(newSec));
183 PetscCall(PetscSectionGetStorageSize(newSec, &newN));
184 PetscCall(VecCreateSeq(PETSC_COMM_SELF, newN, &newVec));
185 PetscCall(VecGetArray(newVec, &newCoordArray));
186 PetscCall(VecGetArray(coords, &coordArray));
187 for (i = 0; i < n / dim; i++) {
188 PetscInt j;
189
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 PetscCall(VecRestoreArray(coords, &coordArray));
194 PetscCall(VecRestoreArray(newVec, &newCoordArray));
195 PetscCall(DMSetCoordinateDim(dm, dimC));
196 PetscCall(DMSetCoordinateSection(dm, dimC, newSec));
197 PetscCall(DMSetCoordinatesLocal(dm, newVec));
198 PetscCall(VecDestroy(&newVec));
199 PetscCall(PetscSectionDestroy(&newSec));
200 PetscCall(DMGetCoordinatesLocal(dm, &coords));
201 PetscCall(DMGetCoordinateDM(dm, &coordDM));
202 if (isFE) {
203 PetscFE fe;
204
205 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dm), dim, dimC, isSimplex ? PETSC_TRUE : PETSC_FALSE, isSimplex ? NULL : "tensor_", PETSC_DEFAULT, &fe));
206 PetscCall(DMSetField(coordDM, 0, NULL, (PetscObject)fe));
207 PetscCall(PetscFEDestroy(&fe));
208 PetscCall(DMCreateDS(coordDM));
209 }
210 PetscCall(DMSetCoordinateDim(coordDM, dimC));
211 PetscCall(VecGetLocalSize(coords, &n));
212 }
213 PetscCall(VecGetArray(coords, &coordArray));
214 for (i = 0; i < n; i++) {
215 PetscCall(PetscRandomGetValueReal(randCtx, &noise));
216 coordArray[i] += noise * perturb;
217 }
218 PetscCall(VecRestoreArray(coords, &coordArray));
219 PetscCall(DMSetCoordinatesLocal(dm, coords));
220 PetscCall(testIdentity(dm, isSimplex ? PETSC_TRUE : PETSC_FALSE, 0, randCtx, numTests, tol));
221 PetscCall(DMDestroy(&dm));
222 }
223 }
224 }
225 }
226 PetscCall(PetscRandomDestroy(&randCtx));
227 PetscCall(PetscFinalize());
228 return 0;
229 }
230
231 /*TEST
232
233 test:
234 suffix: 0
235 args: -petscspace_degree 2 -tensor_petscspace_degree 2
236 output_file: output/empty.out
237
238 TEST*/
239