dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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.

-- 
Anders

Attachment: signature.asc
Description: Digital signature