On Wed, Jan 14, 2009 at 04:17:19PM -0500, Shawn Walker wrote:
>
>
> 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;
>>
>
> 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.
>
> - Shawn2D or 3D will be decided by mesh.geometry().dim() and triangle/tet will be decided by mesh.topology().dim(). -- Anders
signature.asc
Description: Digital signature
_______________________________________________ DOLFIN-dev mailing list [email protected] http://www.fenics.org/mailman/listinfo/dolfin-dev
