Actual source code: spacesubspace.c

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

      3: typedef struct {
      4:   PetscDualSpace dualSubspace;
      5:   PetscSpace     origSpace;
      6:   PetscReal     *x;
      7:   PetscReal     *x_alloc;
      8:   PetscReal     *Jx;
      9:   PetscReal     *Jx_alloc;
     10:   PetscReal     *u;
     11:   PetscReal     *u_alloc;
     12:   PetscReal     *Ju;
     13:   PetscReal     *Ju_alloc;
     14:   PetscReal     *Q;
     15:   PetscInt       Nb;
     16:   PetscBool      setupcalled;
     17: } PetscSpace_Subspace;

     19: static PetscErrorCode PetscSpaceDestroy_Subspace(PetscSpace sp)
     20: {
     21:   PetscSpace_Subspace *subsp;

     23:   PetscFunctionBegin;
     24:   subsp    = (PetscSpace_Subspace *)sp->data;
     25:   subsp->x = NULL;
     26:   PetscCall(PetscFree(subsp->x_alloc));
     27:   subsp->Jx = NULL;
     28:   PetscCall(PetscFree(subsp->Jx_alloc));
     29:   subsp->u = NULL;
     30:   PetscCall(PetscFree(subsp->u_alloc));
     31:   subsp->Ju = NULL;
     32:   PetscCall(PetscFree(subsp->Ju_alloc));
     33:   PetscCall(PetscFree(subsp->Q));
     34:   PetscCall(PetscSpaceDestroy(&subsp->origSpace));
     35:   PetscCall(PetscDualSpaceDestroy(&subsp->dualSubspace));
     36:   PetscCall(PetscFree(subsp));
     37:   sp->data = NULL;
     38:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialGetTensor_C", NULL));
     39:   PetscFunctionReturn(PETSC_SUCCESS);
     40: }

     42: static PetscErrorCode PetscSpaceView_Subspace(PetscSpace sp, PetscViewer viewer)
     43: {
     44:   PetscBool            iascii;
     45:   PetscSpace_Subspace *subsp;

     47:   PetscFunctionBegin;
     48:   subsp = (PetscSpace_Subspace *)sp->data;
     49:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
     50:   if (iascii) {
     51:     PetscInt origDim, subDim, origNc, subNc, o, s;

     53:     PetscCall(PetscSpaceGetNumVariables(subsp->origSpace, &origDim));
     54:     PetscCall(PetscSpaceGetNumComponents(subsp->origSpace, &origNc));
     55:     PetscCall(PetscSpaceGetNumVariables(sp, &subDim));
     56:     PetscCall(PetscSpaceGetNumComponents(sp, &subNc));
     57:     if (subsp->x) {
     58:       PetscCall(PetscViewerASCIIPrintf(viewer, "Subspace-to-space domain shift:\n\n"));
     59:       for (o = 0; o < origDim; o++) PetscCall(PetscViewerASCIIPrintf(viewer, " %g\n", (double)subsp->x[o]));
     60:     }
     61:     if (subsp->Jx) {
     62:       PetscCall(PetscViewerASCIIPrintf(viewer, "Subspace-to-space domain transform:\n\n"));
     63:       for (o = 0; o < origDim; o++) {
     64:         PetscCall(PetscViewerASCIIPrintf(viewer, " %g", (double)subsp->Jx[o * subDim + 0]));
     65:         PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
     66:         for (s = 1; s < subDim; s++) PetscCall(PetscViewerASCIIPrintf(viewer, " %g", (double)subsp->Jx[o * subDim + s]));
     67:         PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
     68:         PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
     69:       }
     70:     }
     71:     if (subsp->u) {
     72:       PetscCall(PetscViewerASCIIPrintf(viewer, "Space-to-subspace range shift:\n\n"));
     73:       for (o = 0; o < origNc; o++) PetscCall(PetscViewerASCIIPrintf(viewer, " %g\n", (double)subsp->u[o]));
     74:     }
     75:     if (subsp->Ju) {
     76:       PetscCall(PetscViewerASCIIPrintf(viewer, "Space-to-subsace domain transform:\n"));
     77:       for (o = 0; o < origNc; o++) {
     78:         PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
     79:         for (s = 0; s < subNc; s++) PetscCall(PetscViewerASCIIPrintf(viewer, " %g", (double)subsp->Ju[o * subNc + s]));
     80:         PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
     81:       }
     82:       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
     83:     }
     84:     PetscCall(PetscViewerASCIIPrintf(viewer, "Original space:\n"));
     85:   }
     86:   PetscCall(PetscViewerASCIIPushTab(viewer));
     87:   PetscCall(PetscSpaceView(subsp->origSpace, viewer));
     88:   PetscCall(PetscViewerASCIIPopTab(viewer));
     89:   PetscFunctionReturn(PETSC_SUCCESS);
     90: }

     92: static PetscErrorCode PetscSpaceEvaluate_Subspace(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
     93: {
     94:   PetscSpace_Subspace *subsp = (PetscSpace_Subspace *)sp->data;
     95:   PetscSpace           origsp;
     96:   PetscInt             origDim, subDim, origNc, subNc, subNb, origNb, i, j, k, l, m, n, o;
     97:   PetscReal           *inpoints, *inB = NULL, *inD = NULL, *inH = NULL;

     99:   PetscFunctionBegin;
    100:   origsp = subsp->origSpace;
    101:   PetscCall(PetscSpaceGetNumVariables(sp, &subDim));
    102:   PetscCall(PetscSpaceGetNumVariables(origsp, &origDim));
    103:   PetscCall(PetscSpaceGetNumComponents(sp, &subNc));
    104:   PetscCall(PetscSpaceGetNumComponents(origsp, &origNc));
    105:   PetscCall(PetscSpaceGetDimension(sp, &subNb));
    106:   PetscCall(PetscSpaceGetDimension(origsp, &origNb));
    107:   PetscCall(DMGetWorkArray(sp->dm, npoints * origDim, MPIU_REAL, &inpoints));
    108:   for (i = 0; i < npoints; i++) {
    109:     if (subsp->x) {
    110:       for (j = 0; j < origDim; j++) inpoints[i * origDim + j] = subsp->x[j];
    111:     } else {
    112:       for (j = 0; j < origDim; j++) inpoints[i * origDim + j] = 0.0;
    113:     }
    114:     if (subsp->Jx) {
    115:       for (j = 0; j < origDim; j++) {
    116:         for (k = 0; k < subDim; k++) inpoints[i * origDim + j] += subsp->Jx[j * subDim + k] * points[i * subDim + k];
    117:       }
    118:     } else {
    119:       for (j = 0; j < PetscMin(subDim, origDim); j++) inpoints[i * origDim + j] += points[i * subDim + j];
    120:     }
    121:   }
    122:   if (B) PetscCall(DMGetWorkArray(sp->dm, npoints * origNb * origNc, MPIU_REAL, &inB));
    123:   if (D) PetscCall(DMGetWorkArray(sp->dm, npoints * origNb * origNc * origDim, MPIU_REAL, &inD));
    124:   if (H) PetscCall(DMGetWorkArray(sp->dm, npoints * origNb * origNc * origDim * origDim, MPIU_REAL, &inH));
    125:   PetscCall(PetscSpaceEvaluate(origsp, npoints, inpoints, inB, inD, inH));
    126:   if (H) {
    127:     PetscReal *phi, *psi;

    129:     PetscCall(DMGetWorkArray(sp->dm, origNc * origDim * origDim, MPIU_REAL, &phi));
    130:     PetscCall(DMGetWorkArray(sp->dm, origNc * subDim * subDim, MPIU_REAL, &psi));
    131:     for (i = 0; i < npoints * subNb * subNc * subDim; i++) D[i] = 0.0;
    132:     for (i = 0; i < subNb; i++) {
    133:       const PetscReal *subq = &subsp->Q[i * origNb];

    135:       for (j = 0; j < npoints; j++) {
    136:         for (k = 0; k < origNc * origDim; k++) phi[k] = 0.;
    137:         for (k = 0; k < origNc * subDim; k++) psi[k] = 0.;
    138:         for (k = 0; k < origNb; k++) {
    139:           for (l = 0; l < origNc * origDim * origDim; l++) phi[l] += inH[(j * origNb + k) * origNc * origDim * origDim + l] * subq[k];
    140:         }
    141:         if (subsp->Jx) {
    142:           for (k = 0; k < subNc; k++) {
    143:             for (l = 0; l < subDim; l++) {
    144:               for (m = 0; m < origDim; m++) {
    145:                 for (n = 0; n < subDim; n++) {
    146:                   for (o = 0; o < origDim; o++) psi[(k * subDim + l) * subDim + n] += subsp->Jx[m * subDim + l] * subsp->Jx[o * subDim + n] * phi[(k * origDim + m) * origDim + o];
    147:                 }
    148:               }
    149:             }
    150:           }
    151:         } else {
    152:           for (k = 0; k < subNc; k++) {
    153:             for (l = 0; l < PetscMin(subDim, origDim); l++) {
    154:               for (m = 0; m < PetscMin(subDim, origDim); m++) psi[(k * subDim + l) * subDim + m] += phi[(k * origDim + l) * origDim + m];
    155:             }
    156:           }
    157:         }
    158:         if (subsp->Ju) {
    159:           for (k = 0; k < subNc; k++) {
    160:             for (l = 0; l < origNc; l++) {
    161:               for (m = 0; m < subDim * subDim; m++) H[((j * subNb + i) * subNc + k) * subDim * subDim + m] += subsp->Ju[k * origNc + l] * psi[l * subDim * subDim + m];
    162:             }
    163:           }
    164:         } else {
    165:           for (k = 0; k < PetscMin(subNc, origNc); k++) {
    166:             for (l = 0; l < subDim * subDim; l++) H[((j * subNb + i) * subNc + k) * subDim * subDim + l] += psi[k * subDim * subDim + l];
    167:           }
    168:         }
    169:       }
    170:     }
    171:     PetscCall(DMRestoreWorkArray(sp->dm, subNc * origDim, MPIU_REAL, &psi));
    172:     PetscCall(DMRestoreWorkArray(sp->dm, origNc * origDim, MPIU_REAL, &phi));
    173:     PetscCall(DMRestoreWorkArray(sp->dm, npoints * origNb * origNc * origDim, MPIU_REAL, &inH));
    174:   }
    175:   if (D) {
    176:     PetscReal *phi, *psi;

    178:     PetscCall(DMGetWorkArray(sp->dm, origNc * origDim, MPIU_REAL, &phi));
    179:     PetscCall(DMGetWorkArray(sp->dm, origNc * subDim, MPIU_REAL, &psi));
    180:     for (i = 0; i < npoints * subNb * subNc * subDim; i++) D[i] = 0.0;
    181:     for (i = 0; i < subNb; i++) {
    182:       const PetscReal *subq = &subsp->Q[i * origNb];

    184:       for (j = 0; j < npoints; j++) {
    185:         for (k = 0; k < origNc * origDim; k++) phi[k] = 0.;
    186:         for (k = 0; k < origNc * subDim; k++) psi[k] = 0.;
    187:         for (k = 0; k < origNb; k++) {
    188:           for (l = 0; l < origNc * origDim; l++) phi[l] += inD[(j * origNb + k) * origNc * origDim + l] * subq[k];
    189:         }
    190:         if (subsp->Jx) {
    191:           for (k = 0; k < subNc; k++) {
    192:             for (l = 0; l < subDim; l++) {
    193:               for (m = 0; m < origDim; m++) psi[k * subDim + l] += subsp->Jx[m * subDim + l] * phi[k * origDim + m];
    194:             }
    195:           }
    196:         } else {
    197:           for (k = 0; k < subNc; k++) {
    198:             for (l = 0; l < PetscMin(subDim, origDim); l++) psi[k * subDim + l] += phi[k * origDim + l];
    199:           }
    200:         }
    201:         if (subsp->Ju) {
    202:           for (k = 0; k < subNc; k++) {
    203:             for (l = 0; l < origNc; l++) {
    204:               for (m = 0; m < subDim; m++) D[((j * subNb + i) * subNc + k) * subDim + m] += subsp->Ju[k * origNc + l] * psi[l * subDim + m];
    205:             }
    206:           }
    207:         } else {
    208:           for (k = 0; k < PetscMin(subNc, origNc); k++) {
    209:             for (l = 0; l < subDim; l++) D[((j * subNb + i) * subNc + k) * subDim + l] += psi[k * subDim + l];
    210:           }
    211:         }
    212:       }
    213:     }
    214:     PetscCall(DMRestoreWorkArray(sp->dm, subNc * origDim, MPIU_REAL, &psi));
    215:     PetscCall(DMRestoreWorkArray(sp->dm, origNc * origDim, MPIU_REAL, &phi));
    216:     PetscCall(DMRestoreWorkArray(sp->dm, npoints * origNb * origNc * origDim, MPIU_REAL, &inD));
    217:   }
    218:   if (B) {
    219:     PetscReal *phi;

    221:     PetscCall(DMGetWorkArray(sp->dm, origNc, MPIU_REAL, &phi));
    222:     if (subsp->u) {
    223:       for (i = 0; i < npoints * subNb; i++) {
    224:         for (j = 0; j < subNc; j++) B[i * subNc + j] = subsp->u[j];
    225:       }
    226:     } else {
    227:       for (i = 0; i < npoints * subNb * subNc; i++) B[i] = 0.0;
    228:     }
    229:     for (i = 0; i < subNb; i++) {
    230:       const PetscReal *subq = &subsp->Q[i * origNb];

    232:       for (j = 0; j < npoints; j++) {
    233:         for (k = 0; k < origNc; k++) phi[k] = 0.;
    234:         for (k = 0; k < origNb; k++) {
    235:           for (l = 0; l < origNc; l++) phi[l] += inB[(j * origNb + k) * origNc + l] * subq[k];
    236:         }
    237:         if (subsp->Ju) {
    238:           for (k = 0; k < subNc; k++) {
    239:             for (l = 0; l < origNc; l++) B[(j * subNb + i) * subNc + k] += subsp->Ju[k * origNc + l] * phi[l];
    240:           }
    241:         } else {
    242:           for (k = 0; k < PetscMin(subNc, origNc); k++) B[(j * subNb + i) * subNc + k] += phi[k];
    243:         }
    244:       }
    245:     }
    246:     PetscCall(DMRestoreWorkArray(sp->dm, origNc, MPIU_REAL, &phi));
    247:     PetscCall(DMRestoreWorkArray(sp->dm, npoints * origNb * origNc, MPIU_REAL, &inB));
    248:   }
    249:   PetscCall(DMRestoreWorkArray(sp->dm, npoints * origDim, MPIU_REAL, &inpoints));
    250:   PetscFunctionReturn(PETSC_SUCCESS);
    251: }

    253: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Subspace(PetscSpace sp)
    254: {
    255:   PetscSpace_Subspace *subsp;

    257:   PetscFunctionBegin;
    258:   PetscCall(PetscNew(&subsp));
    259:   sp->data = (void *)subsp;
    260:   PetscFunctionReturn(PETSC_SUCCESS);
    261: }

    263: static PetscErrorCode PetscSpaceGetDimension_Subspace(PetscSpace sp, PetscInt *dim)
    264: {
    265:   PetscSpace_Subspace *subsp;

    267:   PetscFunctionBegin;
    268:   subsp = (PetscSpace_Subspace *)sp->data;
    269:   *dim  = subsp->Nb;
    270:   PetscFunctionReturn(PETSC_SUCCESS);
    271: }

    273: static PetscErrorCode PetscSpaceSetUp_Subspace(PetscSpace sp)
    274: {
    275:   const PetscReal     *x;
    276:   const PetscReal     *Jx;
    277:   const PetscReal     *u;
    278:   const PetscReal     *Ju;
    279:   PetscDualSpace       dualSubspace;
    280:   PetscSpace           origSpace;
    281:   PetscInt             origDim, subDim, origNc, subNc, origNb, subNb, f, i, j, numPoints, offset;
    282:   PetscReal           *allPoints, *allWeights, *B, *V;
    283:   DM                   dm;
    284:   PetscSpace_Subspace *subsp = (PetscSpace_Subspace *)sp->data;

    286:   PetscFunctionBegin;
    287:   if (subsp->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
    288:   subsp->setupcalled = PETSC_TRUE;
    289:   x                  = subsp->x;
    290:   Jx                 = subsp->Jx;
    291:   u                  = subsp->u;
    292:   Ju                 = subsp->Ju;
    293:   origSpace          = subsp->origSpace;
    294:   dualSubspace       = subsp->dualSubspace;
    295:   PetscCall(PetscSpaceGetNumComponents(origSpace, &origNc));
    296:   PetscCall(PetscSpaceGetNumVariables(origSpace, &origDim));
    297:   PetscCall(PetscDualSpaceGetDM(dualSubspace, &dm));
    298:   PetscCall(DMGetDimension(dm, &subDim));
    299:   PetscCall(PetscSpaceGetDimension(origSpace, &origNb));
    300:   PetscCall(PetscDualSpaceGetDimension(dualSubspace, &subNb));
    301:   PetscCall(PetscDualSpaceGetNumComponents(dualSubspace, &subNc));

    303:   for (f = 0, numPoints = 0; f < subNb; f++) {
    304:     PetscQuadrature q;
    305:     PetscInt        qNp;

    307:     PetscCall(PetscDualSpaceGetFunctional(dualSubspace, f, &q));
    308:     PetscCall(PetscQuadratureGetData(q, NULL, NULL, &qNp, NULL, NULL));
    309:     numPoints += qNp;
    310:   }
    311:   PetscCall(PetscMalloc1(subNb * origNb, &V));
    312:   PetscCall(PetscMalloc3(numPoints * origDim, &allPoints, numPoints * origNc, &allWeights, numPoints * origNb * origNc, &B));
    313:   for (f = 0, offset = 0; f < subNb; f++) {
    314:     PetscQuadrature  q;
    315:     PetscInt         qNp, p;
    316:     const PetscReal *qp;
    317:     const PetscReal *qw;

    319:     PetscCall(PetscDualSpaceGetFunctional(dualSubspace, f, &q));
    320:     PetscCall(PetscQuadratureGetData(q, NULL, NULL, &qNp, &qp, &qw));
    321:     for (p = 0; p < qNp; p++, offset++) {
    322:       if (x) {
    323:         for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = x[i];
    324:       } else {
    325:         for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = 0.0;
    326:       }
    327:       if (Jx) {
    328:         for (i = 0; i < origDim; i++) {
    329:           for (j = 0; j < subDim; j++) allPoints[origDim * offset + i] += Jx[i * subDim + j] * qp[j];
    330:         }
    331:       } else {
    332:         for (i = 0; i < PetscMin(subDim, origDim); i++) allPoints[origDim * offset + i] += qp[i];
    333:       }
    334:       for (i = 0; i < origNc; i++) allWeights[origNc * offset + i] = 0.0;
    335:       if (Ju) {
    336:         for (i = 0; i < origNc; i++) {
    337:           for (j = 0; j < subNc; j++) allWeights[offset * origNc + i] += qw[j] * Ju[j * origNc + i];
    338:         }
    339:       } else {
    340:         for (i = 0; i < PetscMin(subNc, origNc); i++) allWeights[offset * origNc + i] += qw[i];
    341:       }
    342:     }
    343:   }
    344:   PetscCall(PetscSpaceEvaluate(origSpace, numPoints, allPoints, B, NULL, NULL));
    345:   for (f = 0, offset = 0; f < subNb; f++) {
    346:     PetscInt         b, p, s, qNp;
    347:     PetscQuadrature  q;
    348:     const PetscReal *qw;

    350:     PetscCall(PetscDualSpaceGetFunctional(dualSubspace, f, &q));
    351:     PetscCall(PetscQuadratureGetData(q, NULL, NULL, &qNp, NULL, &qw));
    352:     if (u) {
    353:       for (b = 0; b < origNb; b++) {
    354:         for (s = 0; s < subNc; s++) V[f * origNb + b] += qw[s] * u[s];
    355:       }
    356:     } else {
    357:       for (b = 0; b < origNb; b++) V[f * origNb + b] = 0.0;
    358:     }
    359:     for (p = 0; p < qNp; p++, offset++) {
    360:       for (b = 0; b < origNb; b++) {
    361:         for (s = 0; s < origNc; s++) V[f * origNb + b] += B[(offset * origNb + b) * origNc + s] * allWeights[offset * origNc + s];
    362:       }
    363:     }
    364:   }
    365:   /* orthnormalize rows of V */
    366:   for (f = 0; f < subNb; f++) {
    367:     PetscReal rho = 0.0, scal;

    369:     for (i = 0; i < origNb; i++) rho += PetscSqr(V[f * origNb + i]);

    371:     scal = 1. / PetscSqrtReal(rho);

    373:     for (i = 0; i < origNb; i++) V[f * origNb + i] *= scal;
    374:     for (j = f + 1; j < subNb; j++) {
    375:       for (i = 0, scal = 0.; i < origNb; i++) scal += V[f * origNb + i] * V[j * origNb + i];
    376:       for (i = 0; i < origNb; i++) V[j * origNb + i] -= V[f * origNb + i] * scal;
    377:     }
    378:   }
    379:   PetscCall(PetscFree3(allPoints, allWeights, B));
    380:   subsp->Q = V;
    381:   PetscFunctionReturn(PETSC_SUCCESS);
    382: }

    384: static PetscErrorCode PetscSpacePolynomialGetTensor_Subspace(PetscSpace sp, PetscBool *poly)
    385: {
    386:   PetscSpace_Subspace *subsp = (PetscSpace_Subspace *)sp->data;

    388:   PetscFunctionBegin;
    389:   *poly = PETSC_FALSE;
    390:   PetscCall(PetscSpacePolynomialGetTensor(subsp->origSpace, poly));
    391:   if (*poly) {
    392:     if (subsp->Jx) {
    393:       PetscInt subDim, origDim, i, j;
    394:       PetscInt maxnnz;

    396:       PetscCall(PetscSpaceGetNumVariables(subsp->origSpace, &origDim));
    397:       PetscCall(PetscSpaceGetNumVariables(sp, &subDim));
    398:       maxnnz = 0;
    399:       for (i = 0; i < origDim; i++) {
    400:         PetscInt nnz = 0;

    402:         for (j = 0; j < subDim; j++) nnz += (subsp->Jx[i * subDim + j] != 0.);
    403:         maxnnz = PetscMax(maxnnz, nnz);
    404:       }
    405:       for (j = 0; j < subDim; j++) {
    406:         PetscInt nnz = 0;

    408:         for (i = 0; i < origDim; i++) nnz += (subsp->Jx[i * subDim + j] != 0.);
    409:         maxnnz = PetscMax(maxnnz, nnz);
    410:       }
    411:       if (maxnnz > 1) *poly = PETSC_FALSE;
    412:     }
    413:   }
    414:   PetscFunctionReturn(PETSC_SUCCESS);
    415: }

    417: static PetscErrorCode PetscSpaceInitialize_Subspace(PetscSpace sp)
    418: {
    419:   PetscFunctionBegin;
    420:   sp->ops->setup        = PetscSpaceSetUp_Subspace;
    421:   sp->ops->view         = PetscSpaceView_Subspace;
    422:   sp->ops->destroy      = PetscSpaceDestroy_Subspace;
    423:   sp->ops->getdimension = PetscSpaceGetDimension_Subspace;
    424:   sp->ops->evaluate     = PetscSpaceEvaluate_Subspace;
    425:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Subspace));
    426:   PetscFunctionReturn(PETSC_SUCCESS);
    427: }
    428: /*@
    429:   PetscSpaceCreateSubspace - creates a subspace from a an `origSpace` and its dual `dualSubspace`

    431:   Input Parameters:
    432: + origSpace    - the original `PetscSpace`
    433: . dualSubspace - no idea
    434: . x            - no idea
    435: . Jx           - no idea
    436: . u            - no idea
    437: . Ju           - no idea
    438: - copymode     - whether to copy, borrow, or own some of the input arrays I guess

    440:   Output Parameter:
    441: . subspace - the subspace

    443:   Level: advanced

    445: .seealso: `PetscSpace`, `PetscDualSpace`, `PetscCopyMode`, `PetscSpaceType`
    446: @*/
    447: PetscErrorCode PetscSpaceCreateSubspace(PetscSpace origSpace, PetscDualSpace dualSubspace, PetscReal *x, PetscReal *Jx, PetscReal *u, PetscReal *Ju, PetscCopyMode copymode, PetscSpace *subspace)
    448: {
    449:   PetscSpace_Subspace *subsp;
    450:   PetscInt             origDim, subDim, origNc, subNc, subNb;
    451:   PetscInt             order;
    452:   DM                   dm;

    454:   PetscFunctionBegin;
    457:   if (x) PetscAssertPointer(x, 3);
    458:   if (Jx) PetscAssertPointer(Jx, 4);
    459:   if (u) PetscAssertPointer(u, 5);
    460:   if (Ju) PetscAssertPointer(Ju, 6);
    461:   PetscAssertPointer(subspace, 8);
    462:   PetscCall(PetscSpaceGetNumComponents(origSpace, &origNc));
    463:   PetscCall(PetscSpaceGetNumVariables(origSpace, &origDim));
    464:   PetscCall(PetscDualSpaceGetDM(dualSubspace, &dm));
    465:   PetscCall(DMGetDimension(dm, &subDim));
    466:   PetscCall(PetscDualSpaceGetDimension(dualSubspace, &subNb));
    467:   PetscCall(PetscDualSpaceGetNumComponents(dualSubspace, &subNc));
    468:   PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)origSpace), subspace));
    469:   PetscCall(PetscSpaceSetType(*subspace, PETSCSPACESUBSPACE));
    470:   PetscCall(PetscSpaceSetNumVariables(*subspace, subDim));
    471:   PetscCall(PetscSpaceSetNumComponents(*subspace, subNc));
    472:   PetscCall(PetscSpaceGetDegree(origSpace, &order, NULL));
    473:   PetscCall(PetscSpaceSetDegree(*subspace, order, PETSC_DETERMINE));
    474:   subsp     = (PetscSpace_Subspace *)(*subspace)->data;
    475:   subsp->Nb = subNb;
    476:   switch (copymode) {
    477:   case PETSC_OWN_POINTER:
    478:     if (x) subsp->x_alloc = x;
    479:     if (Jx) subsp->Jx_alloc = Jx;
    480:     if (u) subsp->u_alloc = u;
    481:     if (Ju) subsp->Ju_alloc = Ju;
    482:     /* fall through */
    483:   case PETSC_USE_POINTER:
    484:     if (x) subsp->x = x;
    485:     if (Jx) subsp->Jx = Jx;
    486:     if (u) subsp->u = u;
    487:     if (Ju) subsp->Ju = Ju;
    488:     break;
    489:   case PETSC_COPY_VALUES:
    490:     if (x) {
    491:       PetscCall(PetscMalloc1(origDim, &subsp->x_alloc));
    492:       PetscCall(PetscArraycpy(subsp->x_alloc, x, origDim));
    493:       subsp->x = subsp->x_alloc;
    494:     }
    495:     if (Jx) {
    496:       PetscCall(PetscMalloc1(origDim * subDim, &subsp->Jx_alloc));
    497:       PetscCall(PetscArraycpy(subsp->Jx_alloc, Jx, origDim * subDim));
    498:       subsp->Jx = subsp->Jx_alloc;
    499:     }
    500:     if (u) {
    501:       PetscCall(PetscMalloc1(subNc, &subsp->u_alloc));
    502:       PetscCall(PetscArraycpy(subsp->u_alloc, u, subNc));
    503:       subsp->u = subsp->u_alloc;
    504:     }
    505:     if (Ju) {
    506:       PetscCall(PetscMalloc1(origNc * subNc, &subsp->Ju_alloc));
    507:       PetscCall(PetscArraycpy(subsp->Ju_alloc, Ju, origNc * subNc));
    508:       subsp->Ju = subsp->Ju_alloc;
    509:     }
    510:     break;
    511:   default:
    512:     SETERRQ(PetscObjectComm((PetscObject)origSpace), PETSC_ERR_ARG_OUTOFRANGE, "Unknown copy mode");
    513:   }
    514:   PetscCall(PetscObjectReference((PetscObject)origSpace));
    515:   subsp->origSpace = origSpace;
    516:   PetscCall(PetscObjectReference((PetscObject)dualSubspace));
    517:   subsp->dualSubspace = dualSubspace;
    518:   PetscCall(PetscSpaceInitialize_Subspace(*subspace));
    519:   PetscFunctionReturn(PETSC_SUCCESS);
    520: }