Actual source code: ex22.c

  1: const char help[] = "Test DMPlexCoordinatesToReference().\n\n";

  3: #include <petscdm.h>
  4: #include <petscdmplex.h>

  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;

 11:   PetscFunctionBegin;
 12:   PetscCall(DMGetDimension(dm, &dimR));
 13:   PetscCall(DMGetCoordinateDim(dm, &dimC));

 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));

 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)];

 26:         preimage[i * dimR + ((j + 1) % dimR)] = -1. + (y + 1.) * 0.5 * (1. - x);
 27:       }
 28:     }
 29:   }

 31:   PetscCall(DMPlexReferenceToCoordinates(dm, cell, numPoints, preimage, mapped));
 32:   PetscCall(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:       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;

 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:   }

 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: }

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 PETSC_SUCCESS;
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;

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;

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:             void *ctxs[1]                                                                                       = {NULL};

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;

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;

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;

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;

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: }

231: /*TEST

233:   test:
234:     suffix: 0
235:     args: -petscspace_degree 2 -tensor_petscspace_degree 2

237: TEST*/