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