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