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