Subject: [UFC-dev] evaluation of derivatives of basis functions

From: Kristian Oelgaard <[email-masked]>
Date: Wed, 18 Apr 2007 12:40:17 +0200

Hi UFC developers and others,

I've just finished code generation for evaluate_basis() in FFC, and a logical
next step would be to generate code for evaluation of derivatives. For this we
need to agree on an interface and the return values.

My suggestion for ufc/finite_element.py is:

/// Evaluate derivatives of order n of basis function i at given point in cell
virtual void evaluate_basis_derivatives(unsigned int i,
                                        unsigned int n,
                                        double* values,
                                        const double* coordinates,
                                        const ufc::cell& c) const
{
%(evaluate_basis_derivatives)s
}

For the return values, I think the best thing to do, is to return the values of
the derivative with respect to all directions of order 'n' as this will in most
cases be what one wants. So 'n=1' in 2D returns values = [x y], 'n=2' returns
values = [xx xy yx yy] etc.

Otherwise the caller needs to supply an array (or something else) that specifies
the desired direction e.g. 'n=2', direction = [0 1] --> value = xy

In this case one has to call the function 4 times to get all values, and as
a lot of variables inside the function will be the same the overhead could
become quite large.

It could also be possible to only return unique combinations of derivatives e.g.

value = [xx xy yy] in 2D

What do You think?

Kristian

---

From: Anders Logg <[email-masked]>
Date: Wed, 18 Apr 2007 13:47:27 +0200

> It could also be possible to only return unique combinations of derivatives e.g.

I'm all for this.

The reason we didn't include it initially was (at least for me) that
it seemed too much work to implement it, but now that Kristian will
supply the implementation, I see no reason why we should avoid it.

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Wed, 18 Apr 2007 14:27:16 +0200

> supply the implementation, I see no reason why we should avoid it.

Very good.

Only one question: should a finite_element always support this
function for any n? I guess that's the best way. Or do we need a
function to tell the user up to which order it is supported?

martin

---

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

> function to tell the user up to which order it is supported?

Depends on how much code needs to be generated. The compile time might 
really suffer when supporting all n.

Should UFC somehow provide an error message for functions which are not 
implemented? Functions like evaluate_basis can be quite large and will 
only sometimes be required, so it's likely that that programs like FFC 
will generate them only when requested, or only for particular element 
types.

Garth

> martin
> _______________________________________________
> UFC-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufc-dev
>

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Wed, 18 Apr 2007 16:56:47 +0200

> will generate them only when requested, or only for particular element

You could always throw an informative exception.

martin

---

From: Anders Logg <[email-masked]>
Date: Wed, 18 Apr 2007 17:32:29 +0200

> You could always throw an informative exception.

I don't see why this should be a problem. For n > some maximum order,
then the derivative is zero so very little code is needed to compute
that.

Also, at least for the code that FFC generates, taking the derivative
will just be a linear transformation of the expansion coefficients for
the basis functions (in some other basis) so it should not be very
much code.

If it turns out to be a problem, one can always add flags to FFC, like

  -fno-evaluate-basis-derivatives

/Anders

---

From: Kent-Andre Mardal <[email-masked]>
Date: Thu, 19 Apr 2007 09:28:00 +0200

>   -fno-evaluate-basis-derivatives

I think either we should have an integer in the construction of the
elements that spesify the order of the derivative, default behaviour
is derivatives of order 1. 

However, derivatives of order higher than two, who would use them ? 
The elements in which this is useful is a bit special. 
I mean, if we talk about Lagrangian elements of order 2+, then 
the derivative of order 2 is ok inside the elemements, but it would
be the derivatives of dirac functions at the facets of the elements. 
Hence, this feature would not be useful for Lagrangian elements. 

Kent

---

From: Anders Logg <[email-masked]>
Date: Thu, 19 Apr 2007 09:54:45 +0200

> Hence, this feature would not be useful for Lagrangian elements. 

Higher order derivatives might be useful for post processing, for
example, computing the residual when solving Laplace's equation.

/Anders

---

From: Anders Logg <[email-masked]>
Date: Thu, 19 Apr 2007 11:43:18 +0200

On Thu, Apr 19, 2007 at 11:40:29AM +0200, Kent-Andre Mardal wrote:
> tor, 19,.04.2007 kl. 09.54 +0200, skrev Anders Logg:
> Higher than 2 ? 

Ooops. No, not for Poisson but you may have third or forth derivatives
in other equations.

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 28 May 2007 21:56:20 +0200

> virtual void evaluate_basis_derivatives(unsigned int i,
...

> values = [xx xy yx yy] etc.

For vector elements (value_size > 1), how should values be sorted?

(
With
    n = 1
we have
    num_derivatives = num_space_dimensions
)

Let
    k in range(value_size) // subcomponent number
    j in range(num_derivatives),
then the length of values[] is
    value_size * num_derivatives

should the numbering be:

    values[k*num_derivatives + j]
or
    values[j*value_size + k]
?

I vote for
    values[k*num_derivatives + j]
ie packing all derivatives of a single basis function subcomponent together.

martin

---

From: Kristian Oelgaard <[email-masked]>
Date: Tue, 29 May 2007 09:55:55 +0200

> ie packing all derivatives of a single basis function subcomponent together.

Agree, and I think this is how it's implemented.

Kristian