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