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:   subsp    = (PetscSpace_Subspace *)sp->data;
 23:   subsp->x = NULL;
 24:   PetscFree(subsp->x_alloc);
 25:   subsp->Jx = NULL;
 26:   PetscFree(subsp->Jx_alloc);
 27:   subsp->u = NULL;
 28:   PetscFree(subsp->u_alloc);
 29:   subsp->Ju = NULL;
 30:   PetscFree(subsp->Ju_alloc);
 31:   PetscFree(subsp->Q);
 32:   PetscSpaceDestroy(&subsp->origSpace);
 33:   PetscDualSpaceDestroy(&subsp->dualSubspace);
 34:   PetscFree(subsp);
 35:   sp->data = NULL;
 36:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialGetTensor_C", NULL);
 37:   return 0;
 38: }

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

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

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

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

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

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

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

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

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

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

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

249: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Subspace(PetscSpace sp)
250: {
251:   PetscSpace_Subspace *subsp;

253:   PetscNew(&subsp);
254:   sp->data = (void *)subsp;
255:   return 0;
256: }

258: static PetscErrorCode PetscSpaceGetDimension_Subspace(PetscSpace sp, PetscInt *dim)
259: {
260:   PetscSpace_Subspace *subsp;

262:   subsp = (PetscSpace_Subspace *)sp->data;
263:   *dim  = subsp->Nb;
264:   return 0;
265: }

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

280:   subsp        = (PetscSpace_Subspace *)sp->data;
281:   x            = subsp->x;
282:   Jx           = subsp->Jx;
283:   u            = subsp->u;
284:   Ju           = subsp->Ju;
285:   origSpace    = subsp->origSpace;
286:   dualSubspace = subsp->dualSubspace;
287:   PetscSpaceGetNumComponents(origSpace, &origNc);
288:   PetscSpaceGetNumVariables(origSpace, &origDim);
289:   PetscDualSpaceGetDM(dualSubspace, &dm);
290:   DMGetDimension(dm, &subDim);
291:   PetscSpaceGetDimension(origSpace, &origNb);
292:   PetscDualSpaceGetDimension(dualSubspace, &subNb);
293:   PetscDualSpaceGetNumComponents(dualSubspace, &subNc);

295:   for (f = 0, numPoints = 0; f < subNb; f++) {
296:     PetscQuadrature q;
297:     PetscInt        qNp;

299:     PetscDualSpaceGetFunctional(dualSubspace, f, &q);
300:     PetscQuadratureGetData(q, NULL, NULL, &qNp, NULL, NULL);
301:     numPoints += qNp;
302:   }
303:   PetscMalloc1(subNb * origNb, &V);
304:   PetscMalloc3(numPoints * origDim, &allPoints, numPoints * origNc, &allWeights, numPoints * origNb * origNc, &B);
305:   for (f = 0, offset = 0; f < subNb; f++) {
306:     PetscQuadrature  q;
307:     PetscInt         qNp, p;
308:     const PetscReal *qp;
309:     const PetscReal *qw;

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

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

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

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

365:     for (i = 0; i < origNb; i++) V[f * origNb + i] *= scal;
366:     for (j = f + 1; j < subNb; j++) {
367:       for (i = 0, scal = 0.; i < origNb; i++) scal += V[f * origNb + i] * V[j * origNb + i];
368:       for (i = 0; i < origNb; i++) V[j * origNb + i] -= V[f * origNb + i] * scal;
369:     }
370:   }
371:   PetscFree3(allPoints, allWeights, B);
372:   subsp->Q = V;
373:   return 0;
374: }

376: static PetscErrorCode PetscSpacePolynomialGetTensor_Subspace(PetscSpace sp, PetscBool *poly)
377: {
378:   PetscSpace_Subspace *subsp = (PetscSpace_Subspace *)sp->data;

380:   *poly = PETSC_FALSE;
381:   PetscSpacePolynomialGetTensor(subsp->origSpace, poly);
382:   if (*poly) {
383:     if (subsp->Jx) {
384:       PetscInt subDim, origDim, i, j;
385:       PetscInt maxnnz;

387:       PetscSpaceGetNumVariables(subsp->origSpace, &origDim);
388:       PetscSpaceGetNumVariables(sp, &subDim);
389:       maxnnz = 0;
390:       for (i = 0; i < origDim; i++) {
391:         PetscInt nnz = 0;

393:         for (j = 0; j < subDim; j++) nnz += (subsp->Jx[i * subDim + j] != 0.);
394:         maxnnz = PetscMax(maxnnz, nnz);
395:       }
396:       for (j = 0; j < subDim; j++) {
397:         PetscInt nnz = 0;

399:         for (i = 0; i < origDim; i++) nnz += (subsp->Jx[i * subDim + j] != 0.);
400:         maxnnz = PetscMax(maxnnz, nnz);
401:       }
402:       if (maxnnz > 1) *poly = PETSC_FALSE;
403:     }
404:   }
405:   return 0;
406: }

408: static PetscErrorCode PetscSpaceInitialize_Subspace(PetscSpace sp)
409: {
410:   sp->ops->setup        = PetscSpaceSetUp_Subspace;
411:   sp->ops->view         = PetscSpaceView_Subspace;
412:   sp->ops->destroy      = PetscSpaceDestroy_Subspace;
413:   sp->ops->getdimension = PetscSpaceGetDimension_Subspace;
414:   sp->ops->evaluate     = PetscSpaceEvaluate_Subspace;
415:   PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Subspace);
416:   return 0;
417: }

419: PetscErrorCode PetscSpaceCreateSubspace(PetscSpace origSpace, PetscDualSpace dualSubspace, PetscReal *x, PetscReal *Jx, PetscReal *u, PetscReal *Ju, PetscCopyMode copymode, PetscSpace *subspace)
420: {
421:   PetscSpace_Subspace *subsp;
422:   PetscInt             origDim, subDim, origNc, subNc, subNb;
423:   PetscInt             order;
424:   DM                   dm;

433:   PetscSpaceGetNumComponents(origSpace, &origNc);
434:   PetscSpaceGetNumVariables(origSpace, &origDim);
435:   PetscDualSpaceGetDM(dualSubspace, &dm);
436:   DMGetDimension(dm, &subDim);
437:   PetscDualSpaceGetDimension(dualSubspace, &subNb);
438:   PetscDualSpaceGetNumComponents(dualSubspace, &subNc);
439:   PetscSpaceCreate(PetscObjectComm((PetscObject)origSpace), subspace);
440:   PetscSpaceSetType(*subspace, PETSCSPACESUBSPACE);
441:   PetscSpaceSetNumVariables(*subspace, subDim);
442:   PetscSpaceSetNumComponents(*subspace, subNc);
443:   PetscSpaceGetDegree(origSpace, &order, NULL);
444:   PetscSpaceSetDegree(*subspace, order, PETSC_DETERMINE);
445:   subsp     = (PetscSpace_Subspace *)(*subspace)->data;
446:   subsp->Nb = subNb;
447:   switch (copymode) {
448:   case PETSC_OWN_POINTER:
449:     if (x) subsp->x_alloc = x;
450:     if (Jx) subsp->Jx_alloc = Jx;
451:     if (u) subsp->u_alloc = u;
452:     if (Ju) subsp->Ju_alloc = Ju;
453:     /* fall through */
454:   case PETSC_USE_POINTER:
455:     if (x) subsp->x = x;
456:     if (Jx) subsp->Jx = Jx;
457:     if (u) subsp->u = u;
458:     if (Ju) subsp->Ju = Ju;
459:     break;
460:   case PETSC_COPY_VALUES:
461:     if (x) {
462:       PetscMalloc1(origDim, &subsp->x_alloc);
463:       PetscArraycpy(subsp->x_alloc, x, origDim);
464:       subsp->x = subsp->x_alloc;
465:     }
466:     if (Jx) {
467:       PetscMalloc1(origDim * subDim, &subsp->Jx_alloc);
468:       PetscArraycpy(subsp->Jx_alloc, Jx, origDim * subDim);
469:       subsp->Jx = subsp->Jx_alloc;
470:     }
471:     if (u) {
472:       PetscMalloc1(subNc, &subsp->u_alloc);
473:       PetscArraycpy(subsp->u_alloc, u, subNc);
474:       subsp->u = subsp->u_alloc;
475:     }
476:     if (Ju) {
477:       PetscMalloc1(origNc * subNc, &subsp->Ju_alloc);
478:       PetscArraycpy(subsp->Ju_alloc, Ju, origNc * subNc);
479:       subsp->Ju = subsp->Ju_alloc;
480:     }
481:     break;
482:   default:
483:     SETERRQ(PetscObjectComm((PetscObject)origSpace), PETSC_ERR_ARG_OUTOFRANGE, "Unknown copy mode");
484:   }
485:   PetscObjectReference((PetscObject)origSpace);
486:   subsp->origSpace = origSpace;
487:   PetscObjectReference((PetscObject)dualSubspace);
488:   subsp->dualSubspace = dualSubspace;
489:   PetscSpaceInitialize_Subspace(*subspace);
490:   return 0;
491: }