Negative values of local coordiantes

General
1

Hello everybody,

I’ve encountered some unexpected behavior in NGSolve that might be of interest. When meshing global coordinates stored in a NumPy array, the resulting local coordinates can have small negative values (up to approximately -9.9e-5).

From my understanding, these local coordinates serve as barycentric coordinates for the corresponding element, which is indicated by the “nr” attribute. A slightly negative value suggests that, in a strict sense, the point lies just outside the referenced element.

Is this expected behavior? Additionally, is there any way to enforce that the local coordinates correspond exactly to the element in which the point actually is?

Greetings,
Fabian

1 Like
2

Can you post a minimal example?

3

The following code results in the local coordinates
(0.00474244, -6.08618382e-05, 0.00992934).

from netgen.csg import CSGeometry, OrthoBrick, Pnt
from ngsolve import *
from xfem import CutInfo, IF, InterpolateToP1

import numpy as np

# Background mesh
L = 1.5
max_h = 0.375
mesh_refinements = 1

background = CSGeometry()
background.Add(OrthoBrick(Pnt(-L, -L, -L), Pnt(L, L, L)))
mesh = Mesh(background.GenerateMesh(maxh=max_h, quad_dominated=False))

# Level set function
radius = 1.0
order_norm = 2
phi = (x**order_norm + y**order_norm + z**order_norm)**(1/order_norm) - radius
        
# calculating mesh refinements
refine_H1_space = H1(mesh, order=1, autoupdate=True)
refine_levelset = GridFunction(refine_H1_space, autoupdate=True)

for i in range(mesh_refinements):
    InterpolateToP1(phi, refine_levelset)
    cut_info = CutInfo(mesh, refine_levelset).GetElementsOfType(IF)
    mesh.SetRefinementFlags(cut_info)
    mesh.Refine()

global_coordinates = np.array([-0.63795623, -0.47848099, -0.61152876], dtype=np.float64, ndmin=2)
local_coordinates = mesh(global_coordinates[:,0], global_coordinates[:,1], global_coordinates[:,2])
print(local_coordinates)

minimal_working_example.py (1.0 KB)

4

Ah now I see what you mean. Have a look at meshclass.cpp line 5875 there a tolarance value is hard coded.
This allows the Find3dElement function in meshclass.cpp line 15 to find points with this tolerance where barycentric coordinates are allowed 1e-4 off. If you need this to be tightened you’d need to replicate the Find3dElement function in your own c++ library with different/no tolerance. (or self compile netgen with this value changed).

5

here we should introduce a parameter (eps or tol) to provide the tolerance from Python

6

on newest ngsolve master you can now do

local_coordinates = mesh(global_coordinates[:,0], global_coordinates[:,1], global_coordinates[:,2], tol=0)

which results in

[(0.00468158, 0.98532822, 0.00992934, 94766790997952, 0, 4163)]

Best
Christopher