1 const char help[] = "Test DMPlexCoordinatesToReference().\n\n"; 2 3 #include <petscdm.h> 4 #include <petscdmplex.h> 5 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 100 static PetscErrorCode identityEmbedding(PetscInt dim, PetscReal time, const PetscReal *x, PetscInt Nf, PetscScalar *u, void *ctx) 101 { 102 PetscInt i; 103 104 for (i = 0; i < dim; i++) u[i] = x[i]; 105 return PETSC_SUCCESS; 106 } 107 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 void *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 237 TEST*/ 238