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
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> 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
