[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)
a += (nu*InnerProduct(grad(u),grad(v))+InnerProduct(grad(u)u,v)
-div(u)q-div(v)p-lamq-mup)dx
gfu = GridFunction(X)
gfu.components[0].Set(CoefficientFunction((4x(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)
Exporting the results:
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:
Catch any points that might not be in the mesh,
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))
Fill gfuir data with numpy data.
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)
Uncomment to view numpy arrays:
vis = sns.heatmap(np_data[0][0], vmin=np.min(np_data[0][0]), vmax=np.max(np_data[0][0]))
plt.show()
vis = sns.heatmap(np_data[0][1], vmin=np.min(np_data[0][1]), vmax=np.max(np_data[0][1]))
plt.show()
vis = sns.heatmap(np_data[0][2], vmin=np.min(np_data[0][2]), vmax=np.max(np_data[0][2]))
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]