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
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
#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 += grad(u)*grad(v)*dx
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)
Draw (-grad(gfu), mesh, "Flux")
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!!