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