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