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,2); 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 /*@ 144 PetscFEGeomComplete - Calculate derived quntites from base geometry specification 145 146 Input Parameter: 147 . geom - PetscFEGeom object 148 149 Level: intermediate 150 151 .seealso: PetscFEGeomCreate() 152 @*/ 153 PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom) 154 { 155 PetscInt i, j, N, dE; 156 157 PetscFunctionBeginHot; 158 N = geom->numPoints * geom->numCells; 159 dE = geom->dimEmbed; 160 switch (dE) { 161 case 3: 162 for (i = 0; i < N; i++) { 163 DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]); 164 if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]); 165 } 166 break; 167 case 2: 168 for (i = 0; i < N; i++) { 169 DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]); 170 if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]); 171 } 172 break; 173 case 1: 174 for (i = 0; i < N; i++) { 175 geom->detJ[i] = PetscAbsReal(geom->J[i]); 176 if (geom->invJ) geom->invJ[i] = 1. / geom->J[i]; 177 } 178 break; 179 } 180 if (geom->n) { 181 for (i = 0; i < N; i++) { 182 for (j = 0; j < dE; j++) { 183 geom->n[dE*i + j] = geom->J[dE*dE*i + dE*j + dE-1] * ((dE == 2) ? -1. : 1.); 184 } 185 } 186 } 187 PetscFunctionReturn(0); 188 } 189