On Wed, 14 Jan 2009, Anders Logg wrote:
> 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
And, obviously, this set of coordinates would be ordered (0,1,2,3,4,5) for
P2 map on a triangle. Will there need to be a flag or something to say if
the map takes the triangle to 2-D or 3-D? This higher order mapping could
certainly be used to map 2-D triangles into 3-D space.
- Shawn
_______________________________________________
DOLFIN-dev mailing list
[email protected]
http://www.fenics.org/mailman/listinfo/dolfin-dev
