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