1 static char help[] = "Tests for cell geometry\n\n"; 2 3 #include <petscdmplex.h> 4 5 typedef enum {RUN_REFERENCE, RUN_HEX_CURVED, RUN_FILE, RUN_DISPLAY} RunType; 6 7 typedef struct { 8 DM dm; 9 RunType runType; /* Type of mesh to use */ 10 PetscBool transform; /* Use random coordinate transformations */ 11 /* Data for input meshes */ 12 PetscReal *v0, *J, *invJ, *detJ; /* FEM data */ 13 PetscReal *centroid, *normal, *vol; /* FVM data */ 14 } AppCtx; 15 16 static PetscErrorCode ReadMesh(MPI_Comm comm, AppCtx *user, DM *dm) 17 { 18 PetscFunctionBegin; 19 PetscCall(DMCreate(comm, dm)); 20 PetscCall(DMSetType(*dm, DMPLEX)); 21 PetscCall(DMSetFromOptions(*dm)); 22 PetscCall(DMSetApplicationContext(*dm, user)); 23 PetscCall(PetscObjectSetName((PetscObject) *dm, "Input Mesh")); 24 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 25 PetscFunctionReturn(0); 26 } 27 28 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 29 { 30 const char *runTypes[4] = {"reference", "hex_curved", "file", "display"}; 31 PetscInt run; 32 PetscErrorCode ierr; 33 34 PetscFunctionBeginUser; 35 options->runType = RUN_REFERENCE; 36 options->transform = PETSC_FALSE; 37 38 ierr = PetscOptionsBegin(comm, "", "Geometry Test Options", "DMPLEX");PetscCall(ierr); 39 run = options->runType; 40 PetscCall(PetscOptionsEList("-run_type", "The run type", "ex8.c", runTypes, 3, runTypes[options->runType], &run, NULL)); 41 options->runType = (RunType) run; 42 PetscCall(PetscOptionsBool("-transform", "Use random transforms", "ex8.c", options->transform, &options->transform, NULL)); 43 44 if (options->runType == RUN_FILE) { 45 PetscInt dim, cStart, cEnd, numCells, n; 46 PetscBool flg, feFlg; 47 48 PetscCall(ReadMesh(PETSC_COMM_WORLD, options, &options->dm)); 49 PetscCall(DMGetDimension(options->dm, &dim)); 50 PetscCall(DMPlexGetHeightStratum(options->dm, 0, &cStart, &cEnd)); 51 numCells = cEnd-cStart; 52 PetscCall(PetscMalloc4(numCells*dim, &options->v0, numCells*dim*dim, &options->J, numCells*dim*dim, &options->invJ, numCells, &options->detJ)); 53 PetscCall(PetscMalloc1(numCells*dim, &options->centroid)); 54 PetscCall(PetscMalloc1(numCells*dim, &options->normal)); 55 PetscCall(PetscMalloc1(numCells, &options->vol)); 56 n = numCells*dim; 57 PetscCall(PetscOptionsRealArray("-v0", "Input v0 for each cell", "ex8.c", options->v0, &n, &feFlg)); 58 PetscCheckFalse(feFlg && n != numCells*dim,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid size of v0 %D should be %D", n, numCells*dim); 59 n = numCells*dim*dim; 60 PetscCall(PetscOptionsRealArray("-J", "Input Jacobian for each cell", "ex8.c", options->J, &n, &flg)); 61 PetscCheckFalse(flg && n != numCells*dim*dim,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid size of J %D should be %D", n, numCells*dim*dim); 62 n = numCells*dim*dim; 63 PetscCall(PetscOptionsRealArray("-invJ", "Input inverse Jacobian for each cell", "ex8.c", options->invJ, &n, &flg)); 64 PetscCheckFalse(flg && n != numCells*dim*dim,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid size of invJ %D should be %D", n, numCells*dim*dim); 65 n = numCells; 66 PetscCall(PetscOptionsRealArray("-detJ", "Input Jacobian determinant for each cell", "ex8.c", options->detJ, &n, &flg)); 67 PetscCheckFalse(flg && n != numCells,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid size of detJ %D should be %D", n, numCells); 68 n = numCells*dim; 69 if (!feFlg) { 70 PetscCall(PetscFree4(options->v0, options->J, options->invJ, options->detJ)); 71 options->v0 = options->J = options->invJ = options->detJ = NULL; 72 } 73 PetscCall(PetscOptionsRealArray("-centroid", "Input centroid for each cell", "ex8.c", options->centroid, &n, &flg)); 74 PetscCheckFalse(flg && n != numCells*dim,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid size of centroid %D should be %D", n, numCells*dim); 75 if (!flg) { 76 PetscCall(PetscFree(options->centroid)); 77 options->centroid = NULL; 78 } 79 n = numCells*dim; 80 PetscCall(PetscOptionsRealArray("-normal", "Input normal for each cell", "ex8.c", options->normal, &n, &flg)); 81 PetscCheckFalse(flg && n != numCells*dim,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid size of normal %D should be %D", n, numCells*dim); 82 if (!flg) { 83 PetscCall(PetscFree(options->normal)); 84 options->normal = NULL; 85 } 86 n = numCells; 87 PetscCall(PetscOptionsRealArray("-vol", "Input volume for each cell", "ex8.c", options->vol, &n, &flg)); 88 PetscCheckFalse(flg && n != numCells,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid size of vol %D should be %D", n, numCells); 89 if (!flg) { 90 PetscCall(PetscFree(options->vol)); 91 options->vol = NULL; 92 } 93 } else if (options->runType == RUN_DISPLAY) { 94 PetscCall(ReadMesh(PETSC_COMM_WORLD, options, &options->dm)); 95 } 96 ierr = PetscOptionsEnd();PetscCall(ierr); 97 98 if (options->transform) PetscCall(PetscPrintf(comm, "Using random transforms\n")); 99 PetscFunctionReturn(0); 100 } 101 102 static PetscErrorCode ChangeCoordinates(DM dm, PetscInt spaceDim, PetscScalar vertexCoords[]) 103 { 104 PetscSection coordSection; 105 Vec coordinates; 106 PetscScalar *coords; 107 PetscInt vStart, vEnd, v, d, coordSize; 108 109 PetscFunctionBegin; 110 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 111 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 112 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 113 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, spaceDim)); 114 PetscCall(PetscSectionSetChart(coordSection, vStart, vEnd)); 115 for (v = vStart; v < vEnd; ++v) { 116 PetscCall(PetscSectionSetDof(coordSection, v, spaceDim)); 117 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, spaceDim)); 118 } 119 PetscCall(PetscSectionSetUp(coordSection)); 120 PetscCall(PetscSectionGetStorageSize(coordSection, &coordSize)); 121 PetscCall(VecCreate(PETSC_COMM_SELF, &coordinates)); 122 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 123 PetscCall(VecSetSizes(coordinates, coordSize, PETSC_DETERMINE)); 124 PetscCall(VecSetFromOptions(coordinates)); 125 PetscCall(VecGetArray(coordinates, &coords)); 126 for (v = vStart; v < vEnd; ++v) { 127 PetscInt off; 128 129 PetscCall(PetscSectionGetOffset(coordSection, v, &off)); 130 for (d = 0; d < spaceDim; ++d) { 131 coords[off+d] = vertexCoords[(v-vStart)*spaceDim+d]; 132 } 133 } 134 PetscCall(VecRestoreArray(coordinates, &coords)); 135 PetscCall(DMSetCoordinateDim(dm, spaceDim)); 136 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 137 PetscCall(VecDestroy(&coordinates)); 138 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 139 PetscFunctionReturn(0); 140 } 141 142 #define RelativeError(a,b) PetscAbs(a-b)/(1.0+PetscMax(PetscAbs(a),PetscAbs(b))) 143 144 static PetscErrorCode CheckFEMGeometry(DM dm, PetscInt cell, PetscInt spaceDim, PetscReal v0Ex[], PetscReal JEx[], PetscReal invJEx[], PetscReal detJEx) 145 { 146 PetscReal v0[3], J[9], invJ[9], detJ; 147 PetscInt d, i, j; 148 149 PetscFunctionBegin; 150 PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, invJ, &detJ)); 151 for (d = 0; d < spaceDim; ++d) { 152 if (v0[d] != v0Ex[d]) { 153 switch (spaceDim) { 154 case 2: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid v0 (%g, %g) != (%g, %g)", (double)v0[0], (double)v0[1], (double)v0Ex[0], (double)v0Ex[1]); 155 case 3: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid v0 (%g, %g, %g) != (%g, %g, %g)", (double)v0[0], (double)v0[1], (double)v0[2], (double)v0Ex[0], (double)v0Ex[1], (double)v0Ex[2]); 156 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid space dimension %D", spaceDim); 157 } 158 } 159 } 160 for (i = 0; i < spaceDim; ++i) { 161 for (j = 0; j < spaceDim; ++j) { 162 PetscCheckFalse(RelativeError(J[i*spaceDim+j],JEx[i*spaceDim+j]) > 10*PETSC_SMALL,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid J[%D,%D]: %g != %g", i, j, (double)J[i*spaceDim+j], (double)JEx[i*spaceDim+j]); 163 PetscCheckFalse(RelativeError(invJ[i*spaceDim+j],invJEx[i*spaceDim+j]) > 10*PETSC_SMALL,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid invJ[%D,%D]: %g != %g", i, j, (double)invJ[i*spaceDim+j], (double)invJEx[i*spaceDim+j]); 164 } 165 } 166 PetscCheckFalse(RelativeError(detJ,detJEx) > 10*PETSC_SMALL,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid |J| = %g != %g diff %g", (double)detJ, (double)detJEx,(double)(detJ - detJEx)); 167 PetscFunctionReturn(0); 168 } 169 170 static PetscErrorCode CheckFVMGeometry(DM dm, PetscInt cell, PetscInt spaceDim, PetscReal centroidEx[], PetscReal normalEx[], PetscReal volEx) 171 { 172 PetscReal tol = PetscMax(10*PETSC_SMALL, 1e-10); 173 PetscReal centroid[3], normal[3], vol; 174 PetscInt d; 175 176 PetscFunctionBegin; 177 PetscCall(DMPlexComputeCellGeometryFVM(dm, cell, volEx? &vol : NULL, centroidEx? centroid : NULL, normalEx? normal : NULL)); 178 for (d = 0; d < spaceDim; ++d) { 179 if (centroidEx) 180 PetscCheckFalse(RelativeError(centroid[d],centroidEx[d]) > tol,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %D, Invalid centroid[%D]: %g != %g diff %g", cell, d, (double)centroid[d], (double)centroidEx[d],(double)(centroid[d]-centroidEx[d])); 181 if (normalEx) 182 PetscCheckFalse(RelativeError(normal[d],normalEx[d]) > tol,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %D, Invalid normal[%D]: %g != %g", cell, d, (double)normal[d], (double)normalEx[d]); 183 } 184 if (volEx) 185 PetscCheckFalse(RelativeError(volEx,vol) > tol,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %D, Invalid volume = %g != %g diff %g", cell, (double)vol, (double)volEx,(double)(vol - volEx)); 186 PetscFunctionReturn(0); 187 } 188 189 static PetscErrorCode CheckGaussLaw(DM dm, PetscInt cell) 190 { 191 DMPolytopeType ct; 192 PetscReal tol = PetscMax(10*PETSC_SMALL, 1e-10); 193 PetscReal normal[3], integral[3] = {0., 0., 0.}, area; 194 const PetscInt *cone, *ornt; 195 PetscInt coneSize, f, dim, cdim, d; 196 197 PetscFunctionBegin; 198 PetscCall(DMGetDimension(dm, &dim)); 199 PetscCall(DMGetCoordinateDim(dm, &cdim)); 200 if (dim != cdim) PetscFunctionReturn(0); 201 PetscCall(DMPlexGetCellType(dm, cell, &ct)); 202 if (ct == DM_POLYTOPE_TRI_PRISM_TENSOR) PetscFunctionReturn(0); 203 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 204 PetscCall(DMPlexGetCone(dm, cell, &cone)); 205 PetscCall(DMPlexGetConeOrientation(dm, cell, &ornt)); 206 for (f = 0; f < coneSize; ++f) { 207 const PetscInt sgn = dim == 1? (f == 0 ? -1 : 1) : (ornt[f] < 0 ? -1 : 1); 208 209 PetscCall(DMPlexComputeCellGeometryFVM(dm, cone[f], &area, NULL, normal)); 210 for (d = 0; d < cdim; ++d) integral[d] += sgn*area*normal[d]; 211 } 212 for (d = 0; d < cdim; ++d) PetscCheckFalse(PetscAbsReal(integral[d]) > tol,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %D Surface integral for component %D: %g != 0. as it should be for a constant field", cell, d, (double) integral[d]); 213 PetscFunctionReturn(0); 214 } 215 216 static PetscErrorCode CheckCell(DM dm, PetscInt cell, PetscBool transform, PetscReal v0Ex[], PetscReal JEx[], PetscReal invJEx[], PetscReal detJEx, PetscReal centroidEx[], PetscReal normalEx[], PetscReal volEx, PetscReal faceCentroidEx[], PetscReal faceNormalEx[], PetscReal faceVolEx[]) 217 { 218 const PetscInt *cone; 219 PetscInt coneSize, c; 220 PetscInt dim, depth, cdim; 221 222 PetscFunctionBegin; 223 PetscCall(DMPlexGetDepth(dm, &depth)); 224 PetscCall(DMGetDimension(dm, &dim)); 225 PetscCall(DMGetCoordinateDim(dm, &cdim)); 226 if (v0Ex) { 227 PetscCall(CheckFEMGeometry(dm, cell, cdim, v0Ex, JEx, invJEx, detJEx)); 228 } 229 if (dim == depth && centroidEx) { 230 PetscCall(CheckFVMGeometry(dm, cell, cdim, centroidEx, normalEx, volEx)); 231 PetscCall(CheckGaussLaw(dm, cell)); 232 if (faceCentroidEx) { 233 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 234 PetscCall(DMPlexGetCone(dm, cell, &cone)); 235 for (c = 0; c < coneSize; ++c) { 236 PetscCall(CheckFVMGeometry(dm, cone[c], dim, &faceCentroidEx[c*dim], &faceNormalEx[c*dim], faceVolEx[c])); 237 } 238 } 239 } 240 if (transform) { 241 Vec coordinates; 242 PetscSection coordSection; 243 PetscScalar *coords = NULL, *origCoords, *newCoords; 244 PetscReal *v0ExT, *JExT, *invJExT, detJExT=0, *centroidExT, *normalExT, volExT=0; 245 PetscReal *faceCentroidExT, *faceNormalExT, faceVolExT; 246 PetscRandom r, ang, ang2; 247 PetscInt coordSize, numCorners, t; 248 249 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 250 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 251 PetscCall(DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords)); 252 PetscCall(PetscMalloc2(coordSize, &origCoords, coordSize, &newCoords)); 253 PetscCall(PetscMalloc5(cdim, &v0ExT, cdim*cdim, &JExT, cdim*cdim, &invJExT, cdim, ¢roidExT, cdim, &normalExT)); 254 PetscCall(PetscMalloc2(cdim, &faceCentroidExT, cdim, &faceNormalExT)); 255 for (c = 0; c < coordSize; ++c) origCoords[c] = coords[c]; 256 PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords)); 257 numCorners = coordSize/cdim; 258 259 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r)); 260 PetscCall(PetscRandomSetFromOptions(r)); 261 PetscCall(PetscRandomSetInterval(r, 0.0, 10.0)); 262 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &ang)); 263 PetscCall(PetscRandomSetFromOptions(ang)); 264 PetscCall(PetscRandomSetInterval(ang, 0.0, 2*PETSC_PI)); 265 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &ang2)); 266 PetscCall(PetscRandomSetFromOptions(ang2)); 267 PetscCall(PetscRandomSetInterval(ang2, 0.0, PETSC_PI)); 268 for (t = 0; t < 1; ++t) { 269 PetscScalar trans[3]; 270 PetscReal R[9], rot[3], rotM[9]; 271 PetscReal scale, phi, theta, psi = 0.0, norm; 272 PetscInt d, e, f, p; 273 274 for (c = 0; c < coordSize; ++c) newCoords[c] = origCoords[c]; 275 PetscCall(PetscRandomGetValueReal(r, &scale)); 276 PetscCall(PetscRandomGetValueReal(ang, &phi)); 277 PetscCall(PetscRandomGetValueReal(ang2, &theta)); 278 for (d = 0; d < cdim; ++d) { 279 PetscCall(PetscRandomGetValue(r, &trans[d])); 280 } 281 switch (cdim) { 282 case 2: 283 R[0] = PetscCosReal(phi); R[1] = -PetscSinReal(phi); 284 R[2] = PetscSinReal(phi); R[3] = PetscCosReal(phi); 285 break; 286 case 3: 287 { 288 const PetscReal ct = PetscCosReal(theta), st = PetscSinReal(theta); 289 const PetscReal cp = PetscCosReal(phi), sp = PetscSinReal(phi); 290 const PetscReal cs = PetscCosReal(psi), ss = PetscSinReal(psi); 291 R[0] = ct*cs; R[1] = sp*st*cs - cp*ss; R[2] = sp*ss + cp*st*cs; 292 R[3] = ct*ss; R[4] = cp*cs + sp*st*ss; R[5] = cp*st*ss - sp*cs; 293 R[6] = -st; R[7] = sp*ct; R[8] = cp*ct; 294 break; 295 } 296 default: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Invalid coordinate dimension %D", cdim); 297 } 298 if (v0Ex) { 299 detJExT = detJEx; 300 for (d = 0; d < cdim; ++d) { 301 v0ExT[d] = v0Ex[d]; 302 for (e = 0; e < cdim; ++e) { 303 JExT[d*cdim+e] = JEx[d*cdim+e]; 304 invJExT[d*cdim+e] = invJEx[d*cdim+e]; 305 } 306 } 307 for (d = 0; d < cdim; ++d) { 308 v0ExT[d] *= scale; 309 v0ExT[d] += PetscRealPart(trans[d]); 310 /* Only scale dimensions in the manifold */ 311 for (e = 0; e < dim; ++e) { 312 JExT[d*cdim+e] *= scale; 313 invJExT[d*cdim+e] /= scale; 314 } 315 if (d < dim) detJExT *= scale; 316 } 317 /* Do scaling and translation before rotation, so that we can leave out the normal dimension for lower dimensional manifolds */ 318 for (d = 0; d < cdim; ++d) { 319 for (e = 0, rot[d] = 0.0; e < cdim; ++e) { 320 rot[d] += R[d*cdim+e] * v0ExT[e]; 321 } 322 } 323 for (d = 0; d < cdim; ++d) v0ExT[d] = rot[d]; 324 for (d = 0; d < cdim; ++d) { 325 for (e = 0; e < cdim; ++e) { 326 for (f = 0, rotM[d*cdim+e] = 0.0; f < cdim; ++f) { 327 rotM[d*cdim+e] += R[d*cdim+f] * JExT[f*cdim+e]; 328 } 329 } 330 } 331 for (d = 0; d < cdim; ++d) { 332 for (e = 0; e < cdim; ++e) { 333 JExT[d*cdim+e] = rotM[d*cdim+e]; 334 } 335 } 336 for (d = 0; d < cdim; ++d) { 337 for (e = 0; e < cdim; ++e) { 338 for (f = 0, rotM[d*cdim+e] = 0.0; f < cdim; ++f) { 339 rotM[d*cdim+e] += invJExT[d*cdim+f] * R[e*cdim+f]; 340 } 341 } 342 } 343 for (d = 0; d < cdim; ++d) { 344 for (e = 0; e < cdim; ++e) { 345 invJExT[d*cdim+e] = rotM[d*cdim+e]; 346 } 347 } 348 } 349 if (centroidEx) { 350 volExT = volEx; 351 for (d = 0; d < cdim; ++d) { 352 centroidExT[d] = centroidEx[d]; 353 normalExT[d] = normalEx[d]; 354 } 355 for (d = 0; d < cdim; ++d) { 356 centroidExT[d] *= scale; 357 centroidExT[d] += PetscRealPart(trans[d]); 358 normalExT[d] /= scale; 359 /* Only scale dimensions in the manifold */ 360 if (d < dim) volExT *= scale; 361 } 362 /* Do scaling and translation before rotation, so that we can leave out the normal dimension for lower dimensional manifolds */ 363 for (d = 0; d < cdim; ++d) { 364 for (e = 0, rot[d] = 0.0; e < cdim; ++e) { 365 rot[d] += R[d*cdim+e] * centroidExT[e]; 366 } 367 } 368 for (d = 0; d < cdim; ++d) centroidExT[d] = rot[d]; 369 for (d = 0; d < cdim; ++d) { 370 for (e = 0, rot[d] = 0.0; e < cdim; ++e) { 371 rot[d] += R[d*cdim+e] * normalExT[e]; 372 } 373 } 374 for (d = 0; d < cdim; ++d) normalExT[d] = rot[d]; 375 for (d = 0, norm = 0.0; d < cdim; ++d) norm += PetscSqr(normalExT[d]); 376 norm = PetscSqrtReal(norm); 377 for (d = 0; d < cdim; ++d) normalExT[d] /= norm; 378 } 379 for (d = 0; d < cdim; ++d) { 380 for (p = 0; p < numCorners; ++p) { 381 newCoords[p*cdim+d] *= scale; 382 newCoords[p*cdim+d] += trans[d]; 383 } 384 } 385 for (p = 0; p < numCorners; ++p) { 386 for (d = 0; d < cdim; ++d) { 387 for (e = 0, rot[d] = 0.0; e < cdim; ++e) { 388 rot[d] += R[d*cdim+e] * PetscRealPart(newCoords[p*cdim+e]); 389 } 390 } 391 for (d = 0; d < cdim; ++d) newCoords[p*cdim+d] = rot[d]; 392 } 393 394 PetscCall(ChangeCoordinates(dm, cdim, newCoords)); 395 if (v0Ex) { 396 PetscCall(CheckFEMGeometry(dm, 0, cdim, v0ExT, JExT, invJExT, detJExT)); 397 } 398 if (dim == depth && centroidEx) { 399 PetscCall(CheckFVMGeometry(dm, cell, cdim, centroidExT, normalExT, volExT)); 400 PetscCall(CheckGaussLaw(dm, cell)); 401 if (faceCentroidEx) { 402 PetscCall(DMPlexGetConeSize(dm, cell, &coneSize)); 403 PetscCall(DMPlexGetCone(dm, cell, &cone)); 404 for (c = 0; c < coneSize; ++c) { 405 PetscInt off = c*cdim; 406 407 faceVolExT = faceVolEx[c]; 408 for (d = 0; d < cdim; ++d) { 409 faceCentroidExT[d] = faceCentroidEx[off+d]; 410 faceNormalExT[d] = faceNormalEx[off+d]; 411 } 412 for (d = 0; d < cdim; ++d) { 413 faceCentroidExT[d] *= scale; 414 faceCentroidExT[d] += PetscRealPart(trans[d]); 415 faceNormalExT[d] /= scale; 416 /* Only scale dimensions in the manifold */ 417 if (d < dim-1) { 418 faceVolExT *= scale; 419 } 420 } 421 for (d = 0; d < cdim; ++d) { 422 for (e = 0, rot[d] = 0.0; e < cdim; ++e) { 423 rot[d] += R[d*cdim+e] * faceCentroidExT[e]; 424 } 425 } 426 for (d = 0; d < cdim; ++d) faceCentroidExT[d] = rot[d]; 427 for (d = 0; d < cdim; ++d) { 428 for (e = 0, rot[d] = 0.0; e < cdim; ++e) { 429 rot[d] += R[d*cdim+e] * faceNormalExT[e]; 430 } 431 } 432 for (d = 0; d < cdim; ++d) faceNormalExT[d] = rot[d]; 433 for (d = 0, norm = 0.0; d < cdim; ++d) norm += PetscSqr(faceNormalExT[d]); 434 norm = PetscSqrtReal(norm); 435 for (d = 0; d < cdim; ++d) faceNormalExT[d] /= norm; 436 437 PetscCall(CheckFVMGeometry(dm, cone[c], cdim, faceCentroidExT, faceNormalExT, faceVolExT)); 438 } 439 } 440 } 441 } 442 PetscCall(PetscRandomDestroy(&r)); 443 PetscCall(PetscRandomDestroy(&ang)); 444 PetscCall(PetscRandomDestroy(&ang2)); 445 PetscCall(PetscFree2(origCoords, newCoords)); 446 PetscCall(PetscFree5(v0ExT, JExT, invJExT, centroidExT, normalExT)); 447 PetscCall(PetscFree2(faceCentroidExT, faceNormalExT)); 448 } 449 PetscFunctionReturn(0); 450 } 451 452 static PetscErrorCode TestTriangle(MPI_Comm comm, PetscBool transform) 453 { 454 DM dm; 455 456 PetscFunctionBegin; 457 PetscCall(DMPlexCreateReferenceCell(comm, DM_POLYTOPE_TRIANGLE, &dm)); 458 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 459 /* Check reference geometry: determinant is scaled by reference volume (2.0) */ 460 { 461 PetscReal v0Ex[2] = {-1.0, -1.0}; 462 PetscReal JEx[4] = {1.0, 0.0, 0.0, 1.0}; 463 PetscReal invJEx[4] = {1.0, 0.0, 0.0, 1.0}; 464 PetscReal detJEx = 1.0; 465 PetscReal centroidEx[2] = {-((PetscReal)1.)/((PetscReal)3.), -((PetscReal)1.)/((PetscReal)3.)}; 466 PetscReal normalEx[2] = {0.0, 0.0}; 467 PetscReal volEx = 2.0; 468 469 PetscCall(CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL)); 470 } 471 /* Move to 3D: Check reference geometry: determinant is scaled by reference volume (2.0) */ 472 { 473 PetscScalar vertexCoords[9] = {-1.0, -1.0, 0.0, 1.0, -1.0, 0.0, -1.0, 1.0, 0.0}; 474 PetscReal v0Ex[3] = {-1.0, -1.0, 0.0}; 475 PetscReal JEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; 476 PetscReal invJEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; 477 PetscReal detJEx = 1.0; 478 PetscReal centroidEx[3] = {-((PetscReal)1.)/((PetscReal)3.), -((PetscReal)1.)/((PetscReal)3.), 0.0}; 479 PetscReal normalEx[3] = {0.0, 0.0, 1.0}; 480 PetscReal volEx = 2.0; 481 482 PetscCall(ChangeCoordinates(dm, 3, vertexCoords)); 483 PetscCall(CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL)); 484 } 485 /* Cleanup */ 486 PetscCall(DMDestroy(&dm)); 487 PetscFunctionReturn(0); 488 } 489 490 static PetscErrorCode TestQuadrilateral(MPI_Comm comm, PetscBool transform) 491 { 492 DM dm; 493 494 PetscFunctionBegin; 495 PetscCall(DMPlexCreateReferenceCell(comm, DM_POLYTOPE_QUADRILATERAL, &dm)); 496 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 497 /* Check reference geometry: determinant is scaled by reference volume (2.0) */ 498 { 499 PetscReal v0Ex[2] = {-1.0, -1.0}; 500 PetscReal JEx[4] = {1.0, 0.0, 0.0, 1.0}; 501 PetscReal invJEx[4] = {1.0, 0.0, 0.0, 1.0}; 502 PetscReal detJEx = 1.0; 503 PetscReal centroidEx[2] = {0.0, 0.0}; 504 PetscReal normalEx[2] = {0.0, 0.0}; 505 PetscReal volEx = 4.0; 506 507 PetscCall(CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL)); 508 } 509 /* Move to 3D: Check reference geometry: determinant is scaled by reference volume (4.0) */ 510 { 511 PetscScalar vertexCoords[12] = {-1.0, -1.0, 0.0, 1.0, -1.0, 0.0, 1.0, 1.0, 0.0, -1.0, 1.0, 0.0}; 512 PetscReal v0Ex[3] = {-1.0, -1.0, 0.0}; 513 PetscReal JEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; 514 PetscReal invJEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; 515 PetscReal detJEx = 1.0; 516 PetscReal centroidEx[3] = {0.0, 0.0, 0.0}; 517 PetscReal normalEx[3] = {0.0, 0.0, 1.0}; 518 PetscReal volEx = 4.0; 519 520 PetscCall(ChangeCoordinates(dm, 3, vertexCoords)); 521 PetscCall(CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL)); 522 } 523 /* Cleanup */ 524 PetscCall(DMDestroy(&dm)); 525 PetscFunctionReturn(0); 526 } 527 528 static PetscErrorCode TestTetrahedron(MPI_Comm comm, PetscBool transform) 529 { 530 DM dm; 531 532 PetscFunctionBegin; 533 PetscCall(DMPlexCreateReferenceCell(comm, DM_POLYTOPE_TETRAHEDRON, &dm)); 534 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 535 /* Check reference geometry: determinant is scaled by reference volume (4/3) */ 536 { 537 PetscReal v0Ex[3] = {-1.0, -1.0, -1.0}; 538 PetscReal JEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; 539 PetscReal invJEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; 540 PetscReal detJEx = 1.0; 541 PetscReal centroidEx[3] = {-0.5, -0.5, -0.5}; 542 PetscReal normalEx[3] = {0.0, 0.0, 0.0}; 543 PetscReal volEx = (PetscReal)4.0/(PetscReal)3.0; 544 545 PetscCall(CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL)); 546 } 547 /* Cleanup */ 548 PetscCall(DMDestroy(&dm)); 549 PetscFunctionReturn(0); 550 } 551 552 static PetscErrorCode TestHexahedron(MPI_Comm comm, PetscBool transform) 553 { 554 DM dm; 555 556 PetscFunctionBegin; 557 PetscCall(DMPlexCreateReferenceCell(comm, DM_POLYTOPE_HEXAHEDRON, &dm)); 558 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 559 /* Check reference geometry: determinant is scaled by reference volume 8.0 */ 560 { 561 PetscReal v0Ex[3] = {-1.0, -1.0, -1.0}; 562 PetscReal JEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; 563 PetscReal invJEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; 564 PetscReal detJEx = 1.0; 565 PetscReal centroidEx[3] = {0.0, 0.0, 0.0}; 566 PetscReal normalEx[3] = {0.0, 0.0, 0.0}; 567 PetscReal volEx = 8.0; 568 569 PetscCall(CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL)); 570 } 571 /* Cleanup */ 572 PetscCall(DMDestroy(&dm)); 573 PetscFunctionReturn(0); 574 } 575 576 static PetscErrorCode TestHexahedronCurved(MPI_Comm comm) 577 { 578 DM dm; 579 PetscScalar coords[24] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 580 -1.0, -1.0, 1.1, 1.0, -1.0, 1.0, 1.0, 1.0, 1.1, -1.0, 1.0, 1.0}; 581 582 PetscFunctionBegin; 583 PetscCall(DMPlexCreateReferenceCell(comm, DM_POLYTOPE_HEXAHEDRON, &dm)); 584 PetscCall(ChangeCoordinates(dm, 3, coords)); 585 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 586 { 587 PetscReal centroidEx[3] = {0.0, 0.0, 0.016803278688524603}; 588 PetscReal normalEx[3] = {0.0, 0.0, 0.0}; 589 PetscReal volEx = 8.1333333333333346; 590 591 PetscCall(CheckCell(dm, 0, PETSC_FALSE, NULL, NULL, NULL, 0.0, centroidEx, normalEx, volEx, NULL, NULL, NULL)); 592 } 593 PetscCall(DMDestroy(&dm)); 594 PetscFunctionReturn(0); 595 } 596 597 /* This wedge is a tensor product cell, rather than a normal wedge */ 598 static PetscErrorCode TestWedge(MPI_Comm comm, PetscBool transform) 599 { 600 DM dm; 601 602 PetscFunctionBegin; 603 PetscCall(DMPlexCreateReferenceCell(comm, DM_POLYTOPE_TRI_PRISM_TENSOR, &dm)); 604 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 605 /* Check reference geometry: determinant is scaled by reference volume 4.0 */ 606 { 607 #if 0 608 /* FEM geometry is not functional for wedges */ 609 PetscReal v0Ex[3] = {-1.0, -1.0, -1.0}; 610 PetscReal JEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; 611 PetscReal invJEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; 612 PetscReal detJEx = 1.0; 613 #endif 614 615 { 616 PetscReal centroidEx[3] = {-((PetscReal)1.)/((PetscReal)3.), -((PetscReal)1.)/((PetscReal)3.), 0.0}; 617 PetscReal normalEx[3] = {0.0, 0.0, 0.0}; 618 PetscReal volEx = 4.0; 619 PetscReal faceVolEx[5] = {2.0, 2.0, 4.0, PETSC_SQRT2*4.0, 4.0}; 620 PetscReal faceNormalEx[15] = {0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, PETSC_SQRT2/2.0, PETSC_SQRT2/2.0, 0.0, -1.0, 0.0, 0.0}; 621 PetscReal faceCentroidEx[15] = {-((PetscReal)1.)/((PetscReal)3.), -((PetscReal)1.)/((PetscReal)3.), -1.0, 622 -((PetscReal)1.)/((PetscReal)3.), -((PetscReal)1.)/((PetscReal)3.), 1.0, 623 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0}; 624 625 PetscCall(CheckCell(dm, 0, transform, NULL, NULL, NULL, 0.0, centroidEx, normalEx, volEx, faceCentroidEx, faceNormalEx, faceVolEx)); 626 } 627 } 628 /* Cleanup */ 629 PetscCall(DMDestroy(&dm)); 630 PetscFunctionReturn(0); 631 } 632 633 int main(int argc, char **argv) 634 { 635 AppCtx user; 636 637 PetscCall(PetscInitialize(&argc, &argv, NULL,help)); 638 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 639 if (user.runType == RUN_REFERENCE) { 640 PetscCall(TestTriangle(PETSC_COMM_SELF, user.transform)); 641 PetscCall(TestQuadrilateral(PETSC_COMM_SELF, user.transform)); 642 PetscCall(TestTetrahedron(PETSC_COMM_SELF, user.transform)); 643 PetscCall(TestHexahedron(PETSC_COMM_SELF, user.transform)); 644 PetscCall(TestWedge(PETSC_COMM_SELF, user.transform)); 645 } else if (user.runType == RUN_HEX_CURVED) { 646 PetscCall(TestHexahedronCurved(PETSC_COMM_SELF)); 647 } else if (user.runType == RUN_FILE) { 648 PetscInt dim, cStart, cEnd, c; 649 650 PetscCall(DMGetDimension(user.dm, &dim)); 651 PetscCall(DMPlexGetHeightStratum(user.dm, 0, &cStart, &cEnd)); 652 for (c = 0; c < cEnd-cStart; ++c) { 653 PetscReal *v0 = user.v0 ? &user.v0[c*dim] : NULL; 654 PetscReal *J = user.J ? &user.J[c*dim*dim] : NULL; 655 PetscReal *invJ = user.invJ ? &user.invJ[c*dim*dim] : NULL; 656 PetscReal detJ = user.detJ ? user.detJ[c] : 0.0; 657 PetscReal *centroid = user.centroid ? &user.centroid[c*dim] : NULL; 658 PetscReal *normal = user.normal ? &user.normal[c*dim] : NULL; 659 PetscReal vol = user.vol ? user.vol[c] : 0.0; 660 661 PetscCall(CheckCell(user.dm, c+cStart, PETSC_FALSE, v0, J, invJ, detJ, centroid, normal, vol, NULL, NULL, NULL)); 662 } 663 PetscCall(PetscFree4(user.v0,user.J,user.invJ,user.detJ)); 664 PetscCall(PetscFree(user.centroid)); 665 PetscCall(PetscFree(user.normal)); 666 PetscCall(PetscFree(user.vol)); 667 PetscCall(DMDestroy(&user.dm)); 668 } else if (user.runType == RUN_DISPLAY) { 669 DM gdm, dmCell; 670 Vec cellgeom, facegeom; 671 const PetscScalar *cgeom; 672 PetscInt dim, d, cStart, cEnd, cEndInterior, c; 673 674 PetscCall(DMGetCoordinateDim(user.dm, &dim)); 675 PetscCall(DMPlexConstructGhostCells(user.dm, NULL, NULL, &gdm)); 676 if (gdm) { 677 PetscCall(DMDestroy(&user.dm)); 678 user.dm = gdm; 679 } 680 PetscCall(DMPlexComputeGeometryFVM(user.dm, &cellgeom, &facegeom)); 681 PetscCall(DMPlexGetHeightStratum(user.dm, 0, &cStart, &cEnd)); 682 PetscCall(DMPlexGetGhostCellStratum(user.dm, &cEndInterior, NULL)); 683 if (cEndInterior >= 0) cEnd = cEndInterior; 684 PetscCall(VecGetDM(cellgeom, &dmCell)); 685 PetscCall(VecGetArrayRead(cellgeom, &cgeom)); 686 for (c = 0; c < cEnd-cStart; ++c) { 687 PetscFVCellGeom *cg; 688 689 PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg)); 690 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cell %4D: Centroid (", c)); 691 for (d = 0; d < dim; ++d) { 692 if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", ")); 693 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%12.2g", cg->centroid[d])); 694 } 695 PetscCall(PetscPrintf(PETSC_COMM_SELF, ") Vol %12.2g\n", cg->volume)); 696 } 697 PetscCall(VecRestoreArrayRead(cellgeom, &cgeom)); 698 PetscCall(VecDestroy(&cellgeom)); 699 PetscCall(VecDestroy(&facegeom)); 700 PetscCall(DMDestroy(&user.dm)); 701 } 702 PetscCall(PetscFinalize()); 703 return 0; 704 } 705 706 /*TEST 707 708 test: 709 suffix: 1 710 args: -dm_view ascii::ascii_info_detail 711 test: 712 suffix: 2 713 args: -run_type hex_curved 714 test: 715 suffix: 3 716 args: -transform 717 test: 718 suffix: 4 719 requires: exodusii 720 args: -run_type file -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/simpleblock-100.exo -dm_view ascii::ascii_info_detail -v0 -1.5,-0.5,0.5,-0.5,-0.5,0.5,0.5,-0.5,0.5 -J 0.0,0.0,0.5,0.0,0.5,0.0,-0.5,0.0,0.0,0.0,0.0,0.5,0.0,0.5,0.0,-0.5,0.0,0.0,0.0,0.0,0.5,0.0,0.5,0.0,-0.5,0.0,0.0 -invJ 0.0,0.0,-2.0,0.0,2.0,0.0,2.0,0.0,0.0,0.0,0.0,-2.0,0.0,2.0,0.0,2.0,0.0,0.0,0.0,0.0,-2.0,0.0,2.0,0.0,2.0,0.0,0.0 -detJ 0.125,0.125,0.125 -centroid -1.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0 -normal 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 -vol 1.0,1.0,1.0 721 test: 722 suffix: 5 723 args: -run_type file -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 3,1,1 -dm_plex_box_lower -1.5,-0.5,-0.5 -dm_plex_box_upper 1.5,0.5,0.5 -dm_view ascii::ascii_info_detail -centroid -1.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0 -normal 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 -vol 1.0,1.0,1.0 724 test: 725 suffix: 6 726 args: -run_type file -dm_plex_dim 1 -dm_plex_simplex 0 -dm_plex_box_faces 3 -dm_plex_box_lower -1.5 -dm_plex_box_upper 1.5 -dm_view ascii::ascii_info_detail -centroid -1.0,0.0,1.0 -vol 1.0,1.0,1.0 727 TEST*/ 728