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++) {
 57:         PetscViewerASCIIPrintf(viewer," %g\n", (double)subsp->x[o]);
 58:       }
 59:     }
 60:     if (subsp->Jx) {
 61:       PetscViewerASCIIPrintf(viewer,"Subspace-to-space domain transform:\n\n");
 62:       for (o = 0; o < origDim; o++) {
 63:         PetscViewerASCIIPrintf(viewer," %g", (double)subsp->Jx[o * subDim + 0]);
 64:         PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
 65:         for (s = 1; s < subDim; s++) {
 66:           PetscViewerASCIIPrintf(viewer," %g", (double)subsp->Jx[o * subDim + s]);
 67:         }
 68:         PetscViewerASCIIPrintf(viewer,"\n");
 69:         PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
 70:       }
 71:     }
 72:     if (subsp->u) {
 73:       PetscViewerASCIIPrintf(viewer,"Space-to-subspace range shift:\n\n");
 74:       for (o = 0; o < origNc; o++) {
 75:         PetscViewerASCIIPrintf(viewer," %d\n", subsp->u[o]);
 76:       }
 77:     }
 78:     if (subsp->Ju) {
 79:       PetscViewerASCIIPrintf(viewer,"Space-to-subsace domain transform:\n");
 80:       for (o = 0; o < origNc; o++) {
 81:         PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
 82:         for (s = 0; s < subNc; s++) {
 83:           PetscViewerASCIIPrintf(viewer," %d", subsp->Ju[o * subNc + s]);
 84:         }
 85:         PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
 86:       }
 87:       PetscViewerASCIIPrintf(viewer,"\n");
 88:     }
 89:     PetscViewerASCIIPrintf(viewer,"Original space:\n");
 90:   }
 91:   PetscViewerASCIIPushTab(viewer);
 92:   PetscSpaceView(subsp->origSpace,viewer);
 93:   PetscViewerASCIIPopTab(viewer);
 94:   return 0;
 95: }

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

104:   origsp = subsp->origSpace;
105:   PetscSpaceGetNumVariables(sp,&subDim);
106:   PetscSpaceGetNumVariables(origsp,&origDim);
107:   PetscSpaceGetNumComponents(sp,&subNc);
108:   PetscSpaceGetNumComponents(origsp,&origNc);
109:   PetscSpaceGetDimension(sp,&subNb);
110:   PetscSpaceGetDimension(origsp,&origNb);
111:   DMGetWorkArray(sp->dm,npoints*origDim,MPIU_REAL,&inpoints);
112:   for (i = 0; i < npoints; i++) {
113:     if (subsp->x) {
114:       for (j = 0; j < origDim; j++) inpoints[i * origDim + j] = subsp->x[j];
115:     } else {
116:       for (j = 0; j < origDim; j++) inpoints[i * origDim + j] = 0.0;
117:     }
118:     if (subsp->Jx) {
119:       for (j = 0; j < origDim; j++) {
120:         for (k = 0; k < subDim; k++) {
121:           inpoints[i * origDim + j] += subsp->Jx[j * subDim + k] * points[i * subDim + k];
122:         }
123:       }
124:     } else {
125:       for (j = 0; j < PetscMin(subDim, origDim); j++) {
126:         inpoints[i * origDim + j] += points[i * subDim + j];
127:       }
128:     }
129:   }
130:   if (B) {
131:     DMGetWorkArray(sp->dm,npoints*origNb*origNc,MPIU_REAL,&inB);
132:   }
133:   if (D) {
134:     DMGetWorkArray(sp->dm,npoints*origNb*origNc*origDim,MPIU_REAL,&inD);
135:   }
136:   if (H) {
137:     DMGetWorkArray(sp->dm,npoints*origNb*origNc*origDim*origDim,MPIU_REAL,&inH);
138:   }
139:   PetscSpaceEvaluate(origsp,npoints,inpoints,inB,inD,inH);
140:   if (H) {
141:     PetscReal *phi, *psi;

143:     DMGetWorkArray(sp->dm,origNc*origDim*origDim,MPIU_REAL,&phi);
144:     DMGetWorkArray(sp->dm,origNc*subDim*subDim,MPIU_REAL,&psi);
145:     for (i = 0; i < npoints * subNb * subNc * subDim; i++) D[i] = 0.0;
146:     for (i = 0; i < subNb; i++) {
147:       const PetscReal *subq = &subsp->Q[i * origNb];

149:       for (j = 0; j < npoints; j++) {
150:         for (k = 0; k < origNc * origDim; k++) phi[k] = 0.;
151:         for (k = 0; k < origNc * subDim; k++) psi[k] = 0.;
152:         for (k = 0; k < origNb; k++) {
153:           for (l = 0; l < origNc * origDim * origDim; l++) {
154:             phi[l] += inH[(j * origNb + k) * origNc * origDim * origDim + l] * subq[k];
155:           }
156:         }
157:         if (subsp->Jx) {
158:           for (k = 0; k < subNc; k++) {
159:             for (l = 0; l < subDim; l++) {
160:               for (m = 0; m < origDim; m++) {
161:                 for (n = 0; n < subDim; n++) {
162:                   for (o = 0; o < origDim; o++) {
163:                     psi[(k * subDim + l) * subDim + n] += subsp->Jx[m * subDim + l] * subsp->Jx[o * subDim + n] * phi[(k * origDim + m) * origDim + o];
164:                   }
165:                 }
166:               }
167:             }
168:           }
169:         } else {
170:           for (k = 0; k < subNc; k++) {
171:             for (l = 0; l < PetscMin(subDim, origDim); l++) {
172:               for (m = 0; m < PetscMin(subDim, origDim); m++) {
173:                 psi[(k * subDim + l) * subDim + m] += phi[(k * origDim + l) * origDim + m];
174:               }
175:             }
176:           }
177:         }
178:         if (subsp->Ju) {
179:           for (k = 0; k < subNc; k++) {
180:             for (l = 0; l < origNc; l++) {
181:               for (m = 0; m < subDim * subDim; m++) {
182:                 H[((j * subNb + i) * subNc + k) * subDim * subDim + m] += subsp->Ju[k * origNc + l] * psi[l * subDim * subDim + m];
183:               }
184:             }
185:           }
186:         }
187:         else {
188:           for (k = 0; k < PetscMin(subNc, origNc); k++) {
189:             for (l = 0; l < subDim * subDim; l++) {
190:               H[((j * subNb + i) * subNc + k) * subDim * subDim + l] += psi[k * subDim * subDim + l];
191:             }
192:           }
193:         }
194:       }
195:     }
196:     DMRestoreWorkArray(sp->dm,subNc*origDim,MPIU_REAL,&psi);
197:     DMRestoreWorkArray(sp->dm,origNc*origDim,MPIU_REAL,&phi);
198:     DMRestoreWorkArray(sp->dm,npoints*origNb*origNc*origDim,MPIU_REAL,&inH);
199:   }
200:   if (D) {
201:     PetscReal *phi, *psi;

203:     DMGetWorkArray(sp->dm,origNc*origDim,MPIU_REAL,&phi);
204:     DMGetWorkArray(sp->dm,origNc*subDim,MPIU_REAL,&psi);
205:     for (i = 0; i < npoints * subNb * subNc * subDim; i++) D[i] = 0.0;
206:     for (i = 0; i < subNb; i++) {
207:       const PetscReal *subq = &subsp->Q[i * origNb];

209:       for (j = 0; j < npoints; j++) {
210:         for (k = 0; k < origNc * origDim; k++) phi[k] = 0.;
211:         for (k = 0; k < origNc * subDim; k++) psi[k] = 0.;
212:         for (k = 0; k < origNb; k++) {
213:           for (l = 0; l < origNc * origDim; l++) {
214:             phi[l] += inD[(j * origNb + k) * origNc * origDim + l] * subq[k];
215:           }
216:         }
217:         if (subsp->Jx) {
218:           for (k = 0; k < subNc; k++) {
219:             for (l = 0; l < subDim; l++) {
220:               for (m = 0; m < origDim; m++) {
221:                 psi[k * subDim + l] += subsp->Jx[m * subDim + l] * phi[k * origDim + m];
222:               }
223:             }
224:           }
225:         } else {
226:           for (k = 0; k < subNc; k++) {
227:             for (l = 0; l < PetscMin(subDim, origDim); l++) {
228:               psi[k * subDim + l] += phi[k * origDim + l];
229:             }
230:           }
231:         }
232:         if (subsp->Ju) {
233:           for (k = 0; k < subNc; k++) {
234:             for (l = 0; l < origNc; l++) {
235:               for (m = 0; m < subDim; m++) {
236:                 D[((j * subNb + i) * subNc + k) * subDim + m] += subsp->Ju[k * origNc + l] * psi[l * subDim + m];
237:               }
238:             }
239:           }
240:         }
241:         else {
242:           for (k = 0; k < PetscMin(subNc, origNc); k++) {
243:             for (l = 0; l < subDim; l++) {
244:               D[((j * subNb + i) * subNc + k) * subDim + l] += psi[k * subDim + l];
245:             }
246:           }
247:         }
248:       }
249:     }
250:     DMRestoreWorkArray(sp->dm,subNc*origDim,MPIU_REAL,&psi);
251:     DMRestoreWorkArray(sp->dm,origNc*origDim,MPIU_REAL,&phi);
252:     DMRestoreWorkArray(sp->dm,npoints*origNb*origNc*origDim,MPIU_REAL,&inD);
253:   }
254:   if (B) {
255:     PetscReal *phi;

257:     DMGetWorkArray(sp->dm,origNc,MPIU_REAL,&phi);
258:     if (subsp->u) {
259:       for (i = 0; i < npoints * subNb; i++) {
260:         for (j = 0; j < subNc; j++) B[i * subNc + j] = subsp->u[j];
261:       }
262:     } else {
263:       for (i = 0; i < npoints * subNb * subNc; i++) B[i] = 0.0;
264:     }
265:     for (i = 0; i < subNb; i++) {
266:       const PetscReal *subq = &subsp->Q[i * origNb];

268:       for (j = 0; j < npoints; j++) {
269:         for (k = 0; k < origNc; k++) phi[k] = 0.;
270:         for (k = 0; k < origNb; k++) {
271:           for (l = 0; l < origNc; l++) {
272:             phi[l] += inB[(j * origNb + k) * origNc + l] * subq[k];
273:           }
274:         }
275:         if (subsp->Ju) {
276:           for (k = 0; k < subNc; k++) {
277:             for (l = 0; l < origNc; l++) {
278:               B[(j * subNb + i) * subNc + k] += subsp->Ju[k * origNc + l] * phi[l];
279:             }
280:           }
281:         }
282:         else {
283:           for (k = 0; k < PetscMin(subNc, origNc); k++) {
284:             B[(j * subNb + i) * subNc + k] += phi[k];
285:           }
286:         }
287:       }
288:     }
289:     DMRestoreWorkArray(sp->dm,origNc,MPIU_REAL,&phi);
290:     DMRestoreWorkArray(sp->dm,npoints*origNb*origNc,MPIU_REAL,&inB);
291:   }
292:   DMRestoreWorkArray(sp->dm,npoints*origDim,MPIU_REAL,&inpoints);
293:   return 0;
294: }

296: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Subspace(PetscSpace sp)
297: {
298:   PetscSpace_Subspace *subsp;

300:   PetscNewLog(sp,&subsp);
301:   sp->data = (void *) subsp;
302:   return 0;
303: }

305: static PetscErrorCode PetscSpaceGetDimension_Subspace(PetscSpace sp, PetscInt *dim)
306: {
307:   PetscSpace_Subspace *subsp;

309:   subsp = (PetscSpace_Subspace *) sp->data;
310:   *dim = subsp->Nb;
311:   return 0;
312: }

314: static PetscErrorCode PetscSpaceSetUp_Subspace(PetscSpace sp)
315: {
316:   const PetscReal     *x;
317:   const PetscReal     *Jx;
318:   const PetscReal     *u;
319:   const PetscReal     *Ju;
320:   PetscDualSpace      dualSubspace;
321:   PetscSpace          origSpace;
322:   PetscInt            origDim, subDim, origNc, subNc, origNb, subNb, f, i, j, numPoints, offset;
323:   PetscReal           *allPoints, *allWeights, *B, *V;
324:   DM                  dm;
325:   PetscSpace_Subspace *subsp;

327:   subsp = (PetscSpace_Subspace *) sp->data;
328:   x            = subsp->x;
329:   Jx           = subsp->Jx;
330:   u            = subsp->u;
331:   Ju           = subsp->Ju;
332:   origSpace    = subsp->origSpace;
333:   dualSubspace = subsp->dualSubspace;
334:   PetscSpaceGetNumComponents(origSpace,&origNc);
335:   PetscSpaceGetNumVariables(origSpace,&origDim);
336:   PetscDualSpaceGetDM(dualSubspace,&dm);
337:   DMGetDimension(dm,&subDim);
338:   PetscSpaceGetDimension(origSpace,&origNb);
339:   PetscDualSpaceGetDimension(dualSubspace,&subNb);
340:   PetscDualSpaceGetNumComponents(dualSubspace,&subNc);

342:   for (f = 0, numPoints = 0; f < subNb; f++) {
343:     PetscQuadrature q;
344:     PetscInt        qNp;

346:     PetscDualSpaceGetFunctional(dualSubspace,f,&q);
347:     PetscQuadratureGetData(q,NULL,NULL,&qNp,NULL,NULL);
348:     numPoints += qNp;
349:   }
350:   PetscMalloc1(subNb*origNb,&V);
351:   PetscMalloc3(numPoints*origDim,&allPoints,numPoints*origNc,&allWeights,numPoints*origNb*origNc,&B);
352:   for (f = 0, offset = 0; f < subNb; f++) {
353:     PetscQuadrature q;
354:     PetscInt        qNp, p;
355:     const PetscReal *qp;
356:     const PetscReal *qw;

358:     PetscDualSpaceGetFunctional(dualSubspace,f,&q);
359:     PetscQuadratureGetData(q,NULL,NULL,&qNp,&qp,&qw);
360:     for (p = 0; p < qNp; p++, offset++) {
361:       if (x) {
362:         for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = x[i];
363:       } else {
364:         for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = 0.0;
365:       }
366:       if (Jx) {
367:         for (i = 0; i < origDim; i++) {
368:           for (j = 0; j < subDim; j++) {
369:             allPoints[origDim * offset + i] += Jx[i * subDim + j] * qp[j];
370:           }
371:         }
372:       } else {
373:         for (i = 0; i < PetscMin(subDim, origDim); i++) allPoints[origDim * offset + i] += qp[i];
374:       }
375:       for (i = 0; i < origNc; i++) allWeights[origNc * offset + i] = 0.0;
376:       if (Ju) {
377:         for (i = 0; i < origNc; i++) {
378:           for (j = 0; j < subNc; j++) {
379:             allWeights[offset * origNc + i] += qw[j] * Ju[j * origNc + i];
380:           }
381:         }
382:       } else {
383:         for (i = 0; i < PetscMin(subNc, origNc); i++) allWeights[offset * origNc + i] += qw[i];
384:       }
385:     }
386:   }
387:   PetscSpaceEvaluate(origSpace,numPoints,allPoints,B,NULL,NULL);
388:   for (f = 0, offset = 0; f < subNb; f++) {
389:     PetscInt b, p, s, qNp;
390:     PetscQuadrature q;
391:     const PetscReal *qw;

393:     PetscDualSpaceGetFunctional(dualSubspace,f,&q);
394:     PetscQuadratureGetData(q,NULL,NULL,&qNp,NULL,&qw);
395:     if (u) {
396:       for (b = 0; b < origNb; b++) {
397:         for (s = 0; s < subNc; s++) {
398:           V[f * origNb + b] += qw[s] * u[s];
399:         }
400:       }
401:     } else {
402:       for (b = 0; b < origNb; b++) V[f * origNb + b] = 0.0;
403:     }
404:     for (p = 0; p < qNp; p++, offset++) {
405:       for (b = 0; b < origNb; b++) {
406:         for (s = 0; s < origNc; s++) {
407:           V[f * origNb + b] += B[(offset * origNb + b) * origNc + s] * allWeights[offset * origNc + s];
408:         }
409:       }
410:     }
411:   }
412:   /* orthnormalize rows of V */
413:   for (f = 0; f < subNb; f++) {
414:     PetscReal rho = 0.0, scal;

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

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

420:     for (i = 0; i < origNb; i++) V[f * origNb + i] *= scal;
421:     for (j = f + 1; j < subNb; j++) {
422:       for (i = 0, scal = 0.; i < origNb; i++) scal += V[f * origNb + i] * V[j * origNb + i];
423:       for (i = 0; i < origNb; i++) V[j * origNb + i] -= V[f * origNb + i] * scal;
424:     }
425:   }
426:   PetscFree3(allPoints,allWeights,B);
427:   subsp->Q = V;
428:   return 0;
429: }

431: static PetscErrorCode PetscSpacePolynomialGetTensor_Subspace(PetscSpace sp, PetscBool *poly)
432: {
433:   PetscSpace_Subspace *subsp = (PetscSpace_Subspace *) sp->data;

435:   *poly = PETSC_FALSE;
436:   PetscSpacePolynomialGetTensor(subsp->origSpace,poly);
437:   if (*poly) {
438:     if (subsp->Jx) {
439:       PetscInt subDim, origDim, i, j;
440:       PetscInt maxnnz;

442:       PetscSpaceGetNumVariables(subsp->origSpace,&origDim);
443:       PetscSpaceGetNumVariables(sp,&subDim);
444:       maxnnz = 0;
445:       for (i = 0; i < origDim; i++) {
446:         PetscInt nnz = 0;

448:         for (j = 0; j < subDim; j++) nnz += (subsp->Jx[i * subDim + j] != 0.);
449:         maxnnz = PetscMax(maxnnz,nnz);
450:       }
451:       for (j = 0; j < subDim; j++) {
452:         PetscInt nnz = 0;

454:         for (i = 0; i < origDim; i++) nnz += (subsp->Jx[i * subDim + j] != 0.);
455:         maxnnz = PetscMax(maxnnz,nnz);
456:       }
457:       if (maxnnz > 1) *poly = PETSC_FALSE;
458:     }
459:   }
460:   return 0;
461: }

463: static PetscErrorCode PetscSpaceInitialize_Subspace(PetscSpace sp)
464: {
465:   sp->ops->setup = PetscSpaceSetUp_Subspace;
466:   sp->ops->view  = PetscSpaceView_Subspace;
467:   sp->ops->destroy  = PetscSpaceDestroy_Subspace;
468:   sp->ops->getdimension  = PetscSpaceGetDimension_Subspace;
469:   sp->ops->evaluate = PetscSpaceEvaluate_Subspace;
470:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Subspace);
471:   return 0;
472: }

474: PetscErrorCode PetscSpaceCreateSubspace(PetscSpace origSpace, PetscDualSpace dualSubspace, PetscReal *x, PetscReal *Jx, PetscReal *u, PetscReal *Ju, PetscCopyMode copymode, PetscSpace *subspace)
475: {
476:   PetscSpace_Subspace *subsp;
477:   PetscInt            origDim, subDim, origNc, subNc, subNb;
478:   PetscInt            order;
479:   DM                  dm;

488:   PetscSpaceGetNumComponents(origSpace,&origNc);
489:   PetscSpaceGetNumVariables(origSpace,&origDim);
490:   PetscDualSpaceGetDM(dualSubspace,&dm);
491:   DMGetDimension(dm,&subDim);
492:   PetscDualSpaceGetDimension(dualSubspace,&subNb);
493:   PetscDualSpaceGetNumComponents(dualSubspace,&subNc);
494:   PetscSpaceCreate(PetscObjectComm((PetscObject)origSpace),subspace);
495:   PetscSpaceSetType(*subspace,PETSCSPACESUBSPACE);
496:   PetscSpaceSetNumVariables(*subspace,subDim);
497:   PetscSpaceSetNumComponents(*subspace,subNc);
498:   PetscSpaceGetDegree(origSpace,&order,NULL);
499:   PetscSpaceSetDegree(*subspace,order,PETSC_DETERMINE);
500:   subsp = (PetscSpace_Subspace *) (*subspace)->data;
501:   subsp->Nb = subNb;
502:   switch (copymode) {
503:   case PETSC_OWN_POINTER:
504:     if (x) subsp->x_alloc = x;
505:     if (Jx) subsp->Jx_alloc = Jx;
506:     if (u) subsp->u_alloc = u;
507:     if (Ju) subsp->Ju_alloc = Ju;
508:   case PETSC_USE_POINTER:
509:     if (x) subsp->x = x;
510:     if (Jx) subsp->Jx = Jx;
511:     if (u) subsp->u = u;
512:     if (Ju) subsp->Ju = Ju;
513:     break;
514:   case PETSC_COPY_VALUES:
515:     if (x) {
516:       PetscMalloc1(origDim,&subsp->x_alloc);
517:       PetscArraycpy(subsp->x_alloc,x,origDim);
518:       subsp->x = subsp->x_alloc;
519:     }
520:     if (Jx) {
521:       PetscMalloc1(origDim * subDim,&subsp->Jx_alloc);
522:       PetscArraycpy(subsp->Jx_alloc,Jx,origDim * subDim);
523:       subsp->Jx = subsp->Jx_alloc;
524:     }
525:     if (u) {
526:       PetscMalloc1(subNc,&subsp->u_alloc);
527:       PetscArraycpy(subsp->u_alloc,u,subNc);
528:       subsp->u = subsp->u_alloc;
529:     }
530:     if (Ju) {
531:       PetscMalloc1(origNc * subNc,&subsp->Ju_alloc);
532:       PetscArraycpy(subsp->Ju_alloc,Ju,origNc * subNc);
533:       subsp->Ju = subsp->Ju_alloc;
534:     }
535:     break;
536:   default:
537:     SETERRQ(PetscObjectComm((PetscObject)origSpace),PETSC_ERR_ARG_OUTOFRANGE,"Unknown copy mode");
538:   }
539:   PetscObjectReference((PetscObject)origSpace);
540:   subsp->origSpace = origSpace;
541:   PetscObjectReference((PetscObject)dualSubspace);
542:   subsp->dualSubspace = dualSubspace;
543:   PetscSpaceInitialize_Subspace(*subspace);
544:   return 0;
545: }