Accessing Basis Functions (not values) in order to Compute Inner Integral Over Discontinuous Integrand in Libmesh (Read Tutorial, But Completely New to Libmesh)

[Chaztikov]

This isn't so much an issue as it is a question. Apologies, I tried to join both the libmesh-users and libmesh-devel and sent this inquiry originally to libmesh-users@lists.sourceforge.net But for some reason this did not work, so I am posting my inquiry here.

I am very new to libmesh, and in general am looking to enhance my understanding of how to use it as quickly as possible. I read through the tutorial, but my understanding of how to actually use it is very incomplete at this point.

Specifically, I am looking to simulate a solution to a nonlocal diffusion problem which involves a discontinuous kernel in the inner integral of a double integral as follows:

I would like to simulate this in 2D and 3D as well of course, but at the moment I am stuck on how to define an additional quadrature rule to handle the inner integrand. (Numerical) Mathematically, it is clear to me that the outer integral corresponds to one quadrature rule (which I think can be readily implemented in libmesh as in example 3 of the tutorial), while the inner integral requires an additional set of quadrature points (hence also basis function values at these points) which I haven't figured out how to access yet. (to clarify I mean that I haven't figured out how to access the required components in libmesh.)

Is there a way in libmesh to, given a (quadrature) point, return the value of the basis functions at that point? Is there a way to return, given an element, its neighbors? With these, I think I could do what I am trying to achieve.

Thanks for reading, and let me know if any additional explanation is required.

Best, Charlie

PS: As an additional related question, it might be necessary or useful to develop a quadrature rule for this purpose. If anyone knows of a libmesh example related to this or any of the above questions, I'd really appreciate it!

My thoughts on how to achieve part of what I want as outlined above: from fe.h, I see the following:

PPS: /**

• \returns The value of the \f$ i^{th} \f$ shape function at
• point \p p. This method allows you to specify the dimension,
• element type, and order directly. This allows the method to
• be static.

• On a p-refined element, \p o should be the total order of the element. */ static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point & p);

  /**

• \returns The value of the \f$ i^{th} \f$ shape function at
• point \p p. This method allows you to specify the dimension,
• element type, and order directly. This allows the method to
• be static.
• On a p-refined element, \p o should be the base order of the element. / static OutputShape shape(const Elem elem, const Order o, const unsigned int i, const Point & p);

I think this is what I am looking for, as far as the values of the basis functions is concerned. Am I wrong? I still need to find out how to get the neighbors given a pointer to the element elem.

[jwpeterson]

sent this inquiry originally to libmesh-users@lists.sourceforge.net But for some reason this did not work, so I am posting my inquiry here.

I think you have to join both of our email lists before you can actually send mail to them? Most discussion takes place on GitHub these days, so your post is fine.

while the inner integral requires an additional set of quadrature points (hence also basis function values at these points) which I haven't figured out how to access yet. (to clarify I mean that I haven't figured out how to access the required components in libmesh.)

In the 1D case, I suppose you could use an existing 2D quadrature rule... but in 2D that requirement presumably goes up to 4D? Also, I suspect that using "normal" quadrature rules designed for polynomials is going to be very inaccurate for this singular kernel.

Is there a way in libmesh to, given a (quadrature) point, return the value of the basis functions at that point?

You could look at the FEInterface:: functions for this.

Is there a way to return, given an element, its neighbors? With these, I think I could do what I am trying to achieve.

Sure, elem->neighbor_ptr(0), neighbor_ptr(1), etc...

[roystgnr]

In the 1D case, I suppose you could use an existing 2D quadrature rule...

Use it repeatedly, anyway - at some x the y integral is going to extend into multiple elements, and you won't want to use a designed-for-smooth-functions quadrature rule across a point where the shape function derivative has a discontinuous jump, so you'd need to patch together multiple quadrature evaluations. 2 evaluations at most, if you could swear that lambda < hmin... but it looks like lambda is fixed? You'd have arbitrarily many subintegrals if lambda doesn't decrease somehow with mesh refinement.

I suspect that using "normal" quadrature rules designed for polynomials is going to be very inaccurate for this singular kernel.

Oh, definitely. Gaussian quadrature cries and loses convergence rate from as little as a discontinuous derivative; I wouldn't expect it to accomplish anything worthwhile on a singularity. You could use Gaussian quadrature only on neighboring elements, not near the y=x part of the integral.

The exercise of deriving Gaussian quadrature rules might be instructive, though? Maybe there's some similarly standard transformation you could apply to the singular part of the kernel, then solve for the optimal quadrature points/weights for your known basis function set?

[Chaztikov]

Thanks to you both for your responses! I have already garnered some insight from your "elem->neighbor_ptr" suggestion, for instance. Indeed, lambda is fixed.

Right, I don't want to use an existing quadrature rule for smooth integrands for those involving a discontinuous kernel. I think you understand the approach I am taking initially, which is, roughly, to use existing quadrature rules in libmesh in cells of a partition of each element so as to avoid singular (outer-quadrature) points. Roughly, I am trying to determine (outer-)quadrature points for the outer integrand and then use these as the singular points 'x' to partition the element, with sub-elements dependent on inner-quadrature points for the inner integrand, using "elem->neighbor_ptr" as suggested for outer-quadpoints near the element boundary when necessary.

[Chaztikov]

I have thought about this approach, and am revising it, but this is still relevant to the question of how to develop a quadrature rule in the context of libmesh. My guess is that I want to follow the structure of, say, quadrature_gauss.C.

To solve the problem of nonlocal diffusion, I plan to implement the rules for singular kernels developed in http://www.sciencedirect.com/science/article/pii/S0021999116000176 and modify the Gauss-Legendre rule in quadrature_gauss.C to use the above.