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
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> DOLFIN-dev mailing list
> [email protected]
> http://www.fenics.org/mailman/listinfo/dolfin-dev
_______________________________________________
DOLFIN-dev mailing list
[email protected]
http://www.fenics.org/mailman/listinfo/dolfin-dev
