← Back to team overview

dolfin team mailing list archive

  1. dolfin team
  2. Mailing list archive
  3. Message #11715

Re: recap of higher order mesh data

  Thread Previous • Date Previous • • Thread Next 


On Wed, 14 Jan 2009, Anders Logg wrote:


On Wed, Jan 14, 2009 at 03:38:39PM -0500, Shawn Walker wrote:



On Wed, 14 Jan 2009, Anders Logg wrote:


On Wed, Jan 14, 2009 at 08:23:04PM +0000, Garth N. Wells wrote:



Anders Logg wrote:

On Wed, Jan 14, 2009 at 08:12:33PM +0000, Garth N. Wells wrote:


Anders Logg wrote:

On Wed, Jan 14, 2009 at 09:38:23AM -0500, Shawn Walker wrote:

On Wed, 14 Jan 2009, Garth N. Wells wrote:


Anders Logg wrote:

On Wed, Jan 14, 2009 at 08:34:08AM +0000, Garth N. Wells wrote:

Anders Logg wrote:

On Wed, Jan 14, 2009 at 08:10:13AM +0000, Garth N. Wells wrote:

Shawn Walker wrote:

On Tue, 13 Jan 2009, Garth N. Wells wrote:


Shawn Walker wrote:

I cleared out the old email some because the discussion had changed a
little.  See below for a recap of higher order mesh data stuff:

-------------


It will if we want to be able to store a higher-order function space
as a function space with a regular mesh and an additional function
that stores the layout of the coordinates.


Perhaps that is not the best way to do the higher order mesh

coordinates.

If we want the higher order mesh data to be a general Function

(requiring

a FunctionSpace), then I do not see how you can get away from needing

the

FiniteElement signature associated with it, and possibly other things.

Even if you have the vector of data and the DoFmap, that info must

still

be used to create a Function/FunctionSpace in the code.  And in order

for

that to work the DoFmap must be `compatible' with the particular
FiniteElement you will be using.  I probably have this wrong, sorry
for
my confusion.

Another way to do the higher order mesh data would be to keep a little
simpler.  Have a vector of data, a DoFmap, and an indicator about the
degree of polynomial used.  This would be less general but not
bad.  In
case of higher-order mesh data, you will ALWAYS use a continuous

lagrange

finite element.  At least I cannot think of a situation where you
would
use something else.  Would this not be desirable?

If we decide to remove input/output for Functions and FunctionSpaces
(as I've understood is desirable since we then we don't need to rely
on precompiled elements and dofmaps) then how should we read in a
higher-order mesh from file?


Anders wrote:
Here's one option:

  Mesh mesh("mesh");
  LagrangeFunctionSpace V(mesh);
  File file("mesh_coordinate_vector.xml");
  Vector x;
  file >> x;
  V.set_coordinates(x);

That might work, but it's a bit long. There should be room for
improvement.

The discussion on higher-order meshes got a bit confusing for me a
little while back. In summary, exactly what information intended to be
in the mesh file for a high-order mesh?

Garth

-------------------------------------------

Ok, I will try to recap the higher order mesh stuff.

Currently, in a triangulation, there is an implicit assumption on the
form of the map that takes you from the `unit' reference triangle (or
tetrahedron).  The assumption is that the local map is linear.  As
you well know, this makes for various simplifications which can be
used during matrix assembly.

But, for various reasons, it can be more useful (or possibly required
depending on the nature of the FEM method) to have a curved triangle
to better approximate domain boundaries or to better compute higher
order geometric motion!

In this case, one could use a vector quadratic polynomial map and
have a triangle with edges given by a quadratic parametrization.  The
implementation of this only requires a local Lagrange finite element
basis, whose DoFs are just the coordinates of the nodes (for a
quadratic polynomial on a 2-D triangle, this would be 6 nodes per
triangle).  Of course, you will have this for every triangle, and it
makes sense to take the finite element basis to be continuous
lagrange over the whole domain. This continuity is especially
important when deforming the mesh!

So, way back we thought it would be a good idea to have a separate
functionspace to store this `higher order' mesh data.  But that
seemed problematic.


Sounds complicated.


However, in principle, all you need is a DoFmap and a vector of data
containing the node coordinate positions.

This is what I thought. Will we add a field the Mesh xml file to store
this extra data?

Yes.  I don't see why that would be a problem.  And if you don't want to
use the higher order mesh data (that happens to be in a file), then that
should also be fine.


OK, so we won't have the issue that Anders outlined above with respect
to reading in meshes.


And you need a method for

updating the positions (for a deforming mesh) but that isn't a big
deal. Once this information is properly stored, and accessible to the
matrix assembler, THEN...

Then the next step would be to modify FFC to use this higher order
(locally defined) map to compute the local matrices, INSTEAD of the
linear map that is implicitly assumed now.

I realize this will take some time, but we at least need to get a
storage scheme for the higher order mesh data to even proceed!


Kristian is looking at the UFL transition for the FFC quadrature
representation at the moment which will be needed for non-affine maps.

Perhaps a smaller first step in the non-affine direction would be to
support quadrilateral elements.

Garth

Did you mean quadratic elements?  Quadrilaterals are just deformed squares.


I meant quadrilaterals (with just 4 nodes) as a first step in having FFC
generate code for non-affine maps. I expect that quads would require
less initial work on the DOLFIN side, perhaps just an extension of
ufc::cell.


Yes, I agree.  In reality, I cannot forsee the potential difficulties
this will cause.  So, trying to have the full implementation ironed out
before we even put it in may not be helpful.  So, maybe just assuming a
2nd order vector polynomial for the local map may suffice.  This is very
much in line with the current philosophy of implicitly assuming a linear
map.

So, where would the data be stored in the code?  In FunctionSpace by
some extra variable field that contains the vector of coordinate data
and the DoFmap?


Using a FunctionSpace sounds complicated to me. What about letting the
mesh carry this data?

Garth

How would it be represented? We already know how to represent such
fields (by Functions). We would need to reinvent and reimplement
Lagrange elements as part of the Mesh class.


If we're happy with using a Lagrange basis for the map (at least for
now), all the form compiler needs is the locations of all the nodes. I
don't see the need for the complication of a FunctionSpace.

Garth

That would work if we stored the coordinates as a list of
coordinates for each cell:

  cell 0:
     node 0: x y z
     node 1: x y z
     node 2: x y z
     node 3: x y z
     node 4: x y z
     node 5: x y z
  cell 1:
     node 0: x y z
     node 1: x y z
     node 2: x y z
     node 3: x y z
     node 4: x y z
     node 5: x y z
   etc.


During assembly, we need to retrieve the coordinates for all nodes of
the current cell.

The problem is if we want to store a single list of coordinates
(without duplication) and then be able to map ourselves from the local
nodes on each cell to the global coordinate list. For that we need
some kind of dofmap which would be exactly the dofmap that the form
compilers generate for Lagrange elements.


Yes I agree that the dof map needs to be aware of this. What I don't see
is why it's any more complicated than the mesh carrying some extra data
and an appropriate dof map being generated.

Garth

You are all making good points here.  I understand the desire to want to
make this as repeatable and extendable as possible.  But at the same time,
it would be nice to experiment at first and see how it works in a simple
case first.

Perhaps in the long run, it would be better to have a FunctionSpace.  You
may want to interpolate the higher order mesh coordinates in case of a
re-mesh.  Incidentally, whenever you interpolate any function on the
higher order mesh, you will need to account for the higher order mapping.
I don't know how serious that is.  Of course, for now, it may be ok to
cheat a little on this point and just leave the interpolation as it is.

Sorry, I don't know what to say here.  I would just really like this to be
in the code in some way.  :)

- Shawn

How about storing the coordinates in a list and a sort of explicit
dofmap. Something like this:

  coordinates
  -----------
  node 0: x y z
  node 1: x y z
  etc

  cell_nodes
  ----------
  cell 0: 1 3 5 6 7 8
  cell 1: 1 6 7 12 15 16
  etc

So we would eseentially have a dofmap but we store it explicitly so
there's no need to write a function that computes the dofmap.


This is how most codes work and what mesh generators provide. We still
need a dof map to deal with cases other than scalar Lagrange.

Garth


What do you mean? Why would we need a dofmap? The cell_nodes list
above is the dofmap for the coordinates.



If each point has more than one dof.

Garth


I still don't understand. There will be no dofs, only coordinates for
points in the mesh. For a quadratic mapping of a triangular mesh, we
just need to list 6 numbers for each cell, and those numbers indicate
which coordinates are associated with the 6 points (three vertices and
three edge midpoints). Each coordinate (in 2D) is just two double
values.



I think technically, each "higher order vertex" has 3 dofs in 3-D (1 dof
for each scalar coordinate position).  But isn't this already taken into
account by the above data structure?  I guess it depends on how the
coordinate positions are stored.


Yes, we store the coordinates for each node together so we only need
to map one thing.


The above data structure seems reasonable.  It generalizes for higher
order lagrange.

Now, where will this be stored.  This is not a Function.  The coordinates
could be stored as a Vector.  Not sure about the DoFmap.


I think the natural thing would be to store it in MeshGeometry. That
was the original thought (to be able to extend it from just storing
vertex coordinates to higher order mappings).

--
Anders
Yes, I remember that. And this would require a modification of ufc.h and/or UFCCELL.h. I guess that when the mesh cell information is originally read (at the beginning of code), each cell can store it's own DoFmap, which is ONLY used for the higher order mapping. Unless you want to store the DoFmap more globally; that might be more desirable. 

- Shawn

Follow ups

References

  Thread Previous • Date Previous • • Thread Next