Computing Error Norms (C++)

This is the way I did it. It worked for me. I use something different now, because we decided to calculate the error at a different point in our process, but when I was using this it worked fine.

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

uh = Function(element)
u = Function(element)
e = uh - u

M = dot(e, e)*dx

then in the code:

 ENorm::Functional M;

 M.u = Phi0;
 M.uh = Phi1;

 double errorSqr = assemble(M);
 double error = sqrt(errorSqr);

where Phi0 and Phi1 are dolfin::Functions that I had previously solved for. They are both defined over the same function space.

This is obviously least squares. There are built in methods to find error norms if you use python, but as far as I know, the c++ interface hasn't been built yet. However, if you can find the code in the dolfin files, you can use parts of it. That's what I was doing when I was experimenting with different error norms.

I would also like to know if there is a c++ built in way to do this.

I did this a while ago, not sure if the update has any effect on the code.

This isn't a full answer but hopefully it helps.

Thanks Andy R Terrel, that solved my question.

Also, thanks, Phil, your way was more compact than mine.

On Wed, Dec 30, 2009 at 02:00:23PM -0000, Andy R Terrel wrote:
> Question #95601 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/95601
>
> Status: Open => Answered
>
> Andy R Terrel proposed the following answer:
> FFC will take care of interpolation of the computed function for you.
> So you combine your two forms into one:
>
> Error.ufl
> =======
> e3 = FiniteElement("Lagrange", "triangle", 3)
> e10 = FiniteElement("Quadrature", "triangle", 10)
> v = Coefficient(e3)
> w = Coefficient(e10)
> M= ((w - v) * (w - v)) * dx
>
>
> This can be done for any functional (so all your seminorms and what
> not). I don't think this is in dolfin, I think a long time ago
> there were things like this but because the need to support so many
> different combinations of finite elements it got yanked.

Yes, this is something that is not handled automatically in C++
since it would require pregeneration of large amounts of code to cover
different function spaces.

Andy's solution is good, but you might run into stability problems
(round-off errors) for higher order function spaces.

Look at errornorm.py for a better way to compute the error:

  """Compute the error e = u - uh in the given norm. The parameter k
  denotes the degree of accuracy (degree of piecewise polynomials
  approximating u an uh).

  In simple cases, one may just define

    e = u - uh

  and evalute for example the square of the error in the L2 norm by

    e = u - uh
    assemble(e*e*dx, mesh)

  However, this is not stable w.r.t. round-off errors considering
  that the form compiler may expand the expression above to

    u*u*dx + uh*uh*dx - 2*u*uh*dx

  and this might get further expanded into thousands of terms for
  higher order elements. Thus, the error will be evaluated by adding
  a large number of terms which should sum up to something close to
  zero (if the error is small)."""

--
Anders