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