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