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