1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
2
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;
18
PetscSpaceDestroy_Subspace(PetscSpace sp)19 static PetscErrorCode PetscSpaceDestroy_Subspace(PetscSpace sp)
20 {
21 PetscSpace_Subspace *subsp;
22
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 }
41
PetscSpaceView_Subspace(PetscSpace sp,PetscViewer viewer)42 static PetscErrorCode PetscSpaceView_Subspace(PetscSpace sp, PetscViewer viewer)
43 {
44 PetscBool isascii;
45 PetscSpace_Subspace *subsp;
46
47 PetscFunctionBegin;
48 subsp = (PetscSpace_Subspace *)sp->data;
49 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
50 if (isascii) {
51 PetscInt origDim, subDim, origNc, subNc, o, s;
52
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 }
91
PetscSpaceEvaluate_Subspace(PetscSpace sp,PetscInt npoints,const PetscReal points[],PetscReal B[],PetscReal D[],PetscReal H[])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;
98
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;
128
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];
134
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;
177
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];
183
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;
220
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];
231
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 }
252
PetscSpaceCreate_Subspace(PetscSpace sp)253 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Subspace(PetscSpace sp)
254 {
255 PetscSpace_Subspace *subsp;
256
257 PetscFunctionBegin;
258 PetscCall(PetscNew(&subsp));
259 sp->data = (void *)subsp;
260 PetscFunctionReturn(PETSC_SUCCESS);
261 }
262
PetscSpaceGetDimension_Subspace(PetscSpace sp,PetscInt * dim)263 static PetscErrorCode PetscSpaceGetDimension_Subspace(PetscSpace sp, PetscInt *dim)
264 {
265 PetscSpace_Subspace *subsp;
266
267 PetscFunctionBegin;
268 subsp = (PetscSpace_Subspace *)sp->data;
269 *dim = subsp->Nb;
270 PetscFunctionReturn(PETSC_SUCCESS);
271 }
272
PetscSpaceSetUp_Subspace(PetscSpace sp)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;
285
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));
302
303 for (f = 0, numPoints = 0; f < subNb; f++) {
304 PetscQuadrature q;
305 PetscInt qNp;
306
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;
318
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;
349
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;
368
369 for (i = 0; i < origNb; i++) rho += PetscSqr(V[f * origNb + i]);
370
371 scal = 1. / PetscSqrtReal(rho);
372
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 }
383
PetscSpacePolynomialGetTensor_Subspace(PetscSpace sp,PetscBool * poly)384 static PetscErrorCode PetscSpacePolynomialGetTensor_Subspace(PetscSpace sp, PetscBool *poly)
385 {
386 PetscSpace_Subspace *subsp = (PetscSpace_Subspace *)sp->data;
387
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;
395
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;
401
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;
407
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 }
416
PetscSpaceInitialize_Subspace(PetscSpace sp)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`
430
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
439
440 Output Parameter:
441 . subspace - the subspace
442
443 Level: advanced
444
445 .seealso: `PetscSpace`, `PetscDualSpace`, `PetscCopyMode`, `PetscSpaceType`
446 @*/
PetscSpaceCreateSubspace(PetscSpace origSpace,PetscDualSpace dualSubspace,PetscReal * x,PetscReal * Jx,PetscReal * u,PetscReal * Ju,PetscCopyMode copymode,PetscSpace * subspace)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;
453
454 PetscFunctionBegin;
455 PetscValidHeaderSpecific(origSpace, PETSCSPACE_CLASSID, 1);
456 PetscValidHeaderSpecific(dualSubspace, PETSCDUALSPACE_CLASSID, 2);
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 }
521