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