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.
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")