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