# Issues with Gridfunction generated from IntegrationRuleSpace method

Hello,

I would like to use IntegrationRuleSpace to project data from a numpy array back into a Gridfunction. To test this, I have a stationary Navier Stokes example on a unit square, adapted from one of the ngsolve tutorials.

In this code, I am:

1. running Newtons method until I achieve a tolerance of 1e-13.
2. taking that low-tolerance solution, and projecting it to a uniform numpy array
3. Putting this numpy array back into a Gridfunction using the IntegrationRuleSpace method
4. Running Newtons method on this interpolated Gridfunction (I expect the tolerance to be close to 1e-13 and for the method to stop after 1 iteration. However, the tolerance is very high, around 1e-2).

While the results look good visually, the projected gfu and original gfu are very different. Why is the solution so different if I’m projecting from the same, low-tolerance gridfunction?

The attached file is ready to run. The tolerance will be printed out in the terminal for both cases. Additionally, you can uncomment lines 232-239 to visualize the numpy arrays. By default, the code will save two VTU files: The original, <1e-13 Gridfunction, and its projected Gridfunction. These two results look very similar yet produce very different results in the solver.

[code]from netgen import gui
from netgen.geom2d import unit_square
from ngsolve import *
import numpy as np
from ngsolve.comp import IntegrationRuleSpace
from ngsolve.solvers import *
import seaborn as sns
import matplotlib.pyplot as plt
mesh = Mesh(unit_square.GenerateMesh(maxh=0.05))
def SimpleNewtonSolve(gfu,a,tol=1e-13,maxits=25):
res = gfu.vec.CreateVector()
du = gfu.vec.CreateVector()
fes = gfu.space
for it in range(maxits):
print (“Iteration {:3} “.format(it),end=””)
a.Apply(gfu.vec, res)
a.AssembleLinearization(gfu.vec)
du.data = a.mat.Inverse(fes.FreeDofs()) * res
gfu.vec.data -= du
#stopping criteria
stopcritval = sqrt(abs(InnerProduct(du,res)))
print (“<A u”,it,“, A u”,it,">_{-1}^0.5 = ", stopcritval)
if stopcritval < tol:
break
mesh = Mesh(unit_square.GenerateMesh(maxh=0.05)); nu = Parameter(1)
V = VectorH1(mesh,order=3,dirichlet=“bottom|right|top|left”)
Q = H1(mesh,order=2);
N = NumberSpace(mesh);
X = FESpace([V,Q,N])
(u,p,lam), (v,q,mu) = X.TnT()
a = BilinearForm(X)
-div(u)q-div(v)p-lamq-mup)dx
gfu = GridFunction(X)
gfu.components[0].Set(CoefficientFunction((4
x
(1-x),0)),
definedon=mesh.Boundaries(“top”))

SimpleNewtonSolve(gfu,a)
vtk = VTKOutput(ma=mesh,
coefs=[gfu.components[0], gfu.components[1]],
names = [“u”, “p”],
filename=“original_gfu”,
subdivision=3)

vtk.Do()

# Project to uniform Grid

def create_np_data(gfu: GridFunction, mesh: Mesh, n: int, m: int):
“”"
Function to create numpy array from Gridfunction data.
Args:
gfu_prof_components: The original gridfunction to load the values into.
mesh: The mesh used in this simulation.
n: Desired number of discrete cells (elements) in the uniform grid (numpy array)
in the x direction.
m: Desired number of discrete cells (elements) in the uniform grid (numpy array)
in the y direction.
sol_path_str: Optional directory string for the .sol file.
Returns:
Numpy array of field data Ux, Uy, P.
“”"

# Obtain mesh dimensions from vertices.

vertices = mesh.vertices
v1, v3 = vertices[0], vertices[2]
x0, xn, y0, yn = v1.point[0], v3.point[0], v1.point[1], v3.point[1]

# Initialize numpy arrays for field data.

output_Ux = np.zeros((n,m))
output_Uy = np.zeros((n,m))
output_P = np.zeros((n,m))
output_fields = np.empty((1, 3, n, m))

# Create uniform grid points based on mesh dimensions.

x_interp = np.linspace(x0,xn,m)
y_interp = np.linspace(y0,yn,n)
gfu_comp = gfu.components
x_idx = 0
for p1 in x_interp:
y_idx = 0
for p2 in y_interp:
try:
val_Ux = gfu_comp0[0]
val_Uy = gfu_comp[0](mesh(p1, p2))[1]
val_P = gfu_comp[1](mesh(p1, p2))
output_P[y_idx, x_idx] = val_P
output_Ux[y_idx, x_idx] = val_Ux
output_Uy[y_idx, x_idx] = val_Uy
except Exception:

# and set their values to an arbitrary number

output_P[y_idx, x_idx] = -10
output_Ux[y_idx, x_idx] = -10
output_Uy[y_idx, x_idx] = -10
y_idx += 1
x_idx += 1
output_fields[0,0], output_fields[0,1], output_fields[0,2] = output_Ux, output_Uy, output_P
return output_fields

def numpy_to_gfu(np_data: np.array,gfu_init: GridFunction, mesh: Mesh, n: int, m: int):
“”"
Function to create gridfunction from numpy array.
Arguments:
np_data: Numpy array containing projected low-tol solution
gfu_init: The initial gridfunction where the simulation last left off.
mesh: Corresponding mesh
n: Desired number of discrete cells (elements) in the uniform grid (numpy array)
in the x direction.
m: Desired number of discrete cells (elements) in the uniform grid (numpy array)
in the y direction.
“”"
interp_ord = 3
for idx in range(2): # Looping through both components → u and p

# Get FES - same finite element spaces as defined above.

if idx == 0:
X = VectorH1(mesh, order=3, dirichlet=“bottom|right|top|left”)
else:
X = H1(mesh, order=2)
interp_ord -= 1

# Get the dimensions of the problem.

if len(gfu_init.components) > 0:

# Get dimension of field.

dim = gfu_init.components[idx].dim
else:
dim = gfu_init.dim
if dim >= 2:
fesir = IntegrationRuleSpace(mesh=mesh, order=interp_ord) ** dim
fesir_irs = IntegrationRuleSpace(mesh=mesh, order=interp_ord)
else:
fesir = IntegrationRuleSpace(mesh=mesh, order=interp_ord)
fesir_irs = fesir

# Working with 2d data - the fesir coordinate data always has dimension 2.

fesir_coor = IntegrationRuleSpace(mesh=mesh, order=interp_ord) ** 2

# Create IRS Gridfunction for data and corresponding coordinates.

gfuir = GridFunction(fesir)
gfuir_coor = GridFunction(fesir_coor)

# Interpolate point-values in integration points.

if len(gfu_init.components) > 0:
gfuir.Interpolate(gfu_init.components[idx])
else:
gfuir.Interpolate(gfu_init)
gfuir_coor.Interpolate(CF((x, y)))
vertices = mesh.vertices
v1, v3 = vertices[0], vertices[2]
x0, xn, y0, yn = v1.point[0], v3.point[0], v1.point[1], v3.point[1]
########################################################

# Reverse interpolate

########################################################

# Create uniform grid points based on mesh dimensions.

x_interp = np.linspace(x0, xn, m)
y_interp = np.linspace(y0, yn, n)

# Uniform indices.

x_ind = np.arange(m)
y_ind = np.arange(n)

# Coordinates from IntegrationRuleSpace.

coord_x = gfuir_coor.components[0]
coord_y = gfuir_coor.components[1]
for i in range(len(coord_x.vec)):
p1 = coord_x.vec
[i]p2 = coord_y.vec
[i]# Get corresponding indices.
x_idx = int(np.interp(p1, x_interp, x_ind))
y_idx = int(np.interp(p2, y_interp, y_ind))

if dim >= 2:

# Since this is for 2D problems, we are only considering 2 components.

gfuir.components[0].vec[i] = np_data[0][y_idx, x_idx]
gfuir.components[1].vec[i] = np_data[1][y_idx, x_idx]
else:
gfuir.vec[i] = np_data[2][y_idx, x_idx]

# Get integration rules.

irs = fesir_irs.GetIntegrationRules()
fes = X
p, q = fes.TnT()
mass = BilinearForm(p * q * dx).Assemble().mat
invmass = mass.Inverse(inverse=“sparsecholesky”)
rhs = LinearForm(gfuir * q * dx(intrules=irs))
rhs.Assemble()
gfu_init.components[idx].vec.data = invmass * rhs.vec
idx += 1
return gfu_init
n = 256
m = 256

# Create projected numpy arrays containing solutions.

np_data = create_np_data(gfu=gfu, mesh=mesh, n=n, m=m)

# plt.show()

gfu_from_np = numpy_to_gfu(np_data[0], gfu, mesh, n, m)
vtk = VTKOutput(ma=mesh,
coefs=[gfu_from_np.components[0], gfu_from_np.components[1]],
names = [“u”, “p”],
filename=“projected_gfu”,
subdivision=3)

# Exporting the results:

vtk.Do()
print(" NEW SOLVE …")

# Tolerance is still high, even though the data is coming from a low-tolerance solution

print(“Tolerance is still high, even though we’re using the low-tolerance interpolated gfu”)
SimpleNewtonSolve(gfu_from_np,a)[/i][/i][/i][/i][/i][/code]