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.
