dolfin team mailing list archive

[Shawn Walker <walker@xxxxxxxxxxxxxxx>]

---------- Forwarded message ----------
Date: Wed, 20 Aug 2008 17:16:12 -0400 (EDT)
From: Shawn Walker <walker@xxxxxxxxxxxxxxxxxxx>
To: Anders Logg <logg@xxxxxxxxx>
Subject: Re: [DOLFIN-dev] XML format for Higher Order meshes

you mean ufc::mesh class or just mesh class as in Mesh.h?  I was going to 
modify the MeshGeometry class to have the extra information.  But that info 
should also be read into the ufc::cell class so that it would be "ready to use" 
by the tabulate_tensor routine.  This could be done by introducing another 
variable.

I think, no matter what, some code will need to be modified.  It will certainly 
be a little messy, at least not as elegant as the current setup where an affine 
map is assumed everywhere.  But it shouldn't be that bad either.  I think for 
the time being, we should restrict ourselves to maps built from continuous 
lagrange polynomial functions.  I know this won't cover the iso-geometric case 
(ala hughes).  But that will probably take a different implementation anyway. 
And being too general can also be bad.

I'm ok with having a .xml format that has some redundant info.  And it should 
be possible to writing a routine (in ufc maybe) that can interpolate a FEM 
function onto the 6 control points of a P2 triangle. But that can be done 
later.  Right now, I just want to get an .xml format that is backwards 
compatible and will allow me to at least read in the extra "control point" 
(i.e. the 6 vertices) data into ufc so that Tabulate_Tensor can access it. 
Once that is done, then I can see about modifying FFC to generate code for 
assembling curved elements.  At first, it will probably be just quadrature.

I realize I hardly understand this FENICS code at all, being that I am new to 
this project.  I appreciate your patience and time.

- Shawn

On Wed, 20 Aug 2008, Anders Logg wrote:

>  The problem is how to store the extra coordinates in the Mesh, not how
>  to communicate it to UFC. I guess it would be fine to send in 6
>  coordinates to the UFC functions, but we can only store 3 vertices per
>  cell in the Mesh class.
>
>  --
>  Anders
>
>
>  On Wed, Aug 20, 2008 at 02:40:56PM -0400, Shawn Walker wrote:
> >  ok.  so is what I suggest satisfactory?  Just have a separate variable to
> >  store the extra vertices?  And that extra variable would only be used for
> >  the local polynomial map.  Is this fine?
> > 
> >  - Shawn
> > 
> >  On Wed, 20 Aug 2008, Anders Logg wrote:
> > 
> > >  There are many more places in DOLFIN where we may iterate over the
> > >  vertices in a triangle, not just when we send data to the generated
> > >  UFC code. The topology must remain unchanged (it should still be
> > >  triangles) but the geometry may be modified for higher order mappings.
> > > 
> > > 
> > > 
> > >  On Wed, Aug 20, 2008 at 02:32:34PM -0400, Shawn Walker wrote:
> > > >  I understand what you are saying, but the ordering of the vertices is
> > > >  such that the first three vertices (0,1,2) are exactly the same as 
> > > >  they
> > > >  were before.  This would be done no matter what order polynomial map 
> > > >  you
> > > >  used. So, it seems to me it should be fine.
> > > > 
> > > >  However, to be sure, I can modify the ufc stuff such that the 
> > > >  coordinates
> > > >  variable only reads in the first three vertices.  So this would be
> > > >  EXACTLY as it was before.  I would then create a new ufc::cell 
> > > >  variable
> > > >  called `map_coordinates' (if you have a better name, please suggest) 
> > > >  and
> > > >  this would read in the 6 vertices.  And in Tabulate_Tensor, if higher
> > > >  order is desired, then the `map_coordinates' will be used to compute 
> > > >  the
> > > >  FEM matrix.
> > > > 
> > > >  This seems safe, though a little redundant.  Of course, one must be
> > > >  careful when creating the higher order mesh .xml file and ensure that 
> > > >  the
> > > >  first three vertices correspond to the usual triangle vertices.
> > > > 
> > > >  I also want to make sure that the extra `map_coordinates' is available 
> > > >  to
> > > >  be modified with a loop in dolfin.  This would be necessary for an ALE
> > > >  method when the mesh is deforming.
> > > > 
> > > >  - Shawn
> > > > 
> > > >  On Wed, 20 Aug 2008, Anders Logg wrote:
> > > > 
> > > > >  On Wed, Aug 20, 2008 at 07:15:06PM +0100, Garth N. Wells wrote:
> > > > > > 
> > > > > > 
> > > > > >  Anders Logg wrote:
> > > > > > >  I'm not sure this will work. If you attach 6 vertices to a 
> > > > > > >  triangle by
> > > > > > > 
> > > > > > >    <triangle index="0" affine="false" v0="0" v1="1" v2="2" v3="4" 
> > > > > > >    v4="5" v5="6"/>
> > > > > > > 
> > > > > > >  then all sorts of things will break (I imagine). A triangle 
> > > > > > >  always has
> > > > > > >  three vertices.
> > > > > > > 
> > > > > > 
> > > > > >  I noticed this too. It would make some things troublesome.
> > > > > > 
> > > > > > >  The geometry of the triangles is separate.
> > > > > > > 
> > > > > > >  Maybe we could just add extra data which could be "control 
> > > > > > >  points" for
> > > > > > >  the cell facets? For P2 it would be the edge/face midpoints.
> > > > > > > 
> > > > > > 
> > > > > >  I like the idea of defining facet data which would contain the 
> > > > > >  necessary
> > > > > >  info.
> > > > > 
> > > > >  A problem with defining facet data is that the facet numbering is 
> > > > >  not
> > > > >  known a priori. It depends on the algorithm used by DOLFIN to 
> > > > >  compute
> > > > >  the facets from the cells. So we can't store for example a mesh
> > > > >  function over the facets since the facet numbering may change.
> > > > > 
> > > > >  When we read input from VMTK, we need to read facet data (boundary
> > > > >  markers) and these are stored relative to the cell to which the 
> > > > >  facet
> > > > >  belongs and the local number of the facet relative to the cell 
> > > > >  (which
> > > > >  is unique).
> > > > > 
> > > > >  There is an example in data/meshes/aneurysm.xml.gz.
> > > > > 
> > > > > 
> > >