1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 3 PetscErrorCode PetscFEGeomCreate(PetscQuadrature quad, PetscInt numCells, PetscInt dimEmbed, PetscBool faceData, PetscFEGeom **geom) 4 { 5 PetscFEGeom *g; 6 PetscInt dim, Nq, N; 7 const PetscReal *p; 8 PetscErrorCode ierr; 9 10 PetscFunctionBegin; 11 ierr = PetscQuadratureGetData(quad,&dim,NULL,&Nq,&p,NULL);CHKERRQ(ierr); 12 ierr = PetscNew(&g);CHKERRQ(ierr); 13 g->xi = p; 14 g->numCells = numCells; 15 g->numPoints = Nq; 16 g->dim = dim; 17 g->dimEmbed = dimEmbed; 18 N = numCells * Nq; 19 ierr = PetscCalloc3(N * dimEmbed, &g->v, N * dimEmbed * dimEmbed, &g->J, N, &g->detJ);CHKERRQ(ierr); 20 if (faceData) { 21 ierr = PetscCalloc2(numCells, &g->face, N * dimEmbed, &g->n);CHKERRQ(ierr); 22 ierr = PetscCalloc6(N * dimEmbed * dimEmbed, &(g->suppJ[0]), N * dimEmbed * dimEmbed, &(g->suppJ[1]), 23 N * dimEmbed * dimEmbed, &(g->suppInvJ[0]), N * dimEmbed * dimEmbed, &(g->suppInvJ[1]), 24 N, &(g->suppDetJ[0]), N, &(g->suppDetJ[1]));CHKERRQ(ierr); 25 } 26 ierr = PetscCalloc1(N * dimEmbed * dimEmbed, &g->invJ);CHKERRQ(ierr); 27 *geom = g; 28 PetscFunctionReturn(0); 29 } 30 31 PetscErrorCode PetscFEGeomDestroy(PetscFEGeom **geom) 32 { 33 PetscErrorCode ierr; 34 35 PetscFunctionBegin; 36 if (!*geom) PetscFunctionReturn(0); 37 ierr = PetscFree3((*geom)->v,(*geom)->J,(*geom)->detJ);CHKERRQ(ierr); 38 ierr = PetscFree((*geom)->invJ);CHKERRQ(ierr); 39 ierr = PetscFree2((*geom)->face,(*geom)->n);CHKERRQ(ierr); 40 ierr = PetscFree6((*geom)->suppJ[0],(*geom)->suppJ[1],(*geom)->suppInvJ[0],(*geom)->suppInvJ[1],(*geom)->suppDetJ[0],(*geom)->suppDetJ[1]);CHKERRQ(ierr); 41 ierr = PetscFree(*geom);CHKERRQ(ierr); 42 PetscFunctionReturn(0); 43 } 44 45 PetscErrorCode PetscFEGeomGetChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 46 { 47 PetscInt Nq; 48 PetscInt dE; 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 PetscValidPointer(geom,1); 53 PetscValidPointer(chunkGeom,2); 54 if (!(*chunkGeom)) { 55 ierr = PetscNew(chunkGeom);CHKERRQ(ierr); 56 } 57 Nq = geom->numPoints; 58 dE= geom->dimEmbed; 59 (*chunkGeom)->dim = geom->dim; 60 (*chunkGeom)->dimEmbed = geom->dimEmbed; 61 (*chunkGeom)->numPoints = geom->numPoints; 62 (*chunkGeom)->numCells = cEnd - cStart; 63 (*chunkGeom)->xi = geom->xi; 64 (*chunkGeom)->v = &geom->v[Nq*dE*cStart]; 65 (*chunkGeom)->J = &geom->J[Nq*dE*dE*cStart]; 66 (*chunkGeom)->invJ = (geom->invJ) ? &geom->invJ[Nq*dE*dE*cStart] : NULL; 67 (*chunkGeom)->detJ = &geom->detJ[Nq*cStart]; 68 (*chunkGeom)->n = geom->n ? &geom->n[Nq*dE*cStart] : NULL; 69 (*chunkGeom)->face = geom->face ? &geom->face[cStart] : NULL; 70 (*chunkGeom)->suppJ[0] = geom->suppJ[0] ? &geom->suppJ[0][Nq*dE*dE*cStart] : NULL; 71 (*chunkGeom)->suppJ[1] = geom->suppJ[1] ? &geom->suppJ[1][Nq*dE*dE*cStart] : NULL; 72 (*chunkGeom)->suppInvJ[0] = geom->suppInvJ[0] ? &geom->suppInvJ[0][Nq*dE*dE*cStart] : NULL; 73 (*chunkGeom)->suppInvJ[1] = geom->suppInvJ[1] ? &geom->suppInvJ[1][Nq*dE*dE*cStart] : NULL; 74 (*chunkGeom)->suppDetJ[0] = geom->suppDetJ[0] ? &geom->suppDetJ[0][Nq*cStart] : NULL; 75 (*chunkGeom)->suppDetJ[1] = geom->suppDetJ[1] ? &geom->suppDetJ[1][Nq*cStart] : NULL; 76 (*chunkGeom)->isAffine = geom->isAffine; 77 PetscFunctionReturn(0); 78 } 79 80 PetscErrorCode PetscFEGeomRestoreChunk(PetscFEGeom *geom, PetscInt cStart, PetscInt cEnd, PetscFEGeom **chunkGeom) 81 { 82 PetscErrorCode ierr; 83 84 PetscFunctionBegin; 85 ierr = PetscFree(*chunkGeom);CHKERRQ(ierr); 86 PetscFunctionReturn(0); 87 } 88 89 PetscErrorCode PetscFEGeomComplete(PetscFEGeom *geom) 90 { 91 PetscInt i, j, N, dE; 92 93 PetscFunctionBeginHot; 94 N = geom->numPoints * geom->numCells; 95 dE = geom->dimEmbed; 96 switch (dE) { 97 case 3: 98 for (i = 0; i < N; i++) { 99 DMPlex_Det3D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]); 100 if (geom->invJ) DMPlex_Invert3D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]); 101 } 102 break; 103 case 2: 104 for (i = 0; i < N; i++) { 105 DMPlex_Det2D_Internal(&geom->detJ[i], &geom->J[dE*dE*i]); 106 if (geom->invJ) DMPlex_Invert2D_Internal(&geom->invJ[dE*dE*i], &geom->J[dE*dE*i], geom->detJ[i]); 107 } 108 break; 109 case 1: 110 for (i = 0; i < N; i++) { 111 geom->detJ[i] = PetscAbsReal(geom->J[i]); 112 if (geom->invJ) geom->invJ[i] = 1. / geom->J[i]; 113 } 114 break; 115 } 116 if (geom->n) { 117 for (i = 0; i < N; i++) { 118 for (j = 0; j < dE; j++) { 119 geom->n[dE*i + j] = geom->J[dE*dE*i + dE*j + dE-1] * ((dE == 2) ? -1. : 1.); 120 } 121 } 122 } 123 PetscFunctionReturn(0); 124 } 125