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++) { 21 PetscCall(PetscRandomGetValueReal(randCtx, &preimage[i])); 22 } 23 if (dmIsSimplicial && dimR > 1) { 24 for (i = 0; i < numPoints; i++) { 25 for (j = 0; j < dimR; j++) { 26 PetscReal x = preimage[i * dimR + j]; 27 PetscReal y = preimage[i * dimR + ((j + 1) % dimR)]; 28 29 preimage[i * dimR + ((j + 1) % dimR)] = -1. + (y + 1.) * 0.5 * (1. - x); 30 } 31 } 32 } 33 34 PetscCall(DMPlexReferenceToCoordinates(dm,cell,numPoints,preimage,mapped)); 35 PetscCall(DMPlexCoordinatesToReference(dm,cell,numPoints,mapped,inverted)); 36 37 for (i = 0; i < numPoints; i++) { 38 PetscReal max = 0.; 39 40 for (j = 0; j < dimR; j++) { 41 max = PetscMax(max,PetscAbsReal(preimage[i * dimR + j] - inverted[i * dimR + j])); 42 } 43 if (max > tol) { 44 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Bad inversion for cell %" PetscInt_FMT " with error %g (tol %g): (",cell,(double)max,(double)tol)); 45 for (j = 0; j < dimR; j++) { 46 PetscCall(PetscPrintf(PETSC_COMM_SELF,"%+f",(double) preimage[i * dimR + j])); 47 if (j < dimR - 1) { 48 PetscCall(PetscPrintf(PETSC_COMM_SELF,",")); 49 } 50 } 51 PetscCall(PetscPrintf(PETSC_COMM_SELF,") --> (")); 52 for (j = 0; j < dimC; j++) { 53 PetscCall(PetscPrintf(PETSC_COMM_SELF,"%+f",(double) mapped[i * dimC + j])); 54 if (j < dimC - 1) { 55 PetscCall(PetscPrintf(PETSC_COMM_SELF,",")); 56 } 57 } 58 PetscCall(PetscPrintf(PETSC_COMM_SELF,") --> (")); 59 for (j = 0; j < dimR; j++) { 60 PetscCall(PetscPrintf(PETSC_COMM_SELF,"%+f",(double) inverted[i * dimR + j])); 61 if (j < dimR - 1) { 62 PetscCall(PetscPrintf(PETSC_COMM_SELF,",")); 63 } 64 } 65 PetscCall(PetscPrintf(PETSC_COMM_SELF,")\n")); 66 } else { 67 char strBuf[BUFSIZ] = {'\0'}; 68 size_t offset = 0, count; 69 70 PetscCall(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,"Good inversion for cell %" PetscInt_FMT ": (",&count,cell)); 71 offset += count; 72 for (j = 0; j < dimR; j++) { 73 PetscCall(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,"%+f",&count,(double) preimage[i * dimR + j])); 74 offset += count; 75 if (j < dimR - 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 < dimC; j++) { 83 PetscCall(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,"%+f",&count,(double) mapped[i * dimC + j])); 84 offset += count; 85 if (j < dimC - 1) { 86 PetscCall(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,",",&count)); 87 offset += count; 88 } 89 } 90 PetscCall(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,") --> (",&count)); 91 offset += count; 92 for (j = 0; j < dimR; j++) { 93 PetscCall(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,"%+f",&count,(double) inverted[i * dimR + j])); 94 offset += count; 95 if (j < dimR - 1) { 96 PetscCall(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,",",&count)); 97 offset += count; 98 } 99 } 100 PetscCall(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,")\n",&count)); 101 PetscCall(PetscInfo(dm,"%s",strBuf)); 102 } 103 } 104 105 PetscCall(DMRestoreWorkArray(dm,dimR * numPoints,MPIU_REAL,&inverted)); 106 PetscCall(DMRestoreWorkArray(dm,dimC * numPoints,MPIU_REAL,&mapped)); 107 PetscCall(DMRestoreWorkArray(dm,dimR * numPoints,MPIU_REAL,&preimage)); 108 PetscFunctionReturn(0); 109 } 110 111 static PetscErrorCode identityEmbedding(PetscInt dim, PetscReal time, const PetscReal *x, PetscInt Nf, PetscScalar *u, void *ctx) 112 { 113 PetscInt i; 114 115 for (i = 0; i < dim; i++) {u[i] = x[i];}; 116 return 0; 117 } 118 119 int main(int argc, char **argv) 120 { 121 PetscRandom randCtx; 122 PetscInt dim, dimC, isSimplex, isFE, numTests = 10; 123 PetscReal perturb = 0.1, tol = 10. * PETSC_SMALL; 124 125 PetscFunctionBeginUser; 126 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 127 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&randCtx)); 128 PetscCall(PetscRandomSetInterval(randCtx,-1.,1.)); 129 PetscCall(PetscRandomSetFromOptions(randCtx)); 130 PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex21",NULL); 131 PetscCall(PetscOptionsReal("-vertex_perturbation","scale of random vertex distortion",NULL,perturb,&perturb,NULL)); 132 PetscCall(PetscOptionsBoundedInt("-num_test_points","number of points to test",NULL,numTests,&numTests,NULL,0)); 133 PetscOptionsEnd(); 134 for (dim = 1; dim <= 3; dim++) { 135 for (dimC = dim; dimC <= PetscMin(3,dim + 1); dimC++) { 136 for (isSimplex = 0; isSimplex < 2; isSimplex++) { 137 for (isFE = 0; isFE < 2; isFE++) { 138 DM dm; 139 Vec coords; 140 PetscScalar *coordArray; 141 PetscReal noise; 142 PetscInt i, n, order = 1; 143 144 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex ? PETSC_TRUE : PETSC_FALSE), &dm)); 145 if (isFE) { 146 DM dmCoord; 147 PetscSpace sp; 148 PetscFE fe; 149 Vec localCoords; 150 PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *) = {identityEmbedding}; 151 void *ctxs[1] = {NULL}; 152 153 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject) dm),dim,dim,isSimplex ? PETSC_TRUE : PETSC_FALSE,isSimplex ? NULL : "tensor_" ,PETSC_DEFAULT,&fe)); 154 PetscCall(PetscFEGetBasisSpace(fe,&sp)); 155 PetscCall(PetscSpaceGetDegree(sp,&order,NULL)); 156 PetscCall(DMSetField(dm,0,NULL,(PetscObject)fe)); 157 PetscCall(DMCreateDS(dm)); 158 PetscCall(DMCreateLocalVector(dm,&localCoords)); 159 PetscCall(DMProjectFunctionLocal(dm,0,funcs,ctxs,INSERT_VALUES,localCoords)); 160 PetscCall(VecSetDM(localCoords,NULL)); /* This is necessary to prevent a reference loop */ 161 PetscCall(DMClone(dm,&dmCoord)); 162 PetscCall(DMSetField(dmCoord,0,NULL,(PetscObject)fe)); 163 PetscCall(PetscFEDestroy(&fe)); 164 PetscCall(DMCreateDS(dmCoord)); 165 PetscCall(DMSetCoordinateDM(dm,dmCoord)); 166 PetscCall(DMDestroy(&dmCoord)); 167 PetscCall(DMSetCoordinatesLocal(dm,localCoords)); 168 PetscCall(VecDestroy(&localCoords)); 169 } 170 PetscCall(PetscInfo(dm,"Testing %s%" PetscInt_FMT " %" PetscInt_FMT "D mesh embedded in %" PetscInt_FMT "D\n",isSimplex ? "P" : "Q" , order, dim, dimC)); 171 PetscCall(DMGetCoordinatesLocal(dm,&coords)); 172 PetscCall(VecGetLocalSize(coords,&n)); 173 if (dimC > dim) { /* reembed in higher dimension */ 174 PetscSection sec, newSec; 175 PetscInt pStart, pEnd, p, i, newN; 176 Vec newVec; 177 DM coordDM; 178 PetscScalar *newCoordArray; 179 180 PetscCall(DMGetCoordinateSection(dm,&sec)); 181 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)sec),&newSec)); 182 PetscCall(PetscSectionSetNumFields(newSec,1)); 183 PetscCall(PetscSectionGetChart(sec,&pStart,&pEnd)); 184 PetscCall(PetscSectionSetChart(newSec,pStart,pEnd)); 185 for (p = pStart; p < pEnd; p++) { 186 PetscInt nDof; 187 188 PetscCall(PetscSectionGetDof(sec,p,&nDof)); 189 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); 190 PetscCall(PetscSectionSetDof(newSec,p,(nDof/dim)*dimC)); 191 PetscCall(PetscSectionSetFieldDof(newSec,p,0,(nDof/dim)*dimC)); 192 } 193 PetscCall(PetscSectionSetUp(newSec)); 194 PetscCall(PetscSectionGetStorageSize(newSec,&newN)); 195 PetscCall(VecCreateSeq(PETSC_COMM_SELF,newN,&newVec)); 196 PetscCall(VecGetArray(newVec,&newCoordArray)); 197 PetscCall(VecGetArray(coords,&coordArray)); 198 for (i = 0; i < n / dim; i++) { 199 PetscInt j; 200 201 for (j = 0; j < dim; j++) { 202 newCoordArray[i * dimC + j] = coordArray[i * dim + j]; 203 } 204 for (; j < dimC; j++) { 205 newCoordArray[i * dimC + j] = 0.; 206 } 207 } 208 PetscCall(VecRestoreArray(coords,&coordArray)); 209 PetscCall(VecRestoreArray(newVec,&newCoordArray)); 210 PetscCall(DMSetCoordinateDim(dm,dimC)); 211 PetscCall(DMSetCoordinateSection(dm,dimC,newSec)); 212 PetscCall(DMSetCoordinatesLocal(dm,newVec)); 213 PetscCall(VecDestroy(&newVec)); 214 PetscCall(PetscSectionDestroy(&newSec)); 215 PetscCall(DMGetCoordinatesLocal(dm,&coords)); 216 PetscCall(DMGetCoordinateDM(dm,&coordDM)); 217 if (isFE) { 218 PetscFE fe; 219 220 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject) dm),dim,dimC,isSimplex ? PETSC_TRUE : PETSC_FALSE,isSimplex ? NULL : "tensor_",PETSC_DEFAULT,&fe)); 221 PetscCall(DMSetField(coordDM,0,NULL,(PetscObject)fe)); 222 PetscCall(PetscFEDestroy(&fe)); 223 PetscCall(DMCreateDS(coordDM)); 224 } 225 PetscCall(DMSetCoordinateDim(coordDM,dimC)); 226 PetscCall(VecGetLocalSize(coords,&n)); 227 } 228 PetscCall(VecGetArray(coords,&coordArray)); 229 for (i = 0; i < n; i++) { 230 PetscCall(PetscRandomGetValueReal(randCtx,&noise)); 231 coordArray[i] += noise * perturb; 232 } 233 PetscCall(VecRestoreArray(coords,&coordArray)); 234 PetscCall(DMSetCoordinatesLocal(dm, coords)); 235 PetscCall(testIdentity(dm, isSimplex ? PETSC_TRUE : PETSC_FALSE, 0, randCtx, numTests, tol)); 236 PetscCall(DMDestroy(&dm)); 237 } 238 } 239 } 240 } 241 PetscCall(PetscRandomDestroy(&randCtx)); 242 PetscCall(PetscFinalize()); 243 return 0; 244 } 245 246 /*TEST 247 248 test: 249 suffix: 0 250 args: -petscspace_degree 2 -tensor_petscspace_degree 2 251 252 TEST*/ 253