# Solving Poisson problem for a naca airfoil

Hello everyone!
I’m trying to solve the poisson equation for a naca profile inside a rectangle. I defined everything as a 2D geometry and I can mesh it but when solving I get an error.

``````import pandas as pd
import mpmath as mp
import matplotlib.pyplot as plt
import numpy as np
from ngsolve import *
import netgen.geom2d as geom2d

mp.mp.dps = 7 #Set precision to 7 decimal places

geo = geom2d.SplineGeometry()
#To insert a rectangle in the geometry

#To insert the airfoil in the geometry
AF = pd.read_csv("Naca4412",engine = 'python', sep = '\s+', header=0, names=['x', 'y']) #Airfoil data, use Selig format .dat file
#print(AF.shape)
#print(AF)

#To visualize the airfoil given
x = AF["x"]
#print(x)
y = AF["y"]
#print(y)

#plt.plot(x,y)
#plt.show()

#To insert the airfoil 2D geometry
x = AF["x"].to_numpy() #Convert to an array
y = AF["y"].to_numpy()
all_points = list(zip(x, y))

pnts = all_points
p_list = [geo.AppendPoint(*pnt) for pnt in pnts]
#print(p_list)
L = len(p_list)
#print(L)

last_point = len(p_list)-1 #Index of the last point
#curves = [[["line", p_list[i], p_list[i + 1]],10] for i in range(last_point)]
curves = [(p_list[i], p_list[i+1], "airfoil") for i in range(last_point)]
# Avoid index out of bounds for last point
curves.append((p_list[last_point], p_list[0], "airfoil"))
#curves.append([["line", p_list[last_point], p_list[0]],11]) #Add line from last point to point 0
print(curves)

for p0,p1, bn in curves:
#geo.Append(c[0], leftdomain=0, rightdomain=1 , maxh = 0.2)
geo.Append (["line", p0, p1], bc = "airfoil", leftdomain = 0, rightdomain = 1 ,maxh = 0.2)

#To mesh
mesh = Mesh(geo.GenerateMesh(maxh = 0.1))
Draw(mesh)

#Now, to solve the poisson equation -Delta u = f with Dirichlet BC u = 0
# H1-conforming finite element space
fes = H1(mesh, order=3, dirichlet="airfoil|wall|outlet|inlet") # order is the order of the approximation function
#print(fes.ndof)

# define trial- and test-functions
u = fes.TrialFunction()
v = fes.TestFunction()

# the right hand side
f = LinearForm(fes)
f += 32 * (y*(1-y)+x*(1-x)) * v * dx # when we don't know the exact solution this is 0

# the bilinear-form
a = BilinearForm(fes, symmetric=True)

a.Assemble()
f.Assemble()

# the solution field
gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs(), inverse="sparsecholesky") * f.vec
print(gfu)

# plot the solution (netgen-gui only)
Draw (gfu)

exact = 16*x*(1-x)*y*(1-y)
print ("L2-error:", sqrt (Integrate ( (gfu-exact)*(gfu-exact), mesh)))
``````

ERROR:

``````f += 32 * (y*(1-y)+x*(1-x)) * v * dx
TypeError: __iadd__(): incompatible function arguments. The following argument types are supported:
1. (self: ngsolve.comp.LinearForm, lfi: ngsolve.fem.LFI) -> ngsolve.comp.LinearForm
2. (self: ngsolve.comp.LinearForm, lfi: ngsolve.fem.PointEvaluationFunctional) -> ngsolve.comp.LinearForm
3. (self: ngsolve.comp.LinearForm, arg0: ngsolve.comp.SumOfIntegrals) -> ngsolve.comp.LinearForm
``````

Does anyone now how to fix it? I though my boundaries weren’t correctly defined but I defined them and it still doesn’t work.

Thank you!!

You redefined the x and y function to some numpy array. Either use different variable names for that or do something like

``````import ngsolve as ngs
f += 32 * (ngs.y*(1-ngs.y) + ngs.x*(1-ngs.x)) * v * dx
``````

in your integrator to access the ngsolve x and y variables.

best
Christopher