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 12 PetscFunctionBegin; 13 ierr = DMGetDefaultGlobalSection(dm, &gSection);CHKERRQ(ierr); 14 ierr = PetscSectionGetChart(gSection, &pStart, &pEnd);CHKERRQ(ierr); 15 for (p = pStart; p < pEnd; ++p) { 16 PetscInt dof, cdof; 17 18 ierr = PetscSectionGetDof(gSection, p, &dof);CHKERRQ(ierr); 19 ierr = PetscSectionGetConstraintDof(gSection, p, &cdof);CHKERRQ(ierr); 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 ierr = MPI_Allreduce(&blockSize, &bs, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 28 if (blockSize == PETSC_MAX_INT) blockSize = 1; /* Everyone was empty */ 29 ierr = PetscSectionGetConstrainedStorageSize(gSection, &localSize);CHKERRQ(ierr); 30 if (localSize%blockSize) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mismatch between blocksize %d and local storage size %d", blockSize, localSize); 31 ierr = VecCreate(PetscObjectComm((PetscObject)dm), vec);CHKERRQ(ierr); 32 ierr = VecSetSizes(*vec, localSize, PETSC_DETERMINE);CHKERRQ(ierr); 33 ierr = VecSetBlockSize(*vec, bs);CHKERRQ(ierr); 34 ierr = VecSetType(*vec,dm->vectype);CHKERRQ(ierr); 35 ierr = VecSetDM(*vec, dm);CHKERRQ(ierr); 36 /* ierr = VecSetLocalToGlobalMapping(*vec, dm->ltogmap);CHKERRQ(ierr); */ 37 PetscFunctionReturn(0); 38 } 39 40 #undef __FUNCT__ 41 #define __FUNCT__ "DMCreateLocalVector_Section_Private" 42 PetscErrorCode DMCreateLocalVector_Section_Private(DM dm,Vec *vec) 43 { 44 PetscSection section; 45 PetscInt localSize, blockSize = -1, pStart, pEnd, p; 46 PetscErrorCode ierr; 47 48 PetscFunctionBegin; 49 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 50 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 51 for (p = pStart; p < pEnd; ++p) { 52 PetscInt dof; 53 54 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 55 if ((blockSize < 0) && (dof > 0)) blockSize = dof; 56 if ((dof > 0) && (dof != blockSize)) { 57 blockSize = 1; 58 break; 59 } 60 } 61 ierr = PetscSectionGetStorageSize(section, &localSize);CHKERRQ(ierr); 62 ierr = VecCreate(PETSC_COMM_SELF, vec);CHKERRQ(ierr); 63 ierr = VecSetSizes(*vec, localSize, localSize);CHKERRQ(ierr); 64 ierr = VecSetBlockSize(*vec, blockSize);CHKERRQ(ierr); 65 ierr = VecSetType(*vec,dm->vectype);CHKERRQ(ierr); 66 ierr = VecSetDM(*vec, dm);CHKERRQ(ierr); 67 PetscFunctionReturn(0); 68 } 69 70 #undef __FUNCT__ 71 #define __FUNCT__ "DMCreateSubDM_Section_Private" 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; 78 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 if (!numFields) PetscFunctionReturn(0); 82 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 83 ierr = DMGetDefaultGlobalSection(dm, §ionGlobal);CHKERRQ(ierr); 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 ierr = PetscSectionGetNumFields(section, &nF);CHKERRQ(ierr); 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 ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr); 90 for (p = pStart; p < pEnd; ++p) { 91 PetscInt gdof; 92 93 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); 94 if (gdof > 0) { 95 for (f = 0; f < numFields; ++f) { 96 PetscInt fdof, fcdof; 97 98 ierr = PetscSectionGetFieldDof(section, p, fields[f], &fdof);CHKERRQ(ierr); 99 ierr = PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);CHKERRQ(ierr); 100 subSize += fdof-fcdof; 101 } 102 } 103 } 104 ierr = PetscMalloc1(subSize, &subIndices);CHKERRQ(ierr); 105 for (p = pStart; p < pEnd; ++p) { 106 PetscInt gdof, goff; 107 108 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); 109 if (gdof > 0) { 110 ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr); 111 for (f = 0; f < numFields; ++f) { 112 PetscInt fdof, fcdof, fc, f2, poff = 0; 113 114 /* Can get rid of this loop by storing field information in the global section */ 115 for (f2 = 0; f2 < fields[f]; ++f2) { 116 ierr = PetscSectionGetFieldDof(section, p, f2, &fdof);CHKERRQ(ierr); 117 ierr = PetscSectionGetFieldConstraintDof(section, p, f2, &fcdof);CHKERRQ(ierr); 118 poff += fdof-fcdof; 119 } 120 ierr = PetscSectionGetFieldDof(section, p, fields[f], &fdof);CHKERRQ(ierr); 121 ierr = PetscSectionGetFieldConstraintDof(section, p, fields[f], &fcdof);CHKERRQ(ierr); 122 for (fc = 0; fc < fdof-fcdof; ++fc, ++subOff) { 123 subIndices[subOff] = goff+poff+fc; 124 } 125 } 126 } 127 } 128 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm), subSize, subIndices, PETSC_OWN_POINTER, is);CHKERRQ(ierr); 129 } 130 if (subdm) { 131 PetscSection subsection; 132 PetscBool haveNull = PETSC_FALSE; 133 PetscInt f, nf = 0; 134 135 ierr = PetscSectionCreateSubsection(section, numFields, fields, &subsection);CHKERRQ(ierr); 136 ierr = DMSetDefaultSection(*subdm, subsection);CHKERRQ(ierr); 137 ierr = PetscSectionDestroy(&subsection);CHKERRQ(ierr); 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; 147 148 ierr = (*(*subdm)->nullspaceConstructors[nf])(*subdm, nf, &nullSpace);CHKERRQ(ierr); 149 ierr = PetscObjectCompose((PetscObject) *is, "nullspace", (PetscObject) nullSpace);CHKERRQ(ierr); 150 ierr = MatNullSpaceDestroy(&nullSpace);CHKERRQ(ierr); 151 } 152 if (dm->prob) { 153 PetscInt Nf; 154 155 ierr = PetscDSGetNumFields(dm->prob, &Nf);CHKERRQ(ierr); 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 ierr = DMSetNumFields(*subdm, numFields);CHKERRQ(ierr); 158 for (f = 0; f < numFields; ++f) { 159 PetscObject disc; 160 161 ierr = DMGetField(dm, fields[f], &disc);CHKERRQ(ierr); 162 ierr = DMSetField(*subdm, f, disc);CHKERRQ(ierr); 163 } 164 if (numFields == 1 && is) { 165 PetscObject disc, space, pmat; 166 167 ierr = DMGetField(*subdm, 0, &disc);CHKERRQ(ierr); 168 ierr = PetscObjectQuery(disc, "nullspace", &space);CHKERRQ(ierr); 169 if (space) {ierr = PetscObjectCompose((PetscObject) *is, "nullspace", space);CHKERRQ(ierr);} 170 ierr = PetscObjectQuery(disc, "nearnullspace", &space);CHKERRQ(ierr); 171 if (space) {ierr = PetscObjectCompose((PetscObject) *is, "nearnullspace", space);CHKERRQ(ierr);} 172 ierr = PetscObjectQuery(disc, "pmat", &pmat);CHKERRQ(ierr); 173 if (pmat) {ierr = PetscObjectCompose((PetscObject) *is, "pmat", pmat);CHKERRQ(ierr);} 174 } 175 } 176 } 177 PetscFunctionReturn(0); 178 } 179