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++) PetscViewerASCIIPrintf(viewer, " %g\n", (double)subsp->x[o]);
57: }
58: if (subsp->Jx) {
59: PetscViewerASCIIPrintf(viewer, "Subspace-to-space domain transform:\n\n");
60: for (o = 0; o < origDim; o++) {
61: PetscViewerASCIIPrintf(viewer, " %g", (double)subsp->Jx[o * subDim + 0]);
62: PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
63: for (s = 1; s < subDim; s++) PetscViewerASCIIPrintf(viewer, " %g", (double)subsp->Jx[o * subDim + s]);
64: PetscViewerASCIIPrintf(viewer, "\n");
65: PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
66: }
67: }
68: if (subsp->u) {
69: PetscViewerASCIIPrintf(viewer, "Space-to-subspace range shift:\n\n");
70: for (o = 0; o < origNc; o++) PetscViewerASCIIPrintf(viewer, " %g\n", (double)subsp->u[o]);
71: }
72: if (subsp->Ju) {
73: PetscViewerASCIIPrintf(viewer, "Space-to-subsace domain transform:\n");
74: for (o = 0; o < origNc; o++) {
75: PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
76: for (s = 0; s < subNc; s++) PetscViewerASCIIPrintf(viewer, " %g", (double)subsp->Ju[o * subNc + s]);
77: PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
78: }
79: PetscViewerASCIIPrintf(viewer, "\n");
80: }
81: PetscViewerASCIIPrintf(viewer, "Original space:\n");
82: }
83: PetscViewerASCIIPushTab(viewer);
84: PetscSpaceView(subsp->origSpace, viewer);
85: PetscViewerASCIIPopTab(viewer);
86: return 0;
87: }
89: static PetscErrorCode PetscSpaceEvaluate_Subspace(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
90: {
91: PetscSpace_Subspace *subsp = (PetscSpace_Subspace *)sp->data;
92: PetscSpace origsp;
93: PetscInt origDim, subDim, origNc, subNc, subNb, origNb, i, j, k, l, m, n, o;
94: PetscReal *inpoints, *inB = NULL, *inD = NULL, *inH = NULL;
96: origsp = subsp->origSpace;
97: PetscSpaceGetNumVariables(sp, &subDim);
98: PetscSpaceGetNumVariables(origsp, &origDim);
99: PetscSpaceGetNumComponents(sp, &subNc);
100: PetscSpaceGetNumComponents(origsp, &origNc);
101: PetscSpaceGetDimension(sp, &subNb);
102: PetscSpaceGetDimension(origsp, &origNb);
103: DMGetWorkArray(sp->dm, npoints * origDim, MPIU_REAL, &inpoints);
104: for (i = 0; i < npoints; i++) {
105: if (subsp->x) {
106: for (j = 0; j < origDim; j++) inpoints[i * origDim + j] = subsp->x[j];
107: } else {
108: for (j = 0; j < origDim; j++) inpoints[i * origDim + j] = 0.0;
109: }
110: if (subsp->Jx) {
111: for (j = 0; j < origDim; j++) {
112: for (k = 0; k < subDim; k++) inpoints[i * origDim + j] += subsp->Jx[j * subDim + k] * points[i * subDim + k];
113: }
114: } else {
115: for (j = 0; j < PetscMin(subDim, origDim); j++) inpoints[i * origDim + j] += points[i * subDim + j];
116: }
117: }
118: if (B) DMGetWorkArray(sp->dm, npoints * origNb * origNc, MPIU_REAL, &inB);
119: if (D) DMGetWorkArray(sp->dm, npoints * origNb * origNc * origDim, MPIU_REAL, &inD);
120: if (H) DMGetWorkArray(sp->dm, npoints * origNb * origNc * origDim * origDim, MPIU_REAL, &inH);
121: PetscSpaceEvaluate(origsp, npoints, inpoints, inB, inD, inH);
122: if (H) {
123: PetscReal *phi, *psi;
125: DMGetWorkArray(sp->dm, origNc * origDim * origDim, MPIU_REAL, &phi);
126: DMGetWorkArray(sp->dm, origNc * subDim * subDim, MPIU_REAL, &psi);
127: for (i = 0; i < npoints * subNb * subNc * subDim; i++) D[i] = 0.0;
128: for (i = 0; i < subNb; i++) {
129: const PetscReal *subq = &subsp->Q[i * origNb];
131: for (j = 0; j < npoints; j++) {
132: for (k = 0; k < origNc * origDim; k++) phi[k] = 0.;
133: for (k = 0; k < origNc * subDim; k++) psi[k] = 0.;
134: for (k = 0; k < origNb; k++) {
135: for (l = 0; l < origNc * origDim * origDim; l++) phi[l] += inH[(j * origNb + k) * origNc * origDim * origDim + l] * subq[k];
136: }
137: if (subsp->Jx) {
138: for (k = 0; k < subNc; k++) {
139: for (l = 0; l < subDim; l++) {
140: for (m = 0; m < origDim; m++) {
141: for (n = 0; n < subDim; n++) {
142: 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];
143: }
144: }
145: }
146: }
147: } else {
148: for (k = 0; k < subNc; k++) {
149: for (l = 0; l < PetscMin(subDim, origDim); l++) {
150: for (m = 0; m < PetscMin(subDim, origDim); m++) psi[(k * subDim + l) * subDim + m] += phi[(k * origDim + l) * origDim + m];
151: }
152: }
153: }
154: if (subsp->Ju) {
155: for (k = 0; k < subNc; k++) {
156: for (l = 0; l < origNc; l++) {
157: 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];
158: }
159: }
160: } else {
161: for (k = 0; k < PetscMin(subNc, origNc); k++) {
162: for (l = 0; l < subDim * subDim; l++) H[((j * subNb + i) * subNc + k) * subDim * subDim + l] += psi[k * subDim * subDim + l];
163: }
164: }
165: }
166: }
167: DMRestoreWorkArray(sp->dm, subNc * origDim, MPIU_REAL, &psi);
168: DMRestoreWorkArray(sp->dm, origNc * origDim, MPIU_REAL, &phi);
169: DMRestoreWorkArray(sp->dm, npoints * origNb * origNc * origDim, MPIU_REAL, &inH);
170: }
171: if (D) {
172: PetscReal *phi, *psi;
174: DMGetWorkArray(sp->dm, origNc * origDim, MPIU_REAL, &phi);
175: DMGetWorkArray(sp->dm, origNc * subDim, MPIU_REAL, &psi);
176: for (i = 0; i < npoints * subNb * subNc * subDim; i++) D[i] = 0.0;
177: for (i = 0; i < subNb; i++) {
178: const PetscReal *subq = &subsp->Q[i * origNb];
180: for (j = 0; j < npoints; j++) {
181: for (k = 0; k < origNc * origDim; k++) phi[k] = 0.;
182: for (k = 0; k < origNc * subDim; k++) psi[k] = 0.;
183: for (k = 0; k < origNb; k++) {
184: for (l = 0; l < origNc * origDim; l++) phi[l] += inD[(j * origNb + k) * origNc * origDim + l] * subq[k];
185: }
186: if (subsp->Jx) {
187: for (k = 0; k < subNc; k++) {
188: for (l = 0; l < subDim; l++) {
189: for (m = 0; m < origDim; m++) psi[k * subDim + l] += subsp->Jx[m * subDim + l] * phi[k * origDim + m];
190: }
191: }
192: } else {
193: for (k = 0; k < subNc; k++) {
194: for (l = 0; l < PetscMin(subDim, origDim); l++) psi[k * subDim + l] += phi[k * origDim + l];
195: }
196: }
197: if (subsp->Ju) {
198: for (k = 0; k < subNc; k++) {
199: for (l = 0; l < origNc; l++) {
200: for (m = 0; m < subDim; m++) D[((j * subNb + i) * subNc + k) * subDim + m] += subsp->Ju[k * origNc + l] * psi[l * subDim + m];
201: }
202: }
203: } else {
204: for (k = 0; k < PetscMin(subNc, origNc); k++) {
205: for (l = 0; l < subDim; l++) D[((j * subNb + i) * subNc + k) * subDim + l] += psi[k * subDim + l];
206: }
207: }
208: }
209: }
210: DMRestoreWorkArray(sp->dm, subNc * origDim, MPIU_REAL, &psi);
211: DMRestoreWorkArray(sp->dm, origNc * origDim, MPIU_REAL, &phi);
212: DMRestoreWorkArray(sp->dm, npoints * origNb * origNc * origDim, MPIU_REAL, &inD);
213: }
214: if (B) {
215: PetscReal *phi;
217: DMGetWorkArray(sp->dm, origNc, MPIU_REAL, &phi);
218: if (subsp->u) {
219: for (i = 0; i < npoints * subNb; i++) {
220: for (j = 0; j < subNc; j++) B[i * subNc + j] = subsp->u[j];
221: }
222: } else {
223: for (i = 0; i < npoints * subNb * subNc; i++) B[i] = 0.0;
224: }
225: for (i = 0; i < subNb; i++) {
226: const PetscReal *subq = &subsp->Q[i * origNb];
228: for (j = 0; j < npoints; j++) {
229: for (k = 0; k < origNc; k++) phi[k] = 0.;
230: for (k = 0; k < origNb; k++) {
231: for (l = 0; l < origNc; l++) phi[l] += inB[(j * origNb + k) * origNc + l] * subq[k];
232: }
233: if (subsp->Ju) {
234: for (k = 0; k < subNc; k++) {
235: for (l = 0; l < origNc; l++) B[(j * subNb + i) * subNc + k] += subsp->Ju[k * origNc + l] * phi[l];
236: }
237: } else {
238: for (k = 0; k < PetscMin(subNc, origNc); k++) B[(j * subNb + i) * subNc + k] += phi[k];
239: }
240: }
241: }
242: DMRestoreWorkArray(sp->dm, origNc, MPIU_REAL, &phi);
243: DMRestoreWorkArray(sp->dm, npoints * origNb * origNc, MPIU_REAL, &inB);
244: }
245: DMRestoreWorkArray(sp->dm, npoints * origDim, MPIU_REAL, &inpoints);
246: return 0;
247: }
249: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Subspace(PetscSpace sp)
250: {
251: PetscSpace_Subspace *subsp;
253: PetscNew(&subsp);
254: sp->data = (void *)subsp;
255: return 0;
256: }
258: static PetscErrorCode PetscSpaceGetDimension_Subspace(PetscSpace sp, PetscInt *dim)
259: {
260: PetscSpace_Subspace *subsp;
262: subsp = (PetscSpace_Subspace *)sp->data;
263: *dim = subsp->Nb;
264: return 0;
265: }
267: static PetscErrorCode PetscSpaceSetUp_Subspace(PetscSpace sp)
268: {
269: const PetscReal *x;
270: const PetscReal *Jx;
271: const PetscReal *u;
272: const PetscReal *Ju;
273: PetscDualSpace dualSubspace;
274: PetscSpace origSpace;
275: PetscInt origDim, subDim, origNc, subNc, origNb, subNb, f, i, j, numPoints, offset;
276: PetscReal *allPoints, *allWeights, *B, *V;
277: DM dm;
278: PetscSpace_Subspace *subsp;
280: subsp = (PetscSpace_Subspace *)sp->data;
281: x = subsp->x;
282: Jx = subsp->Jx;
283: u = subsp->u;
284: Ju = subsp->Ju;
285: origSpace = subsp->origSpace;
286: dualSubspace = subsp->dualSubspace;
287: PetscSpaceGetNumComponents(origSpace, &origNc);
288: PetscSpaceGetNumVariables(origSpace, &origDim);
289: PetscDualSpaceGetDM(dualSubspace, &dm);
290: DMGetDimension(dm, &subDim);
291: PetscSpaceGetDimension(origSpace, &origNb);
292: PetscDualSpaceGetDimension(dualSubspace, &subNb);
293: PetscDualSpaceGetNumComponents(dualSubspace, &subNc);
295: for (f = 0, numPoints = 0; f < subNb; f++) {
296: PetscQuadrature q;
297: PetscInt qNp;
299: PetscDualSpaceGetFunctional(dualSubspace, f, &q);
300: PetscQuadratureGetData(q, NULL, NULL, &qNp, NULL, NULL);
301: numPoints += qNp;
302: }
303: PetscMalloc1(subNb * origNb, &V);
304: PetscMalloc3(numPoints * origDim, &allPoints, numPoints * origNc, &allWeights, numPoints * origNb * origNc, &B);
305: for (f = 0, offset = 0; f < subNb; f++) {
306: PetscQuadrature q;
307: PetscInt qNp, p;
308: const PetscReal *qp;
309: const PetscReal *qw;
311: PetscDualSpaceGetFunctional(dualSubspace, f, &q);
312: PetscQuadratureGetData(q, NULL, NULL, &qNp, &qp, &qw);
313: for (p = 0; p < qNp; p++, offset++) {
314: if (x) {
315: for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = x[i];
316: } else {
317: for (i = 0; i < origDim; i++) allPoints[origDim * offset + i] = 0.0;
318: }
319: if (Jx) {
320: for (i = 0; i < origDim; i++) {
321: for (j = 0; j < subDim; j++) allPoints[origDim * offset + i] += Jx[i * subDim + j] * qp[j];
322: }
323: } else {
324: for (i = 0; i < PetscMin(subDim, origDim); i++) allPoints[origDim * offset + i] += qp[i];
325: }
326: for (i = 0; i < origNc; i++) allWeights[origNc * offset + i] = 0.0;
327: if (Ju) {
328: for (i = 0; i < origNc; i++) {
329: for (j = 0; j < subNc; j++) allWeights[offset * origNc + i] += qw[j] * Ju[j * origNc + i];
330: }
331: } else {
332: for (i = 0; i < PetscMin(subNc, origNc); i++) allWeights[offset * origNc + i] += qw[i];
333: }
334: }
335: }
336: PetscSpaceEvaluate(origSpace, numPoints, allPoints, B, NULL, NULL);
337: for (f = 0, offset = 0; f < subNb; f++) {
338: PetscInt b, p, s, qNp;
339: PetscQuadrature q;
340: const PetscReal *qw;
342: PetscDualSpaceGetFunctional(dualSubspace, f, &q);
343: PetscQuadratureGetData(q, NULL, NULL, &qNp, NULL, &qw);
344: if (u) {
345: for (b = 0; b < origNb; b++) {
346: for (s = 0; s < subNc; s++) V[f * origNb + b] += qw[s] * u[s];
347: }
348: } else {
349: for (b = 0; b < origNb; b++) V[f * origNb + b] = 0.0;
350: }
351: for (p = 0; p < qNp; p++, offset++) {
352: for (b = 0; b < origNb; b++) {
353: for (s = 0; s < origNc; s++) V[f * origNb + b] += B[(offset * origNb + b) * origNc + s] * allWeights[offset * origNc + s];
354: }
355: }
356: }
357: /* orthnormalize rows of V */
358: for (f = 0; f < subNb; f++) {
359: PetscReal rho = 0.0, scal;
361: for (i = 0; i < origNb; i++) rho += PetscSqr(V[f * origNb + i]);
363: scal = 1. / PetscSqrtReal(rho);
365: for (i = 0; i < origNb; i++) V[f * origNb + i] *= scal;
366: for (j = f + 1; j < subNb; j++) {
367: for (i = 0, scal = 0.; i < origNb; i++) scal += V[f * origNb + i] * V[j * origNb + i];
368: for (i = 0; i < origNb; i++) V[j * origNb + i] -= V[f * origNb + i] * scal;
369: }
370: }
371: PetscFree3(allPoints, allWeights, B);
372: subsp->Q = V;
373: return 0;
374: }
376: static PetscErrorCode PetscSpacePolynomialGetTensor_Subspace(PetscSpace sp, PetscBool *poly)
377: {
378: PetscSpace_Subspace *subsp = (PetscSpace_Subspace *)sp->data;
380: *poly = PETSC_FALSE;
381: PetscSpacePolynomialGetTensor(subsp->origSpace, poly);
382: if (*poly) {
383: if (subsp->Jx) {
384: PetscInt subDim, origDim, i, j;
385: PetscInt maxnnz;
387: PetscSpaceGetNumVariables(subsp->origSpace, &origDim);
388: PetscSpaceGetNumVariables(sp, &subDim);
389: maxnnz = 0;
390: for (i = 0; i < origDim; i++) {
391: PetscInt nnz = 0;
393: for (j = 0; j < subDim; j++) nnz += (subsp->Jx[i * subDim + j] != 0.);
394: maxnnz = PetscMax(maxnnz, nnz);
395: }
396: for (j = 0; j < subDim; j++) {
397: PetscInt nnz = 0;
399: for (i = 0; i < origDim; i++) nnz += (subsp->Jx[i * subDim + j] != 0.);
400: maxnnz = PetscMax(maxnnz, nnz);
401: }
402: if (maxnnz > 1) *poly = PETSC_FALSE;
403: }
404: }
405: return 0;
406: }
408: static PetscErrorCode PetscSpaceInitialize_Subspace(PetscSpace sp)
409: {
410: sp->ops->setup = PetscSpaceSetUp_Subspace;
411: sp->ops->view = PetscSpaceView_Subspace;
412: sp->ops->destroy = PetscSpaceDestroy_Subspace;
413: sp->ops->getdimension = PetscSpaceGetDimension_Subspace;
414: sp->ops->evaluate = PetscSpaceEvaluate_Subspace;
415: PetscObjectComposeFunction((PetscObject)sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Subspace);
416: return 0;
417: }
419: PetscErrorCode PetscSpaceCreateSubspace(PetscSpace origSpace, PetscDualSpace dualSubspace, PetscReal *x, PetscReal *Jx, PetscReal *u, PetscReal *Ju, PetscCopyMode copymode, PetscSpace *subspace)
420: {
421: PetscSpace_Subspace *subsp;
422: PetscInt origDim, subDim, origNc, subNc, subNb;
423: PetscInt order;
424: DM dm;
433: PetscSpaceGetNumComponents(origSpace, &origNc);
434: PetscSpaceGetNumVariables(origSpace, &origDim);
435: PetscDualSpaceGetDM(dualSubspace, &dm);
436: DMGetDimension(dm, &subDim);
437: PetscDualSpaceGetDimension(dualSubspace, &subNb);
438: PetscDualSpaceGetNumComponents(dualSubspace, &subNc);
439: PetscSpaceCreate(PetscObjectComm((PetscObject)origSpace), subspace);
440: PetscSpaceSetType(*subspace, PETSCSPACESUBSPACE);
441: PetscSpaceSetNumVariables(*subspace, subDim);
442: PetscSpaceSetNumComponents(*subspace, subNc);
443: PetscSpaceGetDegree(origSpace, &order, NULL);
444: PetscSpaceSetDegree(*subspace, order, PETSC_DETERMINE);
445: subsp = (PetscSpace_Subspace *)(*subspace)->data;
446: subsp->Nb = subNb;
447: switch (copymode) {
448: case PETSC_OWN_POINTER:
449: if (x) subsp->x_alloc = x;
450: if (Jx) subsp->Jx_alloc = Jx;
451: if (u) subsp->u_alloc = u;
452: if (Ju) subsp->Ju_alloc = Ju;
453: /* fall through */
454: case PETSC_USE_POINTER:
455: if (x) subsp->x = x;
456: if (Jx) subsp->Jx = Jx;
457: if (u) subsp->u = u;
458: if (Ju) subsp->Ju = Ju;
459: break;
460: case PETSC_COPY_VALUES:
461: if (x) {
462: PetscMalloc1(origDim, &subsp->x_alloc);
463: PetscArraycpy(subsp->x_alloc, x, origDim);
464: subsp->x = subsp->x_alloc;
465: }
466: if (Jx) {
467: PetscMalloc1(origDim * subDim, &subsp->Jx_alloc);
468: PetscArraycpy(subsp->Jx_alloc, Jx, origDim * subDim);
469: subsp->Jx = subsp->Jx_alloc;
470: }
471: if (u) {
472: PetscMalloc1(subNc, &subsp->u_alloc);
473: PetscArraycpy(subsp->u_alloc, u, subNc);
474: subsp->u = subsp->u_alloc;
475: }
476: if (Ju) {
477: PetscMalloc1(origNc * subNc, &subsp->Ju_alloc);
478: PetscArraycpy(subsp->Ju_alloc, Ju, origNc * subNc);
479: subsp->Ju = subsp->Ju_alloc;
480: }
481: break;
482: default:
483: SETERRQ(PetscObjectComm((PetscObject)origSpace), PETSC_ERR_ARG_OUTOFRANGE, "Unknown copy mode");
484: }
485: PetscObjectReference((PetscObject)origSpace);
486: subsp->origSpace = origSpace;
487: PetscObjectReference((PetscObject)dualSubspace);
488: subsp->dualSubspace = dualSubspace;
489: PetscSpaceInitialize_Subspace(*subspace);
490: return 0;
491: }