Mapped Integration Points

More
4 years 5 months ago #2766 by philipp
Dear NGSolve community,

I have a very basic question concerning the calculation of mapped integration points like so
Code:
mp = mesh(0.5, 0).
For the mesh generated from the unit_square
Code:
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1)),
I would expect the coordinates
Code:
mp.pnt[0], mp.pnt[1]
of mp to be almost (0.5,0). Instead I get
Code:
(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,
Philipp
More
4 years 5 months ago #2767 by joachim
Replied by joachim on topic Mapped Integration Points
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
Code:
x(mp), y(mp)
to the the global coordinates.

Joachim
The following user(s) said Thank You: philipp
More
4 years 5 months ago - 4 years 5 months ago #2768 by christopher
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:
Code:
from ngsolve.fem import BaseMappedIntegrationPoint mp = mesh(0.5,0.) mapped_mp = BaseMappedIntegrationPoint(mp) print(mapped_mp.point)
I'll let you know when we have fixed this in a nicer way.
Best
Christopher
Last edit: 4 years 5 months ago by christopher.
More
4 years 5 months ago - 4 years 5 months ago #2769 by philipp
Replied by philipp on topic Mapped Integration Points
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 ngsolve.org/docu/latest/i-tutorials/unit...valuate-the-function it says that mesh(...) returns a mapped integration point) but holds the reference coordinates?

Edit:
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,
Philipp
Last edit: 4 years 5 months ago by philipp.
More
4 years 5 months ago - 4 years 5 months ago #2773 by christopher
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):
Code:
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.
Best
Christopher
Last edit: 4 years 5 months ago by christopher.
The following user(s) said Thank You: philipp
More
4 years 5 months ago - 4 years 5 months ago #2776 by philipp
Replied by philipp on topic Mapped Integration Points
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,
Philipp
Last edit: 4 years 5 months ago by philipp.
Time to create page: 0.119 seconds