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:   MPI_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, &section);
 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, &section);
 83:   DMGetDefaultGlobalSection(dm, &sectionGlobal);
 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:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
 90:     for (p = pStart; p < pEnd; ++p) {
 91:       PetscInt gdof;

 93:       PetscSectionGetDof(sectionGlobal, p, &gdof);
 94:       if (gdof > 0) {
 95:         for (f = 0; f < numFields; ++f) {
 96:           PetscInt fdof, fcdof;

 98:           PetscSectionGetFieldDof(section, p, fields[f], &fdof);
 99:           PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);
100:           subSize += fdof-fcdof;
101:         }
102:       }
103:     }
104:     PetscMalloc1(subSize, &subIndices);
105:     for (p = pStart; p < pEnd; ++p) {
106:       PetscInt gdof, goff;

108:       PetscSectionGetDof(sectionGlobal, p, &gdof);
109:       if (gdof > 0) {
110:         PetscSectionGetOffset(sectionGlobal, p, &goff);
111:         for (f = 0; f < numFields; ++f) {
112:           PetscInt fdof, fcdof, fc, f2, poff = 0;

114:           /* Can get rid of this loop by storing field information in the global section */
115:           for (f2 = 0; f2 < fields[f]; ++f2) {
116:             PetscSectionGetFieldDof(section, p, f2, &fdof);
117:             PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof);
118:             poff += fdof-fcdof;
119:           }
120:           PetscSectionGetFieldDof(section, p, fields[f], &fdof);
121:           PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);
122:           for (fc = 0; fc < fdof-fcdof; ++fc, ++subOff) {
123:             subIndices[subOff] = goff+poff+fc;
124:           }
125:         }
126:       }
127:     }
128:     ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is);
129:   }
130:   if (subdm) {
131:     PetscSection subsection;
132:     PetscBool    haveNull = PETSC_FALSE;
133:     PetscInt     f, nf = 0;

135:     PetscSectionCreateSubsection(section, numFields, fields, &subsection);
136:     DMSetDefaultSection(*subdm, subsection);
137:     PetscSectionDestroy(&subsection);
138:     for (f = 0; f < numFields; ++f) {
139:       (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]];
140:       if ((*subdm)->nullspaceConstructors[f]) {
141:         haveNull = PETSC_TRUE;
142:         nf       = f;
143:       }
144:     }
145:     if (haveNull && is) {
146:       MatNullSpace nullSpace;

148:       (*(*subdm)->nullspaceConstructors[nf])(*subdm, nf, &nullSpace);
149:       PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) nullSpace);
150:       MatNullSpaceDestroy(&nullSpace);
151:     }
152:     if (dm->prob) {
153:       PetscInt Nf;

155:       PetscDSGetNumFields(dm->prob, &Nf);
156:       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);
157:       DMSetNumFields(*subdm, numFields);
158:       for (f = 0; f < numFields; ++f) {
159:         PetscObject disc;

161:         DMGetField(dm, fields[f], &disc);
162:         DMSetField(*subdm, f, disc);
163:       }
164:       if (numFields == 1 && is) {
165:         PetscObject disc, space, pmat;

167:         DMGetField(*subdm, 0, &disc);
168:         PetscObjectQuery(disc, "nullspace", &space);
169:         if (space) {PetscObjectCompose((PetscObject) *is, "nullspace", space);}
170:         PetscObjectQuery(disc, "nearnullspace", &space);
171:         if (space) {PetscObjectCompose((PetscObject) *is, "nearnullspace", space);}
172:         PetscObjectQuery(disc, "pmat", &pmat);
173:         if (pmat) {PetscObjectCompose((PetscObject) *is, "pmat", pmat);}
174:       }
175:     }
176:   }
177:   return(0);
178: }