1 2 const char help[] = "Test DMPlexCoordinatesToReference().\n\n"; 3 4 #include <petscdm.h> 5 #include <petscdmplex.h> 6 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; 11 12 PetscFunctionBegin; 13 PetscCall(DMGetDimension(dm, &dimR)); 14 PetscCall(DMGetCoordinateDim(dm, &dimC)); 15 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)); 19 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)]; 26 27 preimage[i * dimR + ((j + 1) % dimR)] = -1. + (y + 1.) * 0.5 * (1. - x); 28 } 29 } 30 } 31 32 PetscCall(DMPlexReferenceToCoordinates(dm, cell, numPoints, preimage, mapped)); 33 PetscCall(DMPlexCoordinatesToReference(dm, cell, numPoints, mapped, inverted)); 34 35 for (i = 0; i < numPoints; i++) { 36 PetscReal max = 0.; 37 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; 59 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 } 94 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 } 100 101 static PetscErrorCode identityEmbedding(PetscInt dim, PetscReal time, const PetscReal *x, PetscInt Nf, PetscScalar *u, void *ctx) 102 { 103 PetscInt i; 104 105 for (i = 0; i < dim; i++) u[i] = x[i]; 106 return PETSC_SUCCESS; 107 } 108 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; 114 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; 133 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}; 142 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; 169 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; 177 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; 190 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; 205 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 } 231 232 /*TEST 233 234 test: 235 suffix: 0 236 args: -petscspace_degree 2 -tensor_petscspace_degree 2 237 238 TEST*/ 239