Defining a CoefficientFunction Using a Mapping Instead of an Explicit Expression

General
1

Hello,

Is it possible to define a CoefficientFunction in NGSolve using a mapping instead of an explicit expression?

Suppose I have a real-valued mapping defined on a domain Ω, meaning that for every point (x, y) ∈ Ω, I know the corresponding value f(x, y) ∈ R.

However, f(x, y) involves complex pointwise operations. For example,
(A) it is a piecewise function, or
(B) it is a composition of a smooth function R^2 → R and a one-dimensional finite element function.

As a result, f(x, y) does not have an explicit formula, so I cannot define it using a mathematical expression. However, I do know the function values at any given point (x, y) ∈ Ω.

In this case, is there a way to define a CoefficientFunction in NGSolve based only on known function values at each point, without requiring an explicit mathematical expression?

Any suggestions would be greatly appreciated. Thank you very much!

2

you can do all that with ngsolve expressions, but maybe it would be easier to showcase in an example? What exactly to you have in mind?
piecewise function per domain or on a grid? How and what kind of composition? where do you need it and how efficient do you need evaluation…

3

Hello,

I’ve been having the same problem. What I found is the function VoxelCoefficient(p1, p2, values) where p1 are the starting coordinates of your quadrature, p2 the ending coordinates, and values ​​the values ​​of the function in your quadrature. I haven’t been able to get it to work, so I’ll stay and read more comments. If you could tell me if it worked for you and how you configured your gridfunction based on your Voxelcoefficient, I would greatly appreciate it.

Best regards and have a nice day.

4

Got it, thank you very much! Let me describe my case more specifically.

Consider the domain Ω = [-1,1]×[-1,1], which is triangulated into a mesh T_h. On Ω, I have a smooth function:

(x,y) ∈ [-1,1]×[-1,1] → arcsin(y) ∈ [-π/2, π/2].

Then, the interval [-π/2, π/2] is also meshed with a 1D mesh of size H, and a finite element space is defined on it:

fes1D = H1(mesh1d, order=order_FEM)
wh = GridFunction(fes1D)  # a one-dimensional finite element solution

Now, I would like to define a CoefficientFunction on Ω that corresponds to the mapping:

(x,y) ∈ [-1,1]×[-1,1] → wh(arcsin(y)) ∈ R

Would it be possible to implement this in NGSolve? Thanks again for your help!

5

Hello,

Thanks for sharing your approach! What I currently know is something that Prof. Joachim Schöberl suggested a few days ago. It is possible to define barycentric coordinates as follows:

fes = Discontinuous(H1(mesh, order=1))
lam1 = GridFunction(fes)
lam2 = GridFunction(fes)
lam3 = GridFunction(fes)

lam1.vec.FV().NumPy()[0::3] = 1
lam2.vec.FV().NumPy()[1::3] = 1
lam3.vec.FV().NumPy()[2::3] = 1

Using these, one can construct linear interpolation via nodal values:

f_h = f(p_1) * lam1 + f(p_2) * lam2 + f(p_3) * lam3

where ( f(p_1), f(p_2), f(p_3) ) are finite element functions in L2(mesh,order=0), representing values at each element’s nodes.

6

This idea can be extended to higher-order approximations. For example,

N1 = lam1**2 - lam1*lam2 - lam1*lam3
N2 = lam2**2 - lam1*lam2 - lam2*lam3
N3 = lam3**2 - lam1*lam3 - lam2*lam3
N12 = 4.0 * lam1 * lam2
N23 = 4.0 * lam2 * lam3
N13 = 4.0 * lam1 * lam3

Then, a higher-order interpolation can be expressed as:

f_h = f(p_1) N_1 + f(p_2) N_2 + f(p_3) N_3 + f(p_{12}) N_{12} + f(p_{13}) N_{13} + f(p_{23}) N_{23}

However, this is still an interpolation approximation, which introduces interpolation errors. Ideally, I would like it to be as accurate as possible.

Moreover, using this approach, writing out higher-order nodal basis functions becomes quite complex. Additionally, I am not sure whether this method would lead to instability for arbitrary higher-order approximations.

7

There were some small things missing for the creation of occ 1d meshes which are now on the current master pushed. With this this code is running and I think does what you were aiming for:


from netgen.occ import *
from ngsolve import *
from ngsolve.fem import CoordinateTrafo

r = Rectangle(2,2).Face().Move((-1,-1,0))
r.edges.Min(X).name = "left"
r.edges.Max(X).name = "right"
r.edges.Min(Y).name = "bottom"
r.edges.Max(Y).name = "top"

mesh2d = Mesh(OCCGeometry(r, dim=2).GenerateMesh(maxh=0.1))

edge = Segment((-pi/2,0,0), (pi/2,0,0))
edge.name = "edge"

mesh1d = Mesh(OCCGeometry(edge, dim=1).GenerateMesh(maxh=0.1))
fes2d = H1(mesh2d, order=2)
fes1d = H1(mesh1d, order=2)

gf1d = GridFunction(fes1d)
gf1d.Set(cos(x))

gf2d = GridFunction(fes2d)

# to cut of numeric values > 1 that lead to nans...
asin_no_nan = IfPos(y**2-1, 0, asin(y))
f = gf1d(CoordinateTrafo(CF((asin_no_nan,0,0)), mesh1d.Materials(".*")))
gf2d.Set(f)
Draw(gf2d)

if you created the 1d mesh by hand it should already work with current releases.
best

8

for performance reasons (this evaluate with CoordinateTrafo is a bit expensive) it is most of the times good to interpolate the transformed function into a gf on the mesh before doing further things.

9

Got it, thank you very much! Currently, I created the 1D mesh by hand. In my current version of NGSolve (6.2.2304), running

gf2d.Set(f)

causes the Jupyter kernel to crash.

I will update my NGSolve to the latest version and try this approach again. Thanks again for your guidance and help!

10

Yes this definitely needs new version, I think even 2501 at least.