dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Wed, Jan 14, 2009 at 03:53:30PM -0500, Shawn Walker wrote:
>
>
> 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).
>>
>>
>
> 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

It would be enough to have a function in MeshGeometry that returned
the coordinates for a given cell, something like

  const double* cell_coordinates(uint cell) const;

or 

  const double * const * cell_coordinates(uint cell) const;

-- 
Anders

Attachment: signature.asc
Description: Digital signature