Subject: [FFC-dev] functions

From: "Garth N. Wells" <[email-masked]>
Date: Tue, 10 Apr 2007 18:14:12 +0200

Connected to Kristian's work on quadrature in FFC, we should think about 
how to work with functions which do not come from a finite element 
space. Such functions (like stress for a plasticity model or viscosity 
for a non-Newtonian flow) are evaluated at quadrature points, rather 
than at nodes.

Garth

---

From: Anders Logg <[email-masked]>
Date: Wed, 11 Apr 2007 11:58:22 +0200

> for a non-Newtonian flow) are evaluated at quadrature points, rather 

The current design (including UFC) assumes that all functions can be
interpolated to a finite element basis, but I think it will work for
the quadrature code generation to imagine that you have a finite
element basis where the quadrature points are the same as the nodes.
(But you never need to know the basis functions.)

The array of coefficients (double** w) that comes in to the function
tabulate_tensor() should contain the coefficients, but for the
quadrature these will be the same as the values at the quadrature
points, and this should work out fine since the values that go in to w
are decided by the evaluate_dof() function that also gets
generated.

So the quadrature code generator just needs to make sure that
evaluate_dof() picks the values at the quadrature points.

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Wed, 11 Apr 2007 15:52:11 +0200

> evaluate_dof() picks the values at the quadrature points.

Perhaps we could define a superclass for finite_element which only
contains the functions needed for this, to make this more obvious?

class coefficient_function
{
  ... global_dimension(...) // == elements * quadpoints
  ... local_dimension(...) // == quadpoints
  ... evaluate_dof(f, ...) // == evaluate f in quadpoint
}

class finite_element: public coefficient_function
...

It would require some more changes in the interface though, and it
wouldn't add any functionality.

martin

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 11 Apr 2007 18:15:30 +0200

> evaluate_dof() picks the values at the quadrature points.

The problem with this approach is that it's not possible to create a 
consistent linearisation when dofs and integration points do not 
coincide. Kristian and I went through this in detail a while ago.

Garth

> /Anders
> _______________________________________________
> DOLFIN-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/dolfin-dev
>

---

From: Anders Logg <[email-masked]>
Date: Wed, 11 Apr 2007 19:25:45 +0200

> coincide. Kristian and I went through this in detail a while ago.

Yes, I remember, but my suggestion was that the dofs and the points
should coincide. Then the linearization should be ok?

/Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 12 Apr 2007 13:32:43 +0200

> should coincide. Then the linearization should be ok?

I don't see how this can work. The locations of the dofs and the Gauss 
quadrature points do not in general coincide.

Garth

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Thu, 12 Apr 2007 17:20:04 +0200

> quadrature points do not in general coincide.

The idea is to construct an "artificial finite element" where the dofs
_are_ values in the Gauss quadrature points. By generating a
corresponding implementation of ufc::finite_element where
evaluate_dof(i, f, ...) evaluates the function f in dof (quadpoint) i,
and local_dimension() returns the number of quadrature points, this
should work.

Some functions like evaluate_basis can have empty implementations,
perhaps raising an exception, if no sensible implementation is
possible.

martin

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 12 Apr 2007 17:58:51 +0200

> and local_dimension() returns the number of quadrature points, this

Isn't this getting overly-complicated? It would also be necessary to 
evaluate the functions and derivatives of functions on the "real" 
element at the dofs of the artificial element (to compute a stress 
update in a plasticity problem, for example).

Garth

-- 

Dr. Garth N. Wells
Faculty of Civil Engineering and Geosciences
Delft University of Technology
Stevinweg 1
2628 CN Delft
Netherlands

tel.   +31 15 278 7922
fax.   +31 15 278 5767
e-mail [email-masked]
url    http://www.mechanics.citg.tudelft.nl/~garth

---

From: Anders Logg <[email-masked]>
Date: Thu, 12 Apr 2007 18:35:24 +0200

> update in a plasticity problem, for example).

I thought you wanted to have the quadrature points coincide with the
nodes?

My understanding of this problem is the following: You have a linear
form L(v) that you assemble into a vector b, and you want to compute
the derivative db_i/db_j by assembling a matrix from some suitable
bilinear form.

The problem is then that the nonlinearity in L is computed by
operating directly on the dofs of the function you want to solve
for. For example, you have

    L(v) = \int v sin(u) dx

and you compute sin(u) in the integral by taking sin() of all the dofs
for u (which makes sense for a Lagrange basis).

The result is that when you differentiate L to get a bilinear form a and
assemble it into a matrix, you will get something that is
approximately the the derivative of b, but not exactly, since you have
approximated the function sin(u).

But (correct me if I'm wrong here) if you take sin() of u at all the
quadrature points you use to compute the integral, you will
approximate both L and a, but in such a way that the matrix A you get
from a will be the exact derivative of b.

So would it be enough to get the values of u at the quadrature points
and then apply the nonlinearity to those values?

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Thu, 12 Apr 2007 20:18:20 +0200

> update in a plasticity problem, for example).

Those coefficient functions should then be passed to tabulate_tensor
in w[][] by their dofs in the basis of the "real element", and
evaluated in the quadrature points within tabulate_tensor.

If you have a discrete function with an attached finite element basis,
you should be able to pass that as a ufc::function to
my_quad_finite_element::evaluate_dof if a precomputation of values in
quadrature points is what you want.

Speaking of which, evaluating derivatives of a finite element basis is
missing from the UFC interface, and can only be done within
tabulate_tensor (since the form compiler knows the basis) or weakly by
defining a separate form and inverting a mass matrix.

The user will have to choose in the FFC input file wether to get a
coefficient in a finite element basis or a "quadrature basis". Instead
of (I don't remember FFC syntax exactly):
  w = Lagrange(...)
  f = Function(w)
it could be:
  w = QuadRule(...) # or some better name...
  f = Function(w)

martin