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*/