High Accuracy Interpolation from Coefficient Function to GridFunction

Main Forum Old-Forum - do not post here
1

I have a vector field which I wish to split up into two fields: 1) unit length direction vector field, 2) vector magnitude scalar field.

I am able to calculate both of these values and get the expected result as a CoefficientFunction, confirmed numerical and visually when graphing the results.

My problem is that I need to then use the vector field for further calculations, and need to turn it back into a GridFunction in order to perform operations on it (e.g. div() and curl()).

I have attached a simplified, and self-contained, version of the code to this post. I have tried gfu.Set(), gfu.Interpolate(), and solving a system manually and setting it to gfu.vec.data.

Neither produce a result as accurate as the CoefficientFunction itself. gfu.Interpolate() does the best job, but I haven’t figured out how to increase the accuracy of the interpolation.

Any help would be great.

Alex

from ngsolve import *
from netgen.csg import unit_cube


def integrate(gfu_solution, gfu_target, fes):
    u, v = fes.TnT()
    mass = BilinearForm(u * v * dx, symmetric=True).Assemble().mat
    invmass = mass.Inverse(inverse="sparsecholesky")

    # important: integrate right hand side in points of fesir
    rhs = LinearForm(gfu_target * v * dx)

    rhs.Assemble()
    gfu_solution.vec.data = invmass * rhs.vec


# Create mesh
mesh = Mesh(unit_cube.GenerateMesh(maxh=0.3))

# Create function space
fes_v = VectorH1(mesh,  order=2)

# Create velocity field
v = GridFunction(fes_v)
v.Set((x*(1-x), y*(1-y), z*(1-z)))

# Create zero vector to use for n if v_mag is approx zero
v_zero = GridFunction(fes_v)

# Calculate magnitude of velocity
v_mag = sqrt(InnerProduct(v, v))

# Create orientation vector
# Using IfPost to make sure we're not dividing by zero
n_coef = IfPos(v_mag - 1e-6, v / v_mag, v_zero)

# Turn coefficient function into grid function
# TODO: This introduces errors. what's a better way of doing it?
n1 = GridFunction(fes_v)
n2 = GridFunction(fes_v)
n3 = GridFunction(fes_v)
# Option 1: Worst performance
n1.Set(n_coef)
# Option 2: Better performance than 1, but still not great
n2.Interpolate(n_coef)
# Option 3: Not better than Option 1
integrate(n3, n_coef, fes_v)

# Output direction and magnitude of velocity
VTKOutput(ma=mesh,
          coefs=[v, v_mag, n_coef, n1, n2, n3],
          names=['velocity', 'velocity_magnitude', 'director_coef',
                 'director_Set', 'director_Interpolate', 'director_integrate_soln'],
          filename='director_test',
          subdivision=1).Do()

print("Done")
2

Hi Alex,

if you normalize the function of your example, you get discontinuities at the corners, and you cannot approximate it well by polynomial finite element spaces.

If you wouldn’t have these singularities, you could improve accuracy by higher order spaces for the interpolation.

You can do some calculus to obtain an expression for the divergence of the normalized vector field, from the full gradient of the field itself (away from the singularities).

Joachim

3

Hi Joachim,

Thank you for the reply, I was able to address the issue using the information you provided.

For anyone else who comes across this post. The solution I ended up going with was to change the calculation slightly so that there is a smooth transition (over the length of an element or two) in the direction vector from the zero vector to unit vector.

I still calculate n_coeff as above

n_coef = IfPos(v_mag - 1e-6, v / v_mag, v_zero)

But now I multiply it by a scalar field derived from v_mag, and this new scalar field is smooth and continuous. The result being that my new values of n (n= n_coeff * S) will be smooth once I Interpolate it onto a GridFunction.

Alex