Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

High Accuracy Interpolation from Coefficient Function to GridFunction

More
2 years 3 months ago #4130 by alex.v
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
Code:
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")
More
2 years 3 months ago #4132 by joachim
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
More
2 years 3 months ago #4139 by alex.v
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
Code:
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
Time to create page: 0.119 seconds