1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 #include <petsc/private/dtimpl.h> /*I "petscdt.h" I*/ 3 #include <petsc/private/dmpleximpl.h> /* For CellRefiner */ 4 #include <petscblaslapack.h> 5 6 static PetscErrorCode PetscFEDestroy_Composite(PetscFE fem) 7 { 8 PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data; 9 PetscErrorCode ierr; 10 11 PetscFunctionBegin; 12 ierr = CellRefinerRestoreAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);CHKERRQ(ierr); 13 ierr = PetscFree(cmp->embedding);CHKERRQ(ierr); 14 ierr = PetscFree(cmp);CHKERRQ(ierr); 15 PetscFunctionReturn(0); 16 } 17 18 static PetscErrorCode PetscFESetUp_Composite(PetscFE fem) 19 { 20 PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data; 21 DM K; 22 PetscReal *subpoint; 23 PetscBLASInt *pivots; 24 PetscBLASInt n, info; 25 PetscScalar *work, *invVscalar; 26 PetscInt dim, pdim, spdim, j, s; 27 PetscErrorCode ierr; 28 29 PetscFunctionBegin; 30 /* Get affine mapping from reference cell to each subcell */ 31 ierr = PetscDualSpaceGetDM(fem->dualSpace, &K);CHKERRQ(ierr); 32 ierr = DMGetDimension(K, &dim);CHKERRQ(ierr); 33 ierr = DMPlexGetCellRefiner_Internal(K, &cmp->cellRefiner);CHKERRQ(ierr); 34 ierr = CellRefinerGetAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);CHKERRQ(ierr); 35 /* Determine dof embedding into subelements */ 36 ierr = PetscDualSpaceGetDimension(fem->dualSpace, &pdim);CHKERRQ(ierr); 37 ierr = PetscSpaceGetDimension(fem->basisSpace, &spdim);CHKERRQ(ierr); 38 ierr = PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding);CHKERRQ(ierr); 39 ierr = DMGetWorkArray(K, dim, MPIU_REAL, &subpoint);CHKERRQ(ierr); 40 for (s = 0; s < cmp->numSubelements; ++s) { 41 PetscInt sd = 0; 42 43 for (j = 0; j < pdim; ++j) { 44 PetscBool inside; 45 PetscQuadrature f; 46 PetscInt d, e; 47 48 ierr = PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);CHKERRQ(ierr); 49 /* Apply transform to first point, and check that point is inside subcell */ 50 for (d = 0; d < dim; ++d) { 51 subpoint[d] = -1.0; 52 for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(f->points[e] - cmp->v0[s*dim+e]); 53 } 54 ierr = CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);CHKERRQ(ierr); 55 if (inside) {cmp->embedding[s*spdim+sd++] = j;} 56 } 57 if (sd != spdim) SETERRQ3(PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %d has %d dual basis vectors != %d", s, sd, spdim); 58 } 59 ierr = DMRestoreWorkArray(K, dim, MPIU_REAL, &subpoint);CHKERRQ(ierr); 60 /* Construct the change of basis from prime basis to nodal basis for each subelement */ 61 ierr = PetscMalloc1(cmp->numSubelements*spdim*spdim,&fem->invV);CHKERRQ(ierr); 62 ierr = PetscMalloc2(spdim,&pivots,spdim,&work);CHKERRQ(ierr); 63 #if defined(PETSC_USE_COMPLEX) 64 ierr = PetscMalloc1(cmp->numSubelements*spdim*spdim,&invVscalar);CHKERRQ(ierr); 65 #else 66 invVscalar = fem->invV; 67 #endif 68 for (s = 0; s < cmp->numSubelements; ++s) { 69 for (j = 0; j < spdim; ++j) { 70 PetscReal *Bf; 71 PetscQuadrature f; 72 const PetscReal *points, *weights; 73 PetscInt Nc, Nq, q, k; 74 75 ierr = PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s*spdim+j], &f);CHKERRQ(ierr); 76 ierr = PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);CHKERRQ(ierr); 77 ierr = PetscMalloc1(f->numPoints*spdim*Nc,&Bf);CHKERRQ(ierr); 78 ierr = PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL);CHKERRQ(ierr); 79 for (k = 0; k < spdim; ++k) { 80 /* n_j \cdot \phi_k */ 81 invVscalar[(s*spdim + j)*spdim+k] = 0.0; 82 for (q = 0; q < Nq; ++q) { 83 invVscalar[(s*spdim + j)*spdim+k] += Bf[q*spdim+k]*weights[q]; 84 } 85 } 86 ierr = PetscFree(Bf);CHKERRQ(ierr); 87 } 88 n = spdim; 89 PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &invVscalar[s*spdim*spdim], &n, pivots, &info)); 90 PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &invVscalar[s*spdim*spdim], &n, pivots, work, &n, &info)); 91 } 92 #if defined(PETSC_USE_COMPLEX) 93 for (s = 0; s <cmp->numSubelements*spdim*spdim; s++) fem->invV[s] = PetscRealPart(invVscalar[s]); 94 ierr = PetscFree(invVscalar);CHKERRQ(ierr); 95 #endif 96 ierr = PetscFree2(pivots,work);CHKERRQ(ierr); 97 PetscFunctionReturn(0); 98 } 99 100 static PetscErrorCode PetscFECreateTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T) 101 { 102 PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data; 103 DM dm; 104 PetscInt pdim; /* Dimension of FE space P */ 105 PetscInt spdim; /* Dimension of subelement FE space P */ 106 PetscInt dim; /* Spatial dimension */ 107 PetscInt comp; /* Field components */ 108 PetscInt *subpoints; 109 PetscReal *B = K >= 0 ? T->T[0] : NULL; 110 PetscReal *D = K >= 1 ? T->T[1] : NULL; 111 PetscReal *H = K >= 2 ? T->T[2] : NULL; 112 PetscReal *tmpB = NULL, *tmpD = NULL, *tmpH = NULL, *subpoint; 113 PetscInt p, s, d, e, j, k; 114 PetscErrorCode ierr; 115 116 PetscFunctionBegin; 117 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 118 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 119 ierr = PetscSpaceGetDimension(fem->basisSpace, &spdim);CHKERRQ(ierr); 120 ierr = PetscDualSpaceGetDimension(fem->dualSpace, &pdim);CHKERRQ(ierr); 121 ierr = PetscFEGetNumComponents(fem, &comp);CHKERRQ(ierr); 122 /* Divide points into subelements */ 123 ierr = DMGetWorkArray(dm, npoints, MPIU_INT, &subpoints);CHKERRQ(ierr); 124 ierr = DMGetWorkArray(dm, dim, MPIU_REAL, &subpoint);CHKERRQ(ierr); 125 for (p = 0; p < npoints; ++p) { 126 for (s = 0; s < cmp->numSubelements; ++s) { 127 PetscBool inside; 128 129 /* Apply transform, and check that point is inside cell */ 130 for (d = 0; d < dim; ++d) { 131 subpoint[d] = -1.0; 132 for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]); 133 } 134 ierr = CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);CHKERRQ(ierr); 135 if (inside) {subpoints[p] = s; break;} 136 } 137 if (s >= cmp->numSubelements) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d was not found in any subelement", p); 138 } 139 ierr = DMRestoreWorkArray(dm, dim, MPIU_REAL, &subpoint);CHKERRQ(ierr); 140 /* Evaluate the prime basis functions at all points */ 141 if (K >= 0) {ierr = DMGetWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);CHKERRQ(ierr);} 142 if (K >= 1) {ierr = DMGetWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);} 143 if (K >= 2) {ierr = DMGetWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);} 144 ierr = PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH);CHKERRQ(ierr); 145 /* Translate to the nodal basis */ 146 if (K >= 0) {ierr = PetscArrayzero(B, npoints*pdim*comp);CHKERRQ(ierr);} 147 if (K >= 1) {ierr = PetscArrayzero(D, npoints*pdim*comp*dim);CHKERRQ(ierr);} 148 if (K >= 2) {ierr = PetscArrayzero(H, npoints*pdim*comp*dim*dim);CHKERRQ(ierr);} 149 for (p = 0; p < npoints; ++p) { 150 const PetscInt s = subpoints[p]; 151 152 if (B) { 153 /* Multiply by V^{-1} (spdim x spdim) */ 154 for (j = 0; j < spdim; ++j) { 155 const PetscInt i = (p*pdim + cmp->embedding[s*spdim+j])*comp; 156 157 B[i] = 0.0; 158 for (k = 0; k < spdim; ++k) { 159 B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k]; 160 } 161 } 162 } 163 if (D) { 164 /* Multiply by V^{-1} (spdim x spdim) */ 165 for (j = 0; j < spdim; ++j) { 166 for (d = 0; d < dim; ++d) { 167 const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d; 168 169 D[i] = 0.0; 170 for (k = 0; k < spdim; ++k) { 171 D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d]; 172 } 173 } 174 } 175 } 176 if (H) { 177 /* Multiply by V^{-1} (pdim x pdim) */ 178 for (j = 0; j < spdim; ++j) { 179 for (d = 0; d < dim*dim; ++d) { 180 const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d; 181 182 H[i] = 0.0; 183 for (k = 0; k < spdim; ++k) { 184 H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d]; 185 } 186 } 187 } 188 } 189 } 190 ierr = DMRestoreWorkArray(dm, npoints, MPIU_INT, &subpoints);CHKERRQ(ierr); 191 if (K >= 0) {ierr = DMRestoreWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);CHKERRQ(ierr);} 192 if (K >= 1) {ierr = DMRestoreWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);} 193 if (K >= 2) {ierr = DMRestoreWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);} 194 PetscFunctionReturn(0); 195 } 196 197 static PetscErrorCode PetscFEInitialize_Composite(PetscFE fem) 198 { 199 PetscFunctionBegin; 200 fem->ops->setfromoptions = NULL; 201 fem->ops->setup = PetscFESetUp_Composite; 202 fem->ops->view = NULL; 203 fem->ops->destroy = PetscFEDestroy_Composite; 204 fem->ops->getdimension = PetscFEGetDimension_Basic; 205 fem->ops->createtabulation = PetscFECreateTabulation_Composite; 206 fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic; 207 fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic; 208 fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */; 209 fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic; 210 PetscFunctionReturn(0); 211 } 212 213 /*MC 214 PETSCFECOMPOSITE = "composite" - A PetscFE object that represents a composite element 215 216 Level: intermediate 217 218 .seealso: PetscFEType, PetscFECreate(), PetscFESetType() 219 M*/ 220 PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem) 221 { 222 PetscFE_Composite *cmp; 223 PetscErrorCode ierr; 224 225 PetscFunctionBegin; 226 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 227 ierr = PetscNewLog(fem, &cmp);CHKERRQ(ierr); 228 fem->data = cmp; 229 230 cmp->cellRefiner = REFINER_NOOP; 231 cmp->numSubelements = -1; 232 cmp->v0 = NULL; 233 cmp->jac = NULL; 234 235 ierr = PetscFEInitialize_Composite(fem);CHKERRQ(ierr); 236 PetscFunctionReturn(0); 237 } 238 239 /*@C 240 PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement 241 242 Not collective 243 244 Input Parameter: 245 . fem - The PetscFE object 246 247 Output Parameters: 248 + blockSize - The number of elements in a block 249 . numBlocks - The number of blocks in a batch 250 . batchSize - The number of elements in a batch 251 - numBatches - The number of batches in a chunk 252 253 Level: intermediate 254 255 .seealso: PetscFECreate() 256 @*/ 257 PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[]) 258 { 259 PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data; 260 261 PetscFunctionBegin; 262 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 263 if (numSubelements) {PetscValidPointer(numSubelements, 2); *numSubelements = cmp->numSubelements;} 264 if (v0) {PetscValidPointer(v0, 3); *v0 = cmp->v0;} 265 if (jac) {PetscValidPointer(jac, 4); *jac = cmp->jac;} 266 if (invjac) {PetscValidPointer(invjac, 5); *invjac = cmp->invjac;} 267 PetscFunctionReturn(0); 268 } 269