Interior/Exterior facet coordinate bug?

[Ido Akkerman]

Hello,

I think I found a bug in dolfin, but I might also a mistake on my side.

I am trying to play around with weak enforcement of Dirichlet boundary conditions and for this I need to evaluate an expression on the boundary. It turns out the expression does not get evaluated at the expected coordinates. It seems there is a mixup of Interior/exterior facets and boundary indexing. When replacing the ds term for a dS term things look better, interior facets are visited but not all of them....

Below is a minimal code showing the problem.

If somebody can give me some pointers i'd like to give it a stab myself.

Thank you,

Ido

PS: I tried the mailing list but that bounced.

test.ufl:

# Define function spaces
Q = FiniteElement("Lagrange", triangle,  1)

# Define trial and test functions
w   = TestFunction(Q)
phi = Coefficient(Q)
g   = Constant(triangle)

# Residual
F = w*(g-phi)*ds + w*phi*dx
##F = jump(w)*avg(g)*dS + w*phi*dx # also funny

# TrialSpace
dw = TrialFunction(Q)

# Jacobian
J = derivative(F,phi,dw)

#include <dolfin.h>
#include "test.h"

main.cpp:
using namespace dolfin;

// Define boundary flow field
class BndValue : public Expression
{
public:

 BndValue() {}

 void eval(Array<double>& values, const Array<double>& x) const
 {
   values[0] = x[0]*x[0];
   std::cout<<"x = "<<x[0]<<" "<<x[1]<<std::endl;
 }
};

int main()
{

 // Create mesh
 RectangleMesh mesh (0,0,1.0,1.0,1,1);

 // Create function spaces and functionals
 test::FunctionSpace X(mesh);
 test::JacobianForm  J(X, X);
 test::ResidualForm  F(X);

 // Create functions
 Function phi(X);

 // Create coefficients
 BndValue bv;

 // Set coefficients
 F.g   = bv;
 F.phi = phi;

 // Initial condition
 Constant zero(0.0);
 phi.interpolate(zero);

 // Solve
 solve(F == 0, phi, J);

 return 0;
}

[Ido Akkerman]

I figured out what was the issue.

The expression gets interpolated on the FE space before actually being used by the Form. This is not what I expected, and also what seems undesirable from a mathematical perspective (programming might be a different thing)

[Prof Garth Wells]

You've defined bv to be a constant (see g = Constant(triangle)), but in its definition on the C++ it varies quadratically.

If you define it on a P2 space, it will be represented exactly.

If you want a function to be evaluated at quadrature points, use a "quadrature space".

[Ido Akkerman]

How do I make a quadrature space? By the way the quadratic was just as a test, I would like to evaluate an arbitrary function at the quad points. Is this still possible with the quadrature space?

[Prof Garth Wells]

Take a look in the FEniCS Book (plasticity chapter) for quadrature elements.

In general, you either make a interpolation error (interpolating functions into a space, but the functions can be integrated exactly) or an integration error (evaluating a function only at discrete points).

[Ido Akkerman]

Thanks. I can get it to work for a volume (dx) term. But for a boundary term (ds) I get the error "Points must be equal to coordinates of quadrature points".

I agree with your comment that in a sense you are trading one error for the other. Which might under normal circumstances be of equal order. However, in my case it turned out that the interpolation is on the volume element but the integration is on the boundary, leading to very different results.

[Kristian Ølgaard]

Quadrature elements are not supported on ds. (See https://bugs.launchpad.net/ffc/+bug/488661) It is not clear how to, correctly, represent the same function using one integration scheme for dx and another for ds. Perhaps the best solution is to use two separate functions in the form?

[Ido Akkerman]

Perhaps a quadrature element is not the best solution for my problem. I have and expression which is valid for every point, so the fact that quadrature space reduces this to quad points is not necessary. Isn't there a straightforward way to get an analytical expression in UFL/FFC? It sounds trivial, so perhaps a RTFM question?

[Martin Sandve Alnæs]

You can express your analytical expression in the form file.

cell = triangle
x = cell.x
g = x[0]*x[1]
...

[Ido Akkerman]

Thanks, that works. But I have reservations regarding its useful in generic use cases; this means the actual definition of boundary conditions and force/source functions trickles down to the UFL level.

I would expect a UFL/FFC feature that is similar to the following: In UFL define an Expression,

 g = Expression()

and then in the C++ code specify this,

 L.g = g;

where g was an Expression instance.

Apparently this feature is not there? Should I try to include it?

[Ido Akkerman]

Basically a user defined MathFunction.

[Kristian Ølgaard]

In your test.ufl file do:

Q2 = FiniteElement("Lagrange", triangle,  2)
g   = Coefficient(Q2)

[Ido Akkerman]

But that will only work because my example function is quadratic, right? What if I have 1/x[0]? Or something more elaborate?

[Prof Garth Wells]

This is a thread for the mailing list. It should be moved over.