dolfin team mailing list archive

[Peter Brune <prbrune@xxxxxxxxx>]

Sorry to interrupt this discussion, but after getting Kristian's update to
FFC to fix the optimization problem, dealing with the unavoidable breaking
that occurs in FEniCS when you update one package but not others ( :) ), and
running the tests, P1 test/trial and P2 geometry looks like:

isoparametric assembly      |        1.4663      1.4663     1
regular assembly            |       0.38364     0.38364     1

Which is a 2x improvement.  Considering the fact that this requires
higher-order quadrature (if we are to not be committing variational crimes
above and beyond parametric coordinates), this is pretty much the marked
improvement I had hoped that optimization would give.  P2 vs P2 is improved
as well.

- Peter


On Mon, Mar 28, 2011 at 11:04 AM, Anders Logg <logg@xxxxxxxxx> wrote:

> On Mon, Mar 28, 2011 at 04:11:59PM +0100, Garth N. Wells wrote:
> >
> >
> > On 28/03/11 13:18, Anders Logg wrote:
> > > On Mon, Mar 28, 2011 at 12:50:08PM +0100, Garth N. Wells wrote:
> > >>
> > >>
> > >> On 28/03/11 12:11, Anders Logg wrote:
> > >>> On Mon, Mar 28, 2011 at 10:37:10AM +0100, Garth N. Wells wrote:
> > >>>>
> > >>>>
> > >>>> On 25/03/11 12:58, Anders Logg wrote:
> > >>>>> On Fri, Mar 25, 2011 at 10:28:20AM +0000, Garth N. Wells wrote:
> > >>>>>>
> > >>>>>>
> > >>>>>> On 24/03/11 19:10, walker@xxxxxxxxxxxx wrote:
> > >>>>>>> On Thu, March 24, 2011 1:32 pm, Anders Logg wrote:
> > >>>>>>>> On Thu, Mar 24, 2011 at 05:19:03PM +0000, Garth N. Wells wrote:
> > >>>>>>>>>
> > >>>>>>>>>
> > >>>>>>>>> On 24/03/11 14:28, walker@xxxxxxxxxxxx wrote:
> > >>>>>>>>>> On Thu, March 24, 2011 8:53 am, Peter Brune wrote:
> > >>>>>>>>>>> Ah, this argument, again.  I continue to disagree.  While you
> could
> > >>>>>>>>>>> certainly home-roll a FEM code that does this under
> assumptions like
> > >>>>>>>>>>> constant quadrature order and small meshes,  this is going to
> fall
> > >>>>>>>>> down
> > >>>>>>>>>>> hard
> > >>>>>>>>>>> as soon as you have varying quadrature degree or type, which
> we
> > >>>>>>>>> always do.
> > >>>>>>>>>>> This would also keep around way more state than the current
> assembly
> > >>>>>>>>>>> paradigm allows, not to mention the per-quadrature-point
> storage
> > >>>>>>>>>>> requirements of a full-blown tensor.
> > >>>>>>>>>>
> > >>>>>>>>>> I don't see why constant quadrature order is so terrible.  Is
> it
> > >>>>>>>>> really
> > >>>>>>>>>> necessary to have varying quadrature degree?  What do you mean
> be more
> > >>>>>>>>>> state info?  This is just another locally defined function.  I
> am a
> > >>>>>>>>> little
> > >>>>>>>>>> unclear about what you are saying.
> > >>>>>>>>>>
> > >>>>>>>>>
> > >>>>>>>>> I don't see what the 'argument' is either. What is it exactly?
> > >>>>>>>>>
> > >>>>>>>>>>> The other problem with this approach is that you severely
> restrict
> > >>>>>>>>> the
> > >>>>>>>>>>> spaces your coordinates can live in by having some
> oddly-defined
> > >>>>>>>>> notion of
> > >>>>>>>>>>> "higher-order vertices.  If you have some coefficient that is
> the
> > >>>>>>>>>>> coordinates instead it's automatic.  I think you definitely
> should
> > >>>>>>>>>>> recreate
> > >>>>>>>>>>> the jacobian in tabulate_tensor; it's what we do now; it's
> what we
> > >>>>>>>>> should
> > >>>>>>>>>>> do
> > >>>>>>>>>>> if we don't want to store per-quadrature-point tensors on the
> DOLFIN
> > >>>>>>>>> side
> > >>>>>>>>>>> when the DOLFIN side doesn't even know about quadrature
> points.
> > >>>>>>>>>>
> > >>>>>>>>>> I am not proposing having higher order vertices.  I am saying
> treat
> > >>>>>>>>> the
> > >>>>>>>>>> geometric map like **it has its own finite element space**!
>  If you
> > >>>>>>>>> want
> > >>>>>>>>>> to have a special map for each bilinear form (so you can have
> varying
> > >>>>>>>>> quad
> > >>>>>>>>>> degree) then fine.  But keep in mind that once you have a
> non-linear
> > >>>>>>>>> map
> > >>>>>>>>>> for the geometry, you will (in general) not be able to compute
> the
> > >>>>>>>>> tensors
> > >>>>>>>>>> exactly.  So you will need to fix the quad degree somehow.
> > >>>>>>>>>>
> > >>>>>>>>>
> > >>>>>>>>> I agree with this. Forgetting forms and PDEs, we should be able
> able to
> > >>>>>>>>> have a pure geometric representation. Say, for example, that I
> have a
> > >>>>>>>>> mesh that uses some mapping and I want to know (possibly
> approximately)
> > >>>>>>>>> its volume. The mesh and cells need to know all the geometry.
> Also, say
> > >>>>>>>>> I want to know in which cell a point (in the real coordiantes)
> lies.
> > >>>>>>>>> This is not related to forms.
> > >>>>>>>>
> > >>>>>>>> We've had this argument many times before... The thing that I've
> > >>>>>>>> always found attractive with Peter's approach is that one can
> > >>>>>>>> represent the geometry (the embedding of the mesh) as a field on
> top of
> > >>>>>>>> a standard mesh (with straight edges). So it doesn't require
> adding
> > >>>>>>>> tons of new data to the mesh class.
> > >>>>>>>>
> > >>>>>>
> > >>>>>> How would this work for very complicated meshes, e.g. a car,
> engine,
> > >>>>>> etc? In these cases one doesn't have a straightforward functional
> > >>>>>> description of the geometry, but just points that lie of the
> boundary or
> > >>>>>> a complex collection of spline.
> > >>>>>
> > >>>>> We would then need a function space that can represent those kinds
> of
> > >>>>> fields (NURBS). If we add such functionality (to DOLFIN and FFC),
> we
> > >>>>> could also use the NURBS as basis functions ("isogeometric
> analysis")
> > >>>>> which seems to be popular these days.
> > >>>>>
> > >>>>
> > >>>> This is missing the point. Complex geometry is given, and not
> > >>>> constructed for convenience by the person performing the analysis.
> It
> > >>>> can come in many formats. It is necessary to support geometries that
> are
> > >>>> provided.
> > >>>
> > >>> No one is arguing against that. The discussion is how we should are
> lotsstore
> > >>> the geometry: either in a special MeshGeometry class or as a
> Function.
> > >>>
> > >>
> > >> This is an over simplified perspective. We are not discussing whether
> or
> > >> not to use Function,
> > >
> > > I thought we were: either store the complex geometry as a Function or
> > > as part of (an extended) MeshGeometry class.
> > >
> > > I don't have any strong objections to doing the latter, but I think it
> > > would be cleaner if we could reuse the Function class and keep
> > > MeshGeometry simple.
> > >
> > > I hardly know anything about NURBS (and the like) but my basic
> > > assumption here is that whatever the geometry is, it can be
> > > represented by a mapping from the standard UFC elements. Maybe that's
> > > not the case?
> > >
> >
> > It's not the case.
>
> ok.
>
> > >> which since we all agree of a functional
> > >> representation is an implementation detail. The question is what is
> > >> required to represent a broad range of geometries.  Splines (the types
> > >> of interest) are not defined on simplices, so it's not a simple
> > >> questions of straightforward extension of FIAT. To evaluate splines,
> > >> geometric data is required.
> > >>
> > >> (I'm wishing that NURBS were never mentioned ;). There is more out
> there
> > >> than just NURBS - we need to deal with whatever is available. Part of
> my
> > >> assertion is that FFC + FIAT cannot possibly support everything.)
> > >
> > > Isn't the alternative that UFC supports everything?
> > >
> >
> > That's not a practical alternative on a practical time scale.
>
> What is not a practical alternative? My point was that if it is not
> supported through FFC, it will need to be supported in the UFC
> interface.
>
> > >>> My point is that if we store it as a Function (which among other
> > >>> things requires adding the corresponding basis functions in FIAT)
> then
> > >>> we could not only represent the geometry using NURBS, but we could
> > >>> also use them to compute. Another advantage would be that we wouldn't
> > >>> need to extend the MeshGeometry class or create a new MeshGeometry
> > >>> class to handle the new geometry.
> > >>>
> > >>
> > >> How could we construct NURBS on a domain without some geometric data
> for
> > >> constructing the NURBS? It's like an FE simulation on a domain, but
> > >> without a mesh to define the domain.
> > >
> > > Just as we do today: by passing data through the UFC interface. My
> > > belief is this can be done by representing the complex geometry as a
> > > suitable Function so that the data can be passed as a coefficient to
> > > tabulate_tensor. Is that not possible?
> > >
> >
> > How does passing data through UFC magically create a geometry and define
> > a domain? Data is required to construct the function that represents the
> > geometry (and which possibly provide a function space for a solution
> > basis). It's not a question of just providing a function of (x, y, z).
> > That is only for trivially simple cases. B-splines, for example, involve
> > knot vectors and control points.
>
> There would have to be an agreement between DOLFIN and FFC that one of
> the coefficients passed to tabulate tensor contains the nodal expansion
> coefficients of a field that describes the geometry, and FFC (FIAT)
> would need to support that basis locally.
>
> > If we're going to delve into assembly via tabulate_tensor, it's highly
> > unlikely that the present interface is suitable for non-FE basis
> > functions since there will no 'cells' in the conventional sense.
>
> Yes, perhaps.
>
> --
> Anders
>