Issues with Gridfunction generated from IntegrationRuleSpace method

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-lam*q-mu*p)*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) # 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_comp[0](mesh(p1,p2))[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]