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 = 7 #Set precision to 7 decimal places 

geo = geom2d.SplineGeometry()
#To insert a rectangle in the geometry 
geo.AddRectangle((-0.5,-0.5),(1.5,0.5),bcs=("wall", "outlet", "wall", "inlet"))

#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 

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


#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]
L = len(p_list)

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

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

#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

# 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 += grad(u)*grad(v)*dx


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

# plot the solution (netgen-gui only)
Draw (gfu)
Draw (-grad(gfu), mesh, "Flux")

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


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.