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: {
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,"Cannot change number of subspaces after setup called\n");
159: if (numSumSpaces == Ns) return(0);
160: if (Ns >= 0) {
161: PetscInt s;
162: for (s=0; s<Ns; ++s) {
163: PetscSpaceDestroy(&sum->sumspaces[s]);
164: }
165: PetscFree(sum->sumspaces);
166: }
168: Ns = sum->numSumSpaces = numSumSpaces;
169: PetscCalloc1(Ns,&sum->sumspaces);
170: return(0);
171: }
173: static PetscErrorCode PetscSpaceSumGetConcatenate_Sum(PetscSpace sp,PetscBool *concatenate)
174: {
175: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
178: *concatenate = sum->concatenate;
179: return(0);
180: }
182: static PetscErrorCode PetscSpaceSumSetConcatenate_Sum(PetscSpace sp,PetscBool concatenate)
183: {
184: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
187: if (sum->setupCalled) SETERRQ(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_WRONGSTATE,"Cannot change space concatenation after setup called.\n");
189: sum->concatenate = concatenate;
190: return(0);
191: }
193: static PetscErrorCode PetscSpaceSumGetSubspace_Sum(PetscSpace space,PetscInt s,PetscSpace *subspace)
194: {
195: PetscSpace_Sum *sum = (PetscSpace_Sum*)space->data;
196: PetscInt Ns = sum->numSumSpaces;
199: if (Ns < 0) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSpaceSumSetNumSubspaces() first\n");
200: if (s<0 || s>=Ns) SETERRQ1 (PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_OUTOFRANGE,"Invalid subspace number %D\n",subspace);
202: *subspace = sum->sumspaces[s];
203: return(0);
204: }
206: static PetscErrorCode PetscSpaceSumSetSubspace_Sum(PetscSpace space,PetscInt s,PetscSpace subspace)
207: {
208: PetscSpace_Sum *sum = (PetscSpace_Sum*)space->data;
209: PetscInt Ns = sum->numSumSpaces;
213: if (sum->setupCalled) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Cannot change subspace after setup called\n");
214: if (Ns < 0) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSpaceSumSetNumSubspaces() first\n");
215: if (s < 0 || s >= Ns) SETERRQ1(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_OUTOFRANGE,"Invalid subspace number %D\n",subspace);
217: PetscObjectReference((PetscObject)subspace);
218: PetscSpaceDestroy(&sum->sumspaces[s]);
219: sum->sumspaces[s] = subspace;
220: return(0);
221: }
223: static PetscErrorCode PetscSpaceSetFromOptions_Sum(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
224: {
225: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
226: PetscInt Ns,Nc,Nv,deg,i;
227: PetscBool concatenate = PETSC_TRUE;
228: const char *prefix;
232: PetscSpaceGetNumVariables(sp,&Nv);
233: if (!Nv) return(0);
234: PetscSpaceGetNumComponents(sp,&Nc);
235: PetscSpaceSumGetNumSubspaces(sp,&Ns);
236: PetscSpaceGetDegree(sp,°,NULL);
237: Ns = (Ns == PETSC_DEFAULT) ? 1 : Ns;
239: PetscOptionsHead(PetscOptionsObject,"PetscSpace sum options");
240: PetscOptionsBoundedInt("-petscspace_sum_spaces","The number of subspaces","PetscSpaceSumSetNumSubspaces",Ns,&Ns,NULL,0);
241: PetscOptionsBool("-petscspace_sum_concatenate","Subspaces are concatenated components of the final space","PetscSpaceSumSetFromOptions",
242: concatenate,&concatenate,NULL);
243: PetscOptionsTail();
245: 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);
246: if (Ns != sum->numSumSpaces) {
247: PetscSpaceSumSetNumSubspaces(sp,Ns);
248: }
249: PetscObjectGetOptionsPrefix((PetscObject)sp,&prefix);
250: for (i=0; i<Ns; ++i) {
251: PetscInt sNv;
252: PetscSpace subspace;
254: PetscSpaceSumGetSubspace(sp,i,&subspace);
255: if (!subspace) {
256: char subspacePrefix[256];
258: PetscSpaceCreate(PetscObjectComm((PetscObject)sp),&subspace);
259: PetscObjectSetOptionsPrefix((PetscObject)subspace,prefix);
260: PetscSNPrintf(subspacePrefix,256,"subspace%d_",i);
261: PetscObjectAppendOptionsPrefix((PetscObject)subspace,subspacePrefix);
262: } else {
263: PetscObjectReference((PetscObject)subspace);
264: }
265: PetscSpaceSetFromOptions(subspace);
266: PetscSpaceGetNumVariables(subspace,&sNv);
267: if (!sNv) SETERRQ1(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_WRONGSTATE,"Subspace %D has not been set properly, number of variables is 0.\n",i);
268: PetscSpaceSumSetSubspace(sp,i,subspace);
269: PetscSpaceDestroy(&subspace);
270: }
271: return(0);
272: }
274: static PetscErrorCode PetscSpaceSetUp_Sum(PetscSpace sp)
275: {
276: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
277: PetscBool concatenate = PETSC_TRUE;
278: PetscInt Nv,Ns,Nc,i,sum_Nc = 0,deg = PETSC_MIN_INT,maxDeg = PETSC_MIN_INT;
282: if (sum->setupCalled) return(0);
284: PetscSpaceGetNumVariables(sp,&Nv);
285: PetscSpaceGetNumComponents(sp,&Nc);
286: PetscSpaceSumGetNumSubspaces(sp,&Ns);
287: if (Ns == PETSC_DEFAULT) {
288: Ns = 1;
289: PetscSpaceSumSetNumSubspaces(sp,Ns);
290: }
291: if (!Ns && Nv) SETERRQ(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,"Cannot have zero subspaces\n");
293: for (i=0; i<Ns; ++i) {
294: PetscInt sNv,sNc,iDeg,iMaxDeg;
295: PetscSpace si;
297: PetscSpaceSumGetSubspace(sp,i,&si);
298: PetscSpaceSetUp(si);
299: PetscSpaceGetNumVariables(si,&sNv);
300: if (sNv != Nv) SETERRQ3(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_WRONGSTATE,"Subspace %D has %D variables, space has %D.\n",i,sNv,Nv);
301: PetscSpaceGetNumComponents(si,&sNc);
302: if (i == 0 && sNc == Nc) concatenate = PETSC_FALSE;
303: sum_Nc += sNc;
304: PetscSpaceSumGetSubspace(sp,i,&si);
305: PetscSpaceGetDegree(si,&iDeg,&iMaxDeg);
306: deg = PetscMax(deg,iDeg);
307: maxDeg = PetscMax(maxDeg,iMaxDeg);
308: }
310: if (concatenate) {
311: if (sum_Nc != Nc) {
312: SETERRQ2(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);
313: }
314: } else {
315: if (sum_Nc != Ns*Nc) SETERRQ(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,"Subspaces must have same number of components as the target space.");
316: }
318: sp->degree = deg;
319: sp->maxDegree = maxDeg;
320: sum->concatenate = concatenate;
321: sum->setupCalled = PETSC_TRUE;
322: return(0);
323: }
325: static PetscErrorCode PetscSpaceSumView_Ascii(PetscSpace sp,PetscViewer v)
326: {
327: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
328: PetscBool concatenate = sum->concatenate;
329: PetscInt i,Ns = sum->numSumSpaces;
333: if (concatenate) {
334: PetscViewerASCIIPrintf(v,"Sum space of %D concatenated subspaces\n",Ns);
335: } else {
336: PetscViewerASCIIPrintf(v,"Sum space of %D subspaces\n",Ns);
337: }
338: for (i=0; i<Ns; ++i) {
339: PetscViewerASCIIPushTab(v);
340: PetscSpaceView(sum->sumspaces[i],v);
341: PetscViewerASCIIPopTab(v);
342: }
343: return(0);
344: }
346: static PetscErrorCode PetscSpaceView_Sum(PetscSpace sp,PetscViewer viewer)
347: {
348: PetscBool iascii;
352: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
353: if (iascii) {
354: PetscSpaceSumView_Ascii(sp,viewer);
355: }
356: return(0);
357: }
359: static PetscErrorCode PetscSpaceDestroy_Sum(PetscSpace sp)
360: {
361: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
362: PetscInt i,Ns = sum->numSumSpaces;
366: for (i=0; i<Ns; ++i) {
367: PetscSpaceDestroy(&sum->sumspaces[i]);
368: }
369: PetscFree(sum->sumspaces);
370: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetSubspace_C",NULL);
371: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetSubspace_C",NULL);
372: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetNumSubspaces_C",NULL);
373: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetNumSubspaces_C",NULL);
374: PetscFree(sum);
375: return(0);
376: }
378: static PetscErrorCode PetscSpaceGetDimension_Sum(PetscSpace sp,PetscInt *dim)
379: {
380: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
381: PetscInt i,d = 0,Ns = sum->numSumSpaces;
385: PetscSpaceSetUp(sp);
387: for (i=0; i<Ns; ++i) {
388: PetscInt id;
390: PetscSpaceGetDimension(sum->sumspaces[i],&id);
391: d += id;
392: }
394: *dim = d;
395: return(0);
396: }
398: static PetscErrorCode PetscSpaceEvaluate_Sum(PetscSpace sp,PetscInt npoints,const PetscReal points[],PetscReal B[],PetscReal D[],PetscReal H[])
399: {
400: PetscSpace_Sum *sum = (PetscSpace_Sum*)sp->data;
401: PetscBool concatenate = sum->concatenate;
402: DM dm = sp->dm;
403: PetscInt Nc = sp->Nc,Nv = sp->Nv,Ns = sum->numSumSpaces;
404: PetscInt i,s,offset,ncoffset,pdimfull,numelB,numelD,numelH;
405: PetscReal *sB = NULL,*sD = NULL,*sH = NULL;
409: if (!sum->setupCalled) {
410: PetscSpaceSetUp(sp);
411: }
412: PetscSpaceGetDimension(sp,&pdimfull);
413: numelB = npoints*pdimfull*Nc;
414: numelD = numelB*Nv;
415: numelH = numelD*Nv;
416: if (B || D || H) {
417: DMGetWorkArray(dm,numelB,MPIU_REAL,&sB);
418: }
419: if (D || H) {
420: DMGetWorkArray(dm,numelD,MPIU_REAL,&sD);
421: }
422: if (H) {
423: DMGetWorkArray(dm,numelH,MPIU_REAL,&sH);
424: }
425: if (B)
426: for (i=0; i<numelB; ++i) B[i] = 0.;
427: if (D)
428: for (i=0; i<numelD; ++i) D[i] = 0.;
429: if (H)
430: for (i=0; i<numelH; ++i) H[i] = 0.;
432: for (s=0,offset=0,ncoffset=0; s<Ns; ++s) {
433: PetscInt sNv,spdim,sNc,p;
435: PetscSpaceGetNumVariables(sum->sumspaces[s],&sNv);
436: PetscSpaceGetNumComponents(sum->sumspaces[s],&sNc);
437: PetscSpaceGetDimension(sum->sumspaces[s],&spdim);
438: if (offset + spdim > pdimfull) SETERRQ(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,"Subspace dimensions exceed target space dimension.\n");
439: PetscSpaceEvaluate(sum->sumspaces[s],npoints,points,sB,sD,sH);
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 PetscSpaceInitialize_Sum(PetscSpace sp)
497: {
501: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Sum;
502: sp->ops->setup = PetscSpaceSetUp_Sum;
503: sp->ops->view = PetscSpaceView_Sum;
504: sp->ops->destroy = PetscSpaceDestroy_Sum;
505: sp->ops->getdimension = PetscSpaceGetDimension_Sum;
506: sp->ops->evaluate = PetscSpaceEvaluate_Sum;
508: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetNumSubspaces_C",PetscSpaceSumGetNumSubspaces_Sum);
509: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetNumSubspaces_C",PetscSpaceSumSetNumSubspaces_Sum);
510: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetSubspace_C",PetscSpaceSumGetSubspace_Sum);
511: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetSubspace_C",PetscSpaceSumSetSubspace_Sum);
512: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumGetConcatenate_C",PetscSpaceSumGetConcatenate_Sum);
513: PetscObjectComposeFunction((PetscObject)sp,"PetscSpaceSumSetConcatenate_C",PetscSpaceSumSetConcatenate_Sum);
514: return(0);
515: }
517: /*MC
518: PETSCSPACESUM = "sum" - A PetscSpace object that encapsulates a sum of subspaces.
519: That sum can either be direct or concatenate a concatenation.For example if A and B are spaces each with 2 components,
520: 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
521: same number of variables.
523: Level: intermediate
525: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
526: M*/
527: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Sum(PetscSpace sp)
528: {
529: PetscSpace_Sum *sum;
534: PetscNewLog(sp,&sum);
535: sp->data = sum;
536: PetscSpaceInitialize_Sum(sp);
537: return(0);
538: }
540: PETSC_EXTERN PetscErrorCode PetscSpaceCreateSum(PetscInt numSubspaces,const PetscSpace subspaces[],PetscBool concatenate,PetscSpace *sumSpace)
541: {
542: PetscInt i,Nv,Nc = 0;
546: if (sumSpace) {
547: PetscSpaceDestroy(sumSpace);
548: }
549: PetscSpaceCreate(PetscObjectComm((PetscObject)subspaces[0]),sumSpace);
550: PetscSpaceSetType(*sumSpace,PETSCSPACESUM);
551: PetscSpaceSumSetNumSubspaces(*sumSpace,numSubspaces);
552: PetscSpaceSumSetConcatenate(*sumSpace,concatenate);
553: for (i=0; i<numSubspaces; ++i) {
554: PetscInt sNc;
556: PetscSpaceSumSetSubspace(*sumSpace,i,subspaces[i]);
557: PetscSpaceGetNumComponents(subspaces[i],&sNc);
558: if (concatenate) Nc += sNc;
559: else Nc = sNc;
560: }
561: PetscSpaceGetNumVariables(subspaces[0],&Nv);
562: PetscSpaceSetNumComponents(*sumSpace,Nc);
563: PetscSpaceSetNumVariables(*sumSpace,Nv);
564: PetscSpaceSetUp(*sumSpace);
566: return(0);
567: }