[petsc-dev] PetscSection

[Blaise A Bourdin]

On Nov 10, 2012, at 2:29 PM, Jed Brown mailto:jedbrown 
at mcs.anl.gov>>
 wrote:

On Sat, Nov 10, 2012 at 2:09 PM, Blaise A Bourdin mailto:bourdin at lsu.edu>> wrote:
Nope... and please, don't break that...

Here is an example of a mesh with 4 element blocks. In an analysis code, for 
each block , I would specify the element type, type of physics, coefficients 
etc in a petscBag. That said, I may be doing it wrong...

block 10 is made of tri
block 22 is made of quad
blocks 100 and 102 are made of bars

I'm not suggesting any demotion of this ability. I'm trying to understand what 
is the advantage of having a stratum with unsorted mixed dimension. If the 
blocks were sorted by dimension, you would still access them using a label 
(which gives you an index set, perhaps contiguous). You would still be able to 
do all the normal operations with those points. If you want to have one loop 
over all elements of the various types, you would union or concatenate the 
index sets of each type. Am I missing something?

No, you are absolutely right. I get your point now.

Blaise



--
Department of Mathematics and Center for Computation & Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin







-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/ae12d685/attachment.html>

[petsc-dev] PetscSection

[Blaise A Bourdin]

On Nov 10, 2012, at 1:54 PM, Jed Brown mailto:jedbrown 
at mcs.anl.gov>> wrote:

On Sat, Nov 10, 2012 at 1:44 PM, Matthew Knepley mailto:knepley at gmail.com>> wrote:
Yes, I agree that they are different geometrically, and in PyLith we
make Labels to
distinguish, which I think is the right way.

Cool. Of course you have to create the labels somehow. Do you currently make 
labels for all the "sidesets" in the mesh you load (e.g., from exodus)?

For exodus, I don't read side sets. The exodus interface for side set is really 
crappy, and does not allow an side set _inside_ the domain. Instead, I do not 
assume that all elements have the same codimension.


Labels in their current form seem to do two things (denote sets and associate 
values with points in those sets). Don't you frequently want one or the other?

> Fine, I care much more about conceptual simplicity in the interface. I think
> adoption will be higher if 95% of users don't see the word "stratum".

Note that this will have to be set explicitly somehow, since you cannot "tell"
from the DAG that something has a given dimension/co-dim.

Are your points not sorted by dimension? (Most of the code I've seen in the 
examples assumes homogeneous strata.)

Nope... and please, don't break that...

Here is an example of a mesh with 4 element blocks. In an analysis code, for 
each block , I would specify the element type, type of physics, coefficients 
etc in a petscBag. That said, I may be doing it wrong...

block 10 is made of tri
block 22 is made of quad
blocks 100 and 102 are made of bars
iMac:Exodus blaise$ ./test1 -i Square-mixed2.gen

dm:
Mesh in 2 dimensions:
  0-cells: 9
  2-cells: 12
[0]: Number of vertices in mesh: 9
[0]: Number of cells in mesh: 12
viewing Mesh 'mesh'
Max sizes cone: 4 support: 5
viewing IFSieve 'mesh sieve'
cap --> base:
[0]: 12 > 0
[0]: 12 > 1
[0]: 12 > 4
[0]: 12 > 3
[0]: 12 > 5
[0]: 13 > 6
[0]: 13 > 0
[0]: 13 > 7
[0]: 13 > 5
[0]: 14 > 1
[0]: 14 > 6
[0]: 14 > 0
[0]: 15 > 2
[0]: 15 > 3
[0]: 15 > 1
[0]: 16 > 11
[0]: 16 > 2
[0]: 17 > 10
[0]: 17 > 3
[0]: 17 > 4
[0]: 17 > 11
[0]: 17 > 2
[0]: 18 > 8
[0]: 18 > 10
[0]: 18 > 4
[0]: 19 > 5
[0]: 19 > 9
[0]: 19 > 8
[0]: 19 > 4
[0]: 20 > 7
[0]: 20 > 5
[0]: 20 > 9
base <-- cap:
[0]: 0<12
[0]: 0<13
[0]: 0<14
[0]: 1<12
[0]: 1<14
[0]: 1<15
[0]: 2<15
[0]: 2<16
[0]: 2<17
[0]: 3<15
[0]: 3<17
[0]: 3<12
[0]: 4<12
[0]: 4<17
[0]: 4<18
[0]: 4<19
[0]: 5<19
[0]: 5<20
[0]: 5<13
[0]: 5<12
[0]: 6<13
[0]: 6<14
[0]: 7<20
[0]: 7<13
[0]: 8<18
[0]: 8<19
[0]: 9<19
[0]: 9<20
[0]: 10<17
[0]: 10<18
[0]: 11<16
[0]: 11<17
Orientation:
[0]: 0<12: 0
[0]: 0<13: 0
[0]: 0<14: 0
[0]: 1<12: 0
[0]: 1<14: 0
[0]: 1<15: 0
[0]: 2<15: 0
[0]: 2<16: 0
[0]: 2<17: 0
[0]: 3<15: 0
[0]: 3<17: 0
[0]: 3<12: 0
[0]: 4<12: 0
[0]: 4<17: 0
[0]: 4<18: 0
[0]: 4<19: 0
[0]: 5<19: 0
[0]: 5<20: 0
[0]: 5<13: 0
[0]: 5<12: 0
[0]: 6<13: 0
[0]: 6<14: 0
[0]: 7<20: 0
[0]: 7<13: 0
[0]: 8<18: 0
[0]: 8<19: 0
[0]: 9<19: 0
[0]: 9<20: 0
[0]: 10<17: 0
[0]: 10<18: 0
[0]: 11<16: 0
[0]: 11<17: 0
viewing GeneralSection 'coordinates'
  Fields: 0
[0]:   12 dim 2 offset 0   0 0
[0]:   13 dim 2 offset 2   -0.5 0
[0]:   14 dim 2 offset 4   -0.5 -0.5
[0]:   15 dim 2 offset 6   0 -0.5
[0]:   16 dim 2 offset 8   0.5 -0.5
[0]:   17 dim 2 offset 10   0.5 0
[0]:   18 dim 2 offset 12   0.5 0.5
[0]:   19 dim 2 offset 14   0 0.5
[0]:   20 dim 2 offset 16   -0.5 0.5
viewing LabelSifter: 'Cell sets'
cap --> base:
[0]: 10>0
[0]: 10>1
[0]: 10>2
[0]: 10>3
[0]: 22>4
[0]: 22>5
[0]: 100>6
[0]: 100>7
[0]: 102>8
[0]: 102>9
[0]: 102>10
[0]: 102>11
viewing LabelSifter: 'Vertex sets'
cap --> base:
[0]: 1>16
[0]: 1>17
[0]: 1>18


--
Department of Mathematics and Center for Computation & Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin







-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/41b921c0/attachment-0001.html>
-- next part --
A non-text attachment was scrubbed...
Name: test1.c
Type: application/octet-stream
Size: 2711 bytes
Desc: test1.c
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/41b921c0/attachment-0002.obj>
-- next part --
A non-text attachment was scrubbed...
Name: Square-mixed2.gen
Type: application/octet-stream
Size: 2572 bytes
Desc: Square-mixed2.gen
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/41b921c0/attachment-0003.obj>

[petsc-dev] PetscSection

[Blaise A Bourdin]

On Nov 10, 2012, at 1:44 PM, Matthew Knepley 
 wrote:

> On Sat, Nov 10, 2012 at 2:33 PM, Jed Brown  wrote:
>> On Sat, Nov 10, 2012 at 1:26 PM, Matthew Knepley  
>> wrote:
>>> 
>>> I do not understand. Codim 0 and 1 are handled in an identical way.
>> 
>> 
>> The user doesn't do the same thing with them. Your data structure treats
>> them the same way (though I'm not sure how you distinguish between a tet and
>> a quad if both are "height 0" and you don't store intermediate levels), but
>> that seems like an implementation detail to me.
> 
> Yes, I agree that they are different geometrically, and in PyLith we
> make Labels to
> distinguish, which I think is the right way.

I actually disagree they may or may not be treated different. It all depends on 
what you want to do. If all you are interested in is computing integrals, 
elements are elements are elements. There is no fundamental difference between 
computing integrals over surfaces or volumes. And the user already needs to 
keep track of the polynomial order, interpolation type, integration type etc on 
a per element or per block of elements basis, right? If you cache the values of 
your FE basis functions at integration points and integration weights, 
integrating over a surface or an element is _exactly_ the same thing.

> 
>>> 
>>> I am not against some convenience interface for dim/co-dim, however it
>>> would
>>> just be implemented by making a label, since we have many equivalent
>>> structural
>>> things which would map to it.
>> 
>> 
>> Fine, I care much more about conceptual simplicity in the interface. I think
>> adoption will be higher if 95% of users don't see the word "stratum".
> 
> Note that this will have to be set explicitly somehow, since you cannot "tell"
> from the DAG that something has a given dimension/co-dim.

You can also not say anything about the element type from the DAG either. Is 
this a problem? Of course no. My point is that if a mesh has cells with 
different codimension, this information is something in the description from 
which the DM was created, and was probably imported somewhere in an IS, right?

Blaise 


-- 
Department of Mathematics and Center for Computation & Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-dev] PetscSection

[Matthew Knepley]

On Sat, Nov 10, 2012 at 7:14 PM, Jed Brown  wrote:
> On Sat, Nov 10, 2012 at 5:31 PM, Matthew Knepley  wrote:
>>
>> I was being too pessimistic. You can assign an intermediate node to any
>> set
>> of shared vertices, but I am not sure whether you start from the top or
>> bottom,
>> and the hard part is preallocation.
>
>
> Hmm, I'm not sure how to do hanging nodes, but for a conforming mesh
> starting with cell-to-vertex c2v[...]
>
> # build inverse map
> v2c = defaultdict(lambda v:[])
> for c in cells:
>   for v in c2v[c]:
> v2c[v].append(c)
>
> # find common subsets
> edges = set()
> c2e = defaultdict(lambda c:set())
> for v in vertices:
>   for (v1,c1,v2,c2) in sorted_vertices_appearing_in_exactly_2_cells(v2c[v]);
> edges.add((v1,v2))
> c2e[c1].add((v1,v2))
> c2e[c2].add((v1,v2))
>
> # repeat using edges to bound cell-to-face
> ...
>
>
>
> The problem with this is that you don't automatically get consistent
> orientation. If you need to work with basis functions, you need to orient,
> in which case the ordering produced by the procedure above needs to be
> fixed.
>
> In my opinion, that orientation is critical, such that connectivity without
> topology is of limited value. The topological dimension and number of
> vertices is just not sufficient to identify the topology of an element. For
> example, a hex can be collapsed to 7 vertices via an edge or across the
> middle of a face.

Definitely, that is much harder than the topology, and where I need specific
stuff for non-simplices, but it might be possible to use the initial orientation
of vertices.

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

[petsc-dev] PetscSection

[Matthew Knepley]

On Sat, Nov 10, 2012 at 6:11 PM, Jed Brown  wrote:
> On Sat, Nov 10, 2012 at 5:03 PM, Matthew Knepley  wrote:
>>
>> Interpolation. I do not have a good, graph theoretic way to phrase it.
>
>
> So if you know the topology of the n-cell, it implies a structure to its
> cone-closure. (You could store that as a cache of single-cell-rooted
> complex.) To interpolate the mesh, you graft in the single-cell, matching
> the intermediate levels (fast query with a vertex index). I've written some
> code for this, though not in full generality (I didn't explicitly store the
> reference single-cells).

I was being too pessimistic. You can assign an intermediate node to any set
of shared vertices, but I am not sure whether you start from the top or bottom,
and the hard part is preallocation.

  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

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 5:31 PM, Matthew Knepley  wrote:

> I was being too pessimistic. You can assign an intermediate node to any set
> of shared vertices, but I am not sure whether you start from the top or
> bottom,
> and the hard part is preallocation.
>

Hmm, I'm not sure how to do hanging nodes, but for a conforming mesh
starting with cell-to-vertex c2v[...]

# build inverse map
v2c = defaultdict(lambda v:[])
for c in cells:
  for v in c2v[c]:
v2c[v].append(c)

# find common subsets
edges = set()
c2e = defaultdict(lambda c:set())
for v in vertices:
  for (v1,c1,v2,c2) in sorted_vertices_appearing_in_exactly_2_cells(v2c[v]);
edges.add((v1,v2))
c2e[c1].add((v1,v2))
c2e[c2].add((v1,v2))

# repeat using edges to bound cell-to-face
...



The problem with this is that you don't automatically get consistent
orientation. If you need to work with basis functions, you need to orient,
in which case the ordering produced by the procedure above needs to be
fixed.

In my opinion, that orientation is critical, such that connectivity without
topology is of limited value. The topological dimension and number of
vertices is just not sufficient to identify the topology of an element. For
example, a hex can be collapsed to 7 vertices via an edge or across the
middle of a face.
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/c409f749/attachment.html>

[petsc-dev] PetscSection

[Matthew Knepley]

On Sat, Nov 10, 2012 at 5:55 PM, Jed Brown  wrote:
> On Sat, Nov 10, 2012 at 4:51 PM, Matthew Knepley  wrote:
>>
>>
>> The is one operation that I do not yet know how to do outside of this.
>> This may
>> be ignorance.
>
>
> What operation?

Interpolation. I do not have a good, graph theoretic way to phrase it.

  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

[petsc-dev] PetscSection

[Blaise A Bourdin]

rdin







-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/2f42d601/attachment.html>

[petsc-dev] PetscSection

[Matthew Knepley]

On Sat, Nov 10, 2012 at 5:01 PM, Jed Brown  wrote:
> On Sat, Nov 10, 2012 at 3:49 PM, Matthew Knepley  wrote:
>>
>> I am not sure I understand. Labels just define point sets. If you want
>> to segregate point
>> sets by dimension, that is fine. Nothing in the interface is sensitive
>> to this, nor should it be.
>
>
> I can think of only one way in which numbering points by stratum versus
> dimension makes a tangible difference: whether contiguous ranges can contain
> unsorted mixed-dimensional points. I want to understand whether this is
> actually important and whether there is a different use case that is
> important. If not, then we could entirely discard the concept of stratum,
> sort points by co/dimension, and query based on labels.
>
> Where is stratum as a concept irreplaceable?

I cannot say off the top of my head. It is very convenient to have the
breadth-first
levels of the graph done this way, and it is also nice for loops, but
I know you are
arguing that people more often write loops over dimension.

>>
>> >>
>> >> You don't need coordinates for interpolation?
>> >
>> >
>> > Are we talking about the same interpolation? If I have cell-to-vertex
>> > connectivity, I can create the faces without coordinates, yes.
>>
>> Ah, we are not. Mesh interpolation.
>>
>> The way I have done interpolation, you have to know the mesh dimension,
>> DMComplexGetDimension(), and then it assumes that height 0 stuff is cells
>> of that dimension.
>
>
> Then we're basically identifying stratum with dimension, suggesting that we
> should be able to remove stratum from the API in favor of co/dimension.

The is one operation that I do not yet know how to do outside of this. This may
be ignorance.

   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

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 5:03 PM, Matthew Knepley  wrote:

> Interpolation. I do not have a good, graph theoretic way to phrase it.


So if you know the topology of the n-cell, it implies a structure to its
cone-closure. (You could store that as a cache of single-cell-rooted
complex.) To interpolate the mesh, you graft in the single-cell, matching
the intermediate levels (fast query with a vertex index). I've written some
code for this, though not in full generality (I didn't explicitly store the
reference single-cells).
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/7d5a0c1e/attachment.html>

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 4:51 PM, Matthew Knepley  wrote:

>
> The is one operation that I do not yet know how to do outside of this.
> This may
> be ignorance.
>

What operation?
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/d5dcb958/attachment-0001.html>

[petsc-dev] PetscSection

[Matthew Knepley]

On Sat, Nov 10, 2012 at 4:30 PM, Jed Brown  wrote:
> On Sat, Nov 10, 2012 at 3:25 PM, Matthew Knepley  wrote:
>>
>> Everything is just a point, so I get back mixed dimension things all the
>> time. I
>> usually separate by stratum.
>
>
> Okay, but what is the compelling use case for DMComplexGetStratumIS()
> returning mixed-dimension that cannot be reasonably accomplished by getting
> dimensional ISs and union/concatenating them?

I am not sure I understand. Labels just define point sets. If you want
to segregate point
sets by dimension, that is fine. Nothing in the interface is sensitive
to this, nor should it be.

>>
>> You don't need coordinates for interpolation?
>
>
> Are we talking about the same interpolation? If I have cell-to-vertex
> connectivity, I can create the faces without coordinates, yes.

Ah, we are not. Mesh interpolation.

The way I have done interpolation, you have to know the mesh dimension,
DMComplexGetDimension(), and then it assumes that height 0 stuff is cells
of that dimension.

   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

[petsc-dev] PetscSection

[Blaise A Bourdin]

Sorry, for the empty message earlier. Fat fingers on a small keyboard...


On Fri, Nov 9, 2012 at 10:20 PM, Blaise A Bourdin mailto:bourdin at lsu.edu>> wrote:
A DM does double duty by describing the geometry of the mesh, and the data 
layout associated with the finite element space. I liked the model where the 
mesh geometry and the data layout on the mesh was split in two objects, but I 
understand the convenience of having everything in the DM, and DMClone works 
just fine. Since I may have to handle scalar, vector, 2nd and 4th order tensor 
on 2 different finite element spaces, in an assembly loop, I may end up dealing 
with as many as 8 DM. I stick all these DM's in the user context of each DM 
associated with a unknown (a Vec on which I may have to call SNESSolve or 
TSSolve), hoping that this is not creating some aliasing problem which as  a 
fortran programmer I can not possibly understand.

;-)

Actually, I am really not sure if passing a dm and a pointer to a user context 
containing this dm is legit or not...

Everything is an IS. I suspect that this means that the Set/GetCone,

I think these just return a range or (offset,size).

Set/GetClosure operations

Closure needs to do orientation, so it's not direct access from the Vec. 
Regardless, I don't think the plan would be to create an IS for that 
granularity.
Fine.  It should be done from the DM anyways.



would return an IS that could be used in VecSetValuesIS, but may not even be 
needed if an equivalent of  DMMesh' SectionRealRestrict / RestrictClosure 
/Update / UpdateCLosure is implemented

I think these are DMCompletVecGet/SetClosure.

Great




Setting Values is good, but getting values is needed too!

Also, in addition to the Barry's simple pseudo-code, I need to be able to get 
an IS for the i^th dof of a field at a given point (think of applying a 
Dirichlet boundary conditions to only one component of a field, for instance). 
It's always been the messy part. How would that fit within the proposed model?

I think you would get an index set for those points on the feature you wanted 
to do something special for.

That's pretty much what I do now. With DMMesh, it required quite a bit of work. 
It looks much easier with DMComplex

Finally, just a comment on stratum vs topological dimension. There is no reason 
why all elements in a mesh should have the same topological dimension. I 
actually find it much easier to have a single mesh where some sets of elements 
have codimension 0 and others codimension 1 and dispatching my assembly 
depending on the physics associated with each block, and the codimension of the 
element rather than having to deal with side sets / boundary mesh / whatever 
you want to call it.

I don't mind that, but can't you have an index set describing the codim 0 
elements (maybe all of them) and another index set for the codim 1 elements on 
the features you care about? You can take their union (or concatenate) for your 
assembly loop if you like. Is there something wrong with this approach?

Thats a very good point. In the end it doesn't really matter. As far as I 
remember, the main reason I ended with my current scheme is that DMMesh did not 
play well with partially interpolated meshes. I don't know what the current 
status of DMComplex is.

Blaise
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/6c3a76bf/attachment-0001.html>

[petsc-dev] PetscSection

[Matthew Knepley]

On Sat, Nov 10, 2012 at 4:13 PM, Jed Brown  wrote:
> On Sat, Nov 10, 2012 at 3:02 PM, Matthew Knepley  wrote:
>>
>>
>> > How does one get an IS containing _only_ the labeled points?
>>
>> DMComplexGetStratumIS()
>
> Great. Under what circumstances do you call this routine and get back an IS
> with mixed dimension?

Everything is just a point, so I get back mixed dimension things all the time. I
usually separate by stratum.

>>
>>
>> >>
>> >> > Can you answer how you can distinguish a quad from a tet in a
>> >> > non-interpolated mesh?
>> >>
>> >> You can't, but that is the point. You are not supposed to distinguish
>> >> them.
>> >
>> >
>> > Hmm, so these are the same at the category-theory level, but definitely
>> > not
>> > at the mesh topology level. It doesn't seem right to not distinguish
>> > them at
>> > all. How do you write DMComplexInterpolate for mixed-dim strata?
>>
>> Interpolation is analysis. You need to know what dimension you are in,
>> etc.
>
>
> But it's still topological rather than geometric. You don't need
> coordinates, for example.

You don't need coordinates for interpolation?

   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

[petsc-dev] PetscSection

[Blaise A Bourdin]

Sent from a handheld device.

On Nov 10, 2012, at 12:38 AM, "Jed Brown" mailto:jedbrown at mcs.anl.gov>> wrote:

On Fri, Nov 9, 2012 at 10:20 PM, Blaise A Bourdin mailto:bourdin at lsu.edu>> wrote:
A DM does double duty by describing the geometry of the mesh, and the data 
layout associated with the finite element space. I liked the model where the 
mesh geometry and the data layout on the mesh was split in two objects, but I 
understand the convenience of having everything in the DM, and DMClone works 
just fine. Since I may have to handle scalar, vector, 2nd and 4th order tensor 
on 2 different finite element spaces, in an assembly loop, I may end up dealing 
with as many as 8 DM. I stick all these DM's in the user context of each DM 
associated with a unknown (a Vec on which I may have to call SNESSolve or 
TSSolve), hoping that this is not creating some aliasing problem which as  a 
fortran programmer I can not possibly understand.

;-)

Everything is an IS. I suspect that this means that the Set/GetCone,

I think these just return a range or (offset,size).

Set/GetClosure operations

Closure needs to do orientation, so it's not direct access from the Vec. 
Regardless, I don't think the plan would be to create an IS for that 
granularity.

would return an IS that could be used in VecSetValuesIS, but may not even be 
needed if an equivalent of  DMMesh' SectionRealRestrict / RestrictClosure 
/Update / UpdateCLosure is implemented

I think these are DMCompletVecGet/SetClosure.


Setting Values is good, but getting values is needed too!

Also, in addition to the Barry's simple pseudo-code, I need to be able to get 
an IS for the i^th dof of a field at a given point (think of applying a 
Dirichlet boundary conditions to only one component of a field, for instance). 
It's always been the messy part. How would that fit within the proposed model?

I think you would get an index set for those points on the feature you wanted 
to do something special for.

Finally, just a comment on stratum vs topological dimension. There is no reason 
why all elements in a mesh should have the same topological dimension. I 
actually find it much easier to have a single mesh where some sets of elements 
have codimension 0 and others codimension 1 and dispatching my assembly 
depending on the physics associated with each block, and the codimension of the 
element rather than having to deal with side sets / boundary mesh / whatever 
you want to call it.

I don't mind that, but can't you have an index set describing the codim 0 
elements (maybe all of them) and another index set for the codim 1 elements on 
the features you care about? You can take their union (or concatenate) for your 
assembly loop if you like. Is there something wrong with this approach?
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/fb0588ca/attachment.html>

[petsc-dev] PetscSection

[Matthew Knepley]

On Sat, Nov 10, 2012 at 4:09 PM, Karl Rupp  wrote:
> On 11/10/2012 02:29 PM, Jed Brown wrote:
>>
>> On Sat, Nov 10, 2012 at 2:09 PM, Blaise A Bourdin > > wrote:
>>
>> Nope... and please, don't break that...
>>
>> Here is an example of a mesh with 4 element blocks. In an analysis
>> code, for each block , I would specify the element type, type of
>> physics, coefficients etc in a petscBag. That said, I may be doing
>> it wrong...
>>
>> block 10 is made of tri
>> block 22 is made of quad
>> blocks 100 and 102 are made of bars
>>
>>
>> I'm not suggesting any demotion of this ability. I'm trying to
>> understand what is the advantage of having a stratum with unsorted mixed
>> dimension. If the blocks were sorted by dimension, you would still
>> access them using a label (which gives you an index set, perhaps
>> contiguous). You would still be able to do all the normal operations
>> with those points. If you want to have one loop over all elements of the
>> various types, you would union or concatenate the index sets of each
>> type. Am I missing something?
>
>
> May I add something that is not specific to the current implementation in
> PETSc, yet still backs up the 'arrays of elements of the same type'
> approach:
>  a) When thinking of accelerators, having, say, all tets stored/referenced
> contiguously in one array and all hexahedra stored contiguously in another
> (separate) array, makes it a lot easier to move all tets to a kernel working
> on tets, and to feed all hexes to a second kernel working on hexes.
>
>  b) In a strongly typed setting (my experience here is in C++), one can do
> pretty nice things by keeping the two arrays separate rather than trying to
> unify them in some base class. For simplicity, suppose you have a mesh and
> two containers holding the tets and hexes. A FEM assembly could then look
> like this:

Yes, I agree completely with batching of this kind, which is why my integration
functions look like they do. However, the topology information is not
really needed
in these kernels. Instead you are packing discretization information, like
FEM coefficients, or geometric information like Jacobians, which are
all controlled
by PetscSection, not by the DMComplex.

   Matt

>   for_each(mesh_cells.begin(), mesh_cells.end(), assembly_functor());
>
> where assembly_functor statically dispatches into the cell type:
>   struct assembly_functor{
> void operator()(tetrahedron const & t) { ... }
> void operator()(hexahedron const & t) { ... }
>   }
> and for_each is implemented such that it internally expands to a for_each
> over each of the two containers for tets and hexes. Thus, instead of
> dispatching for each element over and over again, the dispatch is done once
> for the whole array.
>
> Just my 2 cents,
> Karli
>
>



--
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

[petsc-dev] PetscSection

[Matthew Knepley]

On Sat, Nov 10, 2012 at 3:35 PM, Jed Brown  wrote:
> On Sat, Nov 10, 2012 at 2:28 PM, Matthew Knepley  wrote:
>>
>> > point that is a member. What about if a given level happens to be empty?
>>
>> DMComplexGetLabelIdIS(). Its not problem if a level is empty.
>
>
> Ah, I misinterpreted what that routine did. Aren't these trivially short
> arrays (max length of 4 for meshes in 3D) that are expected to be identical
> on every process?

Not identical on every process.

> How does one get an IS containing _only_ the labeled points?

DMComplexGetStratumIS()

>>
>> > Can you answer how you can distinguish a quad from a tet in a
>> > non-interpolated mesh?
>>
>> You can't, but that is the point. You are not supposed to distinguish
>> them.
>
>
> Hmm, so these are the same at the category-theory level, but definitely not
> at the mesh topology level. It doesn't seem right to not distinguish them at
> all. How do you write DMComplexInterpolate for mixed-dim strata?

Interpolation is analysis. You need to know what dimension you are in, etc.

   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

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 3:49 PM, Matthew Knepley  wrote:

> I am not sure I understand. Labels just define point sets. If you want
> to segregate point
> sets by dimension, that is fine. Nothing in the interface is sensitive
> to this, nor should it be.
>

I can think of only one way in which numbering points by stratum versus
dimension makes a tangible difference: whether contiguous ranges can
contain unsorted mixed-dimensional points. I want to understand whether
this is actually important and whether there is a different use case that
is important. If not, then we could entirely discard the concept of
stratum, sort points by co/dimension, and query based on labels.

Where is stratum as a concept irreplaceable?


> >>
> >> You don't need coordinates for interpolation?
> >
> >
> > Are we talking about the same interpolation? If I have cell-to-vertex
> > connectivity, I can create the faces without coordinates, yes.
>
> Ah, we are not. Mesh interpolation.
>
> The way I have done interpolation, you have to know the mesh dimension,
> DMComplexGetDimension(), and then it assumes that height 0 stuff is cells
> of that dimension.
>

Then we're basically identifying stratum with dimension, suggesting that we
should be able to remove stratum from the API in favor of co/dimension.
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/b2c42c63/attachment.html>

[petsc-dev] TSMonitorSPEigCtxCreate calls KSPSetFromOptions

[Barry Smith]

On Nov 10, 2012, at 12:39 PM, Matthew Knepley  wrote:

> On Sat, Nov 10, 2012 at 12:33 PM, Jed Brown  wrote:
>> which eventually calls PetscOptionsEnd, masking all subsequent options from
>> the TS output. Choices:
>> 
>> 1. circumvent the limitation, moving the call outside the TS options
>> 
>> 2. make PetscOptionsBegin() behave like PetscOptionsHead() when called
>> inside a currently active options section
> 
> I prefer having them nest properly.

  Allowing proper nesting would eliminate this kind of problem once and and for 
all.  So maybe PetscOptionsBegin() could be totally bagged.


   Barry

> 
>   Matt
> 
>> others?
> 
> 
> 
> --
> 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

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 3:25 PM, Matthew Knepley  wrote:

> Everything is just a point, so I get back mixed dimension things all the
> time. I
> usually separate by stratum.
>

Okay, but what is the compelling use case for DMComplexGetStratumIS()
returning mixed-dimension that cannot be reasonably accomplished by getting
dimensional ISs and union/concatenating them?


> You don't need coordinates for interpolation?
>

Are we talking about the same interpolation? If I have cell-to-vertex
connectivity, I can create the faces without coordinates, yes.
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/ba9c62d0/attachment.html>

[petsc-dev] PetscSection

[Matthew Knepley]

On Sat, Nov 10, 2012 at 3:19 PM, Jed Brown  wrote:
> On Sat, Nov 10, 2012 at 2:05 PM, Matthew Knepley  wrote:
>>
>> Well, you have to designate the set with something. The easiest thing to
>> use is
>> an integer, so I use that. This is the only value associated with the set,
>> so I
>> would say Labels really only make sets.
>
>
> I was looking at DMComplexGetLabelValue() and didn't realize it's just one
> value per stratum. So you get the value associated with a point by searching
> to find which stratum it resides in. In the implementation, if point is not

Yes.

> in the set specified by the label, it does not set *value. Is that

It gives back -1 (complex.c:2103)

> intentional? And there is no way to query stratum values without getting a
> point that is a member. What about if a given level happens to be empty?

DMComplexGetLabelIdIS(). Its not problem if a level is empty.


> Also, linear search, oh my. ;-)

Was waiting to fix. Fixed.

>>
>>
>> Exactly, they are sorted by strata, not dimension.
>
>
> Okay, then dimensional query would really need to return an IS, provided it
> didn't switch to sorting by dimension.

Yes. Dimensions are just labels, so they can reuse the Label interface.

> Can you answer how you can distinguish a quad from a tet in a
> non-interpolated mesh?

You can't, but that is the point. You are not supposed to distinguish them.

   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

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 3:12 PM, Matthew Knepley  wrote:

> Yes, I agree completely with batching of this kind, which is why my
> integration
> functions look like they do. However, the topology information is not
> really needed
> in these kernels. Instead you are packing discretization information, like
> FEM coefficients, or geometric information like Jacobians, which are
> all controlled
> by PetscSection, not by the DMComplex.
>

At the time you call an element residual kernel, you definitely need that
information. Now the question is how much should bubble up to higher levels.

Because I think topological dimension is an important attribute (I think
the library should be able to distinguish tets from quads), I was only
proposing using topological dimension more explicitly. Karl suggests one
step further in delineating the specific topology. My concern with that is
keeping a contiguous index space. If there is sorting by cell topology, you
no longer have to store size explicitly for each point in the coneSection.
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/ba5b2fd0/attachment.html>

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 3:02 PM, Matthew Knepley  wrote:

>
> > How does one get an IS containing _only_ the labeled points?
>
> DMComplexGetStratumIS()
>

Great. Under what circumstances do you call this routine and get back an IS
with mixed dimension?


>
> >>
> >> > Can you answer how you can distinguish a quad from a tet in a
> >> > non-interpolated mesh?
> >>
> >> You can't, but that is the point. You are not supposed to distinguish
> >> them.
> >
> >
> > Hmm, so these are the same at the category-theory level, but definitely
> not
> > at the mesh topology level. It doesn't seem right to not distinguish
> them at
> > all. How do you write DMComplexInterpolate for mixed-dim strata?
>
> Interpolation is analysis. You need to know what dimension you are in, etc.
>

But it's still topological rather than geometric. You don't need
coordinates, for example.
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/7377d4d8/attachment.html>

[petsc-dev] PetscSection

[Karl Rupp]

On 11/10/2012 02:29 PM, Jed Brown wrote:
> On Sat, Nov 10, 2012 at 2:09 PM, Blaise A Bourdin  > wrote:
>
> Nope... and please, don't break that...
>
> Here is an example of a mesh with 4 element blocks. In an analysis
> code, for each block , I would specify the element type, type of
> physics, coefficients etc in a petscBag. That said, I may be doing
> it wrong...
>
> block 10 is made of tri
> block 22 is made of quad
> blocks 100 and 102 are made of bars
>
>
> I'm not suggesting any demotion of this ability. I'm trying to
> understand what is the advantage of having a stratum with unsorted mixed
> dimension. If the blocks were sorted by dimension, you would still
> access them using a label (which gives you an index set, perhaps
> contiguous). You would still be able to do all the normal operations
> with those points. If you want to have one loop over all elements of the
> various types, you would union or concatenate the index sets of each
> type. Am I missing something?

May I add something that is not specific to the current implementation 
in PETSc, yet still backs up the 'arrays of elements of the same type' 
approach:
  a) When thinking of accelerators, having, say, all tets 
stored/referenced contiguously in one array and all hexahedra stored 
contiguously in another (separate) array, makes it a lot easier to move 
all tets to a kernel working on tets, and to feed all hexes to a second 
kernel working on hexes.

  b) In a strongly typed setting (my experience here is in C++), one can 
do pretty nice things by keeping the two arrays separate rather than 
trying to unify them in some base class. For simplicity, suppose you 
have a mesh and two containers holding the tets and hexes. A FEM 
assembly could then look like this:

   for_each(mesh_cells.begin(), mesh_cells.end(), assembly_functor());

where assembly_functor statically dispatches into the cell type:
   struct assembly_functor{
 void operator()(tetrahedron const & t) { ... }
 void operator()(hexahedron const & t) { ... }
   }
and for_each is implemented such that it internally expands to a 
for_each over each of the two containers for tets and hexes. Thus, 
instead of dispatching for each element over and over again, the 
dispatch is done once for the whole array.

Just my 2 cents,
Karli

[petsc-dev] PetscSection

[Matthew Knepley]

On Sat, Nov 10, 2012 at 2:54 PM, Jed Brown  wrote:
> On Sat, Nov 10, 2012 at 1:44 PM, Matthew Knepley  wrote:
>>
>> Yes, I agree that they are different geometrically, and in PyLith we
>> make Labels to
>> distinguish, which I think is the right way.
>
>
> Cool. Of course you have to create the labels somehow. Do you currently make
> labels for all the "sidesets" in the mesh you load (e.g., from exodus)?

Yep.

> Labels in their current form seem to do two things (denote sets and
> associate values with points in those sets). Don't you frequently want one
> or the other?

Well, you have to designate the set with something. The easiest thing to use is
an integer, so I use that. This is the only value associated with the set, so I
would say Labels really only make sets.

>>
>> > Fine, I care much more about conceptual simplicity in the interface. I
>> > think
>> > adoption will be higher if 95% of users don't see the word "stratum".
>>
>> Note that this will have to be set explicitly somehow, since you cannot
>> "tell"
>> from the DAG that something has a given dimension/co-dim.
>
>
> Are your points not sorted by dimension? (Most of the code I've seen in the
> examples assumes homogeneous strata.)

Exactly, they are sorted by strata, not dimension.

   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

[petsc-dev] PetscSection

[Matthew Knepley]

On Sat, Nov 10, 2012 at 2:33 PM, Jed Brown  wrote:
> On Sat, Nov 10, 2012 at 1:26 PM, Matthew Knepley  wrote:
>>
>> I do not understand. Codim 0 and 1 are handled in an identical way.
>
>
> The user doesn't do the same thing with them. Your data structure treats
> them the same way (though I'm not sure how you distinguish between a tet and
> a quad if both are "height 0" and you don't store intermediate levels), but
> that seems like an implementation detail to me.

Yes, I agree that they are different geometrically, and in PyLith we
make Labels to
distinguish, which I think is the right way.

>>
>> I am not against some convenience interface for dim/co-dim, however it
>> would
>> just be implemented by making a label, since we have many equivalent
>> structural
>> things which would map to it.
>
>
> Fine, I care much more about conceptual simplicity in the interface. I think
> adoption will be higher if 95% of users don't see the word "stratum".

Note that this will have to be set explicitly somehow, since you cannot "tell"
from the DAG that something has a given dimension/co-dim.

  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

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 2:28 PM, Matthew Knepley  wrote:

> > point that is a member. What about if a given level happens to be empty?
>
> DMComplexGetLabelIdIS(). Its not problem if a level is empty.
>

Ah, I misinterpreted what that routine did. Aren't these trivially short
arrays (max length of 4 for meshes in 3D) that are expected to be identical
on every process?

How does one get an IS containing _only_ the labeled points?


> > Can you answer how you can distinguish a quad from a tet in a
> > non-interpolated mesh?
>
> You can't, but that is the point. You are not supposed to distinguish them.


Hmm, so these are the same at the category-theory level, but definitely not
at the mesh topology level. It doesn't seem right to not distinguish them
at all. How do you write DMComplexInterpolate for mixed-dim strata?
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/21e5e09f/attachment-0001.html>

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 2:09 PM, Blaise A Bourdin  wrote:

> Nope... and please, don't break that...
>
>  Here is an example of a mesh with 4 element blocks. In an analysis code,
> for each block , I would specify the element type, type of physics,
> coefficients etc in a petscBag. That said, I may be doing it wrong...
>
>  block 10 is made of tri
> block 22 is made of quad
> blocks 100 and 102 are made of bars
>

I'm not suggesting any demotion of this ability. I'm trying to understand
what is the advantage of having a stratum with unsorted mixed dimension. If
the blocks were sorted by dimension, you would still access them using a
label (which gives you an index set, perhaps contiguous). You would still
be able to do all the normal operations with those points. If you want to
have one loop over all elements of the various types, you would union or
concatenate the index sets of each type. Am I missing something?
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/4bc9dffd/attachment.html>

[petsc-dev] PetscSection

[Matthew Knepley]

On Sat, Nov 10, 2012 at 2:16 PM, Jed Brown  wrote:
> On Sat, Nov 10, 2012 at 12:45 PM, Matthew Knepley  
> wrote:
>>
>> It just depends on what you want your assumptions to be. I am all for
>> experimenting. The
>> foregoing strategy is good because closure continues to work as
>> advertised, and all our
>> integration code is still dim-independent, but the assumption that
>> height 0 things are all
>> cells breaks down. This is not so bad, since we usually want to group
>> cells by material
>> anyway, using a Label, so I am fine with this.
>
>
> Okay, my complaint is that "stratum" is an additional concept that has no
> established use in the current lexicon and the generality of which doesn't
> really help the user because they have to do different things for codim 0
> and codim 1 anyway. If we are not causing extreme hardship or irreparably
> crippling the flexibility of their code, I think the interface should give
> access to by codim instead of by stratum.

I do not understand. Codim 0 and 1 are handled in an identical way.

I am not against some convenience interface for dim/co-dim, however it would
just be implemented by making a label, since we have many equivalent structural
things which would map to it.

   Matt

> Now in the P1 case, the faces do not store any data, therefore
> closure(some_cell) returns the same thing going directly to vertices as if
> the intermediate dimensions were populated. I think that by definition, this
> sort of mesh does not support getting faces of a cell, therefore it's
> correct to not store the cell -> boundary-face relation. The user can ask to
> interpolate the mesh if they want.
>
> I'm not sure about the multi-domain case where the user wants a high-order
> discretization in one domain and a P1 with non-interpolated mesh in another.
> The benefit of non-interpolated doesn't seem so clear in this case.



--
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

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 2:05 PM, Matthew Knepley  wrote:

> Well, you have to designate the set with something. The easiest thing to
> use is
> an integer, so I use that. This is the only value associated with the set,
> so I
> would say Labels really only make sets.
>

I was looking at DMComplexGetLabelValue() and didn't realize it's just one
value per stratum. So you get the value associated with a point by
searching to find which stratum it resides in. In the implementation, if
point is not in the set specified by the label, it does not set *value. Is
that intentional? And there is no way to query stratum values without
getting a point that is a member. What about if a given level happens to be
empty?

Also, linear search, oh my. ;-)


>
> Exactly, they are sorted by strata, not dimension.
>

Okay, then dimensional query would really need to return an IS, provided it
didn't switch to sorting by dimension.


Can you answer how you can distinguish a quad from a tet in a
non-interpolated mesh?
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/4439964f/attachment.html>

[petsc-dev] Fwd: [mumps-dev] support for distributed right-hand vectors?

[Alexander Grayver]

Garth,

At the time I was tested PaStiX it failed for my problem:
https://lists.mcs.anl.gov/mailman/htdig/petsc-dev/2011-December/006887.html

Since then PaStiX has been updated with several critical bug fixes, so I 
should consider testing new version.

The memory scalability of the MUMPS is not nice, that is true.
Running MUMPS with default parameters on large amount of cores is often 
not optimal. I don't how much you spent tweaking parameters.
MUMPS is among the most robust distributed solvers nowadays and it is 
still being developed and hopefully will improve.

/*To petsc developers:* /are there plans to update PaStiX supplied with 
PETSc? The current version is 5.2 from 2012-06-08 and PETSc-3.3-p3 uses 
5.1.8 from 2011-02-23.

Here is changelog:
https://gforge.inria.fr/frs/shownotes.php?group_id=186&release_id=7096

On 09.11.2012 19:40, Garth N. Wells wrote:
> I've only just joined the petsc-dev list, but I'm hoping with this
> subject line my email will join the right thread . . . . (related to
> MUMPS).
>
> I've been experimenting over the past year with MUMPS and PaStiX for
> parallel LU, and found MUMPS pretty much useless because it uses so
> much memory. PaStiX was vastly superior performance-wise and it
> supports hybrid threads-MPI, which I think is essential for parallel
> LU solvers to make good use of typical multi-socket multi-core compute
> nodes. The interface, build and documentation are a bit clunky (I put
> the last point down to developer language issues), but the performance
> is good and the developers are responsive. I benchmarked PaStiX for P1
> and P2 3D linear elastic finite element problems against a leading
> commercial offering, and PaStiX was marginally faster for P1 and
> marginally slower for P2 (PaStiX performance does depend heavily on
> BLAS). I couldn't even compute the test problems with MUMPS because it
> would blow out the memory. For reference, I tested systems up to 27M
> dofs with PaStiX.
>
> Based on my experience and tests, I'd be happy to see PETSc drop MUMPS
> and focus/enhance/fix support for PaStiX.
>
> Garth


-- 
Regards,
Alexander

-- next part ------
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/e2898c6d/attachment.html>

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 1:44 PM, Matthew Knepley  wrote:

> Yes, I agree that they are different geometrically, and in PyLith we
> make Labels to
> distinguish, which I think is the right way.
>

Cool. Of course you have to create the labels somehow. Do you currently
make labels for all the "sidesets" in the mesh you load (e.g., from exodus)?

Labels in their current form seem to do two things (denote sets and
associate values with points in those sets). Don't you frequently want one
or the other?


> > Fine, I care much more about conceptual simplicity in the interface. I
> think
> > adoption will be higher if 95% of users don't see the word "stratum".
>
> Note that this will have to be set explicitly somehow, since you cannot
> "tell"
> from the DAG that something has a given dimension/co-dim.
>

Are your points not sorted by dimension? (Most of the code I've seen in the
examples assumes homogeneous strata.)
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/10d841ca/attachment.html>

[petsc-dev] PetscSection

[Matthew Knepley]

On Sat, Nov 10, 2012 at 12:04 PM, Jed Brown  wrote:
> On Sat, Nov 10, 2012 at 10:40 AM, Blaise A Bourdin  wrote:
>>
>> On Fri, Nov 9, 2012 at 10:20 PM, Blaise A Bourdin  wrote:
>>>
>>> A DM does double duty by describing the geometry of the mesh, and the
>>> data layout associated with the finite element space. I liked the model
>>> where the mesh geometry and the data layout on the mesh was split in two
>>> objects, but I understand the convenience of having everything in the DM,
>>> and DMClone works just fine. Since I may have to handle scalar, vector, 2nd
>>> and 4th order tensor on 2 different finite element spaces, in an assembly
>>> loop, I may end up dealing with as many as 8 DM. I stick all these DM's in
>>> the user context of each DM associated with a unknown (a Vec on which I may
>>> have to call SNESSolve or TSSolve), hoping that this is not creating some
>>> aliasing problem which as  a fortran programmer I can not possibly
>>> understand.
>>
>>
>> ;-)
>>
>>
>> Actually, I am really not sure if passing a dm and a pointer to a user
>> context containing this dm is legit or not...
>
>
> Ah, okay. To confirm, you have a DM that you are solving for, and in its
> user context, you have several other DMs, each with a Vec, describing the
> "problem data" like coefficients, forcing terms, and internal
> discontinuities? That is completely fine, and not "aliasing", but it does
> not play well with geometric multigrid because coarse grids reference the
> same application context. We have a system of hooks for managing such
> resolution-dependent data, though only with a C interface so far. (We needed
> this to get geometric multigrid and FAS to work with TS. Most non-toy
> applications need it too.)
>
> I'm not sure if there is a way to make this easier. We have been using
> PetscObjectCompose() to attach things to the DM on different levels. We
> could have a slightly friendlier "user" interface for that.
>
> So keeping those things in the app context is just fine, but if you want to
> use geometric multigrid, you'll have to take them out of the app context and
> put them in a different structure attached to the DM that is not
> transparently propagated under coarsening and refinement. If you think you
> might do this eventually, I recommend commenting/organizing your app context
> so that resolution-dependent stuff is easily identifiable.
>
>>
>> I don't mind that, but can't you have an index set describing the codim 0
>> elements (maybe all of them) and another index set for the codim 1 elements
>> on the features you care about? You can take their union (or concatenate)
>> for your assembly loop if you like. Is there something wrong with this
>> approach?
>>
>>
>> Thats a very good point. In the end it doesn't really matter. As far as I
>> remember, the main reason I ended with my current scheme is that DMMesh did
>> not play well with partially interpolated meshes. I don't know what the
>> current status of DMComplex is.
>
>
> Okay, I think it's important to eventually support partially interpolated
> meshes to avoid using a lot of memory when used with low-order
> discretizations. I see no reason why there can't also be a direct cache for
> closure. For a P1 basis, that amounts to a point range
>
> [cells, boundary faces, vertices]
>
> closure: [cells -> vertices, faces -> vertices]
>
> So cell -> face need not be stored anywhere. Presumably there is a reason
> why Matt didn't write it this way. Is it just uniformity of data structure?

It just depends on what you want your assumptions to be. I am all for
experimenting. The
foregoing strategy is good because closure continues to work as
advertised, and all our
integration code is still dim-independent, but the assumption that
height 0 things are all
cells breaks down. This is not so bad, since we usually want to group
cells by material
anyway, using a Label, so I am fine with this.

   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

[petsc-dev] TSMonitorSPEigCtxCreate calls KSPSetFromOptions

[Matthew Knepley]

On Sat, Nov 10, 2012 at 12:33 PM, Jed Brown  wrote:
> which eventually calls PetscOptionsEnd, masking all subsequent options from
> the TS output. Choices:
>
> 1. circumvent the limitation, moving the call outside the TS options
>
> 2. make PetscOptionsBegin() behave like PetscOptionsHead() when called
> inside a currently active options section

I prefer having them nest properly.

   Matt

> others?



--
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

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 1:26 PM, Matthew Knepley  wrote:

> I do not understand. Codim 0 and 1 are handled in an identical way.
>

The user doesn't do the same thing with them. Your data structure treats
them the same way (though I'm not sure how you distinguish between a tet
and a quad if both are "height 0" and you don't store intermediate levels),
but that seems like an implementation detail to me.


> I am not against some convenience interface for dim/co-dim, however it
> would
> just be implemented by making a label, since we have many equivalent
> structural
> things which would map to it.
>

Fine, I care much more about conceptual simplicity in the interface. I
think adoption will be higher if 95% of users don't see the word "stratum".
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/2a7eeeb0/attachment.html>

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 12:45 PM, Matthew Knepley  wrote:

> It just depends on what you want your assumptions to be. I am all for
> experimenting. The
> foregoing strategy is good because closure continues to work as
> advertised, and all our
> integration code is still dim-independent, but the assumption that
> height 0 things are all
> cells breaks down. This is not so bad, since we usually want to group
> cells by material
> anyway, using a Label, so I am fine with this.
>

Okay, my complaint is that "stratum" is an additional concept that has no
established use in the current lexicon and the generality of which doesn't
really help the user because they have to do different things for codim 0
and codim 1 anyway. If we are not causing extreme hardship or irreparably
crippling the flexibility of their code, I think the interface should give
access to by codim instead of by stratum.

Now in the P1 case, the faces do not store any data, therefore
closure(some_cell) returns the same thing going directly to vertices as if
the intermediate dimensions were populated. I think that by definition,
this sort of mesh does not support getting faces of a cell, therefore it's
correct to not store the cell -> boundary-face relation. The user can ask
to interpolate the mesh if they want.

I'm not sure about the multi-domain case where the user wants a high-order
discretization in one domain and a P1 with non-interpolated mesh in
another. The benefit of non-interpolated doesn't seem so clear in this case.
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/19d3e759/attachment.html>

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 11:57 AM, Blaise A Bourdin  wrote:

> Yes, but because of the structure of my problem, the application context
> also contains a link to the original DM.
> Think solving F(u,v) = 0 where and v don't have the same layout, or the
> same dimensionality, without field split.
> My iterate is
> u_{n+1} = argmin_w F(w,v_n)
>v_{n+1} = argmin_w F(u_{n+1},w)
> The DM associated with u is DMu and the one with v is DMv
> For the application context, I can either a create AppCtx with components
> (pointer to u, pointer to v, pointed to DMu, pointer to DMu), or
> AppCtxu with components (pointer to v, pointer to DMv) and AppCtxv with
> components (pointer to u, pointer to DMu)
> In this setting, is calling SNESSolve with Vec u and Application Context
> AppCtx legal, or do I have to use AppCtxu / AppCtxv?
>

Either is okay (and neither is okay for FAS).


>  Can you point me to an example? Are interfaces C only because nobody has
> ever needed fortran versions, or is there a technical reason.
>

SNES ex48 does this, also look at any use of DMCoarsenHookAdd(). There is
no fundamental issue, it just requires custom function pointer juggling
which is messy and really needs to be tested.


>
>
 I had not thought about it, but it is quite feasible. I store all
> coefficients that are constant per block of cells or Dof in PetscBags, and
> everything that has variation at the scale of the finite element space in
> Vecs.  How would the creation of the coarse DMs be handled, though? The
> geometry part is trivial to propagate using DMClone, but you may need to
> user feedback for the data layout part, unless you have a scheme that
> describes it (i.e. for this IS of cells, n dof at the vertices, p at the
> faces etc)
>

You get to implement these two callbacks to create and update data for the
coarse grid problem:

http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCoarsenHookAdd.html


> Do we mean the same by partially interpolated mesh? What I mean is a mesh
> the only faces that are explicitly stored are the one we are interested in
> (typically the boundary mesh). For P1 elements, we need only to know of
> [cells, some faces, vertices].
>

Yup, same thing.


>  I tried to manually build partially interpolated sieves with only some of
> the faces, and distribution would fail. That's how I ended up with a mix of
> cells of co-dimension 0 and 1. If one wants access to ALL faces / edges or
> none of them, there is no problem in the current implementation.
>

Right, I think of this as having a bunch of codim 0 and a few codim 1. As a
user, we call Closure, not Cone, so we go directly to the points from
either. But we would get the points to loop over by requesting codim 0 and
codim 1 separately. I'd expect that in principle, you likely have multiple
groups of codim 1 things that get treated differently. Of course you can
lump them all in one bigger index set and have an if statement, though it's
not as vectorization-friendly so Matt doesn't like that. I would not mind
having an interface for requesting points of multiple dimensions at once
(codim in {0,1}).
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/44352272/attachment.html>

[petsc-dev] TSMonitorSPEigCtxCreate calls KSPSetFromOptions

[Jed Brown]

which eventually calls PetscOptionsEnd, masking all subsequent options from
the TS output. Choices:

1. circumvent the limitation, moving the call outside the TS options

2. make PetscOptionsBegin() behave like PetscOptionsHead() when called
inside a currently active options section

others?
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/5588e375/attachment.html>

[petsc-dev] PetscSection

[Jed Brown]

On Sat, Nov 10, 2012 at 10:40 AM, Blaise A Bourdin  wrote:

> On Fri, Nov 9, 2012 at 10:20 PM, Blaise A Bourdin  wrote:
>
>> A DM does double duty by describing the geometry of the mesh, and the
>> data layout associated with the finite element space. I liked the model
>> where the mesh geometry and the data layout on the mesh was split in two
>> objects, but I understand the convenience of having everything in the DM,
>> and DMClone works just fine. Since I may have to handle scalar, vector, 2nd
>> and 4th order tensor on 2 different finite element spaces, in an assembly
>> loop, I may end up dealing with as many as 8 DM. I stick all these DM's in
>> the user context of each DM associated with a unknown (a Vec on which I may
>> have to call SNESSolve or TSSolve), hoping that this is not creating some
>> aliasing problem which as  a fortran programmer I can not possibly
>> understand.
>>
>
>  ;-)
>
>
>  Actually, I am really not sure if passing a dm and a pointer to a user
> context containing this dm is legit or not...
>

Ah, okay. To confirm, you have a DM that you are solving for, and in its
user context, you have several other DMs, each with a Vec, describing the
"problem data" like coefficients, forcing terms, and internal
discontinuities? That is completely fine, and not "aliasing", but it does
not play well with *geometric* multigrid because coarse grids reference the
same application context. We have a system of hooks for managing such
resolution-dependent data, though only with a C interface so far. (We
needed this to get geometric multigrid and FAS to work with TS. Most
non-toy applications need it too.)

I'm not sure if there is a way to make this easier. We have been using
PetscObjectCompose() to attach things to the DM on different levels. We
could have a slightly friendlier "user" interface for that.

So keeping those things in the app context is just fine, but if you want to
use geometric multigrid, you'll have to take them out of the app context
and put them in a different structure attached to the DM that is not
transparently propagated under coarsening and refinement. If you think you
might do this eventually, I recommend commenting/organizing your app
context so that resolution-dependent stuff is easily identifiable.


> I don't mind that, but can't you have an index set describing the codim 0
> elements (maybe all of them) and another index set for the codim 1 elements
> on the features you care about? You can take their union (or concatenate)
> for your assembly loop if you like. Is there something wrong with this
> approach?
>
>
>   Thats a very good point. In the end it doesn't really matter. As far as
> I remember, the main reason I ended with my current scheme is that DMMesh
> did not play well with partially interpolated meshes. I don't know what the
> current status of DMComplex is.
>

Okay, I think it's important to eventually support partially interpolated
meshes to avoid using a lot of memory when used with low-order
discretizations. I see no reason why there can't also be a direct cache for
closure. For a P1 basis, that amounts to a point range

[cells, boundary faces, vertices]

closure: [cells -> vertices, faces -> vertices]

So cell -> face need not be stored anywhere. Presumably there is a reason
why Matt didn't write it this way. Is it just uniformity of data structure?
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/7f052425/attachment.html>

[petsc-dev] Fwd: [mumps-dev] support for distributed right-hand vectors?

[Satish Balay]

On Sat, 10 Nov 2012, Alexander Grayver wrote:

> /*To petsc developers:* /are there plans to update PaStiX supplied with PETSc?
> The current version is 5.2 from 2012-06-08 and PETSc-3.3-p3 uses 5.1.8 from
> 2011-02-23.
> 
> Here is changelog:
> https://gforge.inria.fr/frs/shownotes.php?group_id=186&release_id=7096

Usually we update petsc-dev with newer releases of externalpackages.
The build appears to go through fine - so I pushed this to petsc-dev
http://petsc.cs.iit.edu/petsc/petsc-dev/rev/b7bdde2395c4


You can try the following with petsc-3.3:

--download-pastix=https://gforge.inria.fr/frs/download.php/30920/pastix_release_3725.tar.bz2

Satish

[petsc-dev] PetscSection

[Blaise A Bourdin]

Hi,

Let me make sure I understand the consensus, since I have a vested interest in 
the unstructured mesh business:

A DM does double duty by describing the geometry of the mesh, and the data 
layout associated with the finite element space. I liked the model where the 
mesh geometry and the data layout on the mesh was split in two objects, but I 
understand the convenience of having everything in the DM, and DMClone works 
just fine. Since I may have to handle scalar, vector, 2nd and 4th order tensor 
on 2 different finite element spaces, in an assembly loop, I may end up dealing 
with as many as 8 DM. I stick all these DM's in the user context of each DM 
associated with a unknown (a Vec on which I may have to call SNESSolve or 
TSSolve), hoping that this is not creating some aliasing problem which as  a 
fortran programmer I can not possibly understand. 

Everything is an IS. I suspect that this means that the Set/GetCone, 
Set/GetClosure operations would return an IS that could be used in 
VecSetValuesIS, but may not even be needed if an equivalent of  DMMesh' 
SectionRealRestrict / RestrictClosure /Update / UpdateCLosure is implemented

Setting Values is good, but getting values is needed too! 

Also, in addition to the Barry's simple pseudo-code, I need to be able to get 
an IS for the i^th dof of a field at a given point (think of applying a 
Dirichlet boundary conditions to only one component of a field, for instance). 
It's always been the messy part. How would that fit within the proposed model?

Finally, just a comment on stratum vs topological dimension. There is no reason 
why all elements in a mesh should have the same topological dimension. I 
actually find it much easier to have a single mesh where some sets of elements 
have codimension 0 and others codimension 1 and dispatching my assembly 
depending on the physics associated with each block, and the codimension of the 
element rather than having to deal with side sets / boundary mesh / whatever 
you want to call it.

Blaise

On Nov 9, 2012, at 9:02 PM, Barry Smith 
 wrote:

> 
>  We have a users manual?
> 
> On Nov 9, 2012, at 8:56 PM, Matthew Knepley  wrote:
> 
>> On Fri, Nov 9, 2012 at 9:46 PM, Barry Smith  wrote:
>>> 
>>> On Nov 9, 2012, at 8:42 PM, Jed Brown  wrote:
>>> 
 On Fri, Nov 9, 2012 at 8:39 PM, Matthew Knepley  
 wrote:
 You are both right :) I don't care where we put them, but I want BCs
 and fields. However, this is no problem.
 Right now I do BCs and fields as PetscSections held inside the primary
 section. If IS can do what PetscSection
 can do, I can jsut use another IS for each one.
 
 Convergence!
 
>>> 
>>>  Write it up!   :-)
>> 
>> In the manual section on Unstructured Grids :)
>> 
>>  Matt
>> 
 
 I actually like it where it is. Stick with DMComplexVecGetClosure()
 with no section argument, and then
 the DM holds an IS++.
 
 Cool, now the "normal" user doesn't even see the ++IS.
>>> 
>> 
>> 
>> 
>> --
>> 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
> 
> 

-- 
Department of Mathematics and Center for Computation & Technology
Louisiana State University, Baton Rouge, LA 70803, USA
Tel. +1 (225) 578 1612, Fax  +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

[petsc-dev] PetscSection

[Jed Brown]

On Fri, Nov 9, 2012 at 10:20 PM, Blaise A Bourdin  wrote:

> A DM does double duty by describing the geometry of the mesh, and the data
> layout associated with the finite element space. I liked the model where
> the mesh geometry and the data layout on the mesh was split in two objects,
> but I understand the convenience of having everything in the DM, and
> DMClone works just fine. Since I may have to handle scalar, vector, 2nd and
> 4th order tensor on 2 different finite element spaces, in an assembly loop,
> I may end up dealing with as many as 8 DM. I stick all these DM's in the
> user context of each DM associated with a unknown (a Vec on which I may
> have to call SNESSolve or TSSolve), hoping that this is not creating some
> aliasing problem which as  a fortran programmer I can not possibly
> understand.
>

;-)

Everything is an IS. I suspect that this means that the Set/GetCone,
>

I think these just return a range or (offset,size).


> Set/GetClosure operations
>

Closure needs to do orientation, so it's not direct access from the Vec.
Regardless, I don't think the plan would be to create an IS for that
granularity.


> would return an IS that could be used in VecSetValuesIS, but may not even
> be needed if an equivalent of  DMMesh' SectionRealRestrict /
> RestrictClosure /Update / UpdateCLosure is implemented
>

I think these are DMCompletVecGet/SetClosure.


>
> Setting Values is good, but getting values is needed too!
>
> Also, in addition to the Barry's simple pseudo-code, I need to be able to
> get an IS for the i^th dof of a field at a given point (think of applying a
> Dirichlet boundary conditions to only one component of a field, for
> instance). It's always been the messy part. How would that fit within the
> proposed model?
>

I think you would get an index set for those points on the feature you
wanted to do something special for.


> Finally, just a comment on stratum vs topological dimension. There is no
> reason why all elements in a mesh should have the same topological
> dimension. I actually find it much easier to have a single mesh where some
> sets of elements have codimension 0 and others codimension 1 and
> dispatching my assembly depending on the physics associated with each
> block, and the codimension of the element rather than having to deal with
> side sets / boundary mesh / whatever you want to call it.
>

I don't mind that, but can't you have an index set describing the codim 0
elements (maybe all of them) and another index set for the codim 1 elements
on the features you care about? You can take their union (or concatenate)
for your assembly loop if you like. Is there something wrong with this
approach?
-- next part --
An HTML attachment was scrubbed...
URL: 
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20121110/d04a3fd5/attachment.html>