Actual source code: fecomposite.c

      1: #include <petsc/private/petscfeimpl.h>
      2: #include <petsc/private/dtimpl.h>
      3: #include <petscblaslapack.h>
      4: #include <petscdmplextransform.h>

      6: static PetscErrorCode PetscFEDestroy_Composite(PetscFE fem)
      7: {
      8:   PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data;

     10:   PetscFunctionBegin;
     11:   PetscCall(PetscFree(cmp->embedding));
     12:   PetscCall(PetscFree(cmp));
     13:   PetscFunctionReturn(PETSC_SUCCESS);
     14: }

     16: static PetscErrorCode PetscFESetUp_Composite(PetscFE fem)
     17: {
     18:   PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data;
     19:   DM                 K;
     20:   DMPolytopeType     ct;
     21:   DMPlexTransform    tr;
     22:   PetscReal         *subpoint;
     23:   PetscBLASInt      *pivots;
     24:   PetscBLASInt       n, info;
     25:   PetscScalar       *work, *invVscalar;
     26:   PetscInt           dim, pdim, spdim, j, s;
     27:   PetscSection       section;

     29:   PetscFunctionBegin;
     30:   /* Get affine mapping from reference cell to each subcell */
     31:   PetscCall(PetscDualSpaceGetDM(fem->dualSpace, &K));
     32:   PetscCall(DMGetDimension(K, &dim));
     33:   PetscCall(DMPlexGetCellType(K, 0, &ct));
     34:   PetscCall(DMPlexTransformCreate(PETSC_COMM_SELF, &tr));
     35:   PetscCall(DMPlexTransformSetType(tr, DMPLEXREFINEREGULAR));
     36:   PetscCall(DMPlexRefineRegularGetAffineTransforms(tr, ct, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac));
     37:   PetscCall(DMPlexTransformDestroy(&tr));
     38:   /* Determine dof embedding into subelements */
     39:   PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, &pdim));
     40:   PetscCall(PetscSpaceGetDimension(fem->basisSpace, &spdim));
     41:   PetscCall(PetscMalloc1(cmp->numSubelements * spdim, &cmp->embedding));
     42:   PetscCall(DMGetWorkArray(K, dim, MPIU_REAL, &subpoint));
     43:   PetscCall(PetscDualSpaceGetSection(fem->dualSpace, &section));
     44:   for (s = 0; s < cmp->numSubelements; ++s) {
     45:     PetscInt  sd = 0;
     46:     PetscInt  closureSize;
     47:     PetscInt *closure = NULL;

     49:     PetscCall(DMPlexGetTransitiveClosure(K, s, PETSC_TRUE, &closureSize, &closure));
     50:     for (j = 0; j < closureSize; j++) {
     51:       PetscInt point = closure[2 * j];
     52:       PetscInt dof, off, k;

     54:       PetscCall(PetscSectionGetDof(section, point, &dof));
     55:       PetscCall(PetscSectionGetOffset(section, point, &off));
     56:       for (k = 0; k < dof; k++) cmp->embedding[s * spdim + sd++] = off + k;
     57:     }
     58:     PetscCall(DMPlexRestoreTransitiveClosure(K, s, PETSC_TRUE, &closureSize, &closure));
     59:     PetscCheck(sd == spdim, PetscObjectComm((PetscObject)fem), PETSC_ERR_PLIB, "Subelement %" PetscInt_FMT " has %" PetscInt_FMT " dual basis vectors != %" PetscInt_FMT, s, sd, spdim);
     60:   }
     61:   PetscCall(DMRestoreWorkArray(K, dim, MPIU_REAL, &subpoint));
     62:   /* Construct the change of basis from prime basis to nodal basis for each subelement */
     63:   PetscCall(PetscMalloc1(cmp->numSubelements * spdim * spdim, &fem->invV));
     64:   PetscCall(PetscMalloc2(spdim, &pivots, spdim, &work));
     65: #if defined(PETSC_USE_COMPLEX)
     66:   PetscCall(PetscMalloc1(cmp->numSubelements * spdim * spdim, &invVscalar));
     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;

     77:       PetscCall(PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s * spdim + j], &f));
     78:       PetscCall(PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights));
     79:       PetscCall(PetscMalloc1(f->numPoints * spdim * Nc, &Bf));
     80:       PetscCall(PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL));
     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) invVscalar[(s * spdim + j) * spdim + k] += Bf[q * spdim + k] * weights[q];
     85:       }
     86:       PetscCall(PetscFree(Bf));
     87:     }
     88:     PetscCall(PetscBLASIntCast(spdim, &n));
     89:     PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &invVscalar[s * spdim * spdim], &n, pivots, &info));
     90:     PetscCallBLAS("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:   PetscCall(PetscFree(invVscalar));
     95: #endif
     96:   PetscCall(PetscFree2(pivots, work));
     97:   PetscFunctionReturn(PETSC_SUCCESS);
     98: }

    100: static PetscErrorCode PetscFEComputeTabulation_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:   DMPolytopeType     ct;
    105:   PetscInt           pdim;  /* Dimension of FE space P */
    106:   PetscInt           spdim; /* Dimension of subelement FE space P */
    107:   PetscInt           dim;   /* Spatial dimension */
    108:   PetscInt           comp;  /* Field components */
    109:   PetscInt          *subpoints;
    110:   PetscReal         *B    = K >= 0 ? T->T[0] : NULL;
    111:   PetscReal         *D    = K >= 1 ? T->T[1] : NULL;
    112:   PetscReal         *H    = K >= 2 ? T->T[2] : NULL;
    113:   PetscReal         *tmpB = NULL, *tmpD = NULL, *tmpH = NULL, *subpoint;
    114:   PetscInt           p, s, d, e, j, k;

    116:   PetscFunctionBegin;
    117:   PetscCall(PetscDualSpaceGetDM(fem->dualSpace, &dm));
    118:   PetscCall(DMGetDimension(dm, &dim));
    119:   PetscCall(DMPlexGetCellType(dm, 0, &ct));
    120:   PetscCall(PetscSpaceGetDimension(fem->basisSpace, &spdim));
    121:   PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, &pdim));
    122:   PetscCall(PetscFEGetNumComponents(fem, &comp));
    123:   /* Divide points into subelements */
    124:   PetscCall(DMGetWorkArray(dm, npoints, MPIU_INT, &subpoints));
    125:   PetscCall(DMGetWorkArray(dm, dim, MPIU_REAL, &subpoint));
    126:   for (p = 0; p < npoints; ++p) {
    127:     for (s = 0; s < cmp->numSubelements; ++s) {
    128:       PetscBool inside;

    130:       /* Apply transform, and check that point is inside cell */
    131:       for (d = 0; d < dim; ++d) {
    132:         subpoint[d] = -1.0;
    133:         for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s * dim + d) * dim + e] * (points[p * dim + e] - cmp->v0[s * dim + e]);
    134:       }
    135:       PetscCall(DMPolytopeInCellTest(ct, subpoint, &inside));
    136:       if (inside) {
    137:         subpoints[p] = s;
    138:         break;
    139:       }
    140:     }
    141:     PetscCheck(s < cmp->numSubelements, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " was not found in any subelement", p);
    142:   }
    143:   PetscCall(DMRestoreWorkArray(dm, dim, MPIU_REAL, &subpoint));
    144:   /* Evaluate the prime basis functions at all points */
    145:   if (K >= 0) PetscCall(DMGetWorkArray(dm, npoints * spdim, MPIU_REAL, &tmpB));
    146:   if (K >= 1) PetscCall(DMGetWorkArray(dm, npoints * spdim * dim, MPIU_REAL, &tmpD));
    147:   if (K >= 2) PetscCall(DMGetWorkArray(dm, npoints * spdim * dim * dim, MPIU_REAL, &tmpH));
    148:   PetscCall(PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH));
    149:   /* Translate to the nodal basis */
    150:   if (K >= 0) PetscCall(PetscArrayzero(B, npoints * pdim * comp));
    151:   if (K >= 1) PetscCall(PetscArrayzero(D, npoints * pdim * comp * dim));
    152:   if (K >= 2) PetscCall(PetscArrayzero(H, npoints * pdim * comp * dim * dim));
    153:   for (p = 0; p < npoints; ++p) {
    154:     const PetscInt s = subpoints[p];

    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;

    161:         B[i] = 0.0;
    162:         for (k = 0; k < spdim; ++k) B[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpB[p * spdim + k];
    163:       }
    164:     }
    165:     if (D) {
    166:       /* Multiply by V^{-1} (spdim x spdim) */
    167:       for (j = 0; j < spdim; ++j) {
    168:         for (d = 0; d < dim; ++d) {
    169:           const PetscInt i = ((p * pdim + cmp->embedding[s * spdim + j]) * comp + 0) * dim + d;

    171:           D[i] = 0.0;
    172:           for (k = 0; k < spdim; ++k) D[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpD[(p * spdim + k) * dim + d];
    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;

    182:           H[i] = 0.0;
    183:           for (k = 0; k < spdim; ++k) H[i] += fem->invV[(s * spdim + k) * spdim + j] * tmpH[(p * spdim + k) * dim * dim + d];
    184:         }
    185:       }
    186:     }
    187:   }
    188:   PetscCall(DMRestoreWorkArray(dm, npoints, MPIU_INT, &subpoints));
    189:   if (K >= 0) PetscCall(DMRestoreWorkArray(dm, npoints * spdim, MPIU_REAL, &tmpB));
    190:   if (K >= 1) PetscCall(DMRestoreWorkArray(dm, npoints * spdim * dim, MPIU_REAL, &tmpD));
    191:   if (K >= 2) PetscCall(DMRestoreWorkArray(dm, npoints * spdim * dim * dim, MPIU_REAL, &tmpH));
    192:   PetscFunctionReturn(PETSC_SUCCESS);
    193: }

    195: static PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
    196: {
    197:   PetscFunctionBegin;
    198:   fem->ops->setfromoptions          = NULL;
    199:   fem->ops->setup                   = PetscFESetUp_Composite;
    200:   fem->ops->view                    = NULL;
    201:   fem->ops->destroy                 = PetscFEDestroy_Composite;
    202:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
    203:   fem->ops->computetabulation       = PetscFEComputeTabulation_Composite;
    204:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
    205:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
    206:   fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */;
    207:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
    208:   PetscFunctionReturn(PETSC_SUCCESS);
    209: }

    211: /*MC
    212:   PETSCFECOMPOSITE = "composite" - A `PetscFEType` that represents a composite element

    214:   Level: intermediate

    216: .seealso: `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
    217: M*/
    218: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
    219: {
    220:   PetscFE_Composite *cmp;

    222:   PetscFunctionBegin;
    224:   PetscCall(PetscNew(&cmp));
    225:   fem->data = cmp;

    227:   cmp->numSubelements = -1;
    228:   cmp->v0             = NULL;
    229:   cmp->jac            = NULL;

    231:   PetscCall(PetscFEInitialize_Composite(fem));
    232:   PetscFunctionReturn(PETSC_SUCCESS);
    233: }

    235: /*@C
    236:   PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement

    238:   Not Collective

    240:   Input Parameter:
    241: . fem - The `PetscFE` object

    243:   Output Parameters:
    244: + numSubelements - The number of sub elements
    245: . v0             - The affine transformation for each element, an array of length $dim * Nc$. Pass `NULL` to ignore.
    246: . jac            - The Jacobian for each element, an array of length $dim^2 * Nc$. Pass `NULL` to ignore.
    247: - invjac         - The inverse of the Jacobian, an array of length $dim^2 * Nc$. Pass `NULL` to ignore.

    249:   Level: intermediate

    251:   Note:
    252:   Do not free the output arrays.

    254: .seealso: `PetscFE`, `PetscFECreate()`
    255: @*/
    256: PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PeOp PetscInt *numSubelements, PeOp const PetscReal *v0[], PeOp const PetscReal *jac[], PeOp const PetscReal *invjac[])
    257: {
    258:   PetscFE_Composite *cmp = (PetscFE_Composite *)fem->data;

    260:   PetscFunctionBegin;
    262:   if (numSubelements) {
    263:     PetscAssertPointer(numSubelements, 2);
    264:     *numSubelements = cmp->numSubelements;
    265:   }
    266:   if (v0) {
    267:     PetscAssertPointer(v0, 3);
    268:     *v0 = cmp->v0;
    269:   }
    270:   if (jac) {
    271:     PetscAssertPointer(jac, 4);
    272:     *jac = cmp->jac;
    273:   }
    274:   if (invjac) {
    275:     PetscAssertPointer(invjac, 5);
    276:     *invjac = cmp->invjac;
    277:   }
    278:   PetscFunctionReturn(PETSC_SUCCESS);
    279: }