dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

The main problem here seems to be how to even represent a higher-order
mesh. Are there any existing conventions for this?

Doing it locally is fine. For each cell, we have a convention for how
we number vertices, edges and faces (it's in the UFC manual) but we
don't have a convention for how to number mesh entities globally.

We could use the dof map (for Lagrange elements) that FFC computes
based on the numbering of mesh entities, but that computation is based
on a global numbering of edges and faces, and that numbering depends
on the algorithm for computing these entities (in TopologyComputation).

So, the only solution I see is to store the coordinates cell-wise,
which would mean duplicating the data, both in the XML file and for
storing the coordinates in MeshGeometry.

-- 
Anders


On Fri, Aug 22, 2008 at 09:22:50AM -0400, Shawn Walker wrote:
> I completely agree with not messing with the topology.  This is fine.  
> The only thing I have modified so far is the read in of the .xml file, 
> but only minimally.  I have no problem modifying MeshGeometry; in fact I  
> already put in another variable to hold the extra coordinates.  If this 
> is all we have to modify, then fine.  I was under the impression that the 
> Tabulate_Tensor routine would need to access this information through  
> ufc::cell, and so would require modifying that also.  But maybe that is  
> not so.
>
> But the extra control point data must be stored so that dolfin knows that 
> certain midpoint vertices are used by two different triangles.  This is  
> necessary if the mesh ever gets updated (i.e. moved), which I would need.
> I think the easiest thing to do here, and the safest as far as backwards  
> compatibility, is to store an additional list of vertex coordinates and  
> triangle connectivity data that has the P2 (or higher) ordering already  
> firgured out.  I would view this as a mesh generation problem, which is  
> completely separate from what I am trying to do.
>
> BTW: having this extra MeshGeometry variable would be useful for having  
> 1-D and 2-D elements in 3-D.  This was also on the FENICS to do list, and 
> would be useful to me.  However, this would come later.
>
> We could do a hack and only have P2 implemented, but I would prefer not  
> to.  It should not be that hard to have general order implemented, at  
> least for reading in the coordinates!  obviously, FFC would eventually be 
> modified to do P2 first, then P3, and so on.
>
> Here is a slightly different example for mesh.xml when the triangulation  
> is just affine:
>
> <dolfin xmlns:dolfin="http://www.fenics.org/dolfin/";>
>   <mesh celltype="triangle" dim="2">
>     <vertices size="4" topology="true">
>       <vertex index="0" x="0.0" y="0.0"/>
>       <vertex index="1" x="1.0" y="0.0"/>
>       <vertex index="2" x="0.0" y="1.0"/>
>       <vertex index="3" x="-1.0" y="0.0"/>
>     </vertices>
>     <cells size="2" topology="true">
>       <triangle index="0" v0="0" v1="1" v2="2"/>
>       <triangle index="1" v0="3" v1="0" v2="2"/>
>     </cells>
>   </mesh>
> </dolfin>
>
> Here is an example when it is P2:
>
> <dolfin xmlns:dolfin="http://www.fenics.org/dolfin/";>
>   <mesh celltype="triangle" dim="2">
>     <vertices size="4" topology="true">
>       <vertex index="0" x="0.0" y="0.0"/>
>       <vertex index="1" x="1.0" y="0.0"/>
>       <vertex index="2" x="0.0" y="1.0"/>
>       <vertex index="3" x="-1.0" y="0.0"/>
>     </vertices>
>     <cells size="2" topology="true">
>       <triangle index="0" v0="0" v1="1" v2="2"/>
>       <triangle index="1" v0="3" v1="0" v2="2"/>
>     </cells>
>
>     <vertices size="9" geometry="P2">
>       <vertex index="0" x="0.0"  y="0.0"/>
>       <vertex index="1" x="1.0"  y="0.0"/>
>       <vertex index="2" x="0.0"  y="1.0"/>
>       <vertex index="3" x="-1.0" y="0.0"/>
>
>       <vertex index="4" x="0.5"  y="0.05"/>
>       <vertex index="5" x="0.6"  y="0.6"/>
>       <vertex index="6" x="0.0"  y="0.5"/>
>       <vertex index="7" x="-0.5" y="0.5"/>
>       <vertex index="8" x="-0.5" y="0.0"/>
>     </vertices>
>
>     <cells size="2" geometry="P2">
>       <triangle index="0" affine="false" v0="0" v1="1" v2="2" v3="4"  
> v4="5" v5="6"/>
>       <triangle index="1" affine="true"  v0="3" v1="0" v2="2" v3="8"  
> v4="6" v5="7"/>
>     </cells>
>   </mesh>
> </dolfin>
>
> - Shawn
>
> On Fri, 22 Aug 2008, Anders Logg wrote:
>
>> On Fri, Aug 22, 2008 at 10:18:11AM +0200, Kent-Andre Mardal wrote:
>>> On to., 2008-08-21 at 22:51 -0400, Shawn Walker wrote:
>>>> Hello, all.  I've been thinking about this higher order mesh stuff, and I
>>>> need to be more clear about what I want to do.   I'm not sure it would be
>>>> good to make an ultra-general higher order mesh format.  It seems to be
>>>> more trouble than it is worth (right now).
>>>>
>>>> What I would like to have is a format for higher order meshes based on
>>>> polynomial mappings (ONLY!).  Currently, everything in FENICS assumes that
>>>> a P1 map is used for each element in the mesh.  Of course, this only
>>>> requires knowledge of the vertex coordinates of the element.  If a higher
>>>> order polynomial map is used, all you need are the positions of the extra
>>>> control points.  And this would extend the P1 map quite naturally.  For
>>>> the .xml format, I think the data for the control points should contain
>>>> the usual vertex coordinates also (not just the extra control points).
>>>> This would be redundant, but a little simpler to read in from the XML
>>>> file.  Note that storing the coordinates for each element separately would
>>>> not be good, because even the edge mid-point positions would be used by 2
>>>> (or more) elements, so you need to keep this connectivity information.
>>>>
>>>> One advantage of using a polynomial map (besides for free boundary
>>>> problems) is that the element tensor can be computed exactly in some cases
>>>> (but not for a stiffness matrix).  This isn't bad.
>>>>
>>>> Anyway, I know this is probably boring some of you.  And I know this is
>>>> not a high priority, but I am willing to put in time implementing this
>>>> (mostly because it would be useful to me).  I will repost another format
>>>> for the mesh.xml file and see what people think.
>>>
>>> This is definitly not boring to me, I would very much like a
>>> higher-order mesh in Dolfin! And I would consider it a great improvement
>>> of Dolfin.
>>>
>>> But as you have noticed there are some nice features concerning eg. the
>>> topology that probably should not be messed with. It might even be that
>>> the mesh should be left alone and that an additional structure can
>>> store the extra coordinates and their relation to mesh in therms of
>>> facets etc.
>>>
>>> Anyway, it seems that every concern is getting on the table and that it
>>> is a matter of sorting it out before implementing, hopefully.
>>>
>>> Finally, there is another option as well. You may write your own
>>> assembly routine that employs another mesh than the Dolfin mesh. UFC was
>>> originally designed such that it should be "easy" to use another mesh.
>>> But we have not really tested it.
>>>
>>> Kent
>>
>> As I mentioned before, what we need to do is to extend the
>> MeshGeometry class so that it can store additional geometry data.
>>
>> A simple solution that would work for P2 mappings would be to just
>> store an extra list of coordinates for edge midpoints.

Attachment: signature.asc
Description: Digital signature