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