Actual source code: spacesubspace.c
petsc-3.13.6 2020-09-29
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;
307: PetscNewLog(sp,&subsp);
308: sp->data = (void *) subsp;
309: return(0);
310: }
312: static PetscErrorCode PetscSpaceGetDimension_Subspace(PetscSpace sp, PetscInt *dim)
313: {
314: PetscSpace_Subspace *subsp;
317: subsp = (PetscSpace_Subspace *) sp->data;
318: *dim = subsp->Nb;
319: return(0);
320: }
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;
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: PetscSpaceGetNumComponents(origSpace,&origNc);
345: PetscSpaceGetNumVariables(origSpace,&origDim);
346: PetscDualSpaceGetDM(dualSubspace,&dm);
347: DMGetDimension(dm,&subDim);
348: PetscSpaceGetDimension(origSpace,&origNb);
349: PetscDualSpaceGetDimension(dualSubspace,&subNb);
350: PetscDualSpaceGetNumComponents(dualSubspace,&subNc);
352: for (f = 0, numPoints = 0; f < subNb; f++) {
353: PetscQuadrature q;
354: PetscInt qNp;
356: PetscDualSpaceGetFunctional(dualSubspace,f,&q);
357: PetscQuadratureGetData(q,NULL,NULL,&qNp,NULL,NULL);
358: numPoints += qNp;
359: }
360: PetscMalloc1(subNb*origNb,&V);
361: PetscMalloc3(numPoints*origDim,&allPoints,numPoints*origNc,&allWeights,numPoints*origNb*origNc,&B);
362: for (f = 0, offset = 0; f < subNb; f++) {
363: PetscQuadrature q;
364: PetscInt qNp, p;
365: const PetscReal *qp;
366: const PetscReal *qw;
368: PetscDualSpaceGetFunctional(dualSubspace,f,&q);
369: PetscQuadratureGetData(q,NULL,NULL,&qNp,&qp,&qw);
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: PetscSpaceEvaluate(origSpace,numPoints,allPoints,B,NULL,NULL);
398: for (f = 0, offset = 0; f < subNb; f++) {
399: PetscInt b, p, s, qNp;
400: PetscQuadrature q;
401: const PetscReal *qw;
403: PetscDualSpaceGetFunctional(dualSubspace,f,&q);
404: PetscQuadratureGetData(q,NULL,NULL,&qNp,NULL,&qw);
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;
426: for (i = 0; i < origNb; i++) rho += PetscSqr(V[f * origNb + i]);
428: scal = 1. / PetscSqrtReal(rho);
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: PetscFree3(allPoints,allWeights,B);
437: subsp->Q = V;
438: return(0);
439: }
441: static PetscErrorCode PetscSpacePolynomialGetTensor_Subspace(PetscSpace sp, PetscBool *poly)
442: {
443: PetscSpace_Subspace *subsp = (PetscSpace_Subspace *) sp->data;
447: *poly = PETSC_FALSE;
448: PetscSpacePolynomialGetTensor(subsp->origSpace,poly);
449: if (*poly) {
450: if (subsp->Jx) {
451: PetscInt subDim, origDim, i, j;
452: PetscInt maxnnz;
454: PetscSpaceGetNumVariables(subsp->origSpace,&origDim);
455: PetscSpaceGetNumVariables(sp,&subDim);
456: maxnnz = 0;
457: for (i = 0; i < origDim; i++) {
458: PetscInt nnz = 0;
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;
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: return(0);
473: }
475: static PetscErrorCode PetscSpaceInitialize_Subspace(PetscSpace sp)
476: {
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: PetscObjectComposeFunction((PetscObject) sp, "PetscSpacePolynomialGetTensor_C", PetscSpacePolynomialGetTensor_Subspace);
486: return(0);
487: }
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;
505: PetscSpaceGetNumComponents(origSpace,&origNc);
506: PetscSpaceGetNumVariables(origSpace,&origDim);
507: PetscDualSpaceGetDM(dualSubspace,&dm);
508: DMGetDimension(dm,&subDim);
509: PetscDualSpaceGetDimension(dualSubspace,&subNb);
510: PetscDualSpaceGetNumComponents(dualSubspace,&subNc);
511: PetscSpaceCreate(PetscObjectComm((PetscObject)origSpace),subspace);
512: PetscSpaceSetType(*subspace,PETSCSPACESUBSPACE);
513: PetscSpaceSetNumVariables(*subspace,subDim);
514: PetscSpaceSetNumComponents(*subspace,subNc);
515: PetscSpaceGetDegree(origSpace,&order,NULL);
516: PetscSpaceSetDegree(*subspace,order,PETSC_DETERMINE);
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: PetscMalloc1(origDim,&subsp->x_alloc);
534: PetscArraycpy(subsp->x_alloc,x,origDim);
535: subsp->x = subsp->x_alloc;
536: }
537: if (Jx) {
538: PetscMalloc1(origDim * subDim,&subsp->Jx_alloc);
539: PetscArraycpy(subsp->Jx_alloc,Jx,origDim * subDim);
540: subsp->Jx = subsp->Jx_alloc;
541: }
542: if (u) {
543: PetscMalloc1(subNc,&subsp->u_alloc);
544: PetscArraycpy(subsp->u_alloc,u,subNc);
545: subsp->u = subsp->u_alloc;
546: }
547: if (Ju) {
548: PetscMalloc1(origNc * subNc,&subsp->Ju_alloc);
549: PetscArraycpy(subsp->Ju_alloc,Ju,origNc * subNc);
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: PetscObjectReference((PetscObject)origSpace);
557: subsp->origSpace = origSpace;
558: PetscObjectReference((PetscObject)dualSubspace);
559: subsp->dualSubspace = dualSubspace;
560: PetscSpaceInitialize_Subspace(*subspace);
561: return(0);
562: }