Actual source code: spacesubspace.c

petsc-master 2019-12-10
Report Typos and Errors
  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;
 21:   PetscErrorCode      ierr;

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

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

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

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

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

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

149:     DMGetWorkArray(sp->dm,origNc*origDim*origDim,MPIU_REAL,&phi);
150:     DMGetWorkArray(sp->dm,origNc*subDim*subDim,MPIU_REAL,&psi);
151:     for (i = 0; i < npoints * subNb * subNc * subDim; i++) D[i] = 0.0;
152:     for (i = 0; i < subNb; i++) {
153:       const PetscReal *subq = &subsp->Q[i * origNb];

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

209:     DMGetWorkArray(sp->dm,origNc*origDim,MPIU_REAL,&phi);
210:     DMGetWorkArray(sp->dm,origNc*subDim,MPIU_REAL,&psi);
211:     for (i = 0; i < npoints * subNb * subNc * subDim; i++) D[i] = 0.0;
212:     for (i = 0; i < subNb; i++) {
213:       const PetscReal *subq = &subsp->Q[i * origNb];

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

263:     DMGetWorkArray(sp->dm,origNc,MPIU_REAL,&phi);
264:     if (subsp->u) {
265:       for (i = 0; i < npoints * subNb; i++) {
266:         for (j = 0; j < subNc; j++) B[i * subNc + j] = subsp->u[j];
267:       }
268:     } else {
269:       for (i = 0; i < npoints * subNb * subNc; i++) B[i] = 0.0;
270:     }
271:     for (i = 0; i < subNb; i++) {
272:       const PetscReal *subq = &subsp->Q[i * origNb];

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

302: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Subspace(PetscSpace sp)
303: {
304:   PetscSpace_Subspace *subsp;

307:   PetscNewLog(sp,&subsp);
308:   sp->data = (void *) subsp;
309:   return(0);
310: }

312: static PetscErrorCode PetscSpaceGetDimension_Subspace(PetscSpace sp, PetscInt *dim)
313: {
314:   PetscSpace_Subspace *subsp;

317:   subsp = (PetscSpace_Subspace *) sp->data;
318:   *dim = subsp->Nb;
319:   return(0);
320: }

322: static PetscErrorCode PetscSpaceSetUp_Subspace(PetscSpace sp)
323: {
324:   const PetscReal     *x;
325:   const PetscReal     *Jx;
326:   const PetscReal     *u;
327:   const PetscReal     *Ju;
328:   PetscDualSpace      dualSubspace;
329:   PetscSpace          origSpace;
330:   PetscInt            origDim, subDim, origNc, subNc, origNb, subNb, f, i, j, numPoints, offset;
331:   PetscReal           *allPoints, *allWeights, *B, *V;
332:   DM                  dm;
333:   PetscSpace_Subspace *subsp;
334:   PetscErrorCode      ierr;

337:   subsp = (PetscSpace_Subspace *) sp->data;
338:   x            = subsp->x;
339:   Jx           = subsp->Jx;
340:   u            = subsp->u;
341:   Ju           = subsp->Ju;
342:   origSpace    = subsp->origSpace;
343:   dualSubspace = subsp->dualSubspace;
344:   PetscSpaceGetNumComponents(origSpace,&origNc);
345:   PetscSpaceGetNumVariables(origSpace,&origDim);
346:   PetscDualSpaceGetDM(dualSubspace,&dm);
347:   DMGetDimension(dm,&subDim);
348:   PetscSpaceGetDimension(origSpace,&origNb);
349:   PetscDualSpaceGetDimension(dualSubspace,&subNb);
350:   PetscDualSpaceGetNumComponents(dualSubspace,&subNc);

352:   for (f = 0, numPoints = 0; f < subNb; f++) {
353:     PetscQuadrature q;
354:     PetscInt        qNp;

356:     PetscDualSpaceGetFunctional(dualSubspace,f,&q);
357:     PetscQuadratureGetData(q,NULL,NULL,&qNp,NULL,NULL);
358:     numPoints += qNp;
359:   }
360:   PetscMalloc1(subNb*origNb,&V);
361:   PetscMalloc3(numPoints*origDim,&allPoints,numPoints*origNc,&allWeights,numPoints*origNb*origNc,&B);
362:   for (f = 0, offset = 0; f < subNb; f++) {
363:     PetscQuadrature q;
364:     PetscInt        qNp, p;
365:     const PetscReal *qp;
366:     const PetscReal *qw;

368:     PetscDualSpaceGetFunctional(dualSubspace,f,&q);
369:     PetscQuadratureGetData(q,NULL,NULL,&qNp,&qp,&qw);
370:     for (p = 0; p < qNp; p++, offset++) {
371:       if (x) {
372:         for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = x[i];
373:       } else {
374:         for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = 0.0;
375:       }
376:       if (Jx) {
377:         for (i = 0; i < origDim; i++) {
378:           for (j = 0; j < subDim; j++) {
379:             allPoints[origDim * offset + i] += Jx[i * subDim + j] * qp[j];
380:           }
381:         }
382:       } else {
383:         for (i = 0; i < PetscMin(subDim, origDim); i++) allPoints[origDim * offset + i] += qp[i];
384:       }
385:       for (i = 0; i < origNc; i++) allWeights[origNc * offset + i] = 0.0;
386:       if (Ju) {
387:         for (i = 0; i < origNc; i++) {
388:           for (j = 0; j < subNc; j++) {
389:             allWeights[offset * origNc + i] += qw[j] * Ju[j * origNc + i];
390:           }
391:         }
392:       } else {
393:         for (i = 0; i < PetscMin(subNc, origNc); i++) allWeights[offset * origNc + i] += qw[i];
394:       }
395:     }
396:   }
397:   PetscSpaceEvaluate(origSpace,numPoints,allPoints,B,NULL,NULL);
398:   for (f = 0, offset = 0; f < subNb; f++) {
399:     PetscInt b, p, s, qNp;
400:     PetscQuadrature q;
401:     const PetscReal *qw;

403:     PetscDualSpaceGetFunctional(dualSubspace,f,&q);
404:     PetscQuadratureGetData(q,NULL,NULL,&qNp,NULL,&qw);
405:     if (u) {
406:       for (b = 0; b < origNb; b++) {
407:         for (s = 0; s < subNc; s++) {
408:           V[f * origNb + b] += qw[s] * u[s];
409:         }
410:       }
411:     } else {
412:       for (b = 0; b < origNb; b++) V[f * origNb + b] = 0.0;
413:     }
414:     for (p = 0; p < qNp; p++, offset++) {
415:       for (b = 0; b < origNb; b++) {
416:         for (s = 0; s < origNc; s++) {
417:           V[f * origNb + b] += B[(offset * origNb + b) * origNc + s] * allWeights[offset * origNc + s];
418:         }
419:       }
420:     }
421:   }
422:   /* orthnormalize rows of V */
423:   for (f = 0; f < subNb; f++) {
424:     PetscReal rho = 0.0, scal;

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

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

430:     for (i = 0; i < origNb; i++) V[f * origNb + i] *= scal;
431:     for (j = f + 1; j < subNb; j++) {
432:       for (i = 0, scal = 0.; i < origNb; i++) scal += V[f * origNb + i] * V[j * origNb + i];
433:       for (i = 0; i < origNb; i++) V[j * origNb + i] -= V[f * origNb + i] * scal;
434:     }
435:   }
436:   PetscFree3(allPoints,allWeights,B);
437:   subsp->Q = V;
438:   return(0);
439: }

441: static PetscErrorCode PetscSpacePolynomialGetTensor_Subspace(PetscSpace sp, PetscBool *poly)
442: {
443:   PetscSpace_Subspace *subsp = (PetscSpace_Subspace *) sp->data;

447:   *poly = PETSC_FALSE;
448:   PetscSpacePolynomialGetTensor(subsp->origSpace,poly);
449:   if (*poly) {
450:     if (subsp->Jx) {
451:       PetscInt subDim, origDim, i, j;
452:       PetscInt maxnnz;

454:       PetscSpaceGetNumVariables(subsp->origSpace,&origDim);
455:       PetscSpaceGetNumVariables(sp,&subDim);
456:       maxnnz = 0;
457:       for (i = 0; i < origDim; i++) {
458:         PetscInt nnz = 0;

460:         for (j = 0; j < subDim; j++) nnz += (subsp->Jx[i * subDim + j] != 0.);
461:         maxnnz = PetscMax(maxnnz,nnz);
462:       }
463:       for (j = 0; j < subDim; j++) {
464:         PetscInt nnz = 0;

466:         for (i = 0; i < origDim; i++) nnz += (subsp->Jx[i * subDim + j] != 0.);
467:         maxnnz = PetscMax(maxnnz,nnz);
468:       }
469:       if (maxnnz > 1) *poly = PETSC_FALSE;
470:     }
471:   }
472:   return(0);
473: }

475: static PetscErrorCode PetscSpaceInitialize_Subspace(PetscSpace sp)
476: {

480:   sp->ops->setup = PetscSpaceSetUp_Subspace;
481:   sp->ops->view  = PetscSpaceView_Subspace;
482:   sp->ops->destroy  = PetscSpaceDestroy_Subspace;
483:   sp->ops->getdimension  = PetscSpaceGetDimension_Subspace;
484:   sp->ops->evaluate = PetscSpaceEvaluate_Subspace;
485:   PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Subspace);
486:   return(0);
487: }

489: PetscErrorCode PetscSpaceCreateSubspace(PetscSpace origSpace, PetscDualSpace dualSubspace, PetscReal *x, PetscReal *Jx, PetscReal *u, PetscReal *Ju, PetscCopyMode copymode, PetscSpace *subspace)
490: {
491:   PetscSpace_Subspace *subsp;
492:   PetscInt            origDim, subDim, origNc, subNc, subNb;
493:   PetscInt            order;
494:   DM                  dm;
495:   PetscErrorCode      ierr;

505:   PetscSpaceGetNumComponents(origSpace,&origNc);
506:   PetscSpaceGetNumVariables(origSpace,&origDim);
507:   PetscDualSpaceGetDM(dualSubspace,&dm);
508:   DMGetDimension(dm,&subDim);
509:   PetscDualSpaceGetDimension(dualSubspace,&subNb);
510:   PetscDualSpaceGetNumComponents(dualSubspace,&subNc);
511:   PetscSpaceCreate(PetscObjectComm((PetscObject)origSpace),subspace);
512:   PetscSpaceSetType(*subspace,PETSCSPACESUBSPACE);
513:   PetscSpaceSetNumVariables(*subspace,subDim);
514:   PetscSpaceSetNumComponents(*subspace,subNc);
515:   PetscSpaceGetDegree(origSpace,&order,NULL);
516:   PetscSpaceSetDegree(*subspace,order,PETSC_DETERMINE);
517:   subsp = (PetscSpace_Subspace *) (*subspace)->data;
518:   subsp->Nb = subNb;
519:   switch (copymode) {
520:   case PETSC_OWN_POINTER:
521:     if (x) subsp->x_alloc = x;
522:     if (Jx) subsp->Jx_alloc = Jx;
523:     if (u) subsp->u_alloc = u;
524:     if (Ju) subsp->Ju_alloc = Ju;
525:   case PETSC_USE_POINTER:
526:     if (x) subsp->x = x;
527:     if (Jx) subsp->Jx = Jx;
528:     if (u) subsp->u = u;
529:     if (Ju) subsp->Ju = Ju;
530:     break;
531:   case PETSC_COPY_VALUES:
532:     if (x) {
533:       PetscMalloc1(origDim,&subsp->x_alloc);
534:       PetscArraycpy(subsp->x_alloc,x,origDim);
535:       subsp->x = subsp->x_alloc;
536:     }
537:     if (Jx) {
538:       PetscMalloc1(origDim * subDim,&subsp->Jx_alloc);
539:       PetscArraycpy(subsp->Jx_alloc,Jx,origDim * subDim);
540:       subsp->Jx = subsp->Jx_alloc;
541:     }
542:     if (u) {
543:       PetscMalloc1(subNc,&subsp->u_alloc);
544:       PetscArraycpy(subsp->u_alloc,u,subNc);
545:       subsp->u = subsp->u_alloc;
546:     }
547:     if (Ju) {
548:       PetscMalloc1(origNc * subNc,&subsp->Ju_alloc);
549:       PetscArraycpy(subsp->Ju_alloc,Ju,origNc * subNc);
550:       subsp->Ju = subsp->Ju_alloc;
551:     }
552:     break;
553:   default:
554:     SETERRQ(PetscObjectComm((PetscObject)origSpace),PETSC_ERR_ARG_OUTOFRANGE,"Unknown copy mode");
555:   }
556:   PetscObjectReference((PetscObject)origSpace);
557:   subsp->origSpace = origSpace;
558:   PetscObjectReference((PetscObject)dualSubspace);
559:   subsp->dualSubspace = dualSubspace;
560:   PetscSpaceInitialize_Subspace(*subspace);
561:   return(0);
562: }