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,
So why do we need a dofmap is there are 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.
>
This will all be in the input - a list of points + cell connectivity -
and the mesh can carry this information. I don't see why we need to deal
with the mid-side nodes any differently to the vertices.
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
