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