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: } PetscSpace_Subspace;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

285:   PetscFunctionBegin;
286:   subsp        = (PetscSpace_Subspace *)sp->data;
287:   x            = subsp->x;
288:   Jx           = subsp->Jx;
289:   u            = subsp->u;
290:   Ju           = subsp->Ju;
291:   origSpace    = subsp->origSpace;
292:   dualSubspace = subsp->dualSubspace;
293:   PetscCall(PetscSpaceGetNumComponents(origSpace, &origNc));
294:   PetscCall(PetscSpaceGetNumVariables(origSpace, &origDim));
295:   PetscCall(PetscDualSpaceGetDM(dualSubspace, &dm));
296:   PetscCall(DMGetDimension(dm, &subDim));
297:   PetscCall(PetscSpaceGetDimension(origSpace, &origNb));
298:   PetscCall(PetscDualSpaceGetDimension(dualSubspace, &subNb));
299:   PetscCall(PetscDualSpaceGetNumComponents(dualSubspace, &subNc));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

438:     Output Parameter:
439: .   subspace - the subspace

441:     Level: advanced

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

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