Actual source code: spacesum.c
petsc-3.14.6 2021-03-30
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: {
22: PetscTryMethod(sp,"PetscSpaceSumGetNumSubspaces_C",(PetscSpace,PetscInt*),(sp,numSumSpaces));
23: return(0);
24: }
26: /*@
27: PetscSpaceSumSetNumSubspaces - Set the number of spaces in the sum
29: Input Parameters:
30: + sp - the function space object
31: - numSumSpaces - the number of spaces
33: Level: intermediate
35: .seealso: PetscSpaceSumGetNumSubspaces(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
36: @*/
37: PetscErrorCode PetscSpaceSumSetNumSubspaces(PetscSpace sp,PetscInt numSumSpaces)
38: {
43: PetscTryMethod(sp,"PetscSpaceSumSetNumSubspaces_C",(PetscSpace,PetscInt),(sp,numSumSpaces));
44: return(0);
45: }
47: /*@
48: PetscSpaceSumGetConcatenate - Get the concatenate flag for this space.
49: A concatenated sum space will have number of components equal to the sum of the number of components of all subspaces.A non-concatenated,
50: or direct sum space will have the same number of components as its subspaces .
52: Input Parameters:
53: . sp - the function space object
55: Output Parameters:
56: . concatenate - flag indicating whether subspaces are concatenated.
58: Level: intermediate
60: .seealso: PetscSpaceSumSetConcatenate()
61: @*/
62: PetscErrorCode PetscSpaceSumGetConcatenate(PetscSpace sp,PetscBool *concatenate)
63: {
68: PetscTryMethod(sp,"PetscSpaceSumGetConcatenate_C",(PetscSpace,PetscBool*),(sp,concatenate));
69: return(0);
70: }
72: /*@
73: PetscSpaceSumSetConcatenate - Sets the concatenate flag for this space.
74: A concatenated sum space will have number of components equal to the sum of the number of components of all subspaces.A non-concatenated,
75: or direct sum space will have the same number of components as its subspaces .
77: Input Parameters:
78: + sp - the function space object
79: - concatenate - are subspaces concatenated components (true) or direct summands (false)
81: Level: intermediate
82: .seealso: PetscSpaceSumGetConcatenate()
83: @*/
84: PetscErrorCode PetscSpaceSumSetConcatenate(PetscSpace sp,PetscBool concatenate)
85: {
90: PetscTryMethod(sp,"PetscSpaceSumSetConcatenate_C",(PetscSpace,PetscBool),(sp,concatenate));
91: return(0);
92: }
94: /*@
95: PetscSpaceSumGetSubspace - Get a space in the sum
97: Input Parameters:
98: + sp - the function space object
99: - s - The space number
101: Output Parameter:
102: . subsp - the PetscSpace
104: Level: intermediate
106: .seealso: PetscSpaceSumSetSubspace(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
107: @*/
108: PetscErrorCode PetscSpaceSumGetSubspace(PetscSpace sp,PetscInt s,PetscSpace *subsp)
109: {
115: PetscTryMethod(sp,"PetscSpaceSumGetSubspace_C",(PetscSpace,PetscInt,PetscSpace*),(sp,s,subsp));
116: return(0);
117: }
119: /*@
120: PetscSpaceSumSetSubspace - Set a space in the sum
122: Input Parameters:
123: + sp - the function space object
124: . s - The space number
125: - subsp - the number of spaces
127: Level: intermediate
129: .seealso: PetscSpaceSumGetSubspace(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
130: @*/
131: PetscErrorCode PetscSpaceSumSetSubspace(PetscSpace sp,PetscInt s,PetscSpace subsp)
132: {
138: PetscTryMethod(sp,"PetscSpaceSumSetSubspace_C",(PetscSpace,PetscInt,PetscSpace),(sp,s,subsp));
139: return(0);
140: }
142: static PetscErrorCode PetscSpaceSumGetNumSubspaces_Sum(PetscSpace space,PetscInt *numSumSpaces)
143: {
144: PetscSpace_Sum *sum = (PetscSpace_Sum*)space->data;
147: *numSumSpaces = sum->numSumSpaces;
148: return(0);
149: }
151: static PetscErrorCode PetscSpaceSumSetNumSubspaces_Sum(PetscSpace space,PetscInt numSumSpaces)
152: {
153: PetscSpace_Sum *sum = (PetscSpace_Sum*)space->data;
154: PetscInt Ns = sum->numSumSpaces;
158: if (sum->setupCalled) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,
159: "Cannot change number of subspaces after setup called\n");
160: if (numSumSpaces == Ns) return(0);
161: if (Ns >= 0) {
162: PetscInt s;
163: for (s=0; s<Ns; ++s) {
164: PetscSpaceDestroy(&sum->sumspaces[s]);
165: }
166: PetscFree(sum->sumspaces);
167: }
169: Ns = sum->numSumSpaces = numSumSpaces;
170: PetscCalloc1(Ns,&sum->sumspaces);
171: return(0);
172: }
174: static PetscErrorCode PetscSpaceSumGetConcatenate_Sum(PetscSpace sp,PetscBool *concatenate)
175: {
176: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
179: *concatenate = sum->concatenate;
180: return(0);
181: }
183: static PetscErrorCode PetscSpaceSumSetConcatenate_Sum(PetscSpace sp,PetscBool concatenate)
184: {
185: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
188: if (sum->setupCalled) SETERRQ(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_WRONGSTATE,"Cannot change space concatenation after setup called.\n");
190: sum->concatenate = concatenate;
191: return(0);
192: }
194: static PetscErrorCode PetscSpaceSumGetSubspace_Sum(PetscSpace space,PetscInt s,PetscSpace *subspace)
195: {
196: PetscSpace_Sum *sum = (PetscSpace_Sum*)space->data;
197: PetscInt Ns = sum->numSumSpaces;
200: if (Ns < 0) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSpaceSumSetNumSubspaces() first\n");
201: if (s<0 || s>=Ns) SETERRQ1 (PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_OUTOFRANGE,"Invalid subspace number %D\n",subspace);
203: *subspace = sum->sumspaces[s];
204: return(0);
205: }
207: static PetscErrorCode PetscSpaceSumSetSubspace_Sum(PetscSpace space,PetscInt s,PetscSpace subspace)
208: {
209: PetscSpace_Sum *sum = (PetscSpace_Sum*)space->data;
210: PetscInt Ns = sum->numSumSpaces;
214: if (sum->setupCalled) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Cannot change subspace after setup called\n");
215: if (Ns < 0) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSpaceSumSetNumSubspaces() first\n");
216: if (s < 0 || s >= Ns) SETERRQ1(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_OUTOFRANGE,"Invalid subspace number %D\n",subspace);
218: PetscObjectReference((PetscObject)subspace);
219: PetscSpaceDestroy(&sum->sumspaces[s]);
220: sum->sumspaces[s] = subspace;
221: return(0);
222: }
224: static PetscErrorCode PetscSpaceSetFromOptions_Sum(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
225: {
226: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
227: PetscInt Ns,Nc,Nv,deg,i;
228: PetscBool concatenate = PETSC_TRUE;
229: const char *prefix;
233: PetscSpaceGetNumVariables(sp,&Nv);
234: if (!Nv) return(0);
235: PetscSpaceGetNumComponents(sp,&Nc);
236: PetscSpaceSumGetNumSubspaces(sp,&Ns);
237: PetscSpaceGetDegree(sp,°,NULL);
238: Ns = (Ns == PETSC_DEFAULT) ? 1 : Ns;
240: PetscOptionsHead(PetscOptionsObject,"PetscSpace sum options");
241: PetscOptionsBoundedInt("-petscspace_sum_spaces","The number of subspaces","PetscSpaceSumSetNumSubspaces",Ns,&Ns,NULL,0);
242: PetscOptionsBool("-petscspace_sum_concatenate","Subspaces are concatenated components of the final space","PetscSpaceSumSetFromOptions",
243: concatenate,&concatenate,NULL);
244: PetscOptionsTail();
246: if (Ns < 0 || (Nv > 0 && Ns == 0)) SETERRQ1(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,"Cannot have a sum space of %D spaces\n",Ns);
247: if (Ns != sum->numSumSpaces) {
248: PetscSpaceSumSetNumSubspaces(sp,Ns);
249: }
250: PetscObjectGetOptionsPrefix((PetscObject)sp,&prefix);
251: for (i=0; i<Ns; ++i) {
252: PetscInt sNv;
253: PetscSpace subspace;
255: PetscSpaceSumGetSubspace(sp,i,&subspace);
256: if (!subspace) {
257: char subspacePrefix[256];
259: PetscSpaceCreate(PetscObjectComm((PetscObject)sp),&subspace);
260: PetscObjectSetOptionsPrefix((PetscObject)subspace,prefix);
261: PetscSNPrintf(subspacePrefix,256,"subspace%d_",i);
262: PetscObjectAppendOptionsPrefix((PetscObject)subspace,subspacePrefix);
263: } else {
264: PetscObjectReference((PetscObject)subspace);
265: }
266: PetscSpaceSetFromOptions(subspace);
267: PetscSpaceGetNumVariables(subspace,&sNv);
268: if (!sNv) SETERRQ1(PetscObjectComm(
269: (PetscObject)sp),PETSC_ERR_ARG_WRONGSTATE,"Subspace %D has not been set properly, number of variables is 0.\n",i);
270: PetscSpaceSumSetSubspace(sp,i,subspace);
271: PetscSpaceDestroy(&subspace);
272: }
273: return(0);
274: }
276: static PetscErrorCode PetscSpaceSetUp_Sum(PetscSpace sp)
277: {
278: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
279: PetscBool concatenate = PETSC_TRUE;
280: PetscInt Nv,Ns,Nc,i,sum_Nc = 0,deg = PETSC_MIN_INT,maxDeg = PETSC_MIN_INT;
284: if (sum->setupCalled) return(0);
286: PetscSpaceGetNumVariables(sp,&Nv);
287: PetscSpaceGetNumComponents(sp,&Nc);
288: PetscSpaceSumGetNumSubspaces(sp,&Ns);
289: if (Ns == PETSC_DEFAULT) {
290: Ns = 1;
291: PetscSpaceSumSetNumSubspaces(sp,Ns);
292: }
293: if (!Ns && Nv) SETERRQ(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,"Cannot have zero subspaces\n");
295: for (i=0; i<Ns; ++i) {
296: PetscInt sNv,sNc,iDeg,iMaxDeg;
297: PetscSpace si;
299: PetscSpaceSumGetSubspace(sp,i,&si);
300: PetscSpaceSetUp(si);
301: PetscSpaceGetNumVariables(si,&sNv);
302: if (sNv != Nv) SETERRQ3(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_WRONGSTATE,"Subspace %D has %D variables, space has %D.\n",i,sNv,Nv);
303: PetscSpaceGetNumComponents(si,&sNc);
304: if (i == 0 && sNc == Nc) concatenate = PETSC_FALSE;
305: sum_Nc += sNc;
306: PetscSpaceSumGetSubspace(sp,i,&si);
307: PetscSpaceGetDegree(si,&iDeg,&iMaxDeg);
308: deg = PetscMax(deg,iDeg);
309: maxDeg = PetscMax(maxDeg,iMaxDeg);
310: }
312: if (concatenate) {
313: if (sum_Nc != Nc) {
314: SETERRQ2(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,
315: "Total number of subspace components (%D) does not match number of target space components (%D).",sum_Nc,Nc);
316: }
317: } else {
318: if (sum_Nc != Ns*Nc) SETERRQ(PetscObjectComm(
319: (PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,"Subspaces must have same number of components as the target space.");
320: }
322: sp->degree = deg;
323: sp->maxDegree = maxDeg;
324: sum->concatenate = concatenate;
325: sum->setupCalled = PETSC_TRUE;
326: return(0);
327: }
329: static PetscErrorCode PetscSpaceSumView_Ascii(PetscSpace sp,PetscViewer v)
330: {
331: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
332: PetscBool concatenate = sum->concatenate;
333: PetscInt i,Ns = sum->numSumSpaces;
337: if (concatenate) {
338: PetscViewerASCIIPrintf(v,"Sum space of %D concatenated subspaces\n",Ns);
339: } else {
340: PetscViewerASCIIPrintf(v,"Sum space of %D subspaces\n",Ns);
341: }
342: for (i=0; i<Ns; ++i) {
343: PetscViewerASCIIPushTab(v);
344: PetscSpaceView(sum->sumspaces[i],v);
345: PetscViewerASCIIPopTab(v);
346: }
347: return(0);
348: }
350: static PetscErrorCode PetscSpaceView_Sum(PetscSpace sp,PetscViewer viewer)
351: {
352: PetscBool iascii;
356: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
357: if (iascii) {
358: PetscSpaceSumView_Ascii(sp,viewer);
359: }
360: return(0);
361: }
363: static PetscErrorCode PetscSpaceDestroy_Sum(PetscSpace sp)
364: {
365: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
366: PetscInt i,Ns = sum->numSumSpaces;
370: for (i=0; i<Ns; ++i) {
371: PetscSpaceDestroy(&sum->sumspaces[i]);
372: }
373: PetscFree(sum->sumspaces);
374: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetSubspace_C",NULL);
375: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetSubspace_C",NULL);
376: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetNumSubspaces_C",NULL);
377: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetNumSubspaces_C",NULL);
378: PetscFree(sum);
379: return(0);
380: }
382: static PetscErrorCode PetscSpaceGetDimension_Sum(PetscSpace sp,PetscInt *dim)
383: {
384: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
385: PetscInt i,d = 0,Ns = sum->numSumSpaces;
389: PetscSpaceSetUp(sp);
391: for (i=0; i<Ns; ++i) {
392: PetscInt id;
394: PetscSpaceGetDimension(sum->sumspaces[i],&id);
395: d += id;
396: }
398: *dim = d;
399: return(0);
400: }
402: static PetscErrorCode PetscSpaceEvaluate_Sum(PetscSpace sp,PetscInt npoints,const PetscReal points[],PetscReal B[],PetscReal D[],PetscReal H[])
403: {
404: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
405: PetscBool concatenate = sum->concatenate;
406: DM dm = sp->dm;
407: PetscInt Nc = sp->Nc,Nv = sp->Nv,Ns = sum->numSumSpaces;
408: PetscInt i,s,offset,ncoffset,pdimfull,numelB,numelD,numelH;
409: PetscReal *sB = NULL,*sD = NULL,*sH = NULL;
413: if (!sum->setupCalled) {
414: PetscSpaceSetUp(sp);
415: }
416: PetscSpaceGetDimension(sp,&pdimfull);
417: numelB = npoints*pdimfull*Nc;
418: numelD = numelB*Nv;
419: numelH = numelD*Nv;
420: if (B || D || H) {
421: DMGetWorkArray(dm,numelB,MPIU_REAL,&sB);
422: }
423: if (D || H) {
424: DMGetWorkArray(dm,numelD,MPIU_REAL,&sD);
425: }
426: if (H) {
427: DMGetWorkArray(dm,numelH,MPIU_REAL,&sH);
428: }
429: if (B)
430: for (i=0; i<numelB; ++i) B[i] = 0.;
431: if (D)
432: for (i=0; i<numelD; ++i) D[i] = 0.;
433: if (H)
434: for (i=0; i<numelH; ++i) H[i] = 0.;
436: for (s=0,offset=0,ncoffset=0; s<Ns; ++s) {
437: PetscInt sNv,spdim,sNc,p;
439: PetscSpaceGetNumVariables(sum->sumspaces[s],&sNv);
440: PetscSpaceGetNumComponents(sum->sumspaces[s],&sNc);
441: PetscSpaceGetDimension(sum->sumspaces[s],&spdim);
442: if (offset + spdim > pdimfull) SETERRQ(PetscObjectComm(
443: (PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,"Subspace dimensions exceed target space dimension.\n");
444: PetscSpaceEvaluate(sum->sumspaces[s],npoints,points,sB,sD,sH);
445: if (B || D || H) {
446: for (p=0; p<npoints; ++p) {
447: PetscInt j;
449: for (j=0; j<spdim; ++j) {
450: PetscInt c;
452: for (c=0; c<sNc; ++c) {
453: PetscInt compoffset,BInd,sBInd;
455: compoffset = concatenate ? c+ncoffset : c;
456: BInd = (p*pdimfull + j + offset)*Nc + compoffset;
457: sBInd = (p*spdim + j)*sNc + c;
458: if (B) B[BInd] += sB[sBInd];
459: if (D || H) {
460: PetscInt v;
462: for (v=0; v<Nv; ++v) {
463: PetscInt DInd,sDInd;
465: DInd = BInd*Nv + v;
466: sDInd = sBInd*Nv + v;
467: if (D) D[DInd] +=sD[sDInd];
468: if (H) {
469: PetscInt v2;
471: for (v2=0; v2<Nv; ++v2) {
472: PetscInt HInd,sHInd;
474: HInd = DInd*Nv + v2;
475: sHInd = sDInd*Nv + v2;
476: H[HInd] += sH[sHInd];
477: }
478: }
479: }
480: }
481: }
482: }
483: }
484: }
485: offset += spdim;
486: ncoffset += sNc;
487: }
489: if (H) {
490: DMRestoreWorkArray(dm,numelH,MPIU_REAL,&sH);
491: }
492: if (D || H) {
493: DMRestoreWorkArray(dm,numelD,MPIU_REAL,&sD);
494: }
495: if (B || D || H) {
496: DMRestoreWorkArray(dm,numelB,MPIU_REAL,&sB);
497: }
498: return(0);
499: }
501: static PetscErrorCode PetscSpaceInitialize_Sum(PetscSpace sp)
502: {
506: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Sum;
507: sp->ops->setup = PetscSpaceSetUp_Sum;
508: sp->ops->view = PetscSpaceView_Sum;
509: sp->ops->destroy = PetscSpaceDestroy_Sum;
510: sp->ops->getdimension = PetscSpaceGetDimension_Sum;
511: sp->ops->evaluate = PetscSpaceEvaluate_Sum;
513: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetNumSubspaces_C",PetscSpaceSumGetNumSubspaces_Sum);
514: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetNumSubspaces_C",PetscSpaceSumSetNumSubspaces_Sum);
515: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetSubspace_C",PetscSpaceSumGetSubspace_Sum);
516: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetSubspace_C",PetscSpaceSumSetSubspace_Sum);
517: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetConcatenate_C",PetscSpaceSumGetConcatenate_Sum);
518: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetConcatenate_C",PetscSpaceSumSetConcatenate_Sum);
519: return(0);
520: }
522: /*MC
523: PETSCSPACESUM = "sum" - A PetscSpace object that encapsulates a sum of subspaces.
524: That sum can either be direct or concatenate a concatenation.For example if A and B are spaces each with 2 components,
525: 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
526: same number of variables.
528: Level: intermediate
530: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
531: M*/
532: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Sum(PetscSpace sp)
533: {
534: PetscSpace_Sum *sum;
539: PetscNewLog(sp,&sum);
540: sp->data = sum;
541: PetscSpaceInitialize_Sum(sp);
542: return(0);
543: }
545: PETSC_EXTERN PetscErrorCode PetscSpaceCreateSum(PetscInt numSubspaces,const PetscSpace subspaces[],PetscBool concatenate,PetscSpace *sumSpace)
546: {
547: PetscInt i,Nv,Nc = 0;
551: if (sumSpace) {
552: PetscSpaceDestroy(sumSpace);
553: }
554: PetscSpaceCreate(PetscObjectComm((PetscObject)subspaces[0]),sumSpace);
555: PetscSpaceSetType(*sumSpace,PETSCSPACESUM);
556: PetscSpaceSumSetNumSubspaces(*sumSpace,numSubspaces);
557: PetscSpaceSumSetConcatenate(*sumSpace,concatenate);
558: for (i=0; i<numSubspaces; ++i) {
559: PetscInt sNc;
561: PetscSpaceSumSetSubspace(*sumSpace,i,subspaces[i]);
562: PetscSpaceGetNumComponents(subspaces[i],&sNc);
563: if (concatenate) Nc += sNc;
564: else Nc = sNc;
565: }
566: PetscSpaceGetNumVariables(subspaces[0],&Nv);
567: PetscSpaceSetNumComponents(*sumSpace,Nc);
568: PetscSpaceSetNumVariables(*sumSpace,Nv);
569: PetscSpaceSetUp(*sumSpace);
571: return(0);
572: }