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