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 CHKERRQ(DMGetDimension(dm,&dimR)); 14 CHKERRQ(DMGetCoordinateDim(dm,&dimC)); 15 16 CHKERRQ(DMGetWorkArray(dm,dimR * numPoints,MPIU_REAL,&preimage)); 17 CHKERRQ(DMGetWorkArray(dm,dimC * numPoints,MPIU_REAL,&mapped)); 18 CHKERRQ(DMGetWorkArray(dm,dimR * numPoints,MPIU_REAL,&inverted)); 19 20 for (i = 0; i < dimR * numPoints; i++) { 21 CHKERRQ(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 CHKERRQ(DMPlexReferenceToCoordinates(dm,cell,numPoints,preimage,mapped)); 35 CHKERRQ(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 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Bad inversion for cell %D with error %g (tol %g): (",cell,(double)max,(double)tol)); 45 for (j = 0; j < dimR; j++) { 46 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"%+f",(double) preimage[i * dimR + j])); 47 if (j < dimR - 1) { 48 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,",")); 49 } 50 } 51 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,") --> (")); 52 for (j = 0; j < dimC; j++) { 53 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"%+f",(double) mapped[i * dimC + j])); 54 if (j < dimC - 1) { 55 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,",")); 56 } 57 } 58 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,") --> (")); 59 for (j = 0; j < dimR; j++) { 60 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"%+f",(double) inverted[i * dimR + j])); 61 if (j < dimR - 1) { 62 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,",")); 63 } 64 } 65 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,")\n")); 66 } else { 67 char strBuf[BUFSIZ] = {'\0'}; 68 size_t offset = 0, count; 69 70 CHKERRQ(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,"Good inversion for cell %D: (",&count,cell)); 71 offset += count; 72 for (j = 0; j < dimR; j++) { 73 CHKERRQ(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,"%+f",&count,(double) preimage[i * dimR + j])); 74 offset += count; 75 if (j < dimR - 1) { 76 CHKERRQ(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,",",&count)); 77 offset += count; 78 } 79 } 80 CHKERRQ(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,") --> (",&count)); 81 offset += count; 82 for (j = 0; j < dimC; j++) { 83 CHKERRQ(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,"%+f",&count,(double) mapped[i * dimC + j])); 84 offset += count; 85 if (j < dimC - 1) { 86 CHKERRQ(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,",",&count)); 87 offset += count; 88 } 89 } 90 CHKERRQ(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,") --> (",&count)); 91 offset += count; 92 for (j = 0; j < dimR; j++) { 93 CHKERRQ(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,"%+f",&count,(double) inverted[i * dimR + j])); 94 offset += count; 95 if (j < dimR - 1) { 96 CHKERRQ(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,",",&count)); 97 offset += count; 98 } 99 } 100 CHKERRQ(PetscSNPrintfCount(strBuf + offset,BUFSIZ-offset,")\n",&count)); 101 CHKERRQ(PetscInfo(dm,"%s",strBuf)); 102 } 103 } 104 105 CHKERRQ(DMRestoreWorkArray(dm,dimR * numPoints,MPIU_REAL,&inverted)); 106 CHKERRQ(DMRestoreWorkArray(dm,dimC * numPoints,MPIU_REAL,&mapped)); 107 CHKERRQ(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 PetscErrorCode ierr; 125 126 CHKERRQ(PetscInitialize(&argc,&argv,NULL,help)); 127 CHKERRQ(PetscRandomCreate(PETSC_COMM_SELF,&randCtx)); 128 CHKERRQ(PetscRandomSetInterval(randCtx,-1.,1.)); 129 CHKERRQ(PetscRandomSetFromOptions(randCtx)); 130 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex21",NULL);CHKERRQ(ierr); 131 CHKERRQ(PetscOptionsReal("-vertex_perturbation","scale of random vertex distortion",NULL,perturb,&perturb,NULL)); 132 CHKERRQ(PetscOptionsBoundedInt("-num_test_points","number of points to test",NULL,numTests,&numTests,NULL,0)); 133 ierr = PetscOptionsEnd();CHKERRQ(ierr); 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 CHKERRQ(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 CHKERRQ(PetscFECreateDefault(PetscObjectComm((PetscObject) dm),dim,dim,isSimplex ? PETSC_TRUE : PETSC_FALSE,isSimplex ? NULL : "tensor_" ,PETSC_DEFAULT,&fe)); 154 CHKERRQ(PetscFEGetBasisSpace(fe,&sp)); 155 CHKERRQ(PetscSpaceGetDegree(sp,&order,NULL)); 156 CHKERRQ(DMSetField(dm,0,NULL,(PetscObject)fe)); 157 CHKERRQ(DMCreateDS(dm)); 158 CHKERRQ(DMCreateLocalVector(dm,&localCoords)); 159 CHKERRQ(DMProjectFunctionLocal(dm,0,funcs,ctxs,INSERT_VALUES,localCoords)); 160 CHKERRQ(VecSetDM(localCoords,NULL)); /* This is necessary to prevent a reference loop */ 161 CHKERRQ(DMClone(dm,&dmCoord)); 162 CHKERRQ(DMSetField(dmCoord,0,NULL,(PetscObject)fe)); 163 CHKERRQ(PetscFEDestroy(&fe)); 164 CHKERRQ(DMCreateDS(dmCoord)); 165 CHKERRQ(DMSetCoordinateDM(dm,dmCoord)); 166 CHKERRQ(DMDestroy(&dmCoord)); 167 CHKERRQ(DMSetCoordinatesLocal(dm,localCoords)); 168 CHKERRQ(VecDestroy(&localCoords)); 169 } 170 CHKERRQ(PetscInfo(dm,"Testing %s%D %DD mesh embedded in %DD\n",isSimplex ? "P" : "Q" , order, dim, dimC)); 171 CHKERRQ(DMGetCoordinatesLocal(dm,&coords)); 172 CHKERRQ(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 CHKERRQ(DMGetCoordinateSection(dm,&sec)); 181 CHKERRQ(PetscSectionCreate(PetscObjectComm((PetscObject)sec),&newSec)); 182 CHKERRQ(PetscSectionSetNumFields(newSec,1)); 183 CHKERRQ(PetscSectionGetChart(sec,&pStart,&pEnd)); 184 CHKERRQ(PetscSectionSetChart(newSec,pStart,pEnd)); 185 for (p = pStart; p < pEnd; p++) { 186 PetscInt nDof; 187 188 CHKERRQ(PetscSectionGetDof(sec,p,&nDof)); 189 PetscCheckFalse(nDof % dim,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Coordinate section point %D has %D dofs != 0 mod %D",p,nDof,dim); 190 CHKERRQ(PetscSectionSetDof(newSec,p,(nDof/dim)*dimC)); 191 CHKERRQ(PetscSectionSetFieldDof(newSec,p,0,(nDof/dim)*dimC)); 192 } 193 CHKERRQ(PetscSectionSetUp(newSec)); 194 CHKERRQ(PetscSectionGetStorageSize(newSec,&newN)); 195 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,newN,&newVec)); 196 CHKERRQ(VecGetArray(newVec,&newCoordArray)); 197 CHKERRQ(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 CHKERRQ(VecRestoreArray(coords,&coordArray)); 209 CHKERRQ(VecRestoreArray(newVec,&newCoordArray)); 210 CHKERRQ(DMSetCoordinateDim(dm,dimC)); 211 CHKERRQ(DMSetCoordinateSection(dm,dimC,newSec)); 212 CHKERRQ(DMSetCoordinatesLocal(dm,newVec)); 213 CHKERRQ(VecDestroy(&newVec)); 214 CHKERRQ(PetscSectionDestroy(&newSec)); 215 CHKERRQ(DMGetCoordinatesLocal(dm,&coords)); 216 CHKERRQ(DMGetCoordinateDM(dm,&coordDM)); 217 if (isFE) { 218 PetscFE fe; 219 220 CHKERRQ(PetscFECreateDefault(PetscObjectComm((PetscObject) dm),dim,dimC,isSimplex ? PETSC_TRUE : PETSC_FALSE,isSimplex ? NULL : "tensor_",PETSC_DEFAULT,&fe)); 221 CHKERRQ(DMSetField(coordDM,0,NULL,(PetscObject)fe)); 222 CHKERRQ(PetscFEDestroy(&fe)); 223 CHKERRQ(DMCreateDS(coordDM)); 224 } 225 CHKERRQ(DMSetCoordinateDim(coordDM,dimC)); 226 CHKERRQ(VecGetLocalSize(coords,&n)); 227 } 228 CHKERRQ(VecGetArray(coords,&coordArray)); 229 for (i = 0; i < n; i++) { 230 CHKERRQ(PetscRandomGetValueReal(randCtx,&noise)); 231 coordArray[i] += noise * perturb; 232 } 233 CHKERRQ(VecRestoreArray(coords,&coordArray)); 234 CHKERRQ(DMSetCoordinatesLocal(dm, coords)); 235 CHKERRQ(testIdentity(dm, isSimplex ? PETSC_TRUE : PETSC_FALSE, 0, randCtx, numTests, tol)); 236 CHKERRQ(DMDestroy(&dm)); 237 } 238 } 239 } 240 } 241 CHKERRQ(PetscRandomDestroy(&randCtx)); 242 CHKERRQ(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