#700 Remove C0LocalBasis base class and traits

Metadata

Property	Value
Reported by	Carsten Gräser (graeser@math.fu-berlin.de)
Reported at	Jan 2, 2010 12:10
Type	Feature Request
Version	Git (pre2.4) [autotools]
Operating System	Unspecified / All
Last edited by	Carsten Gräser (graeser@math.fu-berlin.de)
Last edited at	Jan 19, 2010 09:27
Closed by	Carsten Gräser (graeser@math.fu-berlin.de)
Closed at	Jan 19, 2010 09:27
Closed in version	2.0
Resolution	Implemented
Comment

Description

Currently we have a C0- and a C1LocalBasis base interface and corresponding C0- and C1LocalBasisTraits. This complicates the interface for no reason, thus I suggest to remove this differentiation:

Either C0/C1 means piecewise C0/C1. Then there will hardly be relevant shapefunctions that are piecewise C0 and not C1. In the esoteric case of everywhere continuous and nowhere differentiable functions we could throw an exception.

If the regularity is not meant in a piecewise sense then this differentiation is much to coarse. E.g. the RefinedP1LocalBasis is piecewise but not globally C1. Since one needs the gradients the C0 interface is not sufficient. Furthermore RefinedP0 is not even continuous but one needs the function values.

The same is true for CkLocalBasisTraits. The only usefull information of diffOrder is the maximal order of implemented derivatives. And passing a higher order than implemented could also result in an exception without knowing the max order statically.

I always understood that piecewise-C0 was meant when C0 was written in dune-localfunctions, likewise for C1 and Ck. The different C*LocalBasisInterface classes define different interfaces: C0LocalBasisInterface supports only evaluateFunction() and evaluateCoeffs(). C1LocalBasisInterface adds evaluateJacobian() and friends which are simply not present in C0LocalBasisInterface. CkLocalBasisInterface add methods for arbitrary derivatives, which are present neither in C0LocalBasisInterface nor in C1LocalBasisInterface. So the difference between these classes is more a matter of API than of mathematics.

C1LocalBasisTraits adds the JacobianType typedef to C0LocalBasisTraits, so this is again a matter of API.

So the question is: are the any algorithms which need to do different things at compile time based on whether a LocalBasis supports Jacobian or general differentiation? Do we require Jacobian and/or general differentian methods for all local bases?

In my opinion we should require all these methods since they make a mathematical sense for almost all FEs.

If we want to represent the missing implementation in some way an exception would be nice.

Currently only MimeticLocalFiniteElement uses the C0LocalBasisTraits - and throws exceptions for all methods that are not implemented since.

Sven Marnach (sven@pantoffel-wg.de) commented on the issue:

C0 and C1 could also be understood _element_wise, not piecewise. Piecewise differentiable functions might not be differentiable on the whole element, and there exist some macro elements with these properties. I have no opinion on whether the disctinction between C0 and C1 is useful in the code.

With Refined* some macro elements are already present as (piecewise) C1.

If anyone knows a good reason to keep C0 it would be nice to put a comment. Otherwise we should (IMHO) quickly remove it before the release of dune-localfunctions.

I think the question is not so much between C0/C1 (implementing a pointwise a.e. jacobian should not be that difficult). The question is more how to get higher derivatives into the interface. At the moment something standard (e.g. implemented in our computer coursers) like adaptive FE for laplace with residual based error indicator, cannot be implemented based on the interface. This would require either a hessian method or an evaluate method for higher order derivative. A Ck interface was present but has now been removed. I did not implement the virtualization for Ck because I first wanted to wait for the results of the discussion in Berlin. That does not mean that I think we do not need it; adaptive FE is a standard application which I think has to work with dune-localfunctions. So how should a Ck interface look like?

The idea of throwing exception at run time because a derivative has not been implemented is not a way I would like to use. Compile time errors are always preferable to run time errors.

I think we need an interface of the form template class VirtualLocalBasisInterface; We do not actually need C0/C1.

(Martin and Andreas)

If we keep the present interface for higher order derivatives (it's still used in the monomials) and also want an additional method for the jacobian the

template class VirtualLocalBasisInterface;

approach would still require to explicitly distinguish between maxDerivOrder=0 and >0 or to use exceptions.

How would you statically deal with the case that e.g. only mixed 2nd derivatives are not implemented. Of course they should be, but that's true for all derivatives.

My suggestion was to say that we expect jacobian to always be implemented - so that really is part of the interface and if a DUNE_NOT_IMPLEMENTED can then be used. VirtualLocalBasisInterface<0> would then be equal to VirtualLocalBasisInterface<0> (or <0> is just not meaningful). The case that only mixed derivatives are not implemented would not be covered by this approach and that seems to me to be a bit to fine granularity. In that case I would say the class should be derived from VirtualLocalBasisInterface<1> and then a non-interface method laplace can be added for example. Then anybody can use this function but knows that non-interface stuff is used.

If we keep the current evaluate() semantics the maxDerivOrder would not influence the static interface in any way. It would be the maximal size of the derivative vector in the virtual interface.

Another option would be to make the type of the derivative vector and not it's size a template parameter. Then one could use std::vector in the virtual interface - loosing the compile time fixed maxOrder.

Regarding <0>: Would you redict to evaluateFunction() or just leave the evaluate() method out ?

Oliver and I would also prefer Andreas suggestion as summerized below.

• There is one 'template class VirtualLocalBasisInterface'.
• It contains 'evaluateFunction(...)' AND 'evaluateJacobian(...)' for all maxDerivOrder.
• It contains 'evaluate(const Dune::array<int,k>&,...)' for k=0,...,maxDerivOrder.
• 'VirtualLocalBasisInterface' derives from 'VirtualLocalBasisInterface'
• There is only one LocalBasisTraits template

There are still the following open questions:

1. Should 'evaluate(...)' methods for k=0 always exist ? If yes we must implement them in all FEs.

2. Should maxDerivOrder be part of LocalBasisTraits or handed to the interface seperately ?

3. How should the maxDerivOrder enum be called ?

I'd suggest: @1) Yes, but provide a default default implementation in the interface using 'evaluateFunction(...)' as temp workaround @2) It should be part of the Traits, i.e. CkLocalBasisTraits becomes the new LocalBasisTraits @3) As the jacobian is always there 'maxPartialDerivativeOrder' would be more precise but the current 'diffOrder' would be nicer

I doubt this will work. If 'VirtualLocalBasisInterface' derives from 'VirtualLocalBasisInterface' and defines evaluate, it shadows all previous versions (i.e. all lower order versions) of evaluate.

With proper using statements it should work. See the attached example. The key is to have virtual inheritance and using in the interface and the imp.

Did I miss something in the example?

Attachments

• interface.cc

I've also tried this for the real interface. The recursive inheritance of interfaces and wrappers does indeed work fine and all evaluate() interface/implementation methods are visible (see the test). The attached patch contains:

1)Adapted virtual interfaces and wrappers as suggested 2)The new LocalBasisTraits (without C* prefix) 3)P1LocalBasis using LocalBasisTraits<...,2> WITH implementation of corresponding evaluate() 4)P0LocalBasis using LocalBasisTraits<...,0> WITHOUT implementation of corresponding evaluate() 5)evaluate() for order 0 redirects to evaluateFunction() (needed for P0LocalBasis) 6)Adapeted test case

Are there any objections against this change? Otherwise I'll commit it next week and adapt the LocalBasisTraits in the other FEs.

Attachments

• localfunctions.patch

If it work, then this is exactly the way to go.

Just to point out the necessary changes outside the virtual interface:

A)C*LocalBasisTraits<...> must be changed to LocalBasisTraits<...,order> B)LocalBasisTraits::diffOrder might be zero even if evaluateJacobian() is always there C)All LocalBasis implementations will need (medium-term, see 5) above) at least evaluate<0>(...)

IMHO all these changes should be done since they are reasonable and not to invasive. However, we could also avoid ... A) at the price of needlessly complicated code ... B) by changing the enum 'diffOrder' e.g. to 'maxImplementedPiecewisePartialDiffOrder' ;-) ... C) by saying evaluate<0>(...) should not exist

I commited the changes. Now even double virtualization does work.

If you did not use the virtual interface only the following changes might effect you:

1. The C*LocalBasisTraits are renamed to LocalBasisTraits If you only extract the traits from the basis and do not do template specialization for it everything should work as before.

2. LocalBasisTraits::diffOrder is now the maximal order implemented for evaluate() While diffOrder was 1 for almost all LFEs before since they implemented evaluateJacobian() it is now 0 for almost all LFEs since they do not implement evaluate().

Status changed to closed