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 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 126 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&randCtx)); 127 PetscCall(PetscRandomSetInterval(randCtx,-1.,1.)); 128 PetscCall(PetscRandomSetFromOptions(randCtx)); 129 PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex21",NULL); 130 PetscCall(PetscOptionsReal("-vertex_perturbation","scale of random vertex distortion",NULL,perturb,&perturb,NULL)); 131 PetscCall(PetscOptionsBoundedInt("-num_test_points","number of points to test",NULL,numTests,&numTests,NULL,0)); 132 PetscOptionsEnd(); 133 for (dim = 1; dim <= 3; dim++) { 134 for (dimC = dim; dimC <= PetscMin(3,dim + 1); dimC++) { 135 for (isSimplex = 0; isSimplex < 2; isSimplex++) { 136 for (isFE = 0; isFE < 2; isFE++) { 137 DM dm; 138 Vec coords; 139 PetscScalar *coordArray; 140 PetscReal noise; 141 PetscInt i, n, order = 1; 142 143 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex ? PETSC_TRUE : PETSC_FALSE), &dm)); 144 if (isFE) { 145 DM dmCoord; 146 PetscSpace sp; 147 PetscFE fe; 148 Vec localCoords; 149 PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *) = {identityEmbedding}; 150 void *ctxs[1] = {NULL}; 151 152 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject) dm),dim,dim,isSimplex ? PETSC_TRUE : PETSC_FALSE,isSimplex ? NULL : "tensor_" ,PETSC_DEFAULT,&fe)); 153 PetscCall(PetscFEGetBasisSpace(fe,&sp)); 154 PetscCall(PetscSpaceGetDegree(sp,&order,NULL)); 155 PetscCall(DMSetField(dm,0,NULL,(PetscObject)fe)); 156 PetscCall(DMCreateDS(dm)); 157 PetscCall(DMCreateLocalVector(dm,&localCoords)); 158 PetscCall(DMProjectFunctionLocal(dm,0,funcs,ctxs,INSERT_VALUES,localCoords)); 159 PetscCall(VecSetDM(localCoords,NULL)); /* This is necessary to prevent a reference loop */ 160 PetscCall(DMClone(dm,&dmCoord)); 161 PetscCall(DMSetField(dmCoord,0,NULL,(PetscObject)fe)); 162 PetscCall(PetscFEDestroy(&fe)); 163 PetscCall(DMCreateDS(dmCoord)); 164 PetscCall(DMSetCoordinateDM(dm,dmCoord)); 165 PetscCall(DMDestroy(&dmCoord)); 166 PetscCall(DMSetCoordinatesLocal(dm,localCoords)); 167 PetscCall(VecDestroy(&localCoords)); 168 } 169 PetscCall(PetscInfo(dm,"Testing %s%" PetscInt_FMT " %" PetscInt_FMT "D mesh embedded in %" PetscInt_FMT "D\n",isSimplex ? "P" : "Q" , order, dim, dimC)); 170 PetscCall(DMGetCoordinatesLocal(dm,&coords)); 171 PetscCall(VecGetLocalSize(coords,&n)); 172 if (dimC > dim) { /* reembed in higher dimension */ 173 PetscSection sec, newSec; 174 PetscInt pStart, pEnd, p, i, newN; 175 Vec newVec; 176 DM coordDM; 177 PetscScalar *newCoordArray; 178 179 PetscCall(DMGetCoordinateSection(dm,&sec)); 180 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)sec),&newSec)); 181 PetscCall(PetscSectionSetNumFields(newSec,1)); 182 PetscCall(PetscSectionGetChart(sec,&pStart,&pEnd)); 183 PetscCall(PetscSectionSetChart(newSec,pStart,pEnd)); 184 for (p = pStart; p < pEnd; p++) { 185 PetscInt nDof; 186 187 PetscCall(PetscSectionGetDof(sec,p,&nDof)); 188 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); 189 PetscCall(PetscSectionSetDof(newSec,p,(nDof/dim)*dimC)); 190 PetscCall(PetscSectionSetFieldDof(newSec,p,0,(nDof/dim)*dimC)); 191 } 192 PetscCall(PetscSectionSetUp(newSec)); 193 PetscCall(PetscSectionGetStorageSize(newSec,&newN)); 194 PetscCall(VecCreateSeq(PETSC_COMM_SELF,newN,&newVec)); 195 PetscCall(VecGetArray(newVec,&newCoordArray)); 196 PetscCall(VecGetArray(coords,&coordArray)); 197 for (i = 0; i < n / dim; i++) { 198 PetscInt j; 199 200 for (j = 0; j < dim; j++) { 201 newCoordArray[i * dimC + j] = coordArray[i * dim + j]; 202 } 203 for (; j < dimC; j++) { 204 newCoordArray[i * dimC + j] = 0.; 205 } 206 } 207 PetscCall(VecRestoreArray(coords,&coordArray)); 208 PetscCall(VecRestoreArray(newVec,&newCoordArray)); 209 PetscCall(DMSetCoordinateDim(dm,dimC)); 210 PetscCall(DMSetCoordinateSection(dm,dimC,newSec)); 211 PetscCall(DMSetCoordinatesLocal(dm,newVec)); 212 PetscCall(VecDestroy(&newVec)); 213 PetscCall(PetscSectionDestroy(&newSec)); 214 PetscCall(DMGetCoordinatesLocal(dm,&coords)); 215 PetscCall(DMGetCoordinateDM(dm,&coordDM)); 216 if (isFE) { 217 PetscFE fe; 218 219 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject) dm),dim,dimC,isSimplex ? PETSC_TRUE : PETSC_FALSE,isSimplex ? NULL : "tensor_",PETSC_DEFAULT,&fe)); 220 PetscCall(DMSetField(coordDM,0,NULL,(PetscObject)fe)); 221 PetscCall(PetscFEDestroy(&fe)); 222 PetscCall(DMCreateDS(coordDM)); 223 } 224 PetscCall(DMSetCoordinateDim(coordDM,dimC)); 225 PetscCall(VecGetLocalSize(coords,&n)); 226 } 227 PetscCall(VecGetArray(coords,&coordArray)); 228 for (i = 0; i < n; i++) { 229 PetscCall(PetscRandomGetValueReal(randCtx,&noise)); 230 coordArray[i] += noise * perturb; 231 } 232 PetscCall(VecRestoreArray(coords,&coordArray)); 233 PetscCall(DMSetCoordinatesLocal(dm, coords)); 234 PetscCall(testIdentity(dm, isSimplex ? PETSC_TRUE : PETSC_FALSE, 0, randCtx, numTests, tol)); 235 PetscCall(DMDestroy(&dm)); 236 } 237 } 238 } 239 } 240 PetscCall(PetscRandomDestroy(&randCtx)); 241 PetscCall(PetscFinalize()); 242 return 0; 243 } 244 245 /*TEST 246 247 test: 248 suffix: 0 249 args: -petscspace_degree 2 -tensor_petscspace_degree 2 250 251 TEST*/ 252