1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 3 /*@C 4 PetscFEGeomCreate - Create a PetscFEGeom object to manage geometry for a group of cells 5 6 Input Parameters: 7 + quad - A PetscQuadrature determining the tabulation 8 . numCells - The number of cells in the group 9 . dimEmbed - The coordinate dimension 10 - faceData - Flag to construct geometry data for the faces 11 12 Output Parameter: 13 . geom - The PetscFEGeom object 14 15 Level: beginner 16 17 .seealso: `PetscFEGeomDestroy()`, `PetscFEGeomComplete()` 18 @*/ 19 PetscErrorCode PetscFEGeomCreate(PetscQuadrature quad, PetscInt numCells, PetscInt dimEmbed, PetscBool faceData, PetscFEGeom **geom) { 20 PetscFEGeom *g; 21 PetscInt dim, Nq, N; 22 const PetscReal *p; 23 24 PetscFunctionBegin; 25 PetscCall(PetscQuadratureGetData(quad, &dim, NULL, &Nq, &p, NULL)); 26 PetscCall(PetscNew(&g)); 27 g->xi = p; 28 g->numCells = numCells; 29 g->numPoints = Nq; 30 g->dim = dim; 31 g->dimEmbed = dimEmbed; 32 g->isCohesive = PETSC_FALSE; 33 N = numCells * Nq; 34 PetscCall(PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ)); 35 if (faceData) { 36 PetscCall(PetscCalloc2(numCells, &g->face, N * dimEmbed, &g->n)); 37 PetscCall(PetscCalloc6(N * dimEmbed * dimEmbed, &(g->suppJ[0]), N * dimEmbed * dimEmbed, &(g->suppJ[1]), N * dimEmbed * dimEmbed, &(g->suppInvJ[0]), N * dimEmbed * dimEmbed, &(g->suppInvJ[1]), N, &(g->suppDetJ[0]), N, &(g->suppDetJ[1]))); 38 } 39 PetscCall(PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ)); 40 *geom = g; 41 PetscFunctionReturn(0); 42 } 43 44 /*@C 45 PetscFEGeomDestroy - Destroy a PetscFEGeom object 46 47 Input Parameter: 48 . geom - PetscFEGeom object 49 50 Level: beginner 51 52 .seealso: `PetscFEGeomCreate()` 53 @*/ 54 PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom) { 55 PetscFunctionBegin; 56 if (!*geom) PetscFunctionReturn(0); 57 PetscCall(PetscFree3((*geom)->v, (*geom)->J, (*geom)->detJ)); 58 PetscCall(PetscFree((*geom)->invJ)); 59 PetscCall(PetscFree2((*geom)->face, (*geom)->n)); 60 PetscCall(PetscFree6((*geom)->suppJ[0], (*geom)->suppJ[1], (*geom)->suppInvJ[0], (*geom)->suppInvJ[1], (*geom)->suppDetJ[0], (*geom)->suppDetJ[1])); 61 PetscCall(PetscFree(*geom)); 62 PetscFunctionReturn(0); 63 } 64 65 /*@C 66 PetscFEGeomGetChunk - Get a chunk of cells in the group as a PetscFEGeom 67 68 Input Parameters: 69 + geom - PetscFEGeom object 70 . cStart - The first cell in the chunk 71 - cEnd - The first cell not in the chunk 72 73 Output Parameter: 74 . chunkGeom - The chunk of cells 75 76 Level: intermediate 77 78 .seealso: `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()` 79 @*/ 80 PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) { 81 PetscInt Nq; 82 PetscInt dE; 83 84 PetscFunctionBegin; 85 PetscValidPointer(geom, 1); 86 PetscValidPointer(chunkGeom, 4); 87 if (!(*chunkGeom)) { PetscCall(PetscNew(chunkGeom)); } 88 Nq = geom->numPoints; 89 dE = geom->dimEmbed; 90 (*chunkGeom)->dim = geom->dim; 91 (*chunkGeom)->dimEmbed = geom->dimEmbed; 92 (*chunkGeom)->numPoints = geom->numPoints; 93 (*chunkGeom)->numCells = cEnd - cStart; 94 (*chunkGeom)->xi = geom->xi; 95 (*chunkGeom)->v = &geom->v[Nq * dE * cStart]; 96 (*chunkGeom)->J = &geom->J[Nq * dE * dE * cStart]; 97 (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq * dE * dE * cStart] : NULL; 98 (*chunkGeom)->detJ = &geom->detJ[Nq * cStart]; 99 (*chunkGeom)->n = geom->n ? &geom->n[Nq * dE * cStart] : NULL; 100 (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL; 101 (*chunkGeom)->suppJ[0] = geom->suppJ[0] ? &geom->suppJ[0][Nq * dE * dE * cStart] : NULL; 102 (*chunkGeom)->suppJ[1] = geom->suppJ[1] ? &geom->suppJ[1][Nq * dE * dE * cStart] : NULL; 103 (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq * dE * dE * cStart] : NULL; 104 (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq * dE * dE * cStart] : NULL; 105 (*chunkGeom)->suppDetJ[0] = geom->suppDetJ[0] ? &geom->suppDetJ[0][Nq * cStart] : NULL; 106 (*chunkGeom)->suppDetJ[1] = geom->suppDetJ[1] ? &geom->suppDetJ[1][Nq * cStart] : NULL; 107 (*chunkGeom)->isAffine = geom->isAffine; 108 PetscFunctionReturn(0); 109 } 110 111 /*@C 112 PetscFEGeomRestoreChunk - Restore the chunk 113 114 Input Parameters: 115 + geom - PetscFEGeom object 116 . cStart - The first cell in the chunk 117 . cEnd - The first cell not in the chunk 118 - chunkGeom - The chunk of cells 119 120 Level: intermediate 121 122 .seealso: `PetscFEGeomGetChunk()`, `PetscFEGeomCreate()` 123 @*/ 124 PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) { 125 PetscFunctionBegin; 126 PetscCall(PetscFree(*chunkGeom)); 127 PetscFunctionReturn(0); 128 } 129 130 /*@C 131 PetscFEGeomGetPoint - Get the geometry for cell c at point p as a PetscFEGeom 132 133 Input Parameters: 134 + geom - PetscFEGeom object 135 . c - The cell 136 . p - The point 137 - pcoords - The reference coordinates of point p, or NULL 138 139 Output Parameter: 140 . pgeom - The geometry of cell c at point p 141 142 Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom, 143 nothing needs to be done with it afterwards. 144 145 In the affine case, pgeom must have storage for the integration point coordinates in pgeom->v if pcoords is passed in. 146 147 Level: intermediate 148 149 .seealso: `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()` 150 @*/ 151 PetscErrorCode PetscFEGeomGetPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, const PetscReal pcoords[], PetscFEGeom *pgeom) { 152 const PetscInt dim = geom->dim; 153 const PetscInt dE = geom->dimEmbed; 154 const PetscInt Np = geom->numPoints; 155 156 PetscFunctionBeginHot; 157 pgeom->dim = dim; 158 pgeom->dimEmbed = dE; 159 //pgeom->isAffine = geom->isAffine; 160 if (geom->isAffine) { 161 if (!p) { 162 pgeom->xi = geom->xi; 163 pgeom->J = &geom->J[c * Np * dE * dE]; 164 pgeom->invJ = &geom->invJ[c * Np * dE * dE]; 165 pgeom->detJ = &geom->detJ[c * Np]; 166 pgeom->n = geom->n ? &geom->n[c * Np * dE] : NULL; 167 } 168 if (pcoords) { CoordinatesRefToReal(dE, dim, pgeom->xi, &geom->v[c * Np * dE], pgeom->J, pcoords, pgeom->v); } 169 } else { 170 pgeom->v = &geom->v[(c * Np + p) * dE]; 171 pgeom->J = &geom->J[(c * Np + p) * dE * dE]; 172 pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE]; 173 pgeom->detJ = &geom->detJ[c * Np + p]; 174 pgeom->n = geom->n ? &geom->n[(c * Np + p) * dE] : NULL; 175 } 176 PetscFunctionReturn(0); 177 } 178 179 /*@C 180 PetscFEGeomGetCellPoint - Get the cell geometry for face f at point p as a PetscFEGeom 181 182 Input Parameters: 183 + geom - PetscFEGeom object 184 . f - The face 185 - p - The point 186 187 Output Parameter: 188 . pgeom - The cell geometry of face f at point p 189 190 Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom, 191 nothing needs to be done with it afterwards. 192 193 Level: intermediate 194 195 .seealso: `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()` 196 @*/ 197 PetscErrorCode PetscFEGeomGetCellPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, PetscFEGeom *pgeom) { 198 const PetscBool bd = geom->dimEmbed > geom->dim && !geom->isCohesive ? PETSC_TRUE : PETSC_FALSE; 199 const PetscInt dim = bd ? geom->dimEmbed : geom->dim; 200 const PetscInt dE = geom->dimEmbed; 201 const PetscInt Np = geom->numPoints; 202 203 PetscFunctionBeginHot; 204 pgeom->dim = dim; 205 pgeom->dimEmbed = dE; 206 //pgeom->isAffine = geom->isAffine; 207 if (geom->isAffine) { 208 if (!p) { 209 if (bd) { 210 pgeom->J = &geom->suppJ[0][c * Np * dE * dE]; 211 pgeom->invJ = &geom->suppInvJ[0][c * Np * dE * dE]; 212 pgeom->detJ = &geom->suppDetJ[0][c * Np]; 213 } else { 214 pgeom->J = &geom->J[c * Np * dE * dE]; 215 pgeom->invJ = &geom->invJ[c * Np * dE * dE]; 216 pgeom->detJ = &geom->detJ[c * Np]; 217 } 218 } 219 } else { 220 if (bd) { 221 pgeom->J = &geom->suppJ[0][(c * Np + p) * dE * dE]; 222 pgeom->invJ = &geom->suppInvJ[0][(c * Np + p) * dE * dE]; 223 pgeom->detJ = &geom->suppDetJ[0][c * Np + p]; 224 } else { 225 pgeom->J = &geom->J[(c * Np + p) * dE * dE]; 226 pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE]; 227 pgeom->detJ = &geom->detJ[c * Np + p]; 228 } 229 } 230 PetscFunctionReturn(0); 231 } 232 233 /*@ 234 PetscFEGeomComplete - Calculate derived quntites from base geometry specification 235 236 Input Parameter: 237 . geom - PetscFEGeom object 238 239 Level: intermediate 240 241 .seealso: `PetscFEGeomCreate()` 242 @*/ 243 PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom) { 244 PetscInt i, j, N, dE; 245 246 PetscFunctionBeginHot; 247 N = geom->numPoints * geom->numCells; 248 dE = geom->dimEmbed; 249 switch (dE) { 250 case 3: 251 for (i = 0; i < N; i++) { 252 DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]); 253 if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]); 254 } 255 break; 256 case 2: 257 for (i = 0; i < N; i++) { 258 DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]); 259 if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]); 260 } 261 break; 262 case 1: 263 for (i = 0; i < N; i++) { 264 geom->detJ[i] = PetscAbsReal(geom->J[i]); 265 if (geom->invJ) geom->invJ[i] = 1. / geom->J[i]; 266 } 267 break; 268 } 269 if (geom->n) { 270 for (i = 0; i < N; i++) { 271 for (j = 0; j < dE; j++) { geom->n[dE * i + j] = geom->J[dE * dE * i + dE * j + dE - 1] * ((dE == 2) ? -1. : 1.); } 272 } 273 } 274 PetscFunctionReturn(0); 275 } 276