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 { 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(0); 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: `PetscFEGeomCreate()` 54 @*/ 55 PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom) 56 { 57 PetscFunctionBegin; 58 if (!*geom) PetscFunctionReturn(0); 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(0); 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 - The chunk of cells 77 78 Level: intermediate 79 80 .seealso: `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()` 81 @*/ 82 PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 83 { 84 PetscInt Nq; 85 PetscInt dE; 86 87 PetscFunctionBegin; 88 PetscValidPointer(geom, 1); 89 PetscValidPointer(chunkGeom, 4); 90 if (!(*chunkGeom)) PetscCall(PetscNew(chunkGeom)); 91 Nq = geom->numPoints; 92 dE = geom->dimEmbed; 93 (*chunkGeom)->dim = geom->dim; 94 (*chunkGeom)->dimEmbed = geom->dimEmbed; 95 (*chunkGeom)->numPoints = geom->numPoints; 96 (*chunkGeom)->numCells = cEnd - cStart; 97 (*chunkGeom)->xi = geom->xi; 98 (*chunkGeom)->v = &geom->v[Nq * dE * cStart]; 99 (*chunkGeom)->J = &geom->J[Nq * dE * dE * cStart]; 100 (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq * dE * dE * cStart] : NULL; 101 (*chunkGeom)->detJ = &geom->detJ[Nq * cStart]; 102 (*chunkGeom)->n = geom->n ? &geom->n[Nq * dE * cStart] : NULL; 103 (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL; 104 (*chunkGeom)->suppJ[0] = geom->suppJ[0] ? &geom->suppJ[0][Nq * dE * dE * cStart] : NULL; 105 (*chunkGeom)->suppJ[1] = geom->suppJ[1] ? &geom->suppJ[1][Nq * dE * dE * cStart] : NULL; 106 (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq * dE * dE * cStart] : NULL; 107 (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq * dE * dE * cStart] : NULL; 108 (*chunkGeom)->suppDetJ[0] = geom->suppDetJ[0] ? &geom->suppDetJ[0][Nq * cStart] : NULL; 109 (*chunkGeom)->suppDetJ[1] = geom->suppDetJ[1] ? &geom->suppDetJ[1][Nq * cStart] : NULL; 110 (*chunkGeom)->isAffine = geom->isAffine; 111 PetscFunctionReturn(0); 112 } 113 114 /*@C 115 PetscFEGeomRestoreChunk - Restore the chunk 116 117 Input Parameters: 118 + geom - PetscFEGeom object 119 . cStart - The first cell in the chunk 120 . cEnd - The first cell not in the chunk 121 - chunkGeom - The chunk of cells 122 123 Level: intermediate 124 125 .seealso: `PetscFEGeomGetChunk()`, `PetscFEGeomCreate()` 126 @*/ 127 PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 128 { 129 PetscFunctionBegin; 130 PetscCall(PetscFree(*chunkGeom)); 131 PetscFunctionReturn(0); 132 } 133 134 /*@C 135 PetscFEGeomGetPoint - Get the geometry for cell c at point p as a PetscFEGeom 136 137 Input Parameters: 138 + geom - PetscFEGeom object 139 . c - The cell 140 . p - The point 141 - pcoords - The reference coordinates of point p, or NULL 142 143 Output Parameter: 144 . pgeom - The geometry of cell c at point p 145 146 Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom, 147 nothing needs to be done with it afterwards. 148 149 In the affine case, pgeom must have storage for the integration point coordinates in pgeom->v if pcoords is passed in. 150 151 Level: intermediate 152 153 .seealso: `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()` 154 @*/ 155 PetscErrorCode PetscFEGeomGetPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, const PetscReal pcoords[], PetscFEGeom *pgeom) 156 { 157 const PetscInt dim = geom->dim; 158 const PetscInt dE = geom->dimEmbed; 159 const PetscInt Np = geom->numPoints; 160 161 PetscFunctionBeginHot; 162 pgeom->dim = dim; 163 pgeom->dimEmbed = dE; 164 //pgeom->isAffine = geom->isAffine; 165 if (geom->isAffine) { 166 if (!p) { 167 pgeom->xi = geom->xi; 168 pgeom->J = &geom->J[c * Np * dE * dE]; 169 pgeom->invJ = &geom->invJ[c * Np * dE * dE]; 170 pgeom->detJ = &geom->detJ[c * Np]; 171 pgeom->n = geom->n ? &geom->n[c * Np * dE] : NULL; 172 } 173 if (pcoords) CoordinatesRefToReal(dE, dim, pgeom->xi, &geom->v[c * Np * dE], pgeom->J, pcoords, pgeom->v); 174 } else { 175 pgeom->v = &geom->v[(c * Np + p) * dE]; 176 pgeom->J = &geom->J[(c * Np + p) * dE * dE]; 177 pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE]; 178 pgeom->detJ = &geom->detJ[c * Np + p]; 179 pgeom->n = geom->n ? &geom->n[(c * Np + p) * dE] : NULL; 180 } 181 PetscFunctionReturn(0); 182 } 183 184 /*@C 185 PetscFEGeomGetCellPoint - Get the cell geometry for face f at point p as a PetscFEGeom 186 187 Input Parameters: 188 + geom - PetscFEGeom object 189 . f - The face 190 - p - The point 191 192 Output Parameter: 193 . pgeom - The cell geometry of face f at point p 194 195 Note: For affine geometries, this only copies to pgeom at point 0. Since we copy pointers into pgeom, 196 nothing needs to be done with it afterwards. 197 198 Level: intermediate 199 200 .seealso: `PetscFEGeomRestoreChunk()`, `PetscFEGeomCreate()` 201 @*/ 202 PetscErrorCode PetscFEGeomGetCellPoint(PetscFEGeom *geom, PetscInt c, PetscInt p, PetscFEGeom *pgeom) 203 { 204 const PetscBool bd = geom->dimEmbed > geom->dim && !geom->isCohesive ? PETSC_TRUE : PETSC_FALSE; 205 const PetscInt dim = bd ? geom->dimEmbed : geom->dim; 206 const PetscInt dE = geom->dimEmbed; 207 const PetscInt Np = geom->numPoints; 208 209 PetscFunctionBeginHot; 210 pgeom->dim = dim; 211 pgeom->dimEmbed = dE; 212 //pgeom->isAffine = geom->isAffine; 213 if (geom->isAffine) { 214 if (!p) { 215 if (bd) { 216 pgeom->J = &geom->suppJ[0][c * Np * dE * dE]; 217 pgeom->invJ = &geom->suppInvJ[0][c * Np * dE * dE]; 218 pgeom->detJ = &geom->suppDetJ[0][c * Np]; 219 } else { 220 pgeom->J = &geom->J[c * Np * dE * dE]; 221 pgeom->invJ = &geom->invJ[c * Np * dE * dE]; 222 pgeom->detJ = &geom->detJ[c * Np]; 223 } 224 } 225 } else { 226 if (bd) { 227 pgeom->J = &geom->suppJ[0][(c * Np + p) * dE * dE]; 228 pgeom->invJ = &geom->suppInvJ[0][(c * Np + p) * dE * dE]; 229 pgeom->detJ = &geom->suppDetJ[0][c * Np + p]; 230 } else { 231 pgeom->J = &geom->J[(c * Np + p) * dE * dE]; 232 pgeom->invJ = &geom->invJ[(c * Np + p) * dE * dE]; 233 pgeom->detJ = &geom->detJ[c * Np + p]; 234 } 235 } 236 PetscFunctionReturn(0); 237 } 238 239 /*@ 240 PetscFEGeomComplete - Calculate derived quntites from base geometry specification 241 242 Input Parameter: 243 . geom - PetscFEGeom object 244 245 Level: intermediate 246 247 .seealso: `PetscFEGeomCreate()` 248 @*/ 249 PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom) 250 { 251 PetscInt i, j, N, dE; 252 253 PetscFunctionBeginHot; 254 N = geom->numPoints * geom->numCells; 255 dE = geom->dimEmbed; 256 switch (dE) { 257 case 3: 258 for (i = 0; i < N; i++) { 259 DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]); 260 if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]); 261 } 262 break; 263 case 2: 264 for (i = 0; i < N; i++) { 265 DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE * dE * i]); 266 if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE * dE * i], &geom->J[dE * dE * i], geom->detJ[i]); 267 } 268 break; 269 case 1: 270 for (i = 0; i < N; i++) { 271 geom->detJ[i] = PetscAbsReal(geom->J[i]); 272 if (geom->invJ) geom->invJ[i] = 1. / geom->J[i]; 273 } 274 break; 275 } 276 if (geom->n) { 277 for (i = 0; i < N; i++) { 278 for (j = 0; j < dE; j++) geom->n[dE * i + j] = geom->J[dE * dE * i + dE * j + dE - 1] * ((dE == 2) ? -1. : 1.); 279 } 280 } 281 PetscFunctionReturn(0); 282 } 283