Actual source code: dmi.c
petsc-3.7.7 2017-09-25
1: #include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/
2: #include <petscds.h>
6: PetscErrorCode DMCreateGlobalVector_Section_Private(DM dm,Vec *vec)
7: {
8: PetscSection gSection;
9: PetscInt localSize, bs, blockSize = -1, pStart, pEnd, p;
11: PetscInt in[2],out[2];
14: DMGetDefaultGlobalSection(dm, &gSection);
15: PetscSectionGetChart(gSection, &pStart, &pEnd);
16: for (p = pStart; p < pEnd; ++p) {
17: PetscInt dof, cdof;
19: PetscSectionGetDof(gSection, p, &dof);
20: PetscSectionGetConstraintDof(gSection, p, &cdof);
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: }
34: in[0] = -blockSize;
35: in[1] = blockSize;
36: MPIU_Allreduce(in,out,2,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
37: /* -out[0] = min(blockSize), out[1] = max(blockSize) */
38: if (-out[0] == out[1]) {
39: bs = out[1];
40: } else bs = 1;
42: if (bs < 0) { /* Everyone was empty */
43: blockSize = 1;
44: bs = 1;
45: }
47: PetscSectionGetConstrainedStorageSize(gSection, &localSize);
48: if (localSize%blockSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mismatch between blocksize %d and local storage size %d", blockSize, localSize);
49: VecCreate(PetscObjectComm((PetscObject)dm), vec);
50: VecSetSizes(*vec, localSize, PETSC_DETERMINE);
51: VecSetBlockSize(*vec, bs);
52: VecSetType(*vec,dm->vectype);
53: VecSetDM(*vec, dm);
54: /* VecSetLocalToGlobalMapping(*vec, dm->ltogmap); */
55: return(0);
56: }
60: PetscErrorCode DMCreateLocalVector_Section_Private(DM dm,Vec *vec)
61: {
62: PetscSection section;
63: PetscInt localSize, blockSize = -1, pStart, pEnd, p;
67: DMGetDefaultSection(dm, §ion);
68: PetscSectionGetChart(section, &pStart, &pEnd);
69: for (p = pStart; p < pEnd; ++p) {
70: PetscInt dof;
72: PetscSectionGetDof(section, p, &dof);
73: if ((blockSize < 0) && (dof > 0)) blockSize = dof;
74: if ((dof > 0) && (dof != blockSize)) {
75: blockSize = 1;
76: break;
77: }
78: }
79: PetscSectionGetStorageSize(section, &localSize);
80: VecCreate(PETSC_COMM_SELF, vec);
81: VecSetSizes(*vec, localSize, localSize);
82: VecSetBlockSize(*vec, blockSize);
83: VecSetType(*vec,dm->vectype);
84: VecSetDM(*vec, dm);
85: return(0);
86: }
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;
99: if (!numFields) return(0);
100: DMGetDefaultSection(dm, §ion);
101: DMGetDefaultGlobalSection(dm, §ionGlobal);
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: PetscSectionGetNumFields(section, &nF);
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: PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
109: for (p = pStart; p < pEnd; ++p) {
110: PetscInt gdof, pSubSize = 0;
112: PetscSectionGetDof(sectionGlobal, p, &gdof);
113: if (gdof > 0) {
114: for (f = 0; f < numFields; ++f) {
115: PetscInt fdof, fcdof;
117: PetscSectionGetFieldDof(section, p, fields[f], &fdof);
118: PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);
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: MPIU_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
135: bsLocal = bs < 0 ? bsMax : bs;
136: MPIU_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
137: if (bsMin != bsMax) {
138: bs = 1;
139: } else {
140: bs = bsMax;
141: }
142: PetscMalloc1(subSize, &subIndices);
143: for (p = pStart; p < pEnd; ++p) {
144: PetscInt gdof, goff;
146: PetscSectionGetDof(sectionGlobal, p, &gdof);
147: if (gdof > 0) {
148: PetscSectionGetOffset(sectionGlobal, p, &goff);
149: for (f = 0; f < numFields; ++f) {
150: PetscInt fdof, fcdof, fc, f2, poff = 0;
152: /* Can get rid of this loop by storing field information in the global section */
153: for (f2 = 0; f2 < fields[f]; ++f2) {
154: PetscSectionGetFieldDof(section, p, f2, &fdof);
155: PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof);
156: poff += fdof-fcdof;
157: }
158: PetscSectionGetFieldDof(section, p, fields[f], &fdof);
159: PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);
160: for (fc = 0; fc < fdof-fcdof; ++fc, ++subOff) {
161: subIndices[subOff] = goff+poff+fc;
162: }
163: }
164: }
165: }
166: ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is);
167: ISSetBlockSize(*is, bs);
168: }
169: if (subdm) {
170: PetscSection subsection;
171: PetscBool haveNull = PETSC_FALSE;
172: PetscInt f, nf = 0;
174: PetscSectionCreateSubsection(section, numFields, fields, &subsection);
175: DMSetDefaultSection(*subdm, subsection);
176: PetscSectionDestroy(&subsection);
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;
187: (*(*subdm)->nullspaceConstructors[nf])(*subdm, nf, &nullSpace);
188: PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) nullSpace);
189: MatNullSpaceDestroy(&nullSpace);
190: }
191: if (dm->prob) {
192: PetscInt Nf;
194: PetscDSGetNumFields(dm->prob, &Nf);
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: DMSetNumFields(*subdm, numFields);
197: for (f = 0; f < numFields; ++f) {
198: PetscObject disc;
200: DMGetField(dm, fields[f], &disc);
201: DMSetField(*subdm, f, disc);
202: }
203: if (numFields == 1 && is) {
204: PetscObject disc, space, pmat;
206: DMGetField(*subdm, 0, &disc);
207: PetscObjectQuery(disc, "nullspace", &space);
208: if (space) {PetscObjectCompose((PetscObject) *is, "nullspace", space);}
209: PetscObjectQuery(disc, "nearnullspace", &space);
210: if (space) {PetscObjectCompose((PetscObject) *is, "nearnullspace", space);}
211: PetscObjectQuery(disc, "pmat", &pmat);
212: if (pmat) {PetscObjectCompose((PetscObject) *is, "pmat", pmat);}
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;
224: PetscSectionCreateSubsection(dm->originalSection, numFields, fields, &(*subdm)->originalSection);
225: DMPlexCreateGlobalToNaturalPetscSF(*subdm, &sfNatural);
226: DMPlexSetGlobalToNaturalPetscSF(*subdm, sfNatural);
227: }
228: #endif
229: return(0);
230: }