Actual source code: spacetensor.c

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

      3: static PetscErrorCode PetscSpaceTensorCreateSubspace(PetscSpace space, PetscInt Nvs, PetscInt Ncs, PetscSpace *subspace)
      4: {
      5:   PetscInt    degree;
      6:   const char *prefix;
      7:   const char *name;
      8:   char        subname[PETSC_MAX_PATH_LEN];

     10:   PetscFunctionBegin;
     11:   PetscCall(PetscSpaceGetDegree(space, &degree, NULL));
     12:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)space, &prefix));
     13:   PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)space), subspace));
     14:   PetscCall(PetscSpaceSetType(*subspace, PETSCSPACEPOLYNOMIAL));
     15:   PetscCall(PetscSpaceSetNumVariables(*subspace, Nvs));
     16:   PetscCall(PetscSpaceSetNumComponents(*subspace, Ncs));
     17:   PetscCall(PetscSpaceSetDegree(*subspace, degree, PETSC_DETERMINE));
     18:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*subspace, prefix));
     19:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)*subspace, "tensorcomp_"));
     20:   if (((PetscObject)space)->name) {
     21:     PetscCall(PetscObjectGetName((PetscObject)space, &name));
     22:     PetscCall(PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s tensor component", name));
     23:     PetscCall(PetscObjectSetName((PetscObject)*subspace, subname));
     24:   } else PetscCall(PetscObjectSetName((PetscObject)*subspace, "tensor component"));
     25:   PetscFunctionReturn(PETSC_SUCCESS);
     26: }

     28: static PetscErrorCode PetscSpaceSetFromOptions_Tensor(PetscSpace sp, PetscOptionItems *PetscOptionsObject)
     29: {
     30:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
     31:   PetscInt           Ns, Nc, i, Nv, deg;
     32:   PetscBool          uniform = PETSC_TRUE;

     34:   PetscFunctionBegin;
     35:   PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
     36:   if (!Nv) PetscFunctionReturn(PETSC_SUCCESS);
     37:   PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
     38:   PetscCall(PetscSpaceTensorGetNumSubspaces(sp, &Ns));
     39:   PetscCall(PetscSpaceGetDegree(sp, &deg, NULL));
     40:   if (Ns > 1) {
     41:     PetscSpace s0;

     43:     PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &s0));
     44:     for (i = 1; i < Ns; i++) {
     45:       PetscSpace si;

     47:       PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si));
     48:       if (si != s0) {
     49:         uniform = PETSC_FALSE;
     50:         break;
     51:       }
     52:     }
     53:   }
     54:   Ns = (Ns == PETSC_DEFAULT) ? PetscMax(Nv, 1) : Ns;
     55:   PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace tensor options");
     56:   PetscCall(PetscOptionsBoundedInt("-petscspace_tensor_spaces", "The number of subspaces", "PetscSpaceTensorSetNumSubspaces", Ns, &Ns, NULL, 0));
     57:   PetscCall(PetscOptionsBool("-petscspace_tensor_uniform", "Subspaces are identical", "PetscSpaceTensorSetFromOptions", uniform, &uniform, NULL));
     58:   PetscOptionsHeadEnd();
     59:   PetscCheck(Ns >= 0 && (Nv <= 0 || Ns != 0), PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a tensor space made up of %" PetscInt_FMT " spaces", Ns);
     60:   PetscCheck(Nv <= 0 || Ns <= Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a tensor space with %" PetscInt_FMT " subspaces over %" PetscInt_FMT " variables", Ns, Nv);
     61:   if (Ns != tens->numTensSpaces) PetscCall(PetscSpaceTensorSetNumSubspaces(sp, Ns));
     62:   if (uniform) {
     63:     PetscInt   Nvs = Nv / Ns;
     64:     PetscInt   Ncs;
     65:     PetscSpace subspace;

     67:     PetscCheck(Nv % Ns == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " variable space", Ns, Nv);
     68:     Ncs = (PetscInt)PetscPowReal((PetscReal)Nc, 1. / Ns);
     69:     PetscCheck(Nc % PetscPowInt(Ncs, Ns) == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " component space", Ns, Nc);
     70:     PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &subspace));
     71:     if (!subspace) PetscCall(PetscSpaceTensorCreateSubspace(sp, Nvs, Ncs, &subspace));
     72:     else PetscCall(PetscObjectReference((PetscObject)subspace));
     73:     PetscCall(PetscSpaceSetFromOptions(subspace));
     74:     for (i = 0; i < Ns; i++) PetscCall(PetscSpaceTensorSetSubspace(sp, i, subspace));
     75:     PetscCall(PetscSpaceDestroy(&subspace));
     76:   } else {
     77:     for (i = 0; i < Ns; i++) {
     78:       PetscSpace subspace;

     80:       PetscCall(PetscSpaceTensorGetSubspace(sp, i, &subspace));
     81:       if (!subspace) {
     82:         char tprefix[128];

     84:         PetscCall(PetscSpaceTensorCreateSubspace(sp, 1, 1, &subspace));
     85:         PetscCall(PetscSNPrintf(tprefix, 128, "%d_", (int)i));
     86:         PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subspace, tprefix));
     87:       } else PetscCall(PetscObjectReference((PetscObject)subspace));
     88:       PetscCall(PetscSpaceSetFromOptions(subspace));
     89:       PetscCall(PetscSpaceTensorSetSubspace(sp, i, subspace));
     90:       PetscCall(PetscSpaceDestroy(&subspace));
     91:     }
     92:   }
     93:   PetscFunctionReturn(PETSC_SUCCESS);
     94: }

     96: static PetscErrorCode PetscSpaceTensorView_Ascii(PetscSpace sp, PetscViewer v)
     97: {
     98:   PetscSpace_Tensor *tens    = (PetscSpace_Tensor *)sp->data;
     99:   PetscBool          uniform = PETSC_TRUE;
    100:   PetscInt           Ns      = tens->numTensSpaces, i, n;

    102:   PetscFunctionBegin;
    103:   for (i = 1; i < Ns; i++) {
    104:     if (tens->tensspaces[i] != tens->tensspaces[0]) {
    105:       uniform = PETSC_FALSE;
    106:       break;
    107:     }
    108:   }
    109:   if (uniform) PetscCall(PetscViewerASCIIPrintf(v, "Tensor space of %" PetscInt_FMT " subspaces (all identical)\n", Ns));
    110:   else PetscCall(PetscViewerASCIIPrintf(v, "Tensor space of %" PetscInt_FMT " subspaces\n", Ns));
    111:   n = uniform ? 1 : Ns;
    112:   for (i = 0; i < n; i++) {
    113:     PetscCall(PetscViewerASCIIPushTab(v));
    114:     PetscCall(PetscSpaceView(tens->tensspaces[i], v));
    115:     PetscCall(PetscViewerASCIIPopTab(v));
    116:   }
    117:   PetscFunctionReturn(PETSC_SUCCESS);
    118: }

    120: static PetscErrorCode PetscSpaceView_Tensor(PetscSpace sp, PetscViewer viewer)
    121: {
    122:   PetscBool iascii;

    124:   PetscFunctionBegin;
    125:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    126:   if (iascii) PetscCall(PetscSpaceTensorView_Ascii(sp, viewer));
    127:   PetscFunctionReturn(PETSC_SUCCESS);
    128: }

    130: static PetscErrorCode PetscSpaceSetUp_Tensor(PetscSpace sp)
    131: {
    132:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
    133:   PetscInt           Nc, Nv, Ns;
    134:   PetscBool          uniform = PETSC_TRUE;
    135:   PetscInt           deg, maxDeg;
    136:   PetscInt           Ncprod;

    138:   PetscFunctionBegin;
    139:   if (tens->setupCalled) PetscFunctionReturn(PETSC_SUCCESS);
    140:   PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
    141:   PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
    142:   PetscCall(PetscSpaceTensorGetNumSubspaces(sp, &Ns));
    143:   if (Ns == PETSC_DEFAULT) {
    144:     Ns = Nv;
    145:     PetscCall(PetscSpaceTensorSetNumSubspaces(sp, Ns));
    146:   }
    147:   if (!Ns) {
    148:     SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zero subspaces");
    149:   } else {
    150:     PetscSpace s0 = NULL;

    152:     PetscCheck(Nv <= 0 || Ns <= Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a tensor space with %" PetscInt_FMT " subspaces over %" PetscInt_FMT " variables", Ns, Nv);
    153:     PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &s0));
    154:     for (PetscInt i = 1; i < Ns; i++) {
    155:       PetscSpace si;

    157:       PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si));
    158:       if (si != s0) {
    159:         uniform = PETSC_FALSE;
    160:         break;
    161:       }
    162:     }
    163:     if (uniform) {
    164:       PetscInt Nvs = Nv / Ns;
    165:       PetscInt Ncs;

    167:       PetscCheck(Nv % Ns == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " variable space", Ns, Nv);
    168:       Ncs = (PetscInt)(PetscPowReal((PetscReal)Nc, 1. / Ns));
    169:       PetscCheck(Nc % PetscPowInt(Ncs, Ns) == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " component space", Ns, Nc);
    170:       if (!s0) PetscCall(PetscSpaceTensorCreateSubspace(sp, Nvs, Ncs, &s0));
    171:       else PetscCall(PetscObjectReference((PetscObject)s0));
    172:       PetscCall(PetscSpaceSetUp(s0));
    173:       for (PetscInt i = 0; i < Ns; i++) PetscCall(PetscSpaceTensorSetSubspace(sp, i, s0));
    174:       PetscCall(PetscSpaceDestroy(&s0));
    175:       Ncprod = PetscPowInt(Ncs, Ns);
    176:     } else {
    177:       PetscInt Nvsum = 0;

    179:       Ncprod = 1;
    180:       for (PetscInt i = 0; i < Ns; i++) {
    181:         PetscInt   Nvs, Ncs;
    182:         PetscSpace si = NULL;

    184:         PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si));
    185:         if (!si) PetscCall(PetscSpaceTensorCreateSubspace(sp, 1, 1, &si));
    186:         else PetscCall(PetscObjectReference((PetscObject)si));
    187:         PetscCall(PetscSpaceSetUp(si));
    188:         PetscCall(PetscSpaceTensorSetSubspace(sp, i, si));
    189:         PetscCall(PetscSpaceGetNumVariables(si, &Nvs));
    190:         PetscCall(PetscSpaceGetNumComponents(si, &Ncs));
    191:         PetscCall(PetscSpaceDestroy(&si));

    193:         Nvsum += Nvs;
    194:         Ncprod *= Ncs;
    195:       }

    197:       PetscCheck(Nvsum == Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Sum of subspace variables %" PetscInt_FMT " does not equal the number of variables %" PetscInt_FMT, Nvsum, Nv);
    198:       PetscCheck(Nc % Ncprod == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Product of subspace components %" PetscInt_FMT " does not divide the number of components %" PetscInt_FMT, Ncprod, Nc);
    199:     }
    200:     if (Ncprod != Nc) {
    201:       PetscInt    Ncopies = Nc / Ncprod;
    202:       PetscInt    Nv      = sp->Nv;
    203:       const char *prefix;
    204:       const char *name;
    205:       char        subname[PETSC_MAX_PATH_LEN];
    206:       PetscSpace  subsp;

    208:       PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp));
    209:       PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix));
    210:       PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix));
    211:       PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_"));
    212:       if (((PetscObject)sp)->name) {
    213:         PetscCall(PetscObjectGetName((PetscObject)sp, &name));
    214:         PetscCall(PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s sum component", name));
    215:         PetscCall(PetscObjectSetName((PetscObject)subsp, subname));
    216:       } else PetscCall(PetscObjectSetName((PetscObject)subsp, "sum component"));
    217:       PetscCall(PetscSpaceSetType(subsp, PETSCSPACETENSOR));
    218:       PetscCall(PetscSpaceSetNumVariables(subsp, Nv));
    219:       PetscCall(PetscSpaceSetNumComponents(subsp, Ncprod));
    220:       PetscCall(PetscSpaceTensorSetNumSubspaces(subsp, Ns));
    221:       for (PetscInt i = 0; i < Ns; i++) {
    222:         PetscSpace ssp;

    224:         PetscCall(PetscSpaceTensorGetSubspace(sp, i, &ssp));
    225:         PetscCall(PetscSpaceTensorSetSubspace(subsp, i, ssp));
    226:       }
    227:       PetscCall(PetscSpaceSetUp(subsp));
    228:       PetscCall(PetscSpaceSetType(sp, PETSCSPACESUM));
    229:       PetscCall(PetscSpaceSumSetNumSubspaces(sp, Ncopies));
    230:       for (PetscInt i = 0; i < Ncopies; i++) PetscCall(PetscSpaceSumSetSubspace(sp, i, subsp));
    231:       PetscCall(PetscSpaceDestroy(&subsp));
    232:       PetscCall(PetscSpaceSetUp(sp));
    233:       PetscFunctionReturn(PETSC_SUCCESS);
    234:     }
    235:   }
    236:   deg    = PETSC_INT_MAX;
    237:   maxDeg = 0;
    238:   for (PetscInt i = 0; i < Ns; i++) {
    239:     PetscSpace si;
    240:     PetscInt   iDeg, iMaxDeg;

    242:     PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si));
    243:     PetscCall(PetscSpaceGetDegree(si, &iDeg, &iMaxDeg));
    244:     deg = PetscMin(deg, iDeg);
    245:     maxDeg += iMaxDeg;
    246:   }
    247:   sp->degree        = deg;
    248:   sp->maxDegree     = maxDeg;
    249:   tens->uniform     = uniform;
    250:   tens->setupCalled = PETSC_TRUE;
    251:   PetscFunctionReturn(PETSC_SUCCESS);
    252: }

    254: static PetscErrorCode PetscSpaceDestroy_Tensor(PetscSpace sp)
    255: {
    256:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
    257:   PetscInt           Ns, i;

    259:   PetscFunctionBegin;
    260:   Ns = tens->numTensSpaces;
    261:   if (tens->heightsubspaces) {
    262:     PetscInt d;

    264:     /* sp->Nv is the spatial dimension, so it is equal to the number
    265:      * of subspaces on higher co-dimension points */
    266:     for (d = 0; d < sp->Nv; ++d) PetscCall(PetscSpaceDestroy(&tens->heightsubspaces[d]));
    267:   }
    268:   PetscCall(PetscFree(tens->heightsubspaces));
    269:   for (i = 0; i < Ns; i++) PetscCall(PetscSpaceDestroy(&tens->tensspaces[i]));
    270:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", NULL));
    271:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", NULL));
    272:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", NULL));
    273:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", NULL));
    274:   PetscCall(PetscFree(tens->tensspaces));
    275:   PetscCall(PetscFree(tens));
    276:   PetscFunctionReturn(PETSC_SUCCESS);
    277: }

    279: static PetscErrorCode PetscSpaceGetDimension_Tensor(PetscSpace sp, PetscInt *dim)
    280: {
    281:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
    282:   PetscInt           i, Ns, d;

    284:   PetscFunctionBegin;
    285:   PetscCall(PetscSpaceSetUp(sp));
    286:   Ns = tens->numTensSpaces;
    287:   d  = 1;
    288:   for (i = 0; i < Ns; i++) {
    289:     PetscInt id;

    291:     PetscCall(PetscSpaceGetDimension(tens->tensspaces[i], &id));
    292:     d *= id;
    293:   }
    294:   *dim = d;
    295:   PetscFunctionReturn(PETSC_SUCCESS);
    296: }

    298: static PetscErrorCode PetscSpaceEvaluate_Tensor(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
    299: {
    300:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
    301:   DM                 dm   = sp->dm;
    302:   PetscInt           Nc   = sp->Nc;
    303:   PetscInt           Nv   = sp->Nv;
    304:   PetscInt           Ns;
    305:   PetscReal         *lpoints, *sB = NULL, *sD = NULL, *sH = NULL;
    306:   PetscInt           pdim;

    308:   PetscFunctionBegin;
    309:   if (!tens->setupCalled) {
    310:     PetscCall(PetscSpaceSetUp(sp));
    311:     PetscCall(PetscSpaceEvaluate(sp, npoints, points, B, D, H));
    312:     PetscFunctionReturn(PETSC_SUCCESS);
    313:   }
    314:   Ns = tens->numTensSpaces;
    315:   PetscCall(PetscSpaceGetDimension(sp, &pdim));
    316:   PetscCall(DMGetWorkArray(dm, npoints * Nv, MPIU_REAL, &lpoints));
    317:   if (B || D || H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &sB));
    318:   if (D || H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * Nv, MPIU_REAL, &sD));
    319:   if (H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * Nv * Nv, MPIU_REAL, &sH));
    320:   if (B) {
    321:     for (PetscInt i = 0; i < npoints * pdim * Nc; i++) B[i] = 1.;
    322:   }
    323:   if (D) {
    324:     for (PetscInt i = 0; i < npoints * pdim * Nc * Nv; i++) D[i] = 1.;
    325:   }
    326:   if (H) {
    327:     for (PetscInt i = 0; i < npoints * pdim * Nc * Nv * Nv; i++) H[i] = 1.;
    328:   }
    329:   for (PetscInt s = 0, d = 0, vstep = 1, cstep = 1; s < Ns; s++) {
    330:     PetscInt sNv, sNc, spdim;
    331:     PetscInt vskip, cskip;

    333:     PetscCall(PetscSpaceGetNumVariables(tens->tensspaces[s], &sNv));
    334:     PetscCall(PetscSpaceGetNumComponents(tens->tensspaces[s], &sNc));
    335:     PetscCall(PetscSpaceGetDimension(tens->tensspaces[s], &spdim));
    336:     PetscCheck(!(pdim % vstep) && !(pdim % spdim), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad tensor loop: Nv %" PetscInt_FMT ", Ns %" PetscInt_FMT ", pdim %" PetscInt_FMT ", s %" PetscInt_FMT ", vstep %" PetscInt_FMT ", spdim %" PetscInt_FMT, Nv, Ns, pdim, s, vstep, spdim);
    337:     PetscCheck(!(Nc % cstep) && !(Nc % sNc), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad tensor loop: Nv %" PetscInt_FMT ", Ns %" PetscInt_FMT ", Nc %" PetscInt_FMT ", s %" PetscInt_FMT ", cstep %" PetscInt_FMT ", sNc %" PetscInt_FMT, Nv, Ns, Nc, s, cstep, spdim);
    338:     vskip = pdim / (vstep * spdim);
    339:     cskip = Nc / (cstep * sNc);
    340:     for (PetscInt p = 0; p < npoints; ++p) {
    341:       for (PetscInt i = 0; i < sNv; i++) lpoints[p * sNv + i] = points[p * Nv + d + i];
    342:     }
    343:     PetscCall(PetscSpaceEvaluate(tens->tensspaces[s], npoints, lpoints, sB, sD, sH));
    344:     if (B) {
    345:       for (PetscInt k = 0; k < vskip; k++) {
    346:         for (PetscInt si = 0; si < spdim; si++) {
    347:           for (PetscInt j = 0; j < vstep; j++) {
    348:             PetscInt i = (k * spdim + si) * vstep + j;

    350:             for (PetscInt l = 0; l < cskip; l++) {
    351:               for (PetscInt sc = 0; sc < sNc; sc++) {
    352:                 for (PetscInt m = 0; m < cstep; m++) {
    353:                   PetscInt c = (l * sNc + sc) * cstep + m;

    355:                   for (PetscInt p = 0; p < npoints; p++) B[(pdim * p + i) * Nc + c] *= sB[(spdim * p + si) * sNc + sc];
    356:                 }
    357:               }
    358:             }
    359:           }
    360:         }
    361:       }
    362:     }
    363:     if (D) {
    364:       for (PetscInt k = 0; k < vskip; k++) {
    365:         for (PetscInt si = 0; si < spdim; si++) {
    366:           for (PetscInt j = 0; j < vstep; j++) {
    367:             PetscInt i = (k * spdim + si) * vstep + j;

    369:             for (PetscInt l = 0; l < cskip; l++) {
    370:               for (PetscInt sc = 0; sc < sNc; sc++) {
    371:                 for (PetscInt m = 0; m < cstep; m++) {
    372:                   PetscInt c = (l * sNc + sc) * cstep + m;

    374:                   for (PetscInt der = 0; der < Nv; der++) {
    375:                     if (der >= d && der < d + sNv) {
    376:                       for (PetscInt p = 0; p < npoints; p++) D[((pdim * p + i) * Nc + c) * Nv + der] *= sD[((spdim * p + si) * sNc + sc) * sNv + der - d];
    377:                     } else {
    378:                       for (PetscInt p = 0; p < npoints; p++) D[((pdim * p + i) * Nc + c) * Nv + der] *= sB[(spdim * p + si) * sNc + sc];
    379:                     }
    380:                   }
    381:                 }
    382:               }
    383:             }
    384:           }
    385:         }
    386:       }
    387:     }
    388:     if (H) {
    389:       for (PetscInt k = 0; k < vskip; k++) {
    390:         for (PetscInt si = 0; si < spdim; si++) {
    391:           for (PetscInt j = 0; j < vstep; j++) {
    392:             PetscInt i = (k * spdim + si) * vstep + j;

    394:             for (PetscInt l = 0; l < cskip; l++) {
    395:               for (PetscInt sc = 0; sc < sNc; sc++) {
    396:                 for (PetscInt m = 0; m < cstep; m++) {
    397:                   PetscInt c = (l * sNc + sc) * cstep + m;

    399:                   for (PetscInt der = 0; der < Nv; der++) {
    400:                     for (PetscInt der2 = 0; der2 < Nv; der2++) {
    401:                       if (der >= d && der < d + sNv && der2 >= d && der2 < d + sNv) {
    402:                         for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sH[(((spdim * p + si) * sNc + sc) * sNv + der - d) * sNv + der2 - d];
    403:                       } else if (der >= d && der < d + sNv) {
    404:                         for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sD[((spdim * p + si) * sNc + sc) * sNv + der - d];
    405:                       } else if (der2 >= d && der2 < d + sNv) {
    406:                         for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sD[((spdim * p + si) * sNc + sc) * sNv + der2 - d];
    407:                       } else {
    408:                         for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sB[(spdim * p + si) * sNc + sc];
    409:                       }
    410:                     }
    411:                   }
    412:                 }
    413:               }
    414:             }
    415:           }
    416:         }
    417:       }
    418:     }
    419:     d += sNv;
    420:     vstep *= spdim;
    421:     cstep *= sNc;
    422:   }
    423:   if (H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * Nv * Nv, MPIU_REAL, &sH));
    424:   if (D || H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * Nv, MPIU_REAL, &sD));
    425:   if (B || D || H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &sB));
    426:   PetscCall(DMRestoreWorkArray(dm, npoints * Nv, MPIU_REAL, &lpoints));
    427:   PetscFunctionReturn(PETSC_SUCCESS);
    428: }

    430: /*@
    431:   PetscSpaceTensorSetNumSubspaces - Set the number of spaces in the tensor product space

    433:   Input Parameters:
    434: + sp            - the function space object
    435: - numTensSpaces - the number of spaces

    437:   Level: intermediate

    439:   Note:
    440:   The name NumSubspaces is misleading because it is actually setting the number of defining spaces of the tensor product space, not a number of Subspaces of it

    442: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
    443: @*/
    444: PetscErrorCode PetscSpaceTensorSetNumSubspaces(PetscSpace sp, PetscInt numTensSpaces)
    445: {
    446:   PetscFunctionBegin;
    448:   PetscTryMethod(sp, "PetscSpaceTensorSetNumSubspaces_C", (PetscSpace, PetscInt), (sp, numTensSpaces));
    449:   PetscFunctionReturn(PETSC_SUCCESS);
    450: }

    452: /*@
    453:   PetscSpaceTensorGetNumSubspaces - Get the number of spaces in the tensor product space

    455:   Input Parameter:
    456: . sp - the function space object

    458:   Output Parameter:
    459: . numTensSpaces - the number of spaces

    461:   Level: intermediate

    463:   Note:
    464:   The name NumSubspaces is misleading because it is actually getting the number of defining spaces of the tensor product space, not a number of Subspaces of it

    466: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
    467: @*/
    468: PetscErrorCode PetscSpaceTensorGetNumSubspaces(PetscSpace sp, PetscInt *numTensSpaces)
    469: {
    470:   PetscFunctionBegin;
    472:   PetscAssertPointer(numTensSpaces, 2);
    473:   PetscTryMethod(sp, "PetscSpaceTensorGetNumSubspaces_C", (PetscSpace, PetscInt *), (sp, numTensSpaces));
    474:   PetscFunctionReturn(PETSC_SUCCESS);
    475: }

    477: /*@
    478:   PetscSpaceTensorSetSubspace - Set a space in the tensor product space

    480:   Input Parameters:
    481: + sp    - the function space object
    482: . s     - The space number
    483: - subsp - the number of spaces

    485:   Level: intermediate

    487:   Note:
    488:   The name SetSubspace is misleading because it is actually setting one of the defining spaces of the tensor product space, not a Subspace of it

    490: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
    491: @*/
    492: PetscErrorCode PetscSpaceTensorSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
    493: {
    494:   PetscFunctionBegin;
    497:   PetscTryMethod(sp, "PetscSpaceTensorSetSubspace_C", (PetscSpace, PetscInt, PetscSpace), (sp, s, subsp));
    498:   PetscFunctionReturn(PETSC_SUCCESS);
    499: }

    501: /*@
    502:   PetscSpaceTensorGetSubspace - Get a space in the tensor product space

    504:   Input Parameters:
    505: + sp - the function space object
    506: - s  - The space number

    508:   Output Parameter:
    509: . subsp - the `PetscSpace`

    511:   Level: intermediate

    513:   Note:
    514:   The name GetSubspace is misleading because it is actually getting one of the defining spaces of the tensor product space, not a Subspace of it

    516: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
    517: @*/
    518: PetscErrorCode PetscSpaceTensorGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
    519: {
    520:   PetscFunctionBegin;
    522:   PetscAssertPointer(subsp, 3);
    523:   PetscTryMethod(sp, "PetscSpaceTensorGetSubspace_C", (PetscSpace, PetscInt, PetscSpace *), (sp, s, subsp));
    524:   PetscFunctionReturn(PETSC_SUCCESS);
    525: }

    527: static PetscErrorCode PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space, PetscInt numTensSpaces)
    528: {
    529:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
    530:   PetscInt           Ns;

    532:   PetscFunctionBegin;
    533:   PetscCheck(!tens->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of subspaces after setup called");
    534:   Ns = tens->numTensSpaces;
    535:   if (numTensSpaces == Ns) PetscFunctionReturn(PETSC_SUCCESS);
    536:   if (Ns >= 0) {
    537:     PetscInt s;

    539:     for (s = 0; s < Ns; s++) PetscCall(PetscSpaceDestroy(&tens->tensspaces[s]));
    540:     PetscCall(PetscFree(tens->tensspaces));
    541:   }
    542:   Ns = tens->numTensSpaces = numTensSpaces;
    543:   PetscCall(PetscCalloc1(Ns, &tens->tensspaces));
    544:   PetscFunctionReturn(PETSC_SUCCESS);
    545: }

    547: static PetscErrorCode PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space, PetscInt *numTensSpaces)
    548: {
    549:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;

    551:   PetscFunctionBegin;
    552:   *numTensSpaces = tens->numTensSpaces;
    553:   PetscFunctionReturn(PETSC_SUCCESS);
    554: }

    556: static PetscErrorCode PetscSpaceTensorSetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace subspace)
    557: {
    558:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
    559:   PetscInt           Ns;

    561:   PetscFunctionBegin;
    562:   PetscCheck(!tens->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change subspace after setup called");
    563:   Ns = tens->numTensSpaces;
    564:   PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceTensorSetNumSubspaces() first");
    565:   PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
    566:   PetscCall(PetscObjectReference((PetscObject)subspace));
    567:   PetscCall(PetscSpaceDestroy(&tens->tensspaces[s]));
    568:   tens->tensspaces[s] = subspace;
    569:   PetscFunctionReturn(PETSC_SUCCESS);
    570: }

    572: static PetscErrorCode PetscSpaceGetHeightSubspace_Tensor(PetscSpace sp, PetscInt height, PetscSpace *subsp)
    573: {
    574:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
    575:   PetscInt           Nc, dim, order, i;
    576:   PetscSpace         bsp;

    578:   PetscFunctionBegin;
    579:   PetscCall(PetscSpaceGetNumVariables(sp, &dim));
    580:   PetscCheck(tens->uniform && tens->numTensSpaces == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can only get a generic height subspace of a uniform tensor space of 1d spaces.");
    581:   PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
    582:   PetscCall(PetscSpaceGetDegree(sp, &order, NULL));
    583:   PetscCheck(height <= dim && height >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %" PetscInt_FMT " for dimension %" PetscInt_FMT " space", height, dim);
    584:   if (!tens->heightsubspaces) PetscCall(PetscCalloc1(dim, &tens->heightsubspaces));
    585:   if (height <= dim) {
    586:     if (!tens->heightsubspaces[height - 1]) {
    587:       PetscSpace  sub;
    588:       const char *name;

    590:       PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &bsp));
    591:       PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub));
    592:       PetscCall(PetscObjectGetName((PetscObject)sp, &name));
    593:       PetscCall(PetscObjectSetName((PetscObject)sub, name));
    594:       PetscCall(PetscSpaceSetType(sub, PETSCSPACETENSOR));
    595:       PetscCall(PetscSpaceSetNumComponents(sub, Nc));
    596:       PetscCall(PetscSpaceSetDegree(sub, order, PETSC_DETERMINE));
    597:       PetscCall(PetscSpaceSetNumVariables(sub, dim - height));
    598:       PetscCall(PetscSpaceTensorSetNumSubspaces(sub, dim - height));
    599:       for (i = 0; i < dim - height; i++) PetscCall(PetscSpaceTensorSetSubspace(sub, i, bsp));
    600:       PetscCall(PetscSpaceSetUp(sub));
    601:       tens->heightsubspaces[height - 1] = sub;
    602:     }
    603:     *subsp = tens->heightsubspaces[height - 1];
    604:   } else {
    605:     *subsp = NULL;
    606:   }
    607:   PetscFunctionReturn(PETSC_SUCCESS);
    608: }

    610: static PetscErrorCode PetscSpaceTensorGetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace *subspace)
    611: {
    612:   PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
    613:   PetscInt           Ns;

    615:   PetscFunctionBegin;
    616:   Ns = tens->numTensSpaces;
    617:   PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceTensorSetNumSubspaces() first");
    618:   PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
    619:   *subspace = tens->tensspaces[s];
    620:   PetscFunctionReturn(PETSC_SUCCESS);
    621: }

    623: static PetscErrorCode PetscSpaceInitialize_Tensor(PetscSpace sp)
    624: {
    625:   PetscFunctionBegin;
    626:   sp->ops->setfromoptions    = PetscSpaceSetFromOptions_Tensor;
    627:   sp->ops->setup             = PetscSpaceSetUp_Tensor;
    628:   sp->ops->view              = PetscSpaceView_Tensor;
    629:   sp->ops->destroy           = PetscSpaceDestroy_Tensor;
    630:   sp->ops->getdimension      = PetscSpaceGetDimension_Tensor;
    631:   sp->ops->evaluate          = PetscSpaceEvaluate_Tensor;
    632:   sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Tensor;
    633:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", PetscSpaceTensorGetNumSubspaces_Tensor));
    634:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", PetscSpaceTensorSetNumSubspaces_Tensor));
    635:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", PetscSpaceTensorGetSubspace_Tensor));
    636:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", PetscSpaceTensorSetSubspace_Tensor));
    637:   PetscFunctionReturn(PETSC_SUCCESS);
    638: }

    640: /*MC
    641:   PETSCSPACETENSOR = "tensor" - A `PetscSpace` object that encapsulates a tensor product space.
    642:                      A tensor product is created of the components of the subspaces as well.

    644:   Level: intermediate

    646: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`, `PetscSpaceTensorGetSubspace()`, `PetscSpaceTensorSetSubspace()`,
    647:           `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceTensorSetNumSubspaces()`
    648: M*/

    650: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Tensor(PetscSpace sp)
    651: {
    652:   PetscSpace_Tensor *tens;

    654:   PetscFunctionBegin;
    656:   PetscCall(PetscNew(&tens));
    657:   sp->data = tens;

    659:   tens->numTensSpaces = PETSC_DEFAULT;

    661:   PetscCall(PetscSpaceInitialize_Tensor(sp));
    662:   PetscFunctionReturn(PETSC_SUCCESS);
    663: }