xref: /petsc/src/dm/interface/dmi.c (revision efd4aadf157bf1ba2d80c2be092fcf4247860003) ! 
1 #include <petsc/private/dmimpl.h>     /*I      "petscdm.h"     I*/
2 #include <petscds.h>
3 
4 PetscErrorCode DMCreateGlobalVector_Section_Private(DM dm,Vec *vec)
5 {
6   PetscSection   gSection;
7   PetscInt       localSize, bs, blockSize = -1, pStart, pEnd, p;
8   PetscErrorCode ierr;
9   PetscInt       in[2],out[2];
10 
11   PetscFunctionBegin;
12   ierr = DMGetDefaultGlobalSection(dm, &gSection);CHKERRQ(ierr);
13   ierr = PetscSectionGetChart(gSection, &pStart, &pEnd);CHKERRQ(ierr);
14   for (p = pStart; p < pEnd; ++p) {
15     PetscInt dof, cdof;
16 
17     ierr = PetscSectionGetDof(gSection, p, &dof);CHKERRQ(ierr);
18     ierr = PetscSectionGetConstraintDof(gSection, p, &cdof);CHKERRQ(ierr);
19 
20     if (dof > 0) {
21       if (blockSize < 0 && dof-cdof > 0) {
22         /* set blockSize */
23         blockSize = dof-cdof;
24       } else if (dof-cdof != blockSize) {
25         /* non-identical blockSize, set it as 1 */
26         blockSize = 1;
27         break;
28       }
29     }
30   }
31 
32   in[0] = -blockSize;
33   in[1] = blockSize;
34   ierr = MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
35   /* -out[0] = min(blockSize), out[1] = max(blockSize) */
36   if (-out[0] == out[1]) {
37     bs = out[1];
38   } else bs = 1;
39 
40   if (bs < 0) { /* Everyone was empty */
41     blockSize = 1;
42     bs        = 1;
43   }
44 
45   ierr = PetscSectionGetConstrainedStorageSize(gSection, &localSize);CHKERRQ(ierr);
46   if (localSize%blockSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mismatch between blocksize %d and local storage size %d", blockSize, localSize);
47   ierr = VecCreate(PetscObjectComm((PetscObject)dm), vec);CHKERRQ(ierr);
48   ierr = VecSetSizes(*vec, localSize, PETSC_DETERMINE);CHKERRQ(ierr);
49   ierr = VecSetBlockSize(*vec, bs);CHKERRQ(ierr);
50   ierr = VecSetType(*vec,dm->vectype);CHKERRQ(ierr);
51   ierr = VecSetDM(*vec, dm);CHKERRQ(ierr);
52   /* ierr = VecSetLocalToGlobalMapping(*vec, dm->ltogmap);CHKERRQ(ierr); */
53   PetscFunctionReturn(0);
54 }
55 
56 PetscErrorCode DMCreateLocalVector_Section_Private(DM dm,Vec *vec)
57 {
58   PetscSection   section;
59   PetscInt       localSize, blockSize = -1, pStart, pEnd, p;
60   PetscErrorCode ierr;
61 
62   PetscFunctionBegin;
63   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
64   ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr);
65   for (p = pStart; p < pEnd; ++p) {
66     PetscInt dof;
67 
68     ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr);
69     if ((blockSize < 0) && (dof > 0)) blockSize = dof;
70     if ((dof > 0) && (dof != blockSize)) {
71       blockSize = 1;
72       break;
73     }
74   }
75   ierr = PetscSectionGetStorageSize(section, &localSize);CHKERRQ(ierr);
76   ierr = VecCreate(PETSC_COMM_SELF, vec);CHKERRQ(ierr);
77   ierr = VecSetSizes(*vec, localSize, localSize);CHKERRQ(ierr);
78   ierr = VecSetBlockSize(*vec, blockSize);CHKERRQ(ierr);
79   ierr = VecSetType(*vec,dm->vectype);CHKERRQ(ierr);
80   ierr = VecSetDM(*vec, dm);CHKERRQ(ierr);
81   PetscFunctionReturn(0);
82 }
83 
84 /* This assumes that the DM has been cloned prior to the call */
85 PetscErrorCode DMCreateSubDM_Section_Private(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
86 {
87   PetscSection   section, sectionGlobal;
88   PetscInt      *subIndices;
89   PetscInt       subSize = 0, subOff = 0, nF, f, pStart, pEnd, p;
90   PetscErrorCode ierr;
91 
92   PetscFunctionBegin;
93   if (!numFields) PetscFunctionReturn(0);
94   ierr = DMGetDefaultSection(dm, &section);CHKERRQ(ierr);
95   ierr = DMGetDefaultGlobalSection(dm, &sectionGlobal);CHKERRQ(ierr);
96   if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields");
97   if (!sectionGlobal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields");
98   ierr = PetscSectionGetNumFields(section, &nF);CHKERRQ(ierr);
99   if (numFields > nF) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Number of requested fields %d greater than number of DM fields %d", numFields, nF);
100   if (is) {
101     PetscInt bs = -1, bsLocal, bsMax, bsMin;
102     ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr);
103     for (p = pStart; p < pEnd; ++p) {
104       PetscInt gdof, pSubSize  = 0;
105 
106       ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr);
107       if (gdof > 0) {
108         for (f = 0; f < numFields; ++f) {
109           PetscInt fdof, fcdof;
110 
111           ierr     = PetscSectionGetFieldDof(section, p, fields[f], &fdof);CHKERRQ(ierr);
112           ierr     = PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);CHKERRQ(ierr);
113           pSubSize += fdof-fcdof;
114         }
115         subSize += pSubSize;
116         if (pSubSize) {
117           if (bs < 0) {
118             bs = pSubSize;
119           } else if (bs != pSubSize) {
120             /* Layout does not admit a pointwise block size */
121             bs = 1;
122           }
123         }
124       }
125     }
126     /* Must have same blocksize on all procs (some might have no points) */
127     bsLocal = bs;
128     ierr = MPIU_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
129     bsLocal = bs < 0 ? bsMax : bs;
130     ierr = MPIU_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
131     if (bsMin != bsMax) {
132       bs = 1;
133     } else {
134       bs = bsMax;
135     }
136     ierr = PetscMalloc1(subSize, &subIndices);CHKERRQ(ierr);
137     for (p = pStart; p < pEnd; ++p) {
138       PetscInt gdof, goff;
139 
140       ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr);
141       if (gdof > 0) {
142         ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
143         for (f = 0; f < numFields; ++f) {
144           PetscInt fdof, fcdof, fc, f2, poff = 0;
145 
146           /* Can get rid of this loop by storing field information in the global section */
147           for (f2 = 0; f2 < fields[f]; ++f2) {
148             ierr  = PetscSectionGetFieldDof(section, p, f2, &fdof);CHKERRQ(ierr);
149             ierr  = PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof);CHKERRQ(ierr);
150             poff += fdof-fcdof;
151           }
152           ierr = PetscSectionGetFieldDof(section, p, fields[f], &fdof);CHKERRQ(ierr);
153           ierr = PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);CHKERRQ(ierr);
154           for (fc = 0; fc < fdof-fcdof; ++fc, ++subOff) {
155             subIndices[subOff] = goff+poff+fc;
156           }
157         }
158       }
159     }
160     ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is);CHKERRQ(ierr);
161     ierr = ISSetBlockSize(*is, bs);CHKERRQ(ierr);
162   }
163   if (subdm) {
164     PetscSection subsection;
165     PetscBool    haveNull = PETSC_FALSE;
166     PetscInt     f, nf = 0;
167 
168     ierr = PetscSectionCreateSubsection(section, numFields, fields, &subsection);CHKERRQ(ierr);
169     ierr = DMSetDefaultSection(*subdm, subsection);CHKERRQ(ierr);
170     ierr = PetscSectionDestroy(&subsection);CHKERRQ(ierr);
171     for (f = 0; f < numFields; ++f) {
172       (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]];
173       if ((*subdm)->nullspaceConstructors[f]) {
174         haveNull = PETSC_TRUE;
175         nf       = f;
176       }
177     }
178     if (haveNull && is) {
179       MatNullSpace nullSpace;
180 
181       ierr = (*(*subdm)->nullspaceConstructors[nf])(*subdm, nf, &nullSpace);CHKERRQ(ierr);
182       ierr = PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) nullSpace);CHKERRQ(ierr);
183       ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr);
184     }
185     if (dm->prob) {
186       PetscInt Nf;
187 
188       ierr = PetscDSGetNumFields(dm->prob, &Nf);CHKERRQ(ierr);
189       if (nF != Nf) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "The number of DM fields %d does not match the number of Section fields %d", Nf, nF);
190       ierr = DMSetNumFields(*subdm, numFields);CHKERRQ(ierr);
191       for (f = 0; f < numFields; ++f) {
192         PetscObject disc;
193 
194         ierr = DMGetField(dm, fields[f], &disc);CHKERRQ(ierr);
195         ierr = DMSetField(*subdm, f, disc);CHKERRQ(ierr);
196       }
197       if (numFields == 1 && is) {
198         PetscObject disc, space, pmat;
199 
200         ierr = DMGetField(*subdm, 0, &disc);CHKERRQ(ierr);
201         ierr = PetscObjectQuery(disc, "nullspace", &space);CHKERRQ(ierr);
202         if (space) {ierr = PetscObjectCompose((PetscObject) *is, "nullspace", space);CHKERRQ(ierr);}
203         ierr = PetscObjectQuery(disc, "nearnullspace", &space);CHKERRQ(ierr);
204         if (space) {ierr = PetscObjectCompose((PetscObject) *is, "nearnullspace", space);CHKERRQ(ierr);}
205         ierr = PetscObjectQuery(disc, "pmat", &pmat);CHKERRQ(ierr);
206         if (pmat) {ierr = PetscObjectCompose((PetscObject) *is, "pmat", pmat);CHKERRQ(ierr);}
207       }
208     }
209   }
210 #if 0
211   /* We need a way to filter the original SF for given fields:
212        - Keeping the original section around it too much I think
213        - We could keep the distributed section, and subset it
214    */
215   if (dm->sfNatural) {
216     PetscSF sfNatural;
217 
218     ierr = PetscSectionCreateSubsection(dm->originalSection, numFields, fields, &(*subdm)->originalSection);CHKERRQ(ierr);
219     ierr = DMPlexCreateGlobalToNaturalPetscSF(*subdm, &sfNatural);CHKERRQ(ierr);
220     ierr = DMPlexSetGlobalToNaturalPetscSF(*subdm, sfNatural);CHKERRQ(ierr);
221   }
222 #endif
223   PetscFunctionReturn(0);
224 }
225