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