> I have DMComplexGetDefaultSection(), which is the first thing I call when
> constructing the FieldSplit info.
Found it. So you want PCFieldSplit to troll through those inner loops of
/* Define fields */
PetscSection section, sectionGlobal;
PetscInt *fieldSizes, **fieldIndices;
PetscInt nF, f, pStart, pEnd, p;
ierr = DMComplexGetDefaultSection(pc->dm, §ion);CHKERRQ(ierr);
ierr = DMComplexGetDefaultGlobalSection(pc->dm,
§ionGlobal);CHKERRQ(ierr);
ierr = PetscSectionGetNumFields(section, &nF);CHKERRQ(ierr);
ierr = PetscMalloc2(nF,PetscInt,&fieldSizes,nF,PetscInt
*,&fieldIndices);CHKERRQ(ierr);
ierr = PetscSectionGetChart(sectionGlobal, &pStart,
&pEnd);CHKERRQ(ierr);
for(f = 0; f < nF; ++f) {
fieldSizes[f] = 0;
}
for(p = pStart; p < pEnd; ++p) {
PetscInt gdof;
ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr);
if (gdof > 0) {
for(f = 0; f < nF; ++f) {
PetscInt fcdof;
ierr = PetscSectionGetFieldConstraintDof(section, p, f,
&fcdof);CHKERRQ(ierr);
fieldSizes[f] += fcdof;
}
}
}
for(f = 0; f < nF; ++f) {
ierr = PetscMalloc(fieldSizes[f] * sizeof(PetscInt),
&fieldIndices[f]);CHKERRQ(ierr);
fieldSizes[f] = 0;
}
for(p = pStart; p < pEnd; ++p) {
PetscInt gdof, goff;
ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr);
if (gdof > 0) {
ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr);
for(f = 0; f < nF; ++f) {
PetscInt fcdof, fc;
ierr = PetscSectionGetFieldConstraintDof(section, p, f,
&fcdof);CHKERRQ(ierr);
for(fc = 0; fc < fcdof; ++fc, ++fieldSizes[f]) {
fieldIndices[f][fieldSizes[f]] = goff++;
}
}
}
}
for(f = 0; f < nF; ++f) {
IS field;
const char *fieldname;
ierr = PetscSectionGetFieldName(section, f, &fieldname);CHKERRQ(ierr);
ierr = ISCreateGeneral(((PetscObject) pc)->comm, fieldSizes[f],
fieldIndices[f], PETSC_OWN_POINTER, &field);CHKERRQ(ierr);
ierr = PetscPrintf(((PetscObject) pc)->comm, "Field %s\n",
fieldname);CHKERRQ(ierr);
ierr = ISView(field, PETSC_NULL);CHKERRQ(ierr);
ierr = PCFieldSplitSetIS(pc, fieldname, field);CHKERRQ(ierr);
ierr = ISDestroy(&field);CHKERRQ(ierr);
}
ierr = PetscFree2(fieldSizes,fieldIndices);CHKERRQ(ierr);
}
every time instead of call DMGetFieldsIS(dm&is); That's insane, just take the
crud above and call it DMGetFieldsIS_Complex() and we'll have PCFieldSplit
independent of any DM guts.
I don't understand what we are fighting about.
On Feb 24, 2012, at 3:49 PM, Matthew Knepley wrote:
>
> It is the correct level of granularity, especially fr FEM which PETSc has
> almost none of right now
> because we make it unnaturally hard due to our ridiculous infatuation with
> Finite Difference. That kind
> of thinking colors everything.
I never said that PetscSection doesn't belong anywhere. It just doesn't seem
to belong anywhere in PC/KSP land, it is unnecessarily low for that. Sure it
may be ok for the finite element assembly etc.
Barry
>
> Regardless, say the concept of PetscSection as a way to keep the information
> about fields is fine. PetscSection is a better representation of an array of
> ISs.
>
> 1) Imbedding these constraint thingys in it seems completely wrong.
> Information about the constraint things should just be in another
> PetscSection not inside the PetscSection that tells you the field division.
>
> I disagree since you need this information in every situation where you ask
> about sizes, etc. Again, there is only ONE
> problem in PETSc that does this, written by me. It is often useful to try
> something before trying to figure out what is the
> best way to do it.
>
> 2) One could then have DMGetPetscSection(DM,section) that tells you the
> decomposition, but since the solvers need DMs for the pieces (in
> PCFieldSplit) it would be DMGetSections(DM,section,DM *subdms) and not so
> different from my first email with PetscSection replacing my array of ISs. So
> why not fix IS to be a bit more dynamic and build there is no real need for
> PetscSection except as an array of ISs.
>
> I have DMComplexGetDefaultSection(), which is the first thing I call when
> constructing the FieldSplit info.
>
> Matt
>
> Jed is laughing his pants off at the reaction to his side comment in that
> email :-)
>
>
> Barry
>
> On Feb 24, 2012, at 3:32 PM, Matthew Knepley wrote:
>
> > On Fri, Feb 24, 2012 at 3:26 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> > On Feb 24, 2012, at 3:13 PM, Matthew Knepley wrote:
> >
> > > On Fri, Feb 24, 2012 at 3:05 PM, Barry Smith <bsmith at mcs.anl.gov>
> > > wrote:
> > >
> > > On Feb 24, 2012, at 2:41 PM, Matthew Knepley wrote:
> > >
> > > > On Fri, Feb 24, 2012 at 2:39 PM, Barry Smith <bsmith at mcs.anl.gov>
> > > > wrote:
> > > >
> > > > On Feb 24, 2012, at 2:23 PM, Jed Brown wrote:
> > > >
> > > > > We need a richer field interface. Some way for a client of DM to
> > > > > discover vector and tensor structure.
> > > >
> > > > No idea what a tensor is, but DMGetSubDMS(dm,PetscInt *cnt,IS
> > > > *ises,DM *dms); and DMPushSubDMType(DM,"name of a type of splitting
> > > > like u,v,p or velocity,p"); DMPopSubDMType(); and
> > > > DMGetAvailableSubDMTypes(DM,char **);
> > > >
> > > > Do you need more than that?
> > > >
> > > > Yuck Yuck Yuck. Why are we persisting in shoving all of the index stuff
> > > > into the DM interface. I think this is so
> > > > unwieldy and complicated. PetscSection is small, powerful, and does
> > > > everything suggested and more. I see
> > > > no argument in support of stuff like the above.
> > >
> > > I have no idea what a PetscSection is and
> > >
> > > PetscSection - This is a mapping from DMMESH points to sets of values,
> > > which is
> > > our presentation of a fibre bundle.
> > >
> > > doesn't help explain it. PetscSection is small? It has more methods
> > > than DM, that is not small.
> > >
> > > Don't make me return the count. It has fewer, and definitely fewer
> > > concepts. Most of those are
> > > simple get/set.
> > >
> > > PetscSectionGetFieldComponents() seems to be intuitive, cool. By why does
> > > it return an array of integers? The universal way in PETSc would have it
> > > returning an IS. Wait a minute, it returns one int, the number of
> > > components? Then why is called GetFieldComponents and not
> > > GetNumFieldComponents()
> > >
> > > It could be changed to GetNumFieldComponents.
> > >
> > > And what is all this constraint stuff? PetscSectionGetConstraintDof()
> > > Huh? It looks inside the "section" and gets an array of ints from a
> > > section inside the section?
> > >
> > > Yes. It keeps track of constrained dofs. People, who through laziness or
> > > willful blindness, refuse to acknowledge that some
> > > problems make more sense when constrained dofs ar eliminated.
> > >
> > > Here's my guess: a PetscSection is way of "chopping up" a canonical
> > > vector into a bunch of fields. Fine, but then why does it have a bunch of
> > > other stuff in it (like constraints, what the hell are constraints)? And
> > > why isn't a PetscSection just an array of IS that define the fields?
> > >
> > > Constraints are just that. IS is designed to be immutable. Building on it
> > > internally is not a good idea for this reason.
> > >
> > > The reason to shove all the index stuff into the DM interface is because
> > > generally one doesn't just want to know the entries of a field but one
> > > also wants to create vectors for those field variables only, create
> > > matrices, know something about the mesh for that field so one can
> > > generate interpolations and coarsen for multigrid. One cannot do any of
> > > those things with PetscSection.
> > >
> > > Yes you can.
> >
> > Hmm, I guess I haven't pulled recently. I don't see a
> > PetscSectionGetInterpolation() a PetscSectionCoarsen() a
> > PetscSectionGetGlobalVector() etc etc.
> >
> > >
> > > Please explain in freshman linear algebra what a PetscSection is.
> > >
> > > Exactly what it says in the help. It is a mapping from mesh points
> > > (vertices, edges, faces, cells) to dofs. How can this be easier?
> >
> > Say I want an IS that represents all the cell center dofs, how do I
> > obtain that from a PetscSection? I don't see the functionality for that
> > but it seems the most basic kind of thing needed.
> >
> > DMComplexGetHeight(dm, 0, &cStart, &cEnd); /* Get the range of mesh points
> > which are cells, height is like co-dimension */
> > numCellDof = 0;
> > for(c = cStart; c < cEnd; ++c) {
> > PetscSectionGetDof(section, c, &dof);
> > numCellDof += dof;
> > }
> > PetscMalloc(numCellDof * sizeof(PetscInt), &indices);
> > numCellDof = 0;
> > for(c = cStart; c < cEnd; ++c) {
> > PetscSectionGetDof(section, c, &dof);
> > PetscSectionGetOffset(section, c, &off);
> > for(d = 0; d < dof; ++d, ++numCellDof) {
> > indices[numCellDof] = off+d;
> > }
> > }
> >
> > Matt
> >
> >
> > Barry
> >
> > >
> > > Matt
> > >
> > >
> > > Barry
> > >
> > >
> > >
> > >
> > >
> > >
> > > >
> > > > Matt
> > > >
> > > > Barry
> > > >
> > > > Note that current DMDA's default subdms served up would be a DM for
> > > > each field consisting of one of the degrees of freedom per node and it
> > > > could also serve up DMs that represent any collections of those for
> > > > example dof 0 and 1 per node, then dof 2 per node. We'll need to come
> > > > up with some naming convention for these.
> > > >
> > > >
> > > > > Possibly also (eventually) a way to solve equations of state so that
> > > > > we can formulate reduced problems in non-conservative variables.
> > > > >
> > > > > On Feb 24, 2012 2:19 PM, "Matthew Knepley" <knepley at gmail.com>
> > > > > wrote:
> > > > > On Fri, Feb 24, 2012 at 1:58 PM, Jed Brown <jedbrown at mcs.anl.gov>
> > > > > wrote:
> > > > > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/518ff70e8a0a
> > > > >
> > > > > Matt, I want to get rid of these conditionals, not add more. We
> > > > > should have a DM base interface for getting fields on which to split,
> > > > > that common interface should _replace_ the DMComposite specialization.
> > > > >
> > > > > I agree. However, I won't put anything in until I do it by hand once.
> > > > >
> > > > > Also, I was tempted to just promote DMGetGlobalISes(), but its not
> > > > > quite right. I am planning
> > > > > to put the IS method in PetscSection. I am hoping eventually DMDA
> > > > > uses PetscSection for
> > > > > layout.
> > > > >
> > > > > Matt
> > > > >
> > > > > --
> > > > > What most experimenters take for granted before they begin their
> > > > > experiments is infinitely more interesting than any results to which
> > > > > their experiments lead.
> > > > > -- Norbert Wiener
> > > >
> > > >
> > > >
> > > >
> > > > --
> > > > What most experimenters take for granted before they begin their
> > > > experiments is infinitely more interesting than any results to which
> > > > their experiments lead.
> > > > -- Norbert Wiener
> > >
> > >
> > >
> > >
> > > --
> > > What most experimenters take for granted before they begin their
> > > experiments is infinitely more interesting than any results to which
> > > their experiments lead.
> > > -- Norbert Wiener
> >
> >
> >
> >
> > --
> > What most experimenters take for granted before they begin their
> > experiments is infinitely more interesting than any results to which their
> > experiments lead.
> > -- Norbert Wiener
>
>
>
>
> --
> What most experimenters take for granted before they begin their experiments
> is infinitely more interesting than any results to which their experiments
> lead.
> -- Norbert Wiener
