Actual source code: spacesum.c
1: #include <petsc/private/petscfeimpl.h>
2: /*@
3: PetscSpaceSumGetNumSubspaces - Get the number of spaces in the sum
5: Input Parameter:
6: . sp - the function space object
8: Output Parameter:
9: . numSumSpaces - the number of spaces
11: Level: intermediate
13: .seealso: PetscSpaceSumSetNumSubspaces(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
14: @*/
15: PetscErrorCode PetscSpaceSumGetNumSubspaces(PetscSpace sp,PetscInt *numSumSpaces)
16: {
19: PetscTryMethod(sp,"PetscSpaceSumGetNumSubspaces_C",(PetscSpace,PetscInt*),(sp,numSumSpaces));
20: return 0;
21: }
23: /*@
24: PetscSpaceSumSetNumSubspaces - Set the number of spaces in the sum
26: Input Parameters:
27: + sp - the function space object
28: - numSumSpaces - the number of spaces
30: Level: intermediate
32: .seealso: PetscSpaceSumGetNumSubspaces(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
33: @*/
34: PetscErrorCode PetscSpaceSumSetNumSubspaces(PetscSpace sp,PetscInt numSumSpaces)
35: {
37: PetscTryMethod(sp,"PetscSpaceSumSetNumSubspaces_C",(PetscSpace,PetscInt),(sp,numSumSpaces));
38: return 0;
39: }
41: /*@
42: PetscSpaceSumGetConcatenate - Get the concatenate flag for this space.
43: A concatenated sum space will have number of components equal to the sum of the number of components of all subspaces.A non-concatenated,
44: or direct sum space will have the same number of components as its subspaces .
46: Input Parameters:
47: . sp - the function space object
49: Output Parameters:
50: . concatenate - flag indicating whether subspaces are concatenated.
52: Level: intermediate
54: .seealso: PetscSpaceSumSetConcatenate()
55: @*/
56: PetscErrorCode PetscSpaceSumGetConcatenate(PetscSpace sp,PetscBool *concatenate)
57: {
59: PetscTryMethod(sp,"PetscSpaceSumGetConcatenate_C",(PetscSpace,PetscBool*),(sp,concatenate));
60: return 0;
61: }
63: /*@
64: PetscSpaceSumSetConcatenate - Sets the concatenate flag for this space.
65: A concatenated sum space will have number of components equal to the sum of the number of components of all subspaces.A non-concatenated,
66: or direct sum space will have the same number of components as its subspaces .
68: Input Parameters:
69: + sp - the function space object
70: - concatenate - are subspaces concatenated components (true) or direct summands (false)
72: Level: intermediate
73: .seealso: PetscSpaceSumGetConcatenate()
74: @*/
75: PetscErrorCode PetscSpaceSumSetConcatenate(PetscSpace sp,PetscBool concatenate)
76: {
78: PetscTryMethod(sp,"PetscSpaceSumSetConcatenate_C",(PetscSpace,PetscBool),(sp,concatenate));
79: return 0;
80: }
82: /*@
83: PetscSpaceSumGetSubspace - Get a space in the sum
85: Input Parameters:
86: + sp - the function space object
87: - s - The space number
89: Output Parameter:
90: . subsp - the PetscSpace
92: Level: intermediate
94: .seealso: PetscSpaceSumSetSubspace(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
95: @*/
96: PetscErrorCode PetscSpaceSumGetSubspace(PetscSpace sp,PetscInt s,PetscSpace *subsp)
97: {
100: PetscTryMethod(sp,"PetscSpaceSumGetSubspace_C",(PetscSpace,PetscInt,PetscSpace*),(sp,s,subsp));
101: return 0;
102: }
104: /*@
105: PetscSpaceSumSetSubspace - Set a space in the sum
107: Input Parameters:
108: + sp - the function space object
109: . s - The space number
110: - subsp - the number of spaces
112: Level: intermediate
114: .seealso: PetscSpaceSumGetSubspace(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
115: @*/
116: PetscErrorCode PetscSpaceSumSetSubspace(PetscSpace sp,PetscInt s,PetscSpace subsp)
117: {
120: PetscTryMethod(sp,"PetscSpaceSumSetSubspace_C",(PetscSpace,PetscInt,PetscSpace),(sp,s,subsp));
121: return 0;
122: }
124: static PetscErrorCode PetscSpaceSumGetNumSubspaces_Sum(PetscSpace space,PetscInt *numSumSpaces)
125: {
126: PetscSpace_Sum *sum = (PetscSpace_Sum*)space->data;
128: *numSumSpaces = sum->numSumSpaces;
129: return 0;
130: }
132: static PetscErrorCode PetscSpaceSumSetNumSubspaces_Sum(PetscSpace space,PetscInt numSumSpaces)
133: {
134: PetscSpace_Sum *sum = (PetscSpace_Sum*)space->data;
135: PetscInt Ns = sum->numSumSpaces;
138: if (numSumSpaces == Ns) return 0;
139: if (Ns >= 0) {
140: PetscInt s;
141: for (s=0; s<Ns; ++s) {
142: PetscSpaceDestroy(&sum->sumspaces[s]);
143: }
144: PetscFree(sum->sumspaces);
145: }
147: Ns = sum->numSumSpaces = numSumSpaces;
148: PetscCalloc1(Ns,&sum->sumspaces);
149: return 0;
150: }
152: static PetscErrorCode PetscSpaceSumGetConcatenate_Sum(PetscSpace sp,PetscBool *concatenate)
153: {
154: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
156: *concatenate = sum->concatenate;
157: return 0;
158: }
160: static PetscErrorCode PetscSpaceSumSetConcatenate_Sum(PetscSpace sp,PetscBool concatenate)
161: {
162: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
166: sum->concatenate = concatenate;
167: return 0;
168: }
170: static PetscErrorCode PetscSpaceSumGetSubspace_Sum(PetscSpace space,PetscInt s,PetscSpace *subspace)
171: {
172: PetscSpace_Sum *sum = (PetscSpace_Sum*)space->data;
173: PetscInt Ns = sum->numSumSpaces;
178: *subspace = sum->sumspaces[s];
179: return 0;
180: }
182: static PetscErrorCode PetscSpaceSumSetSubspace_Sum(PetscSpace space,PetscInt s,PetscSpace subspace)
183: {
184: PetscSpace_Sum *sum = (PetscSpace_Sum*)space->data;
185: PetscInt Ns = sum->numSumSpaces;
191: PetscObjectReference((PetscObject)subspace);
192: PetscSpaceDestroy(&sum->sumspaces[s]);
193: sum->sumspaces[s] = subspace;
194: return 0;
195: }
197: static PetscErrorCode PetscSpaceSetFromOptions_Sum(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
198: {
199: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
200: PetscInt Ns,Nc,Nv,deg,i;
201: PetscBool concatenate = PETSC_TRUE;
202: const char *prefix;
204: PetscSpaceGetNumVariables(sp,&Nv);
205: if (!Nv) return 0;
206: PetscSpaceGetNumComponents(sp,&Nc);
207: PetscSpaceSumGetNumSubspaces(sp,&Ns);
208: PetscSpaceGetDegree(sp,°,NULL);
209: Ns = (Ns == PETSC_DEFAULT) ? 1 : Ns;
211: PetscOptionsHead(PetscOptionsObject,"PetscSpace sum options");
212: PetscOptionsBoundedInt("-petscspace_sum_spaces","The number of subspaces","PetscSpaceSumSetNumSubspaces",Ns,&Ns,NULL,0);
213: PetscCall(PetscOptionsBool("-petscspace_sum_concatenate","Subspaces are concatenated components of the final space","PetscSpaceSumSetFromOptions",
214: concatenate,&concatenate,NULL));
215: PetscOptionsTail();
218: if (Ns != sum->numSumSpaces) {
219: PetscSpaceSumSetNumSubspaces(sp,Ns);
220: }
221: PetscObjectGetOptionsPrefix((PetscObject)sp,&prefix);
222: for (i=0; i<Ns; ++i) {
223: PetscInt sNv;
224: PetscSpace subspace;
226: PetscSpaceSumGetSubspace(sp,i,&subspace);
227: if (!subspace) {
228: char subspacePrefix[256];
230: PetscSpaceCreate(PetscObjectComm((PetscObject)sp),&subspace);
231: PetscObjectSetOptionsPrefix((PetscObject)subspace,prefix);
232: PetscSNPrintf(subspacePrefix,256,"sumcomp_%D_",i);
233: PetscObjectAppendOptionsPrefix((PetscObject)subspace,subspacePrefix);
234: } else {
235: PetscObjectReference((PetscObject)subspace);
236: }
237: PetscSpaceSetFromOptions(subspace);
238: PetscSpaceGetNumVariables(subspace,&sNv);
240: PetscSpaceSumSetSubspace(sp,i,subspace);
241: PetscSpaceDestroy(&subspace);
242: }
243: return 0;
244: }
246: static PetscErrorCode PetscSpaceSetUp_Sum(PetscSpace sp)
247: {
248: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
249: PetscBool concatenate = PETSC_TRUE;
250: PetscBool uniform;
251: PetscInt Nv,Ns,Nc,i,sum_Nc = 0,deg = PETSC_MAX_INT,maxDeg = PETSC_MIN_INT;
252: PetscInt minNc,maxNc;
254: if (sum->setupCalled) return 0;
256: PetscSpaceGetNumVariables(sp,&Nv);
257: PetscSpaceGetNumComponents(sp,&Nc);
258: PetscSpaceSumGetNumSubspaces(sp,&Ns);
259: if (Ns == PETSC_DEFAULT) {
260: Ns = 1;
261: PetscSpaceSumSetNumSubspaces(sp,Ns);
262: }
264: uniform = PETSC_TRUE;
265: if (Ns) {
266: PetscSpace s0;
268: PetscSpaceSumGetSubspace(sp,0,&s0);
269: for (PetscInt i = 1; i < Ns; i++) {
270: PetscSpace si;
272: PetscSpaceSumGetSubspace(sp,i,&si);
273: if (si != s0) {
274: uniform = PETSC_FALSE;
275: break;
276: }
277: }
278: }
280: minNc = Nc;
281: maxNc = Nc;
282: for (i=0; i<Ns; ++i) {
283: PetscInt sNv,sNc,iDeg,iMaxDeg;
284: PetscSpace si;
286: PetscSpaceSumGetSubspace(sp,i,&si);
287: PetscSpaceSetUp(si);
288: PetscSpaceGetNumVariables(si,&sNv);
290: PetscSpaceGetNumComponents(si,&sNc);
291: if (i == 0 && sNc == Nc) concatenate = PETSC_FALSE;
292: minNc = PetscMin(minNc, sNc);
293: maxNc = PetscMax(maxNc, sNc);
294: sum_Nc += sNc;
295: PetscSpaceSumGetSubspace(sp,i,&si);
296: PetscSpaceGetDegree(si,&iDeg,&iMaxDeg);
297: deg = PetscMin(deg,iDeg);
298: maxDeg = PetscMax(maxDeg,iMaxDeg);
299: }
301: if (concatenate) {
302: if (sum_Nc != Nc) {
303: SETERRQ(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,"Total number of subspace components (%D) does not match number of target space components (%D).",sum_Nc,Nc);
304: }
305: } else {
307: }
309: sp->degree = deg;
310: sp->maxDegree = maxDeg;
311: sum->concatenate = concatenate;
312: sum->uniform = uniform;
313: sum->setupCalled = PETSC_TRUE;
314: return 0;
315: }
317: static PetscErrorCode PetscSpaceSumView_Ascii(PetscSpace sp,PetscViewer v)
318: {
319: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
320: PetscBool concatenate = sum->concatenate;
321: PetscInt i,Ns = sum->numSumSpaces;
323: if (concatenate) {
324: PetscViewerASCIIPrintf(v,"Sum space of %D concatenated subspaces%s\n",Ns, sum->uniform ? " (all identical)": "");
325: } else {
326: PetscViewerASCIIPrintf(v,"Sum space of %D subspaces%s\n",Ns, sum->uniform ? " (all identical)" : "");
327: }
328: for (i=0; i < (sum->uniform ? (Ns > 0 ? 1 : 0) : Ns); ++i) {
329: PetscViewerASCIIPushTab(v);
330: PetscSpaceView(sum->sumspaces[i],v);
331: PetscViewerASCIIPopTab(v);
332: }
333: return 0;
334: }
336: static PetscErrorCode PetscSpaceView_Sum(PetscSpace sp,PetscViewer viewer)
337: {
338: PetscBool iascii;
340: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
341: if (iascii) {
342: PetscSpaceSumView_Ascii(sp,viewer);
343: }
344: return 0;
345: }
347: static PetscErrorCode PetscSpaceDestroy_Sum(PetscSpace sp)
348: {
349: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
350: PetscInt i,Ns = sum->numSumSpaces;
352: for (i=0; i<Ns; ++i) {
353: PetscSpaceDestroy(&sum->sumspaces[i]);
354: }
355: PetscFree(sum->sumspaces);
356: if (sum->heightsubspaces) {
357: PetscInt d;
359: /* sp->Nv is the spatial dimension, so it is equal to the number
360: * of subspaces on higher co-dimension points */
361: for (d = 0; d < sp->Nv; ++d) {
362: PetscSpaceDestroy(&sum->heightsubspaces[d]);
363: }
364: }
365: PetscFree(sum->heightsubspaces);
366: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetSubspace_C",NULL);
367: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetSubspace_C",NULL);
368: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetNumSubspaces_C",NULL);
369: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetNumSubspaces_C",NULL);
370: PetscFree(sum);
371: return 0;
372: }
374: static PetscErrorCode PetscSpaceGetDimension_Sum(PetscSpace sp,PetscInt *dim)
375: {
376: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
377: PetscInt i,d = 0,Ns = sum->numSumSpaces;
379: if (!sum->setupCalled) {
380: PetscSpaceSetUp(sp);
381: PetscSpaceGetDimension(sp, dim);
382: return 0;
383: }
385: for (i=0; i<Ns; ++i) {
386: PetscInt id;
388: PetscSpaceGetDimension(sum->sumspaces[i],&id);
389: d += id;
390: }
392: *dim = d;
393: return 0;
394: }
396: static PetscErrorCode PetscSpaceEvaluate_Sum(PetscSpace sp,PetscInt npoints,const PetscReal points[],PetscReal B[],PetscReal D[],PetscReal H[])
397: {
398: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
399: PetscBool concatenate = sum->concatenate;
400: DM dm = sp->dm;
401: PetscInt Nc = sp->Nc,Nv = sp->Nv,Ns = sum->numSumSpaces;
402: PetscInt i,s,offset,ncoffset,pdimfull,numelB,numelD,numelH;
403: PetscReal *sB = NULL,*sD = NULL,*sH = NULL;
405: if (!sum->setupCalled) {
406: PetscSpaceSetUp(sp);
407: PetscSpaceEvaluate(sp, npoints, points, B, D, H);
408: return 0;
409: }
410: PetscSpaceGetDimension(sp,&pdimfull);
411: numelB = npoints*pdimfull*Nc;
412: numelD = numelB*Nv;
413: numelH = numelD*Nv;
414: if (B || D || H) {
415: DMGetWorkArray(dm,numelB,MPIU_REAL,&sB);
416: }
417: if (D || H) {
418: DMGetWorkArray(dm,numelD,MPIU_REAL,&sD);
419: }
420: if (H) {
421: DMGetWorkArray(dm,numelH,MPIU_REAL,&sH);
422: }
423: if (B)
424: for (i=0; i<numelB; ++i) B[i] = 0.;
425: if (D)
426: for (i=0; i<numelD; ++i) D[i] = 0.;
427: if (H)
428: for (i=0; i<numelH; ++i) H[i] = 0.;
430: for (s=0,offset=0,ncoffset=0; s<Ns; ++s) {
431: PetscInt sNv,spdim,sNc,p;
433: PetscSpaceGetNumVariables(sum->sumspaces[s],&sNv);
434: PetscSpaceGetNumComponents(sum->sumspaces[s],&sNc);
435: PetscSpaceGetDimension(sum->sumspaces[s],&spdim);
437: if (s == 0 || !sum->uniform) {
438: PetscSpaceEvaluate(sum->sumspaces[s],npoints,points,sB,sD,sH);
439: }
440: if (B || D || H) {
441: for (p=0; p<npoints; ++p) {
442: PetscInt j;
444: for (j=0; j<spdim; ++j) {
445: PetscInt c;
447: for (c=0; c<sNc; ++c) {
448: PetscInt compoffset,BInd,sBInd;
450: compoffset = concatenate ? c+ncoffset : c;
451: BInd = (p*pdimfull + j + offset)*Nc + compoffset;
452: sBInd = (p*spdim + j)*sNc + c;
453: if (B) B[BInd] = sB[sBInd];
454: if (D || H) {
455: PetscInt v;
457: for (v=0; v<Nv; ++v) {
458: PetscInt DInd,sDInd;
460: DInd = BInd*Nv + v;
461: sDInd = sBInd*Nv + v;
462: if (D) D[DInd] = sD[sDInd];
463: if (H) {
464: PetscInt v2;
466: for (v2=0; v2<Nv; ++v2) {
467: PetscInt HInd,sHInd;
469: HInd = DInd*Nv + v2;
470: sHInd = sDInd*Nv + v2;
471: H[HInd] = sH[sHInd];
472: }
473: }
474: }
475: }
476: }
477: }
478: }
479: }
480: offset += spdim;
481: ncoffset += sNc;
482: }
484: if (H) {
485: DMRestoreWorkArray(dm,numelH,MPIU_REAL,&sH);
486: }
487: if (D || H) {
488: DMRestoreWorkArray(dm,numelD,MPIU_REAL,&sD);
489: }
490: if (B || D || H) {
491: DMRestoreWorkArray(dm,numelB,MPIU_REAL,&sB);
492: }
493: return 0;
494: }
496: static PetscErrorCode PetscSpaceGetHeightSubspace_Sum(PetscSpace sp, PetscInt height, PetscSpace *subsp)
497: {
498: PetscSpace_Sum *sum = (PetscSpace_Sum *) sp->data;
499: PetscInt Nc, dim, order;
500: PetscBool tensor;
502: PetscSpaceGetNumComponents(sp, &Nc);
503: PetscSpaceGetNumVariables(sp, &dim);
504: PetscSpaceGetDegree(sp, &order, NULL);
505: PetscSpacePolynomialGetTensor(sp, &tensor);
507: if (!sum->heightsubspaces) PetscCalloc1(dim, &sum->heightsubspaces);
508: if (height <= dim) {
509: if (!sum->heightsubspaces[height-1]) {
510: PetscSpace sub;
511: const char *name;
513: PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);
514: PetscObjectGetName((PetscObject) sp, &name);
515: PetscObjectSetName((PetscObject) sub, name);
516: PetscSpaceSetType(sub, PETSCSPACESUM);
517: PetscSpaceSumSetNumSubspaces(sub, sum->numSumSpaces);
518: PetscSpaceSumSetConcatenate(sub, sum->concatenate);
519: PetscSpaceSetNumComponents(sub, Nc);
520: PetscSpaceSetNumVariables(sub, dim-height);
521: for (PetscInt i = 0; i < sum->numSumSpaces; i++) {
522: PetscSpace subh;
524: PetscSpaceGetHeightSubspace(sum->sumspaces[i], height, &subh);
525: PetscSpaceSumSetSubspace(sub, i, subh);
526: }
527: PetscSpaceSetUp(sub);
528: sum->heightsubspaces[height-1] = sub;
529: }
530: *subsp = sum->heightsubspaces[height-1];
531: } else {
532: *subsp = NULL;
533: }
534: return 0;
535: }
537: static PetscErrorCode PetscSpaceInitialize_Sum(PetscSpace sp)
538: {
539: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Sum;
540: sp->ops->setup = PetscSpaceSetUp_Sum;
541: sp->ops->view = PetscSpaceView_Sum;
542: sp->ops->destroy = PetscSpaceDestroy_Sum;
543: sp->ops->getdimension = PetscSpaceGetDimension_Sum;
544: sp->ops->evaluate = PetscSpaceEvaluate_Sum;
545: sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Sum;
547: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetNumSubspaces_C",PetscSpaceSumGetNumSubspaces_Sum);
548: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetNumSubspaces_C",PetscSpaceSumSetNumSubspaces_Sum);
549: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetSubspace_C",PetscSpaceSumGetSubspace_Sum);
550: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetSubspace_C",PetscSpaceSumSetSubspace_Sum);
551: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetConcatenate_C",PetscSpaceSumGetConcatenate_Sum);
552: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetConcatenate_C",PetscSpaceSumSetConcatenate_Sum);
553: return 0;
554: }
556: /*MC
557: PETSCSPACESUM = "sum" - A PetscSpace object that encapsulates a sum of subspaces.
558: That sum can either be direct or concatenate a concatenation.For example if A and B are spaces each with 2 components,
559: the direct sum of A and B will also have 2 components while the concatenated sum will have 4 components.In both cases A and B must be defined over the
560: same number of variables.
562: Level: intermediate
564: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
565: M*/
566: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Sum(PetscSpace sp)
567: {
568: PetscSpace_Sum *sum;
571: PetscNewLog(sp,&sum);
572: sum->numSumSpaces = PETSC_DEFAULT;
573: sp->data = sum;
574: PetscSpaceInitialize_Sum(sp);
575: return 0;
576: }
578: PETSC_EXTERN PetscErrorCode PetscSpaceCreateSum(PetscInt numSubspaces,const PetscSpace subspaces[],PetscBool concatenate,PetscSpace *sumSpace)
579: {
580: PetscInt i,Nv,Nc = 0;
582: if (sumSpace) {
583: PetscSpaceDestroy(sumSpace);
584: }
585: PetscSpaceCreate(PetscObjectComm((PetscObject)subspaces[0]),sumSpace);
586: PetscSpaceSetType(*sumSpace,PETSCSPACESUM);
587: PetscSpaceSumSetNumSubspaces(*sumSpace,numSubspaces);
588: PetscSpaceSumSetConcatenate(*sumSpace,concatenate);
589: for (i=0; i<numSubspaces; ++i) {
590: PetscInt sNc;
592: PetscSpaceSumSetSubspace(*sumSpace,i,subspaces[i]);
593: PetscSpaceGetNumComponents(subspaces[i],&sNc);
594: if (concatenate) Nc += sNc;
595: else Nc = sNc;
596: }
597: PetscSpaceGetNumVariables(subspaces[0],&Nv);
598: PetscSpaceSetNumComponents(*sumSpace,Nc);
599: PetscSpaceSetNumVariables(*sumSpace,Nv);
600: PetscSpaceSetUp(*sumSpace);
602: return 0;
603: }