Exact numerical integration of right hand side

First of all, thank you all for the cool stuff.

I got the following problem:

I want to show how the regularity of the data enters the order of convergence for higher-order methods.

Consider therefore the classical Poisson problem with a piecewise constant function f as right-hand side.
then f is nearly in H^{1/2}(\Omega). The solution u is then nearly in H^{2.5}. We, therefore, expect in the L^2 norm, for example, the rate h^\min(p+1, 2.5), where p is the polynomial order of the FEM space.

To then seen the effect the line of discontinuity should not be fitted to the mesh. (see [attachment=undefined]mesh.jpeg[/attachment]) Otherwise, the method is only limited by the order of polynomials employed.

The problem now is that in this setting the numerical integration always produces an error of order h. Due to the Strang lemma, this limits the performance to a rate of h... In 1D the problem can be resolved by a setting my own integration rule. However, I want to do it in 2D with smooth boundary and this can't be done by a special integration rule now...

I got the advice to then use the XFEM addon. And do the following:

-) Use ngsolve to setup the bilinear form and the linear form for my problem.
-) Use the ngsxfem to calculate the correct right-hand side.
-) Pick the right entries of the ngsxfem - generated right-hand side and put them at the correct place in the standard formulation.

But to be honest I have no idea how to do the last one... Is there an easier way or can someone help?

Thanks in advance,

Max

Hi Max,

You can use ngsxfem directly to compute the linearform. It's really simple .

I make a simple example that you can easily adapt to your problem. I want to have f as the source term that is 1 for x<0.5 and 0 otherwise (independent of the mesh).



You will need to install the ngsxfem package for this:
github.com/ngsxfem/ngsxfem

Some remarks:
* Importing SymbolicLFI from xfem overwrites (monkey-patches) the usual NGSolve definition. But no worries, the syntax is compatible with standard ngsolve (in fact, standard ngsolve SymbolicLFI is called if no levelset_domain is provided). It's actually just a wrapper.
* NEG is imported as a keyword for describing that you want to integrate only where the level set function is negative. If you want to integrate on other parts of the level set domain you may want to use POS or IF (zero level set) for that.
* For this linear interface the integration in the SymbolicLFI will be geometrically exact.
* For essentially arbitary level set functions the ngsxfem-integration will introduce an additional error.
* This error is of second order unless you do something about it. There are two ways to improve the situation:
1. Either you accept the second order convergence and add more refinements (for the numerical integration). For that you provide the additional key "subdivlvl" with a value > 0.
2. Or you use an isoparametric mapping, see the jupyter- and non-jupyter examples on github for that.

Best,
Christoph

Hi Christoph,

Thanks a lot! Worked out perfectly fine! I just set the additional key "subdivlvl" to deal with the second order error. Now I see the expected convergence rate.

Best,

Max