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;
22: subsp = (PetscSpace_Subspace *) sp->data;
23: subsp->x = NULL;
24: PetscFree(subsp->x_alloc);
25: subsp->Jx = NULL;
26: PetscFree(subsp->Jx_alloc);
27: subsp->u = NULL;
28: PetscFree(subsp->u_alloc);
29: subsp->Ju = NULL;
30: PetscFree(subsp->Ju_alloc);
31: PetscFree(subsp->Q);
32: PetscSpaceDestroy(&subsp->origSpace);
33: PetscDualSpaceDestroy(&subsp->dualSubspace);
34: PetscFree(subsp);
35: sp->data = NULL;
36: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", NULL);
37: return 0;
38: }
40: static PetscErrorCode PetscSpaceView_Subspace(PetscSpace sp, PetscViewer viewer)
41: {
42: PetscBool iascii;
43: PetscSpace_Subspace *subsp;
45: subsp = (PetscSpace_Subspace *) sp->data;
46: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
47: if (iascii) {
48: PetscInt origDim, subDim, origNc, subNc, o, s;
50: PetscSpaceGetNumVariables(subsp->origSpace,&origDim);
51: PetscSpaceGetNumComponents(subsp->origSpace,&origNc);
52: PetscSpaceGetNumVariables(sp,&subDim);
53: PetscSpaceGetNumComponents(sp,&subNc);
54: if (subsp->x) {
55: PetscViewerASCIIPrintf(viewer,"Subspace-to-space domain shift:\n\n");
56: for (o = 0; o < origDim; o++) {
57: PetscViewerASCIIPrintf(viewer," %g\n", (double)subsp->x[o]);
58: }
59: }
60: if (subsp->Jx) {
61: PetscViewerASCIIPrintf(viewer,"Subspace-to-space domain transform:\n\n");
62: for (o = 0; o < origDim; o++) {
63: PetscViewerASCIIPrintf(viewer," %g", (double)subsp->Jx[o * subDim + 0]);
64: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
65: for (s = 1; s < subDim; s++) {
66: PetscViewerASCIIPrintf(viewer," %g", (double)subsp->Jx[o * subDim + s]);
67: }
68: PetscViewerASCIIPrintf(viewer,"\n");
69: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
70: }
71: }
72: if (subsp->u) {
73: PetscViewerASCIIPrintf(viewer,"Space-to-subspace range shift:\n\n");
74: for (o = 0; o < origNc; o++) {
75: PetscViewerASCIIPrintf(viewer," %d\n", subsp->u[o]);
76: }
77: }
78: if (subsp->Ju) {
79: PetscViewerASCIIPrintf(viewer,"Space-to-subsace domain transform:\n");
80: for (o = 0; o < origNc; o++) {
81: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
82: for (s = 0; s < subNc; s++) {
83: PetscViewerASCIIPrintf(viewer," %d", subsp->Ju[o * subNc + s]);
84: }
85: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
86: }
87: PetscViewerASCIIPrintf(viewer,"\n");
88: }
89: PetscViewerASCIIPrintf(viewer,"Original space:\n");
90: }
91: PetscViewerASCIIPushTab(viewer);
92: PetscSpaceView(subsp->origSpace,viewer);
93: PetscViewerASCIIPopTab(viewer);
94: return 0;
95: }
97: static PetscErrorCode PetscSpaceEvaluate_Subspace(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
98: {
99: PetscSpace_Subspace *subsp = (PetscSpace_Subspace *) sp->data;
100: PetscSpace origsp;
101: PetscInt origDim, subDim, origNc, subNc, subNb, origNb, i, j, k, l, m, n, o;
102: PetscReal *inpoints, *inB = NULL, *inD = NULL, *inH = NULL;
104: origsp = subsp->origSpace;
105: PetscSpaceGetNumVariables(sp,&subDim);
106: PetscSpaceGetNumVariables(origsp,&origDim);
107: PetscSpaceGetNumComponents(sp,&subNc);
108: PetscSpaceGetNumComponents(origsp,&origNc);
109: PetscSpaceGetDimension(sp,&subNb);
110: PetscSpaceGetDimension(origsp,&origNb);
111: DMGetWorkArray(sp->dm,npoints*origDim,MPIU_REAL,&inpoints);
112: for (i = 0; i < npoints; i++) {
113: if (subsp->x) {
114: for (j = 0; j < origDim; j++) inpoints[i * origDim + j] = subsp->x[j];
115: } else {
116: for (j = 0; j < origDim; j++) inpoints[i * origDim + j] = 0.0;
117: }
118: if (subsp->Jx) {
119: for (j = 0; j < origDim; j++) {
120: for (k = 0; k < subDim; k++) {
121: inpoints[i * origDim + j] += subsp->Jx[j * subDim + k] * points[i * subDim + k];
122: }
123: }
124: } else {
125: for (j = 0; j < PetscMin(subDim, origDim); j++) {
126: inpoints[i * origDim + j] += points[i * subDim + j];
127: }
128: }
129: }
130: if (B) {
131: DMGetWorkArray(sp->dm,npoints*origNb*origNc,MPIU_REAL,&inB);
132: }
133: if (D) {
134: DMGetWorkArray(sp->dm,npoints*origNb*origNc*origDim,MPIU_REAL,&inD);
135: }
136: if (H) {
137: DMGetWorkArray(sp->dm,npoints*origNb*origNc*origDim*origDim,MPIU_REAL,&inH);
138: }
139: PetscSpaceEvaluate(origsp,npoints,inpoints,inB,inD,inH);
140: if (H) {
141: PetscReal *phi, *psi;
143: DMGetWorkArray(sp->dm,origNc*origDim*origDim,MPIU_REAL,&phi);
144: DMGetWorkArray(sp->dm,origNc*subDim*subDim,MPIU_REAL,&psi);
145: for (i = 0; i < npoints * subNb * subNc * subDim; i++) D[i] = 0.0;
146: for (i = 0; i < subNb; i++) {
147: const PetscReal *subq = &subsp->Q[i * origNb];
149: for (j = 0; j < npoints; j++) {
150: for (k = 0; k < origNc * origDim; k++) phi[k] = 0.;
151: for (k = 0; k < origNc * subDim; k++) psi[k] = 0.;
152: for (k = 0; k < origNb; k++) {
153: for (l = 0; l < origNc * origDim * origDim; l++) {
154: phi[l] += inH[(j * origNb + k) * origNc * origDim * origDim + l] * subq[k];
155: }
156: }
157: if (subsp->Jx) {
158: for (k = 0; k < subNc; k++) {
159: for (l = 0; l < subDim; l++) {
160: for (m = 0; m < origDim; m++) {
161: for (n = 0; n < subDim; n++) {
162: for (o = 0; o < origDim; o++) {
163: psi[(k * subDim + l) * subDim + n] += subsp->Jx[m * subDim + l] * subsp->Jx[o * subDim + n] * phi[(k * origDim + m) * origDim + o];
164: }
165: }
166: }
167: }
168: }
169: } else {
170: for (k = 0; k < subNc; k++) {
171: for (l = 0; l < PetscMin(subDim, origDim); l++) {
172: for (m = 0; m < PetscMin(subDim, origDim); m++) {
173: psi[(k * subDim + l) * subDim + m] += phi[(k * origDim + l) * origDim + m];
174: }
175: }
176: }
177: }
178: if (subsp->Ju) {
179: for (k = 0; k < subNc; k++) {
180: for (l = 0; l < origNc; l++) {
181: for (m = 0; m < subDim * subDim; m++) {
182: H[((j * subNb + i) * subNc + k) * subDim * subDim + m] += subsp->Ju[k * origNc + l] * psi[l * subDim * subDim + m];
183: }
184: }
185: }
186: }
187: else {
188: for (k = 0; k < PetscMin(subNc, origNc); k++) {
189: for (l = 0; l < subDim * subDim; l++) {
190: H[((j * subNb + i) * subNc + k) * subDim * subDim + l] += psi[k * subDim * subDim + l];
191: }
192: }
193: }
194: }
195: }
196: DMRestoreWorkArray(sp->dm,subNc*origDim,MPIU_REAL,&psi);
197: DMRestoreWorkArray(sp->dm,origNc*origDim,MPIU_REAL,&phi);
198: DMRestoreWorkArray(sp->dm,npoints*origNb*origNc*origDim,MPIU_REAL,&inH);
199: }
200: if (D) {
201: PetscReal *phi, *psi;
203: DMGetWorkArray(sp->dm,origNc*origDim,MPIU_REAL,&phi);
204: DMGetWorkArray(sp->dm,origNc*subDim,MPIU_REAL,&psi);
205: for (i = 0; i < npoints * subNb * subNc * subDim; i++) D[i] = 0.0;
206: for (i = 0; i < subNb; i++) {
207: const PetscReal *subq = &subsp->Q[i * origNb];
209: for (j = 0; j < npoints; j++) {
210: for (k = 0; k < origNc * origDim; k++) phi[k] = 0.;
211: for (k = 0; k < origNc * subDim; k++) psi[k] = 0.;
212: for (k = 0; k < origNb; k++) {
213: for (l = 0; l < origNc * origDim; l++) {
214: phi[l] += inD[(j * origNb + k) * origNc * origDim + l] * subq[k];
215: }
216: }
217: if (subsp->Jx) {
218: for (k = 0; k < subNc; k++) {
219: for (l = 0; l < subDim; l++) {
220: for (m = 0; m < origDim; m++) {
221: psi[k * subDim + l] += subsp->Jx[m * subDim + l] * phi[k * origDim + m];
222: }
223: }
224: }
225: } else {
226: for (k = 0; k < subNc; k++) {
227: for (l = 0; l < PetscMin(subDim, origDim); l++) {
228: psi[k * subDim + l] += phi[k * origDim + l];
229: }
230: }
231: }
232: if (subsp->Ju) {
233: for (k = 0; k < subNc; k++) {
234: for (l = 0; l < origNc; l++) {
235: for (m = 0; m < subDim; m++) {
236: D[((j * subNb + i) * subNc + k) * subDim + m] += subsp->Ju[k * origNc + l] * psi[l * subDim + m];
237: }
238: }
239: }
240: }
241: else {
242: for (k = 0; k < PetscMin(subNc, origNc); k++) {
243: for (l = 0; l < subDim; l++) {
244: D[((j * subNb + i) * subNc + k) * subDim + l] += psi[k * subDim + l];
245: }
246: }
247: }
248: }
249: }
250: DMRestoreWorkArray(sp->dm,subNc*origDim,MPIU_REAL,&psi);
251: DMRestoreWorkArray(sp->dm,origNc*origDim,MPIU_REAL,&phi);
252: DMRestoreWorkArray(sp->dm,npoints*origNb*origNc*origDim,MPIU_REAL,&inD);
253: }
254: if (B) {
255: PetscReal *phi;
257: DMGetWorkArray(sp->dm,origNc,MPIU_REAL,&phi);
258: if (subsp->u) {
259: for (i = 0; i < npoints * subNb; i++) {
260: for (j = 0; j < subNc; j++) B[i * subNc + j] = subsp->u[j];
261: }
262: } else {
263: for (i = 0; i < npoints * subNb * subNc; i++) B[i] = 0.0;
264: }
265: for (i = 0; i < subNb; i++) {
266: const PetscReal *subq = &subsp->Q[i * origNb];
268: for (j = 0; j < npoints; j++) {
269: for (k = 0; k < origNc; k++) phi[k] = 0.;
270: for (k = 0; k < origNb; k++) {
271: for (l = 0; l < origNc; l++) {
272: phi[l] += inB[(j * origNb + k) * origNc + l] * subq[k];
273: }
274: }
275: if (subsp->Ju) {
276: for (k = 0; k < subNc; k++) {
277: for (l = 0; l < origNc; l++) {
278: B[(j * subNb + i) * subNc + k] += subsp->Ju[k * origNc + l] * phi[l];
279: }
280: }
281: }
282: else {
283: for (k = 0; k < PetscMin(subNc, origNc); k++) {
284: B[(j * subNb + i) * subNc + k] += phi[k];
285: }
286: }
287: }
288: }
289: DMRestoreWorkArray(sp->dm,origNc,MPIU_REAL,&phi);
290: DMRestoreWorkArray(sp->dm,npoints*origNb*origNc,MPIU_REAL,&inB);
291: }
292: DMRestoreWorkArray(sp->dm,npoints*origDim,MPIU_REAL,&inpoints);
293: return 0;
294: }
296: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Subspace(PetscSpace sp)
297: {
298: PetscSpace_Subspace *subsp;
300: PetscNewLog(sp,&subsp);
301: sp->data = (void *) subsp;
302: return 0;
303: }
305: static PetscErrorCode PetscSpaceGetDimension_Subspace(PetscSpace sp, PetscInt *dim)
306: {
307: PetscSpace_Subspace *subsp;
309: subsp = (PetscSpace_Subspace *) sp->data;
310: *dim = subsp->Nb;
311: return 0;
312: }
314: static PetscErrorCode PetscSpaceSetUp_Subspace(PetscSpace sp)
315: {
316: const PetscReal *x;
317: const PetscReal *Jx;
318: const PetscReal *u;
319: const PetscReal *Ju;
320: PetscDualSpace dualSubspace;
321: PetscSpace origSpace;
322: PetscInt origDim, subDim, origNc, subNc, origNb, subNb, f, i, j, numPoints, offset;
323: PetscReal *allPoints, *allWeights, *B, *V;
324: DM dm;
325: PetscSpace_Subspace *subsp;
327: subsp = (PetscSpace_Subspace *) sp->data;
328: x = subsp->x;
329: Jx = subsp->Jx;
330: u = subsp->u;
331: Ju = subsp->Ju;
332: origSpace = subsp->origSpace;
333: dualSubspace = subsp->dualSubspace;
334: PetscSpaceGetNumComponents(origSpace,&origNc);
335: PetscSpaceGetNumVariables(origSpace,&origDim);
336: PetscDualSpaceGetDM(dualSubspace,&dm);
337: DMGetDimension(dm,&subDim);
338: PetscSpaceGetDimension(origSpace,&origNb);
339: PetscDualSpaceGetDimension(dualSubspace,&subNb);
340: PetscDualSpaceGetNumComponents(dualSubspace,&subNc);
342: for (f = 0, numPoints = 0; f < subNb; f++) {
343: PetscQuadrature q;
344: PetscInt qNp;
346: PetscDualSpaceGetFunctional(dualSubspace,f,&q);
347: PetscQuadratureGetData(q,NULL,NULL,&qNp,NULL,NULL);
348: numPoints += qNp;
349: }
350: PetscMalloc1(subNb*origNb,&V);
351: PetscMalloc3(numPoints*origDim,&allPoints,numPoints*origNc,&allWeights,numPoints*origNb*origNc,&B);
352: for (f = 0, offset = 0; f < subNb; f++) {
353: PetscQuadrature q;
354: PetscInt qNp, p;
355: const PetscReal *qp;
356: const PetscReal *qw;
358: PetscDualSpaceGetFunctional(dualSubspace,f,&q);
359: PetscQuadratureGetData(q,NULL,NULL,&qNp,&qp,&qw);
360: for (p = 0; p < qNp; p++, offset++) {
361: if (x) {
362: for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = x[i];
363: } else {
364: for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = 0.0;
365: }
366: if (Jx) {
367: for (i = 0; i < origDim; i++) {
368: for (j = 0; j < subDim; j++) {
369: allPoints[origDim * offset + i] += Jx[i * subDim + j] * qp[j];
370: }
371: }
372: } else {
373: for (i = 0; i < PetscMin(subDim, origDim); i++) allPoints[origDim * offset + i] += qp[i];
374: }
375: for (i = 0; i < origNc; i++) allWeights[origNc * offset + i] = 0.0;
376: if (Ju) {
377: for (i = 0; i < origNc; i++) {
378: for (j = 0; j < subNc; j++) {
379: allWeights[offset * origNc + i] += qw[j] * Ju[j * origNc + i];
380: }
381: }
382: } else {
383: for (i = 0; i < PetscMin(subNc, origNc); i++) allWeights[offset * origNc + i] += qw[i];
384: }
385: }
386: }
387: PetscSpaceEvaluate(origSpace,numPoints,allPoints,B,NULL,NULL);
388: for (f = 0, offset = 0; f < subNb; f++) {
389: PetscInt b, p, s, qNp;
390: PetscQuadrature q;
391: const PetscReal *qw;
393: PetscDualSpaceGetFunctional(dualSubspace,f,&q);
394: PetscQuadratureGetData(q,NULL,NULL,&qNp,NULL,&qw);
395: if (u) {
396: for (b = 0; b < origNb; b++) {
397: for (s = 0; s < subNc; s++) {
398: V[f * origNb + b] += qw[s] * u[s];
399: }
400: }
401: } else {
402: for (b = 0; b < origNb; b++) V[f * origNb + b] = 0.0;
403: }
404: for (p = 0; p < qNp; p++, offset++) {
405: for (b = 0; b < origNb; b++) {
406: for (s = 0; s < origNc; s++) {
407: V[f * origNb + b] += B[(offset * origNb + b) * origNc + s] * allWeights[offset * origNc + s];
408: }
409: }
410: }
411: }
412: /* orthnormalize rows of V */
413: for (f = 0; f < subNb; f++) {
414: PetscReal rho = 0.0, scal;
416: for (i = 0; i < origNb; i++) rho += PetscSqr(V[f * origNb + i]);
418: scal = 1. / PetscSqrtReal(rho);
420: for (i = 0; i < origNb; i++) V[f * origNb + i] *= scal;
421: for (j = f + 1; j < subNb; j++) {
422: for (i = 0, scal = 0.; i < origNb; i++) scal += V[f * origNb + i] * V[j * origNb + i];
423: for (i = 0; i < origNb; i++) V[j * origNb + i] -= V[f * origNb + i] * scal;
424: }
425: }
426: PetscFree3(allPoints,allWeights,B);
427: subsp->Q = V;
428: return 0;
429: }
431: static PetscErrorCode PetscSpacePolynomialGetTensor_Subspace(PetscSpace sp, PetscBool *poly)
432: {
433: PetscSpace_Subspace *subsp = (PetscSpace_Subspace *) sp->data;
435: *poly = PETSC_FALSE;
436: PetscSpacePolynomialGetTensor(subsp->origSpace,poly);
437: if (*poly) {
438: if (subsp->Jx) {
439: PetscInt subDim, origDim, i, j;
440: PetscInt maxnnz;
442: PetscSpaceGetNumVariables(subsp->origSpace,&origDim);
443: PetscSpaceGetNumVariables(sp,&subDim);
444: maxnnz = 0;
445: for (i = 0; i < origDim; i++) {
446: PetscInt nnz = 0;
448: for (j = 0; j < subDim; j++) nnz += (subsp->Jx[i * subDim + j] != 0.);
449: maxnnz = PetscMax(maxnnz,nnz);
450: }
451: for (j = 0; j < subDim; j++) {
452: PetscInt nnz = 0;
454: for (i = 0; i < origDim; i++) nnz += (subsp->Jx[i * subDim + j] != 0.);
455: maxnnz = PetscMax(maxnnz,nnz);
456: }
457: if (maxnnz > 1) *poly = PETSC_FALSE;
458: }
459: }
460: return 0;
461: }
463: static PetscErrorCode PetscSpaceInitialize_Subspace(PetscSpace sp)
464: {
465: sp->ops->setup = PetscSpaceSetUp_Subspace;
466: sp->ops->view = PetscSpaceView_Subspace;
467: sp->ops->destroy = PetscSpaceDestroy_Subspace;
468: sp->ops->getdimension = PetscSpaceGetDimension_Subspace;
469: sp->ops->evaluate = PetscSpaceEvaluate_Subspace;
470: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Subspace);
471: return 0;
472: }
474: PetscErrorCode PetscSpaceCreateSubspace(PetscSpace origSpace, PetscDualSpace dualSubspace, PetscReal *x, PetscReal *Jx, PetscReal *u, PetscReal *Ju, PetscCopyMode copymode, PetscSpace *subspace)
475: {
476: PetscSpace_Subspace *subsp;
477: PetscInt origDim, subDim, origNc, subNc, subNb;
478: PetscInt order;
479: DM dm;
488: PetscSpaceGetNumComponents(origSpace,&origNc);
489: PetscSpaceGetNumVariables(origSpace,&origDim);
490: PetscDualSpaceGetDM(dualSubspace,&dm);
491: DMGetDimension(dm,&subDim);
492: PetscDualSpaceGetDimension(dualSubspace,&subNb);
493: PetscDualSpaceGetNumComponents(dualSubspace,&subNc);
494: PetscSpaceCreate(PetscObjectComm((PetscObject)origSpace),subspace);
495: PetscSpaceSetType(*subspace,PETSCSPACESUBSPACE);
496: PetscSpaceSetNumVariables(*subspace,subDim);
497: PetscSpaceSetNumComponents(*subspace,subNc);
498: PetscSpaceGetDegree(origSpace,&order,NULL);
499: PetscSpaceSetDegree(*subspace,order,PETSC_DETERMINE);
500: subsp = (PetscSpace_Subspace *) (*subspace)->data;
501: subsp->Nb = subNb;
502: switch (copymode) {
503: case PETSC_OWN_POINTER:
504: if (x) subsp->x_alloc = x;
505: if (Jx) subsp->Jx_alloc = Jx;
506: if (u) subsp->u_alloc = u;
507: if (Ju) subsp->Ju_alloc = Ju;
508: case PETSC_USE_POINTER:
509: if (x) subsp->x = x;
510: if (Jx) subsp->Jx = Jx;
511: if (u) subsp->u = u;
512: if (Ju) subsp->Ju = Ju;
513: break;
514: case PETSC_COPY_VALUES:
515: if (x) {
516: PetscMalloc1(origDim,&subsp->x_alloc);
517: PetscArraycpy(subsp->x_alloc,x,origDim);
518: subsp->x = subsp->x_alloc;
519: }
520: if (Jx) {
521: PetscMalloc1(origDim * subDim,&subsp->Jx_alloc);
522: PetscArraycpy(subsp->Jx_alloc,Jx,origDim * subDim);
523: subsp->Jx = subsp->Jx_alloc;
524: }
525: if (u) {
526: PetscMalloc1(subNc,&subsp->u_alloc);
527: PetscArraycpy(subsp->u_alloc,u,subNc);
528: subsp->u = subsp->u_alloc;
529: }
530: if (Ju) {
531: PetscMalloc1(origNc * subNc,&subsp->Ju_alloc);
532: PetscArraycpy(subsp->Ju_alloc,Ju,origNc * subNc);
533: subsp->Ju = subsp->Ju_alloc;
534: }
535: break;
536: default:
537: SETERRQ(PetscObjectComm((PetscObject)origSpace),PETSC_ERR_ARG_OUTOFRANGE,"Unknown copy mode");
538: }
539: PetscObjectReference((PetscObject)origSpace);
540: subsp->origSpace = origSpace;
541: PetscObjectReference((PetscObject)dualSubspace);
542: subsp->dualSubspace = dualSubspace;
543: PetscSpaceInitialize_Subspace(*subspace);
544: return 0;
545: }