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