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;
21: PetscErrorCode ierr;
24: subsp = (PetscSpace_Subspace *) sp->data;
25: subsp->x = NULL;
26: PetscFree(subsp->x_alloc);
27: subsp->Jx = NULL;
28: PetscFree(subsp->Jx_alloc);
29: subsp->u = NULL;
30: PetscFree(subsp->u_alloc);
31: subsp->Ju = NULL;
32: PetscFree(subsp->Ju_alloc);
33: PetscFree(subsp->Q);
34: PetscSpaceDestroy(&subsp->origSpace);
35: PetscDualSpaceDestroy(&subsp->dualSubspace);
36: PetscFree(subsp);
37: sp->data = NULL;
38: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", NULL);
39: return(0);
40: }
42: static PetscErrorCode PetscSpaceView_Subspace(PetscSpace sp, PetscViewer viewer)
43: {
44: PetscBool iascii;
45: PetscSpace_Subspace *subsp;
46: PetscErrorCode ierr;
49: subsp = (PetscSpace_Subspace *) sp->data;
50: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
51: if (iascii) {
52: PetscInt origDim, subDim, origNc, subNc, o, s;
54: PetscSpaceGetNumVariables(subsp->origSpace,&origDim);
55: PetscSpaceGetNumComponents(subsp->origSpace,&origNc);
56: PetscSpaceGetNumVariables(sp,&subDim);
57: PetscSpaceGetNumComponents(sp,&subNc);
58: if (subsp->x) {
59: PetscViewerASCIIPrintf(viewer,"Subspace-to-space domain shift:\n\n");
60: for (o = 0; o < origDim; o++) {
61: PetscViewerASCIIPrintf(viewer," %g\n", (double)subsp->x[o]);
62: }
63: }
64: if (subsp->Jx) {
65: PetscViewerASCIIPrintf(viewer,"Subspace-to-space domain transform:\n\n");
66: for (o = 0; o < origDim; o++) {
67: PetscViewerASCIIPrintf(viewer," %g", (double)subsp->Jx[o * subDim + 0]);
68: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
69: for (s = 1; s < subDim; s++) {
70: PetscViewerASCIIPrintf(viewer," %g", (double)subsp->Jx[o * subDim + s]);
71: }
72: PetscViewerASCIIPrintf(viewer,"\n");
73: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
74: }
75: }
76: if (subsp->u) {
77: PetscViewerASCIIPrintf(viewer,"Space-to-subspace range shift:\n\n");
78: for (o = 0; o < origNc; o++) {
79: PetscViewerASCIIPrintf(viewer," %d\n", subsp->u[o]);
80: }
81: }
82: if (subsp->Ju) {
83: PetscViewerASCIIPrintf(viewer,"Space-to-subsace domain transform:\n");
84: for (o = 0; o < origNc; o++) {
85: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
86: for (s = 0; s < subNc; s++) {
87: PetscViewerASCIIPrintf(viewer," %d", subsp->Ju[o * subNc + s]);
88: }
89: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
90: }
91: PetscViewerASCIIPrintf(viewer,"\n");
92: }
93: PetscViewerASCIIPrintf(viewer,"Original space:\n");
94: }
95: PetscViewerASCIIPushTab(viewer);
96: PetscSpaceView(subsp->origSpace,viewer);
97: PetscViewerASCIIPopTab(viewer);
98: return(0);
99: }
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;
110: origsp = subsp->origSpace;
111: PetscSpaceGetNumVariables(sp,&subDim);
112: PetscSpaceGetNumVariables(origsp,&origDim);
113: PetscSpaceGetNumComponents(sp,&subNc);
114: PetscSpaceGetNumComponents(origsp,&origNc);
115: PetscSpaceGetDimension(sp,&subNb);
116: PetscSpaceGetDimension(origsp,&origNb);
117: DMGetWorkArray(sp->dm,npoints*origDim,MPIU_REAL,&inpoints);
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: DMGetWorkArray(sp->dm,npoints*origNb*origNc,MPIU_REAL,&inB);
138: }
139: if (D) {
140: DMGetWorkArray(sp->dm,npoints*origNb*origNc*origDim,MPIU_REAL,&inD);
141: }
142: if (H) {
143: DMGetWorkArray(sp->dm,npoints*origNb*origNc*origDim*origDim,MPIU_REAL,&inH);
144: }
145: PetscSpaceEvaluate(origsp,npoints,inpoints,inB,inD,inH);
146: if (H) {
147: PetscReal *phi, *psi;
149: DMGetWorkArray(sp->dm,origNc*origDim*origDim,MPIU_REAL,&phi);
150: DMGetWorkArray(sp->dm,origNc*subDim*subDim,MPIU_REAL,&psi);
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];
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: DMRestoreWorkArray(sp->dm,subNc*origDim,MPIU_REAL,&psi);
203: DMRestoreWorkArray(sp->dm,origNc*origDim,MPIU_REAL,&phi);
204: DMRestoreWorkArray(sp->dm,npoints*origNb*origNc*origDim,MPIU_REAL,&inH);
205: }
206: if (D) {
207: PetscReal *phi, *psi;
209: DMGetWorkArray(sp->dm,origNc*origDim,MPIU_REAL,&phi);
210: DMGetWorkArray(sp->dm,origNc*subDim,MPIU_REAL,&psi);
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];
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: DMRestoreWorkArray(sp->dm,subNc*origDim,MPIU_REAL,&psi);
257: DMRestoreWorkArray(sp->dm,origNc*origDim,MPIU_REAL,&phi);
258: DMRestoreWorkArray(sp->dm,npoints*origNb*origNc*origDim,MPIU_REAL,&inD);
259: }
260: if (B) {
261: PetscReal *phi;
263: DMGetWorkArray(sp->dm,origNc,MPIU_REAL,&phi);
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];
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: DMRestoreWorkArray(sp->dm,origNc,MPIU_REAL,&phi);
296: DMRestoreWorkArray(sp->dm,npoints*origNb*origNc,MPIU_REAL,&inB);
297: }
298: DMRestoreWorkArray(sp->dm,npoints*origDim,MPIU_REAL,&inpoints);
299: return(0);
300: }
302: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Subspace(PetscSpace sp)
303: {
304: PetscSpace_Subspace *subsp;
308: PetscNewLog(sp,&subsp);
309: sp->data = (void *) subsp;
310: return(0);
311: }
313: static PetscErrorCode PetscSpaceGetDimension_Subspace(PetscSpace sp, PetscInt *dim)
314: {
315: PetscSpace_Subspace *subsp;
318: subsp = (PetscSpace_Subspace *) sp->data;
319: *dim = subsp->Nb;
320: return(0);
321: }
323: static PetscErrorCode PetscSpaceSetUp_Subspace(PetscSpace sp)
324: {
325: const PetscReal *x;
326: const PetscReal *Jx;
327: const PetscReal *u;
328: const PetscReal *Ju;
329: PetscDualSpace dualSubspace;
330: PetscSpace origSpace;
331: PetscInt origDim, subDim, origNc, subNc, origNb, subNb, f, i, j, numPoints, offset;
332: PetscReal *allPoints, *allWeights, *B, *V;
333: DM dm;
334: PetscSpace_Subspace *subsp;
335: PetscErrorCode ierr;
338: subsp = (PetscSpace_Subspace *) sp->data;
339: x = subsp->x;
340: Jx = subsp->Jx;
341: u = subsp->u;
342: Ju = subsp->Ju;
343: origSpace = subsp->origSpace;
344: dualSubspace = subsp->dualSubspace;
345: PetscSpaceGetNumComponents(origSpace,&origNc);
346: PetscSpaceGetNumVariables(origSpace,&origDim);
347: PetscDualSpaceGetDM(dualSubspace,&dm);
348: DMGetDimension(dm,&subDim);
349: PetscSpaceGetDimension(origSpace,&origNb);
350: PetscDualSpaceGetDimension(dualSubspace,&subNb);
351: PetscDualSpaceGetNumComponents(dualSubspace,&subNc);
353: for (f = 0, numPoints = 0; f < subNb; f++) {
354: PetscQuadrature q;
355: PetscInt qNp;
357: PetscDualSpaceGetFunctional(dualSubspace,f,&q);
358: PetscQuadratureGetData(q,NULL,NULL,&qNp,NULL,NULL);
359: numPoints += qNp;
360: }
361: PetscMalloc1(subNb*origNb,&V);
362: PetscMalloc3(numPoints*origDim,&allPoints,numPoints*origNc,&allWeights,numPoints*origNb*origNc,&B);
363: for (f = 0, offset = 0; f < subNb; f++) {
364: PetscQuadrature q;
365: PetscInt qNp, p;
366: const PetscReal *qp;
367: const PetscReal *qw;
369: PetscDualSpaceGetFunctional(dualSubspace,f,&q);
370: PetscQuadratureGetData(q,NULL,NULL,&qNp,&qp,&qw);
371: for (p = 0; p < qNp; p++, offset++) {
372: if (x) {
373: for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = x[i];
374: } else {
375: for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = 0.0;
376: }
377: if (Jx) {
378: for (i = 0; i < origDim; i++) {
379: for (j = 0; j < subDim; j++) {
380: allPoints[origDim * offset + i] += Jx[i * subDim + j] * qp[j];
381: }
382: }
383: } else {
384: for (i = 0; i < PetscMin(subDim, origDim); i++) allPoints[origDim * offset + i] += qp[i];
385: }
386: for (i = 0; i < origNc; i++) allWeights[origNc * offset + i] = 0.0;
387: if (Ju) {
388: for (i = 0; i < origNc; i++) {
389: for (j = 0; j < subNc; j++) {
390: allWeights[offset * origNc + i] += qw[j] * Ju[j * origNc + i];
391: }
392: }
393: } else {
394: for (i = 0; i < PetscMin(subNc, origNc); i++) allWeights[offset * origNc + i] += qw[i];
395: }
396: }
397: }
398: PetscSpaceEvaluate(origSpace,numPoints,allPoints,B,NULL,NULL);
399: for (f = 0, offset = 0; f < subNb; f++) {
400: PetscInt b, p, s, qNp;
401: PetscQuadrature q;
402: const PetscReal *qw;
404: PetscDualSpaceGetFunctional(dualSubspace,f,&q);
405: PetscQuadratureGetData(q,NULL,NULL,&qNp,NULL,&qw);
406: if (u) {
407: for (b = 0; b < origNb; b++) {
408: for (s = 0; s < subNc; s++) {
409: V[f * origNb + b] += qw[s] * u[s];
410: }
411: }
412: } else {
413: for (b = 0; b < origNb; b++) V[f * origNb + b] = 0.0;
414: }
415: for (p = 0; p < qNp; p++, offset++) {
416: for (b = 0; b < origNb; b++) {
417: for (s = 0; s < origNc; s++) {
418: V[f * origNb + b] += B[(offset * origNb + b) * origNc + s] * allWeights[offset * origNc + s];
419: }
420: }
421: }
422: }
423: /* orthnormalize rows of V */
424: for (f = 0; f < subNb; f++) {
425: PetscReal rho = 0.0, scal;
427: for (i = 0; i < origNb; i++) rho += PetscSqr(V[f * origNb + i]);
429: scal = 1. / PetscSqrtReal(rho);
431: for (i = 0; i < origNb; i++) V[f * origNb + i] *= scal;
432: for (j = f + 1; j < subNb; j++) {
433: for (i = 0, scal = 0.; i < origNb; i++) scal += V[f * origNb + i] * V[j * origNb + i];
434: for (i = 0; i < origNb; i++) V[j * origNb + i] -= V[f * origNb + i] * scal;
435: }
436: }
437: PetscFree3(allPoints,allWeights,B);
438: subsp->Q = V;
439: return(0);
440: }
442: static PetscErrorCode PetscSpacePolynomialGetTensor_Subspace(PetscSpace sp, PetscBool *poly)
443: {
444: PetscSpace_Subspace *subsp = (PetscSpace_Subspace *) sp->data;
448: *poly = PETSC_FALSE;
449: PetscSpacePolynomialGetTensor(subsp->origSpace,poly);
450: if (*poly) {
451: if (subsp->Jx) {
452: PetscInt subDim, origDim, i, j;
453: PetscInt maxnnz;
455: PetscSpaceGetNumVariables(subsp->origSpace,&origDim);
456: PetscSpaceGetNumVariables(sp,&subDim);
457: maxnnz = 0;
458: for (i = 0; i < origDim; i++) {
459: PetscInt nnz = 0;
461: for (j = 0; j < subDim; j++) nnz += (subsp->Jx[i * subDim + j] != 0.);
462: maxnnz = PetscMax(maxnnz,nnz);
463: }
464: for (j = 0; j < subDim; j++) {
465: PetscInt nnz = 0;
467: for (i = 0; i < origDim; i++) nnz += (subsp->Jx[i * subDim + j] != 0.);
468: maxnnz = PetscMax(maxnnz,nnz);
469: }
470: if (maxnnz > 1) *poly = PETSC_FALSE;
471: }
472: }
473: return(0);
474: }
476: static PetscErrorCode PetscSpaceInitialize_Subspace(PetscSpace sp)
477: {
481: sp->ops->setup = PetscSpaceSetUp_Subspace;
482: sp->ops->view = PetscSpaceView_Subspace;
483: sp->ops->destroy = PetscSpaceDestroy_Subspace;
484: sp->ops->getdimension = PetscSpaceGetDimension_Subspace;
485: sp->ops->evaluate = PetscSpaceEvaluate_Subspace;
486: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Subspace);
487: return(0);
488: }
490: PetscErrorCode PetscSpaceCreateSubspace(PetscSpace origSpace, PetscDualSpace dualSubspace, PetscReal *x, PetscReal *Jx, PetscReal *u, PetscReal *Ju, PetscCopyMode copymode, PetscSpace *subspace)
491: {
492: PetscSpace_Subspace *subsp;
493: PetscInt origDim, subDim, origNc, subNc, subNb;
494: PetscInt order;
495: DM dm;
496: PetscErrorCode ierr;
506: PetscSpaceGetNumComponents(origSpace,&origNc);
507: PetscSpaceGetNumVariables(origSpace,&origDim);
508: PetscDualSpaceGetDM(dualSubspace,&dm);
509: DMGetDimension(dm,&subDim);
510: PetscDualSpaceGetDimension(dualSubspace,&subNb);
511: PetscDualSpaceGetNumComponents(dualSubspace,&subNc);
512: PetscSpaceCreate(PetscObjectComm((PetscObject)origSpace),subspace);
513: PetscSpaceSetType(*subspace,PETSCSPACESUBSPACE);
514: PetscSpaceSetNumVariables(*subspace,subDim);
515: PetscSpaceSetNumComponents(*subspace,subNc);
516: PetscSpaceGetDegree(origSpace,&order,NULL);
517: PetscSpaceSetDegree(*subspace,order,PETSC_DETERMINE);
518: subsp = (PetscSpace_Subspace *) (*subspace)->data;
519: subsp->Nb = subNb;
520: switch (copymode) {
521: case PETSC_OWN_POINTER:
522: if (x) subsp->x_alloc = x;
523: if (Jx) subsp->Jx_alloc = Jx;
524: if (u) subsp->u_alloc = u;
525: if (Ju) subsp->Ju_alloc = Ju;
526: case PETSC_USE_POINTER:
527: if (x) subsp->x = x;
528: if (Jx) subsp->Jx = Jx;
529: if (u) subsp->u = u;
530: if (Ju) subsp->Ju = Ju;
531: break;
532: case PETSC_COPY_VALUES:
533: if (x) {
534: PetscMalloc1(origDim,&subsp->x_alloc);
535: PetscArraycpy(subsp->x_alloc,x,origDim);
536: subsp->x = subsp->x_alloc;
537: }
538: if (Jx) {
539: PetscMalloc1(origDim * subDim,&subsp->Jx_alloc);
540: PetscArraycpy(subsp->Jx_alloc,Jx,origDim * subDim);
541: subsp->Jx = subsp->Jx_alloc;
542: }
543: if (u) {
544: PetscMalloc1(subNc,&subsp->u_alloc);
545: PetscArraycpy(subsp->u_alloc,u,subNc);
546: subsp->u = subsp->u_alloc;
547: }
548: if (Ju) {
549: PetscMalloc1(origNc * subNc,&subsp->Ju_alloc);
550: PetscArraycpy(subsp->Ju_alloc,Ju,origNc * subNc);
551: subsp->Ju = subsp->Ju_alloc;
552: }
553: break;
554: default:
555: SETERRQ(PetscObjectComm((PetscObject)origSpace),PETSC_ERR_ARG_OUTOFRANGE,"Unknown copy mode");
556: }
557: PetscObjectReference((PetscObject)origSpace);
558: subsp->origSpace = origSpace;
559: PetscObjectReference((PetscObject)dualSubspace);
560: subsp->dualSubspace = dualSubspace;
561: PetscSpaceInitialize_Subspace(*subspace);
562: return(0);
563: }