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 } PetscSpace_Subspace;
17
PetscSpaceDestroy_Subspace(PetscSpace sp)18 static PetscErrorCode PetscSpaceDestroy_Subspace(PetscSpace sp)
19 {
20 PetscSpace_Subspace *subsp;
21 PetscErrorCode ierr;
22
23 PetscFunctionBegin;
24 subsp = (PetscSpace_Subspace *) sp->data;
25 subsp->x = NULL;
26 ierr = PetscFree(subsp->x_alloc);CHKERRQ(ierr);
27 subsp->Jx = NULL;
28 ierr = PetscFree(subsp->Jx_alloc);CHKERRQ(ierr);
29 subsp->u = NULL;
30 ierr = PetscFree(subsp->u_alloc);CHKERRQ(ierr);
31 subsp->Ju = NULL;
32 ierr = PetscFree(subsp->Ju_alloc);CHKERRQ(ierr);
33 ierr = PetscFree(subsp->Q);CHKERRQ(ierr);
34 ierr = PetscSpaceDestroy(&subsp->origSpace);CHKERRQ(ierr);
35 ierr = PetscDualSpaceDestroy(&subsp->dualSubspace);CHKERRQ(ierr);
36 ierr = PetscFree(subsp);CHKERRQ(ierr);
37 sp->data = NULL;
38 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", NULL);CHKERRQ(ierr);
39 PetscFunctionReturn(0);
40 }
41
PetscSpaceView_Subspace(PetscSpace sp,PetscViewer viewer)42 static PetscErrorCode PetscSpaceView_Subspace(PetscSpace sp, PetscViewer viewer)
43 {
44 PetscBool iascii;
45 PetscSpace_Subspace *subsp;
46 PetscErrorCode ierr;
47
48 PetscFunctionBegin;
49 subsp = (PetscSpace_Subspace *) sp->data;
50 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
51 if (iascii) {
52 PetscInt origDim, subDim, origNc, subNc, o, s;
53
54 ierr = PetscSpaceGetNumVariables(subsp->origSpace,&origDim);CHKERRQ(ierr);
55 ierr = PetscSpaceGetNumComponents(subsp->origSpace,&origNc);CHKERRQ(ierr);
56 ierr = PetscSpaceGetNumVariables(sp,&subDim);CHKERRQ(ierr);
57 ierr = PetscSpaceGetNumComponents(sp,&subNc);CHKERRQ(ierr);
58 if (subsp->x) {
59 ierr = PetscViewerASCIIPrintf(viewer,"Subspace-to-space domain shift:\n\n");CHKERRQ(ierr);
60 for (o = 0; o < origDim; o++) {
61 ierr = PetscViewerASCIIPrintf(viewer," %g\n", (double)subsp->x[o]);CHKERRQ(ierr);
62 }
63 }
64 if (subsp->Jx) {
65 ierr = PetscViewerASCIIPrintf(viewer,"Subspace-to-space domain transform:\n\n");CHKERRQ(ierr);
66 for (o = 0; o < origDim; o++) {
67 ierr = PetscViewerASCIIPrintf(viewer," %g", (double)subsp->Jx[o * subDim + 0]);CHKERRQ(ierr);
68 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);
69 for (s = 1; s < subDim; s++) {
70 ierr = PetscViewerASCIIPrintf(viewer," %g", (double)subsp->Jx[o * subDim + s]);CHKERRQ(ierr);
71 }
72 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
73 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);
74 }
75 }
76 if (subsp->u) {
77 ierr = PetscViewerASCIIPrintf(viewer,"Space-to-subspace range shift:\n\n");CHKERRQ(ierr);
78 for (o = 0; o < origNc; o++) {
79 ierr = PetscViewerASCIIPrintf(viewer," %d\n", subsp->u[o]);CHKERRQ(ierr);
80 }
81 }
82 if (subsp->Ju) {
83 ierr = PetscViewerASCIIPrintf(viewer,"Space-to-subsace domain transform:\n");CHKERRQ(ierr);
84 for (o = 0; o < origNc; o++) {
85 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr);
86 for (s = 0; s < subNc; s++) {
87 ierr = PetscViewerASCIIPrintf(viewer," %d", subsp->Ju[o * subNc + s]);CHKERRQ(ierr);
88 }
89 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr);
90 }
91 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr);
92 }
93 ierr = PetscViewerASCIIPrintf(viewer,"Original space:\n");CHKERRQ(ierr);
94 }
95 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
96 ierr = PetscSpaceView(subsp->origSpace,viewer);CHKERRQ(ierr);
97 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
98 PetscFunctionReturn(0);
99 }
100
PetscSpaceEvaluate_Subspace(PetscSpace sp,PetscInt npoints,const PetscReal points[],PetscReal B[],PetscReal D[],PetscReal H[])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;
108
109 PetscFunctionBegin;
110 origsp = subsp->origSpace;
111 ierr = PetscSpaceGetNumVariables(sp,&subDim);CHKERRQ(ierr);
112 ierr = PetscSpaceGetNumVariables(origsp,&origDim);CHKERRQ(ierr);
113 ierr = PetscSpaceGetNumComponents(sp,&subNc);CHKERRQ(ierr);
114 ierr = PetscSpaceGetNumComponents(origsp,&origNc);CHKERRQ(ierr);
115 ierr = PetscSpaceGetDimension(sp,&subNb);CHKERRQ(ierr);
116 ierr = PetscSpaceGetDimension(origsp,&origNb);CHKERRQ(ierr);
117 ierr = DMGetWorkArray(sp->dm,npoints*origDim,MPIU_REAL,&inpoints);CHKERRQ(ierr);
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 ierr = DMGetWorkArray(sp->dm,npoints*origNb*origNc,MPIU_REAL,&inB);CHKERRQ(ierr);
138 }
139 if (D) {
140 ierr = DMGetWorkArray(sp->dm,npoints*origNb*origNc*origDim,MPIU_REAL,&inD);CHKERRQ(ierr);
141 }
142 if (H) {
143 ierr = DMGetWorkArray(sp->dm,npoints*origNb*origNc*origDim*origDim,MPIU_REAL,&inH);CHKERRQ(ierr);
144 }
145 ierr = PetscSpaceEvaluate(origsp,npoints,inpoints,inB,inD,inH);CHKERRQ(ierr);
146 if (H) {
147 PetscReal *phi, *psi;
148
149 ierr = DMGetWorkArray(sp->dm,origNc*origDim*origDim,MPIU_REAL,&phi);CHKERRQ(ierr);
150 ierr = DMGetWorkArray(sp->dm,origNc*subDim*subDim,MPIU_REAL,&psi);CHKERRQ(ierr);
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];
154
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 ierr = DMRestoreWorkArray(sp->dm,subNc*origDim,MPIU_REAL,&psi);CHKERRQ(ierr);
203 ierr = DMRestoreWorkArray(sp->dm,origNc*origDim,MPIU_REAL,&phi);CHKERRQ(ierr);
204 ierr = DMRestoreWorkArray(sp->dm,npoints*origNb*origNc*origDim,MPIU_REAL,&inH);CHKERRQ(ierr);
205 }
206 if (D) {
207 PetscReal *phi, *psi;
208
209 ierr = DMGetWorkArray(sp->dm,origNc*origDim,MPIU_REAL,&phi);CHKERRQ(ierr);
210 ierr = DMGetWorkArray(sp->dm,origNc*subDim,MPIU_REAL,&psi);CHKERRQ(ierr);
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];
214
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 ierr = DMRestoreWorkArray(sp->dm,subNc*origDim,MPIU_REAL,&psi);CHKERRQ(ierr);
257 ierr = DMRestoreWorkArray(sp->dm,origNc*origDim,MPIU_REAL,&phi);CHKERRQ(ierr);
258 ierr = DMRestoreWorkArray(sp->dm,npoints*origNb*origNc*origDim,MPIU_REAL,&inD);CHKERRQ(ierr);
259 }
260 if (B) {
261 PetscReal *phi;
262
263 ierr = DMGetWorkArray(sp->dm,origNc,MPIU_REAL,&phi);CHKERRQ(ierr);
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];
273
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 ierr = DMRestoreWorkArray(sp->dm,origNc,MPIU_REAL,&phi);CHKERRQ(ierr);
296 ierr = DMRestoreWorkArray(sp->dm,npoints*origNb*origNc,MPIU_REAL,&inB);CHKERRQ(ierr);
297 }
298 ierr = DMRestoreWorkArray(sp->dm,npoints*origDim,MPIU_REAL,&inpoints);CHKERRQ(ierr);
299 PetscFunctionReturn(0);
300 }
301
PetscSpaceCreate_Subspace(PetscSpace sp)302 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Subspace(PetscSpace sp)
303 {
304 PetscSpace_Subspace *subsp;
305
306 PetscErrorCode ierr;
307 ierr = PetscNewLog(sp,&subsp);CHKERRQ(ierr);
308 sp->data = (void *) subsp;
309 PetscFunctionReturn(0);
310 }
311
PetscSpaceGetDimension_Subspace(PetscSpace sp,PetscInt * dim)312 static PetscErrorCode PetscSpaceGetDimension_Subspace(PetscSpace sp, PetscInt *dim)
313 {
314 PetscSpace_Subspace *subsp;
315
316 PetscFunctionBegin;
317 subsp = (PetscSpace_Subspace *) sp->data;
318 *dim = subsp->Nb;
319 PetscFunctionReturn(0);
320 }
321
PetscSpaceSetUp_Subspace(PetscSpace sp)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;
335
336 PetscFunctionBegin;
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 ierr = PetscSpaceGetNumComponents(origSpace,&origNc);CHKERRQ(ierr);
345 ierr = PetscSpaceGetNumVariables(origSpace,&origDim);CHKERRQ(ierr);
346 ierr = PetscDualSpaceGetDM(dualSubspace,&dm);CHKERRQ(ierr);
347 ierr = DMGetDimension(dm,&subDim);CHKERRQ(ierr);
348 ierr = PetscSpaceGetDimension(origSpace,&origNb);CHKERRQ(ierr);
349 ierr = PetscDualSpaceGetDimension(dualSubspace,&subNb);CHKERRQ(ierr);
350 ierr = PetscDualSpaceGetNumComponents(dualSubspace,&subNc);CHKERRQ(ierr);
351
352 for (f = 0, numPoints = 0; f < subNb; f++) {
353 PetscQuadrature q;
354 PetscInt qNp;
355
356 ierr = PetscDualSpaceGetFunctional(dualSubspace,f,&q);CHKERRQ(ierr);
357 ierr = PetscQuadratureGetData(q,NULL,NULL,&qNp,NULL,NULL);CHKERRQ(ierr);
358 numPoints += qNp;
359 }
360 ierr = PetscMalloc1(subNb*origNb,&V);CHKERRQ(ierr);
361 ierr = PetscMalloc3(numPoints*origDim,&allPoints,numPoints*origNc,&allWeights,numPoints*origNb*origNc,&B);CHKERRQ(ierr);
362 for (f = 0, offset = 0; f < subNb; f++) {
363 PetscQuadrature q;
364 PetscInt qNp, p;
365 const PetscReal *qp;
366 const PetscReal *qw;
367
368 ierr = PetscDualSpaceGetFunctional(dualSubspace,f,&q);CHKERRQ(ierr);
369 ierr = PetscQuadratureGetData(q,NULL,NULL,&qNp,&qp,&qw);CHKERRQ(ierr);
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 ierr = PetscSpaceEvaluate(origSpace,numPoints,allPoints,B,NULL,NULL);CHKERRQ(ierr);
398 for (f = 0, offset = 0; f < subNb; f++) {
399 PetscInt b, p, s, qNp;
400 PetscQuadrature q;
401 const PetscReal *qw;
402
403 ierr = PetscDualSpaceGetFunctional(dualSubspace,f,&q);CHKERRQ(ierr);
404 ierr = PetscQuadratureGetData(q,NULL,NULL,&qNp,NULL,&qw);CHKERRQ(ierr);
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;
425
426 for (i = 0; i < origNb; i++) rho += PetscSqr(V[f * origNb + i]);
427
428 scal = 1. / PetscSqrtReal(rho);
429
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 ierr = PetscFree3(allPoints,allWeights,B);CHKERRQ(ierr);
437 subsp->Q = V;
438 PetscFunctionReturn(0);
439 }
440
PetscSpacePolynomialGetTensor_Subspace(PetscSpace sp,PetscBool * poly)441 static PetscErrorCode PetscSpacePolynomialGetTensor_Subspace(PetscSpace sp, PetscBool *poly)
442 {
443 PetscSpace_Subspace *subsp = (PetscSpace_Subspace *) sp->data;
444 PetscErrorCode ierr;
445
446 PetscFunctionBegin;
447 *poly = PETSC_FALSE;
448 ierr = PetscSpacePolynomialGetTensor(subsp->origSpace,poly);CHKERRQ(ierr);
449 if (*poly) {
450 if (subsp->Jx) {
451 PetscInt subDim, origDim, i, j;
452 PetscInt maxnnz;
453
454 ierr = PetscSpaceGetNumVariables(subsp->origSpace,&origDim);CHKERRQ(ierr);
455 ierr = PetscSpaceGetNumVariables(sp,&subDim);CHKERRQ(ierr);
456 maxnnz = 0;
457 for (i = 0; i < origDim; i++) {
458 PetscInt nnz = 0;
459
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;
465
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 PetscFunctionReturn(0);
473 }
474
PetscSpaceInitialize_Subspace(PetscSpace sp)475 static PetscErrorCode PetscSpaceInitialize_Subspace(PetscSpace sp)
476 {
477 PetscErrorCode ierr;
478
479 PetscFunctionBegin;
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 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Subspace);CHKERRQ(ierr);
486 PetscFunctionReturn(0);
487 }
488
PetscSpaceCreateSubspace(PetscSpace origSpace,PetscDualSpace dualSubspace,PetscReal * x,PetscReal * Jx,PetscReal * u,PetscReal * Ju,PetscCopyMode copymode,PetscSpace * subspace)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;
496
497 PetscFunctionBegin;
498 PetscValidHeaderSpecific(origSpace,PETSCSPACE_CLASSID,1);
499 PetscValidHeaderSpecific(dualSubspace,PETSCDUALSPACE_CLASSID,2);
500 if (x) PetscValidRealPointer(x,3);
501 if (Jx) PetscValidRealPointer(Jx,4);
502 if (u) PetscValidRealPointer(u,5);
503 if (Ju) PetscValidRealPointer(Ju,6);
504 PetscValidPointer(subspace,7);
505 ierr = PetscSpaceGetNumComponents(origSpace,&origNc);CHKERRQ(ierr);
506 ierr = PetscSpaceGetNumVariables(origSpace,&origDim);CHKERRQ(ierr);
507 ierr = PetscDualSpaceGetDM(dualSubspace,&dm);CHKERRQ(ierr);
508 ierr = DMGetDimension(dm,&subDim);CHKERRQ(ierr);
509 ierr = PetscDualSpaceGetDimension(dualSubspace,&subNb);CHKERRQ(ierr);
510 ierr = PetscDualSpaceGetNumComponents(dualSubspace,&subNc);CHKERRQ(ierr);
511 ierr = PetscSpaceCreate(PetscObjectComm((PetscObject)origSpace),subspace);CHKERRQ(ierr);
512 ierr = PetscSpaceSetType(*subspace,PETSCSPACESUBSPACE);CHKERRQ(ierr);
513 ierr = PetscSpaceSetNumVariables(*subspace,subDim);CHKERRQ(ierr);
514 ierr = PetscSpaceSetNumComponents(*subspace,subNc);CHKERRQ(ierr);
515 ierr = PetscSpaceGetDegree(origSpace,&order,NULL);CHKERRQ(ierr);
516 ierr = PetscSpaceSetDegree(*subspace,order,PETSC_DETERMINE);CHKERRQ(ierr);
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 ierr = PetscMalloc1(origDim,&subsp->x_alloc);CHKERRQ(ierr);
534 ierr = PetscArraycpy(subsp->x_alloc,x,origDim);CHKERRQ(ierr);
535 subsp->x = subsp->x_alloc;
536 }
537 if (Jx) {
538 ierr = PetscMalloc1(origDim * subDim,&subsp->Jx_alloc);CHKERRQ(ierr);
539 ierr = PetscArraycpy(subsp->Jx_alloc,Jx,origDim * subDim);CHKERRQ(ierr);
540 subsp->Jx = subsp->Jx_alloc;
541 }
542 if (u) {
543 ierr = PetscMalloc1(subNc,&subsp->u_alloc);CHKERRQ(ierr);
544 ierr = PetscArraycpy(subsp->u_alloc,u,subNc);CHKERRQ(ierr);
545 subsp->u = subsp->u_alloc;
546 }
547 if (Ju) {
548 ierr = PetscMalloc1(origNc * subNc,&subsp->Ju_alloc);CHKERRQ(ierr);
549 ierr = PetscArraycpy(subsp->Ju_alloc,Ju,origNc * subNc);CHKERRQ(ierr);
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 ierr = PetscObjectReference((PetscObject)origSpace);CHKERRQ(ierr);
557 subsp->origSpace = origSpace;
558 ierr = PetscObjectReference((PetscObject)dualSubspace);CHKERRQ(ierr);
559 subsp->dualSubspace = dualSubspace;
560 ierr = PetscSpaceInitialize_Subspace(*subspace);CHKERRQ(ierr);
561 PetscFunctionReturn(0);
562 }
563
564