Compute the functional questions

Question

I work on the Poisson equation. I use functionals a=(u-ue)(u-ue)dx and a=(u-ue)(u-ue)dx + inner(grad(u-ue), grad(u-ue))dx (ufl form) to compute the L2 and H1 error. The situations are:
when I use P1 element, the error orders are perfect. But when I use P2 or P3 element, the H1 error order are good, the L2 error order degenerates.
I compute functional a=(u-ue)(u-ue)dx and a=inner(grad(u-ue), grad(u-ue))dx, then add them together. I found the summation does not equal the value a=(u-ue)(u-ue)dx + inner(grad(u-ue), grad(u-ue))*dx.
This is a troublesome problem!!! Any suggestions are appreciated! Thanks in advance.

code:
1. main.cpp:

include <dolfin.h>

include "Poisson.h"

include "L2error.h"

include "H1error.h"

include "Grad.h"

using namespace dolfin;

// Source term (right-hand side)
class Source : public Expression
{
void eval(Array& values, const Array& x) const
{

values[0] = -x[1]*(1-x[1])*(1-x[0]-0.5*x[0]*x[0])*exp(x[0]+x[1])-x[0]*(1-0.5*x[0])*(-3*x[1]-x[1]*x[1])*exp(x[0]+x[1]);

}
};

// Sub domain for Dirichlet boundary condition
class DirichletBoundary : public SubDomain
{
bool inside(const Array& x, bool on_boundary) const
{
return on_boundary;
}
};

// Exact solution
class ExactSolution : public Expression
{
void eval(Array& values, const Array& x) const
{
values[0] = x[0]x[1](1-x[0]/2)(1-x[1])exp(x[0]+x[1]);
}
};

int main()
{
// Create mesh and function space
int n=64;
UnitSquareMesh mesh(n, n);
Poisson::FunctionSpace V(mesh);

// Define boundary condition
ExactSolution u0;
DirichletBoundary boundary;
DirichletBC bc(V, u0, boundary);

// Define variational forms
Poisson::BilinearForm a(V, V);
Poisson::LinearForm L(V);

Source f;
L.f = f;

// Compute solution
Function u(V);
solve(a == L, u, bc);

ExactSolution ue;

// Compute L2 error
L2error:: Functional L2error_form (mesh);
L2error_form.ComputedSolution=u;
L2error_form.ExactSolution=ue;
double MyL2error = sqrt(assemble(L2error_form));
double error1 = assemble(L2error_form);

// Compute Grad energy
Grad:: Functional Grad_form (mesh);
Grad_form.ComputedSolution=u;
Grad_form.ExactSolution=ue;
double error2=assemble(Grad_form);

// Compute H1 error
H1error:: Functional H1error_form (mesh);
H1error_form.ComputedSolution=u;
H1error_form.ExactSolution=ue;
double MyH1error=sqrt(assemble(H1error_form));
double error3=assemble(H1error_form);

// print error norm
cout<<"L2 error =" <<MyL2error<<endl;
cout<<"H1 error =" <<MyH1error<<endl;
cout <<error1+error2<<endl;
cout<<error3<<endl;

return 0;
}

Poisson.ufl

element = FiniteElement("Lagrange", triangle, 2)

u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)
g = Coefficient(element)

a = inner(grad(u), grad(v))dx
L = fv*dx

Grad.ufl :

order = 4
element0 = FiniteElement("Lagrange", triangle,order)
element1 = FiniteElement("Lagrange", triangle,order+2)

ComputedSolution = Coefficient(element0)
ExactSolution = Coefficient(element1)

a= inner(grad(ComputedSolution-ExactSolution), grad(ComputedSolution-ExactSolution))*dx

H1error.ufl:

order = 4
element0 = FiniteElement("Lagrange", triangle,order)
element1 = FiniteElement("Lagrange", triangle,order+2)

ComputedSolution = Coefficient(element0)
ExactSolution = Coefficient(element1)

a=(ComputedSolution-ExactSolution)(ComputedSolution-ExactSolution)dx + inner(grad(ComputedSolution-ExactSolution), grad(ComputedSolution-ExactSolution))*dx

L2error.ufl:
order = 4
element0 = FiniteElement("Lagrange", triangle,order)
element1 = FiniteElement("Lagrange", triangle,order+2)

ComputedSolution = Coefficient(element0)
ExactSolution = Coefficient(element1)

a=(ComputedSolution-ExactSolution)(ComputedSolution-ExactSolution)dx

The above is my complete code. Thank you for help.

Answer 1

Hello,

is it possible for you to post your code? It is very hard to help without looking at the code.
what, for example, is ue? Is it a Function, an Expression....?

Also, what command do you use to compute the error norm from the form?

But as I said, anything without looking at your code I can't even make educated guesses about your problem...

best regards

Comments

I had pasted my code. Thank you for help!

---

Oh sorry I am not familiar with the c++ interface.... but the code will probably help others in understanding.

Btw, If you want the code to look "nicer" in the forum, you can just indent every line of code by 4 spaces.

I can only make some suggestions:

• Have you tried increasing the integration order?
• according to the docs (http://fenicsproject.org/documentation/dolfin/dev/python/programmers-reference/fem/norms/errornorm.html) you should use builtin functions to compute error norms because using your own functions is not stable w.r.t. roundoff errors. The link refers to the python interface but there is probably something similar for c++.

regards