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 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 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 PetscErrorCode PetscFEGetTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H) 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 *tmpB, *tmpD, *tmpH, *subpoint; 110 PetscInt p, s, d, e, j, k; 111 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr); 115 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 116 ierr = PetscSpaceGetDimension(fem->basisSpace, &spdim);CHKERRQ(ierr); 117 ierr = PetscDualSpaceGetDimension(fem->dualSpace, &pdim);CHKERRQ(ierr); 118 ierr = PetscFEGetNumComponents(fem, &comp);CHKERRQ(ierr); 119 /* Divide points into subelements */ 120 ierr = DMGetWorkArray(dm, npoints, MPIU_INT, &subpoints);CHKERRQ(ierr); 121 ierr = DMGetWorkArray(dm, dim, MPIU_REAL, &subpoint);CHKERRQ(ierr); 122 for (p = 0; p < npoints; ++p) { 123 for (s = 0; s < cmp->numSubelements; ++s) { 124 PetscBool inside; 125 126 /* Apply transform, and check that point is inside cell */ 127 for (d = 0; d < dim; ++d) { 128 subpoint[d] = -1.0; 129 for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]); 130 } 131 ierr = CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);CHKERRQ(ierr); 132 if (inside) {subpoints[p] = s; break;} 133 } 134 if (s >= cmp->numSubelements) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d was not found in any subelement", p); 135 } 136 ierr = DMRestoreWorkArray(dm, dim, MPIU_REAL, &subpoint);CHKERRQ(ierr); 137 /* Evaluate the prime basis functions at all points */ 138 if (B) {ierr = DMGetWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);CHKERRQ(ierr);} 139 if (D) {ierr = DMGetWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);} 140 if (H) {ierr = DMGetWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);} 141 ierr = PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);CHKERRQ(ierr); 142 /* Translate to the nodal basis */ 143 if (B) {ierr = PetscMemzero(B, npoints*pdim*comp * sizeof(PetscReal));CHKERRQ(ierr);} 144 if (D) {ierr = PetscMemzero(D, npoints*pdim*comp*dim * sizeof(PetscReal));CHKERRQ(ierr);} 145 if (H) {ierr = PetscMemzero(H, npoints*pdim*comp*dim*dim * sizeof(PetscReal));CHKERRQ(ierr);} 146 for (p = 0; p < npoints; ++p) { 147 const PetscInt s = subpoints[p]; 148 149 if (B) { 150 /* Multiply by V^{-1} (spdim x spdim) */ 151 for (j = 0; j < spdim; ++j) { 152 const PetscInt i = (p*pdim + cmp->embedding[s*spdim+j])*comp; 153 154 B[i] = 0.0; 155 for (k = 0; k < spdim; ++k) { 156 B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k]; 157 } 158 } 159 } 160 if (D) { 161 /* Multiply by V^{-1} (spdim x spdim) */ 162 for (j = 0; j < spdim; ++j) { 163 for (d = 0; d < dim; ++d) { 164 const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d; 165 166 D[i] = 0.0; 167 for (k = 0; k < spdim; ++k) { 168 D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d]; 169 } 170 } 171 } 172 } 173 if (H) { 174 /* Multiply by V^{-1} (pdim x pdim) */ 175 for (j = 0; j < spdim; ++j) { 176 for (d = 0; d < dim*dim; ++d) { 177 const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d; 178 179 H[i] = 0.0; 180 for (k = 0; k < spdim; ++k) { 181 H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d]; 182 } 183 } 184 } 185 } 186 } 187 ierr = DMRestoreWorkArray(dm, npoints, MPIU_INT, &subpoints);CHKERRQ(ierr); 188 if (B) {ierr = DMRestoreWorkArray(dm, npoints*spdim, MPIU_REAL, &tmpB);CHKERRQ(ierr);} 189 if (D) {ierr = DMRestoreWorkArray(dm, npoints*spdim*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);} 190 if (H) {ierr = DMRestoreWorkArray(dm, npoints*spdim*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);} 191 PetscFunctionReturn(0); 192 } 193 194 PetscErrorCode PetscFEInitialize_Composite(PetscFE fem) 195 { 196 PetscFunctionBegin; 197 fem->ops->setfromoptions = NULL; 198 fem->ops->setup = PetscFESetUp_Composite; 199 fem->ops->view = NULL; 200 fem->ops->destroy = PetscFEDestroy_Composite; 201 fem->ops->getdimension = PetscFEGetDimension_Basic; 202 fem->ops->gettabulation = PetscFEGetTabulation_Composite; 203 fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic; 204 fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic; 205 fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */; 206 fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic; 207 PetscFunctionReturn(0); 208 } 209 210 /*MC 211 PETSCFECOMPOSITE = "composite" - A PetscFE object that represents a composite element 212 213 Level: intermediate 214 215 .seealso: PetscFEType, PetscFECreate(), PetscFESetType() 216 M*/ 217 PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem) 218 { 219 PetscFE_Composite *cmp; 220 PetscErrorCode ierr; 221 222 PetscFunctionBegin; 223 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 224 ierr = PetscNewLog(fem, &cmp);CHKERRQ(ierr); 225 fem->data = cmp; 226 227 cmp->cellRefiner = REFINER_NOOP; 228 cmp->numSubelements = -1; 229 cmp->v0 = NULL; 230 cmp->jac = NULL; 231 232 ierr = PetscFEInitialize_Composite(fem);CHKERRQ(ierr); 233 PetscFunctionReturn(0); 234 } 235 236 /*@C 237 PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement 238 239 Not collective 240 241 Input Parameter: 242 . fem - The PetscFE object 243 244 Output Parameters: 245 + blockSize - The number of elements in a block 246 . numBlocks - The number of blocks in a batch 247 . batchSize - The number of elements in a batch 248 - numBatches - The number of batches in a chunk 249 250 Level: intermediate 251 252 .seealso: PetscFECreate() 253 @*/ 254 PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[]) 255 { 256 PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data; 257 258 PetscFunctionBegin; 259 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1); 260 if (numSubelements) {PetscValidPointer(numSubelements, 2); *numSubelements = cmp->numSubelements;} 261 if (v0) {PetscValidPointer(v0, 3); *v0 = cmp->v0;} 262 if (jac) {PetscValidPointer(jac, 4); *jac = cmp->jac;} 263 if (invjac) {PetscValidPointer(invjac, 5); *invjac = cmp->invjac;} 264 PetscFunctionReturn(0); 265 } 266