Actual source code: fevector.c

      1: #include <petsc/private/petscfeimpl.h>
      2: #include <petsc/private/petscimpl.h>

      4: typedef struct _n_PetscFE_Vec {
      5:   PetscFE   scalar_fe;
      6:   PetscInt  num_copies;
      7:   PetscBool interleave_basis;
      8:   PetscBool interleave_components;
      9: } PetscFE_Vec;

     11: static PetscErrorCode PetscFEDestroy_Vector(PetscFE fe)
     12: {
     13:   PetscFE_Vec *v;

     15:   PetscFunctionBegin;
     16:   v = (PetscFE_Vec *)fe->data;
     17:   PetscCall(PetscFEDestroy(&v->scalar_fe));
     18:   PetscCall(PetscFree(v));
     19:   PetscFunctionReturn(PETSC_SUCCESS);
     20: }

     22: static PetscErrorCode PetscFEView_Vector_Ascii(PetscFE fe, PetscViewer v)
     23: {
     24:   PetscInt          dim, Nc, scalar_Nc;
     25:   PetscSpace        basis = NULL;
     26:   PetscDualSpace    dual  = NULL;
     27:   PetscQuadrature   quad  = NULL;
     28:   PetscFE_Vec      *vec;
     29:   PetscViewerFormat fmt;

     31:   PetscFunctionBegin;
     32:   vec = (PetscFE_Vec *)fe->data;
     33:   PetscCall(PetscFEGetSpatialDimension(fe, &dim));
     34:   PetscCall(PetscFEGetNumComponents(fe, &Nc));
     35:   PetscCall(PetscFEGetNumComponents(vec->scalar_fe, &scalar_Nc));
     36:   PetscCall(PetscFEGetBasisSpace(fe, &basis));
     37:   PetscCall(PetscFEGetDualSpace(fe, &dual));
     38:   PetscCall(PetscFEGetQuadrature(fe, &quad));
     39:   PetscCall(PetscViewerGetFormat(v, &fmt));
     40:   PetscCall(PetscViewerASCIIPushTab(v));
     41:   if (scalar_Nc == 1) {
     42:     PetscCall(PetscViewerASCIIPrintf(v, "Vector Finite Element in %" PetscInt_FMT " dimensions with %" PetscInt_FMT " components\n", dim, Nc));
     43:   } else {
     44:     PetscCall(PetscViewerASCIIPrintf(v, "Vector Finite Element in %" PetscInt_FMT " dimensions with %" PetscInt_FMT " components (%" PetscInt_FMT " copies of finite element with %" PetscInt_FMT " components)\n", dim, Nc, vec->num_copies, scalar_Nc));
     45:   }
     46:   if (fmt == PETSC_VIEWER_ASCII_INFO_DETAIL) {
     47:     PetscCall(PetscViewerASCIIPrintf(v, "Original finite element:\n"));
     48:     PetscCall(PetscViewerASCIIPushTab(v));
     49:     PetscCall(PetscFEView(vec->scalar_fe, v));
     50:     PetscCall(PetscViewerASCIIPopTab(v));
     51:   }
     52:   if (basis) PetscCall(PetscSpaceView(basis, v));
     53:   if (dual) PetscCall(PetscDualSpaceView(dual, v));
     54:   if (quad) PetscCall(PetscQuadratureView(quad, v));
     55:   PetscCall(PetscViewerASCIIPopTab(v));
     56:   PetscFunctionReturn(PETSC_SUCCESS);
     57: }

     59: static PetscErrorCode PetscFEView_Vector(PetscFE fe, PetscViewer v)
     60: {
     61:   PetscBool isascii;

     63:   PetscFunctionBegin;
     64:   PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
     65:   if (isascii) PetscCall(PetscFEView_Vector_Ascii(fe, v));
     66:   PetscFunctionReturn(PETSC_SUCCESS);
     67: }

     69: static PetscErrorCode PetscFESetUp_Vector(PetscFE fe)
     70: {
     71:   PetscFE_Vec        *v = (PetscFE_Vec *)fe->data;
     72:   PetscDualSpace      dsp;
     73:   PetscInt            n, Ncopies = v->num_copies;
     74:   PetscInt            scalar_n;
     75:   PetscInt           *d, *d_mapped;
     76:   PetscDualSpace_Sum *sum;
     77:   PetscBool           is_sum;

     79:   PetscFunctionBegin;
     80:   PetscCall(PetscFESetUp(v->scalar_fe));
     81:   PetscCall(PetscFEGetDimension(v->scalar_fe, &scalar_n));
     82:   PetscCall(PetscFEGetDualSpace(fe, &dsp));
     83:   PetscCall(PetscObjectTypeCompare((PetscObject)dsp, PETSCDUALSPACESUM, &is_sum));
     84:   PetscCheck(is_sum, PetscObjectComm((PetscObject)fe), PETSC_ERR_ARG_INCOMP, "Expected PETSCDUALSPACESUM dual space");
     85:   sum = (PetscDualSpace_Sum *)dsp->data;
     86:   n   = Ncopies * scalar_n;
     87:   PetscCall(PetscCalloc1(n * n, &fe->invV));
     88:   PetscCall(PetscMalloc2(scalar_n, &d, scalar_n, &d_mapped));
     89:   for (PetscInt i = 0; i < scalar_n; i++) d[i] = i;
     90:   for (PetscInt c = 0; c < Ncopies; c++) {
     91:     PetscCall(ISLocalToGlobalMappingApply(sum->all_rows[c], scalar_n, d, d_mapped));
     92:     for (PetscInt i = 0; i < scalar_n; i++) {
     93:       PetscInt         iw      = d_mapped[i];
     94:       PetscReal       *row_w   = &fe->invV[iw * n];
     95:       const PetscReal *row_r   = &v->scalar_fe->invV[i * scalar_n];
     96:       PetscInt         j0      = v->interleave_basis ? c : c * scalar_n;
     97:       PetscInt         jstride = v->interleave_basis ? Ncopies : 1;

     99:       for (PetscInt j = 0; j < scalar_n; j++) row_w[j0 + j * jstride] = row_r[j];
    100:     }
    101:   }
    102:   PetscCall(PetscFree2(d, d_mapped));
    103:   PetscFunctionReturn(PETSC_SUCCESS);
    104: }

    106: static PetscErrorCode PetscFEGetDimension_Vector(PetscFE fe, PetscInt *dim)
    107: {
    108:   PetscFE_Vec *v = (PetscFE_Vec *)fe->data;

    110:   PetscFunctionBegin;
    111:   PetscCall(PetscFEGetDimension(v->scalar_fe, dim));
    112:   *dim *= v->num_copies;
    113:   PetscFunctionReturn(PETSC_SUCCESS);
    114: }

    116: static PetscErrorCode PetscFEVectorInsertTabulation(PetscFE fe, PetscInt npoints, const PetscReal points[], PetscInt k, PetscInt scalar_Nb, PetscInt scalar_point_stride, const PetscReal scalar_Tk[], PetscReal Tk[])
    117: {
    118:   PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
    119:   PetscInt     scalar_Nc;
    120:   PetscInt     Nc;
    121:   PetscInt     cdim;
    122:   PetscInt     dblock;
    123:   PetscInt     block;
    124:   PetscInt     scalar_block;

    126:   PetscFunctionBegin;
    127:   PetscCall(PetscFEGetNumComponents(v->scalar_fe, &scalar_Nc));
    128:   Nc = scalar_Nc * v->num_copies;
    129:   PetscCall(PetscFEGetSpatialDimension(v->scalar_fe, &cdim));
    130:   dblock       = PetscPowInt(cdim, k);
    131:   block        = Nc * dblock;
    132:   scalar_block = scalar_Nc * dblock;
    133:   for (PetscInt p = 0; p < npoints; p++) {
    134:     const PetscReal *scalar_Tp = &scalar_Tk[p * scalar_point_stride];
    135:     PetscReal       *Tp        = &Tk[p * scalar_point_stride * v->num_copies * v->num_copies];

    137:     for (PetscInt j = 0; j < v->num_copies; j++) {
    138:       for (PetscInt scalar_i = 0; scalar_i < scalar_Nb; scalar_i++) {
    139:         PetscInt         i         = v->interleave_basis ? (scalar_i * v->num_copies + j) : (j * scalar_Nb + scalar_i);
    140:         const PetscReal *scalar_Ti = &scalar_Tp[scalar_i * scalar_block];
    141:         PetscReal       *Ti        = &Tp[i * block];

    143:         for (PetscInt scalar_c = 0; scalar_c < scalar_Nc; scalar_c++) {
    144:           PetscInt         c         = v->interleave_components ? (scalar_c * v->num_copies + j) : (j * scalar_Nc + scalar_c);
    145:           const PetscReal *scalar_Tc = &scalar_Ti[scalar_c * dblock];
    146:           PetscReal       *Tc        = &Ti[c * dblock];

    148:           PetscCall(PetscArraycpy(Tc, scalar_Tc, dblock));
    149:         }
    150:       }
    151:     }
    152:   }
    153:   PetscFunctionReturn(PETSC_SUCCESS);
    154: }

    156: static PetscErrorCode PetscFEComputeTabulation_Vector(PetscFE fe, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
    157: {
    158:   PetscFE_Vec    *v = (PetscFE_Vec *)fe->data;
    159:   PetscInt        scalar_Nc;
    160:   PetscInt        scalar_Nb;
    161:   PetscInt        cdim;
    162:   PetscTabulation scalar_T;

    164:   PetscFunctionBegin;
    165:   PetscAssert(npoints == T->Nr * T->Np, PetscObjectComm((PetscObject)fe), PETSC_ERR_PLIB, "Expected to be able decode PetscFECreateTabulation() from PetscTabulation fields");
    166:   PetscCall(PetscFECreateTabulation(v->scalar_fe, T->Nr, T->Np, points, K, &scalar_T));
    167:   PetscCall(PetscFEGetNumComponents(v->scalar_fe, &scalar_Nc));
    168:   PetscCall(PetscFEGetDimension(v->scalar_fe, &scalar_Nb));
    169:   PetscCall(PetscFEGetSpatialDimension(v->scalar_fe, &cdim));
    170:   for (PetscInt k = 0; k <= T->K; k++) {
    171:     PetscReal       *Tk                  = T->T[k];
    172:     const PetscReal *scalar_Tk           = scalar_T->T[k];
    173:     PetscInt         dblock              = PetscPowInt(cdim, k);
    174:     PetscInt         scalar_block        = scalar_Nc * dblock;
    175:     PetscInt         scalar_point_stride = scalar_Nb * scalar_block;

    177:     if (v->interleave_basis) {
    178:       PetscCall(PetscFEVectorInsertTabulation(fe, npoints, points, k, scalar_Nb, scalar_point_stride, scalar_Tk, Tk));
    179:     } else {
    180:       PetscDualSpace scalar_dsp;
    181:       PetscSection   ref_section;
    182:       PetscInt       pStart, pEnd;

    184:       PetscCall(PetscFEGetDualSpace(v->scalar_fe, &scalar_dsp));
    185:       PetscCall(PetscDualSpaceGetSection(scalar_dsp, &ref_section));
    186:       PetscCall(PetscSectionGetChart(ref_section, &pStart, &pEnd));
    187:       for (PetscInt p = pStart; p < pEnd; p++) {
    188:         PetscInt         dof, off;
    189:         PetscReal       *Tp;
    190:         const PetscReal *scalar_Tp;

    192:         PetscCall(PetscSectionGetDof(ref_section, p, &dof));
    193:         PetscCall(PetscSectionGetOffset(ref_section, p, &off));
    194:         scalar_Tp = &scalar_Tk[off * scalar_block];
    195:         Tp        = &Tk[off * scalar_block * v->num_copies * v->num_copies];
    196:         PetscCall(PetscFEVectorInsertTabulation(fe, npoints, points, k, dof, scalar_point_stride, scalar_Tp, Tp));
    197:       }
    198:     }
    199:   }
    200:   PetscCall(PetscTabulationDestroy(&scalar_T));
    201:   PetscFunctionReturn(PETSC_SUCCESS);
    202: }

    204: static PetscErrorCode PetscFECreatePointTrace_Vector(PetscFE fe, PetscInt refPoint, PetscFE *trFE)
    205: {
    206:   PetscFE_Vec *v = (PetscFE_Vec *)fe->data;
    207:   PetscFE      scalar_trFE;
    208:   const char  *name;

    210:   PetscFunctionBegin;
    211:   PetscCall(PetscFECreatePointTrace(v->scalar_fe, refPoint, &scalar_trFE));
    212:   PetscCall(PetscFECreateVector(scalar_trFE, v->num_copies, v->interleave_basis, v->interleave_components, trFE));
    213:   PetscCall(PetscFEDestroy(&scalar_trFE));
    214:   PetscCall(PetscObjectGetName((PetscObject)fe, &name));
    215:   if (name) PetscCall(PetscFESetName(*trFE, name));
    216:   PetscFunctionReturn(PETSC_SUCCESS);
    217: }

    219: PETSC_INTERN PetscErrorCode PetscFEIntegrate_Basic(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
    220: PETSC_INTERN PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS, PetscInt, PetscBdPointFn *, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
    221: PETSC_INTERN PetscErrorCode PetscFEIntegrateHybridResidual_Basic(PetscDS, PetscDS, PetscFormKey, PetscInt, PetscInt, PetscFEGeom *, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
    222: PETSC_INTERN PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS, PetscWeakForm, PetscFEJacobianType, PetscFormKey, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
    223: PETSC_INTERN PetscErrorCode PetscFEIntegrateHybridJacobian_Basic(PetscDS, PetscDS, PetscFEJacobianType, PetscFormKey, PetscInt, PetscInt, PetscFEGeom *, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);

    225: static PetscErrorCode PetscFEInitialize_Vector(PetscFE fe)
    226: {
    227:   PetscFunctionBegin;
    228:   fe->ops->setfromoptions          = NULL;
    229:   fe->ops->setup                   = PetscFESetUp_Vector;
    230:   fe->ops->view                    = PetscFEView_Vector;
    231:   fe->ops->destroy                 = PetscFEDestroy_Vector;
    232:   fe->ops->getdimension            = PetscFEGetDimension_Vector;
    233:   fe->ops->createpointtrace        = PetscFECreatePointTrace_Vector;
    234:   fe->ops->computetabulation       = PetscFEComputeTabulation_Vector;
    235:   fe->ops->integrate               = PetscFEIntegrate_Basic;
    236:   fe->ops->integratebd             = PetscFEIntegrateBd_Basic;
    237:   fe->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
    238:   fe->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
    239:   fe->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
    240:   fe->ops->integratejacobianaction = NULL;
    241:   fe->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
    242:   fe->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
    243:   fe->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
    244:   PetscFunctionReturn(PETSC_SUCCESS);
    245: }

    247: /*MC
    248:   PETSCFEVECTOR = "vector" - A vector-valued `PetscFE` object that is repeated copies
    249:   of the same underlying finite element.

    251:   Level: intermediate

    253: .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`, `PETSCFEBASIC`, `PetscFECreateVector()`
    254: M*/
    255: PETSC_EXTERN PetscErrorCode PetscFECreate_Vector(PetscFE fe)
    256: {
    257:   PetscFE_Vec *v;

    259:   PetscFunctionBegin;
    261:   PetscCall(PetscNew(&v));
    262:   fe->data = v;

    264:   PetscCall(PetscFEInitialize_Vector(fe));
    265:   PetscFunctionReturn(PETSC_SUCCESS);
    266: }

    268: /*@
    269:   PetscFECreateVector - Create a vector-valued `PetscFE` from multiple copies of an underlying
    270:   `PetscFE`.

    272:   Collective

    274:   Input Parameters:
    275: + scalar_fe             - a `PetscFE` finite element
    276: . num_copies            - a positive integer
    277: . interleave_basis      - if `PETSC_TRUE`, the first `num_copies` basis vectors
    278:                           of the output finite element will be copies of the
    279:                           first basis vector of `scalar_fe`, and so on for the
    280:                           other basis vectors; otherwise all of the first-copy
    281:                           basis vectors will come first, followed by all of the
    282:                           second-copy, and so on.
    283: - interleave_components - if `PETSC_TRUE`, the first `num_copies` components
    284:                           of the output finite element will be copies of the
    285:                           first component of `scalar_fe`, and so on for the
    286:                           other components; otherwise all of the first-copy
    287:                           components will come first, followed by all of the
    288:                           second-copy, and so on.

    290:   Output Parameter:
    291: . vector_fe - a `PetscFE` of type `PETSCFEVECTOR` that represent a discretization space with `num_copies` copies of `scalar_fe`

    293:   Level: intermediate

    295: .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`, `PETSCFEBASIC`, `PETSCFEVECTOR`
    296: @*/
    297: PetscErrorCode PetscFECreateVector(PetscFE scalar_fe, PetscInt num_copies, PetscBool interleave_basis, PetscBool interleave_components, PetscFE *vector_fe)
    298: {
    299:   MPI_Comm        comm;
    300:   PetscFE         fe_vec;
    301:   PetscFE_Vec    *v;
    302:   PetscInt        scalar_Nc;
    303:   PetscQuadrature quad;

    305:   PetscFunctionBegin;
    307:   PetscCall(PetscObjectGetComm((PetscObject)scalar_fe, &comm));
    308:   PetscCheck(num_copies >= 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "Expected positive number of copies, got %" PetscInt_FMT, num_copies);
    309:   PetscCall(PetscFECreate(comm, vector_fe));
    310:   fe_vec = *vector_fe;
    311:   PetscCall(PetscFESetType(fe_vec, PETSCFEVECTOR));
    312:   v = (PetscFE_Vec *)fe_vec->data;
    313:   PetscCall(PetscObjectReference((PetscObject)scalar_fe));
    314:   v->scalar_fe             = scalar_fe;
    315:   v->num_copies            = num_copies;
    316:   v->interleave_basis      = interleave_basis;
    317:   v->interleave_components = interleave_components;
    318:   PetscCall(PetscFEGetNumComponents(scalar_fe, &scalar_Nc));
    319:   PetscCall(PetscFESetNumComponents(fe_vec, scalar_Nc * num_copies));
    320:   PetscCall(PetscFEGetQuadrature(scalar_fe, &quad));
    321:   PetscCall(PetscFESetQuadrature(fe_vec, quad));
    322:   PetscCall(PetscFEGetFaceQuadrature(scalar_fe, &quad));
    323:   PetscCall(PetscFESetFaceQuadrature(fe_vec, quad));
    324:   {
    325:     PetscSpace  scalar_sp;
    326:     PetscSpace *copies;
    327:     PetscSpace  sp;

    329:     PetscCall(PetscFEGetBasisSpace(scalar_fe, &scalar_sp));
    330:     PetscCall(PetscMalloc1(num_copies, &copies));
    331:     for (PetscInt i = 0; i < num_copies; i++) {
    332:       PetscCall(PetscObjectReference((PetscObject)scalar_sp));
    333:       copies[i] = scalar_sp;
    334:     }
    335:     PetscCall(PetscSpaceCreateSum(num_copies, copies, PETSC_TRUE, &sp));
    336:     PetscCall(PetscSpaceSumSetInterleave(sp, interleave_basis, interleave_components));
    337:     PetscCall(PetscSpaceSetUp(sp));
    338:     PetscCall(PetscFESetBasisSpace(fe_vec, sp));
    339:     PetscCall(PetscSpaceDestroy(&sp));
    340:     for (PetscInt i = 0; i < num_copies; i++) PetscCall(PetscSpaceDestroy(&copies[i]));
    341:     PetscCall(PetscFree(copies));
    342:   }
    343:   {
    344:     PetscDualSpace  scalar_sp;
    345:     PetscDualSpace *copies;
    346:     PetscDualSpace  sp;

    348:     PetscCall(PetscFEGetDualSpace(scalar_fe, &scalar_sp));
    349:     PetscCall(PetscMalloc1(num_copies, &copies));
    350:     for (PetscInt i = 0; i < num_copies; i++) {
    351:       PetscCall(PetscObjectReference((PetscObject)scalar_sp));
    352:       copies[i] = scalar_sp;
    353:     }
    354:     PetscCall(PetscDualSpaceCreateSum(num_copies, copies, PETSC_TRUE, &sp));
    355:     PetscCall(PetscDualSpaceSumSetInterleave(sp, interleave_basis, interleave_components));
    356:     PetscCall(PetscDualSpaceSetUp(sp));
    357:     PetscCall(PetscFESetDualSpace(fe_vec, sp));
    358:     PetscCall(PetscDualSpaceDestroy(&sp));
    359:     for (PetscInt i = 0; i < num_copies; i++) PetscCall(PetscDualSpaceDestroy(&copies[i]));
    360:     PetscCall(PetscFree(copies));
    361:   }
    362:   PetscFunctionReturn(PETSC_SUCCESS);
    363: }