Mapped Integration Points

Dear NGSolve community,

I have a very basic question concerning the calculation of mapped integration points like so

mp = mesh(0.5, 0).

For the mesh generated from the unit_square

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1)),

I would expect the coordinates

mp.pnt[0], mp.pnt[1]

of mp to be almost (0.5,0). Instead I get

(0.9999999999999445, 5.551115123127694e-14)

I tried to dig into the meshing code but found myself hopelessly lost in it to be honest.
What goes on here? To my understanding, mp differs from the coordinates given to mesh
only in the sense that it is mapped to the next mesh point and transformed if there is any transformation of the domain. Since I have not registered any transformation, I would expect to be working directly on the unit square.

Sorry for the basic question, but I couldn’t find an answer here or in the docs.

All the best,

the MeshPoint knows the element number the point belongs to, and the coordinates on the reference elements (essentially, the barycentric coordinates).

Thats the information you need to evaluate all kind of CoefficientFunctions

You can write

x(mp), y(mp)

to the the global coordinates.


you are right, this is not nice… MeshPoint does only hold the reference coordinates and pnt returns them… but it is neither documented nor described how to get the physical ones. We will have a discussion about this.
For now you can map the reference point using the constructor of the BaseMappedIntegrationPoint which you have to import from ngsolve.fem:

from ngsolve.fem import BaseMappedIntegrationPoint

mp = mesh(0.5,0.)
mapped_mp = BaseMappedIntegrationPoint(mp)

I’ll let you know when we have fixed this in a nicer way.

Hello Joachim, Hello Christopher,

thank you very much for your response.
So do I understand it right that mesh point is a subclass of mapped integration point
(here 1.2 CoefficientFunctions — NGS-Py 6.2.2302 documentation it says that mesh(…) returns a mapped integration point) but holds the reference coordinates?

I looked into the cpp source and found that the MeshPoint struct defined in fem/intrule.hpp (which is then registered as the python class MeshPoint) is not derived from BaseMappedIntegrationPoint.
But nevertheless I can call a coefficient function that takes BaseMappedIntegrationPoint with mp
which is supposed to be a MeshPoint.
Could you please elaborate on the cast that is going on here?

From a pure mathematical point of view: I understand that you want to integrate coefficient functions on a fixed reference domain but evaluation of coefficient functions would rather be wanted in the physical domain, right?

Sorry for all the questions. I am currently implementing my own coefficient function and really want to make sure I understand what I am doing.

All the best,

Ok, so the reason for this confusion thing is that numpy only allows very simple types (POD) to be used in numpy arrays. We wanted to be able to have the call operator of mesh to be able to use numpy arrays and produce mesh points efficiently in numpy arrays without the python overhead. This is why MeshPoint was created. It is a very simple struct with all the needed information of a MappedIntegrationPoint but without all the virtual functions and stuff.
This allows things like this very efficiently (the function evaluation is even parallel with TaskManager):

vals = func(mesh(np.linspace(0,1,100), 0.5))

to work with the downside of not having the full information of MappedIntegrationPoints (like eltrafo,…) right at hand. But they can be created using the BaseMappedIntegrationPoint constructor.
The physical domain values can be created using the coordinate coefficientfunctions as well like Joachim explained.

Hope this explains what is happening.

I’ve added a docstring explaining the behavior in the current master.

Ah, that clears things up. Thanks for the detailed description.
I also understand now how a CoefficientFunction can directly be called
with a MeshPoint although its C++ code only has Evaluate methods
taking BaseMappedIntegrationPoint: This should be due to the
constructor added in the python interface to the python BaseMappedIntegrationPoint
class that takes MeshPoint as argument.

I think the additional hint is a good idea. With this knowledge working with
the integration points is very smooth.

All the best,