Hello everyone,
I’m trying to define a levelset function of an airfoil to use with the xfem add-on.
I have the data points from the airofil and I was able to obtain its curve using parametric spline fitting. Then, I created a function that computes the signed distance of a single point with respect to the airfoil.
Now I want to create the level set function but I always get errors related to the CoefficientFunction definition.
Here’s the code I’m using :
import numpy as np
import pandas as pd
import mpmath as mp
import matplotlib.pyplot as plt
from ngsolve import*
from netgen.geom2d import SplineGeometry
from xfem import *
from xfem.mlset import *
from scipy.interpolate import splprep, splev
from scipy.spatial.distance import cdist
import math
#To insert the airfoil in the geometry
AF = pd.read_csv("Naca0012",engine = 'python', sep = '\s+', header=0, names=['x_airfoil', 'y_airfoil']) #Airfoil data, use Selig format .dat file
x_airfoil = AF["x_airfoil"]
y_airfoil = AF["y_airfoil"]
x_airfoil = AF["x_airfoil"].to_numpy() #Convert to an array
y_airfoil = AF["y_airfoil"].to_numpy()
all_points = list(zip(x_airfoil, y_airfoil))
airfoil_points = all_points
#print(airfoil_points)
# Ensure the shape is closed by repeating the first point at the end
x_airfoil = np.append(x_airfoil, x_airfoil[0])
y_airfoil = np.append(y_airfoil, y_airfoil[0])
#print(x_airfoil)
# Parametric spline fitting
tck, u = splprep([x_airfoil, y_airfoil], s=0, per=True)
# Evaluate the spline at more points for a smooth curve
u_fine = np.linspace(0, 1, 1000)
x_fit, y_fit = splev(u_fine, tck)
boundary_points = np.column_stack((x_fit, y_fit))
#To plot
plt.plot(x_airfoil, y_airfoil, 'o', x_fit, y_fit)
labels = ['Data points', 'Spline fit']
plt.legend(labels)
plt.show()
# Print the spline representation (tck)
#print("Knots:", tck[0]) #the knot vector
#print("Coefficients:", tck[1]) #the coefficients (control points)
#print("Degree:", tck[2]) #the degree of the spline
#Background mesh and geometry
geo = SplineGeometry()
ll = (-0.2, -0.2)
ur = (1.2, 1.2)
geo.AddRectangle(ll, ur, bcs=("wall", "outlet", "wall", "inlet "))
mesh = Mesh(geo.GenerateMesh(maxh=0.1))
# Define the level set function
#First, define a function that given an (x,y) input determines the signed distance to
#the boundaries of the airfoil.
def signed_distance (x_given,y_given):
#Signed distance function for a single point wrt the airfoil
distances = cdist([[x_given,y_given]], boundary_points, metric='euclidean') #For the euclidean distance between 2 points
min_distance_index = np.argmin(distances)
closest_point = boundary_points[min_distance_index]
print("Closest point is: ", closest_point)
y_point = closest_point[1]
min_distance = np.min(distances)
#Determine the sign
if x_given < 0 or x_given > 1 or abs(y_given) > abs(y_point):
sign = 1
elif x_given > 0 and x_given < 1 and abs(y_given) < abs(y_point):
sign = -1
else:
sign = 0 #The point is on the boundary
return sign * min_distance
print("Minimum distance is: ", signed_distance( 0.177579 , 0.055886))
#Create a CoefficientFunction for the levelset
#help(CoefficientFunction)
def levelset():
return signed_distance(x,y)
# Visualization
Draw(levelset, mesh, "levelset")
File "/home/angela/Thesis/XFEM/Poisson/AirfoilPoisson.py", line 97, in <module>
Draw(levelset, mesh, "levelset")
TypeError: Draw(): incompatible function arguments. The following argument types are supported:
1. (cf: ngsolve.fem.CoefficientFunction, mesh: ngsolve.comp.Mesh, name: str, sd: int = 2, autoscale: bool = True, min: float = 0.0, max: float = 1.0, draw_vol: bool = True, draw_surf: bool = True, reset: bool = False, title: str = '', number_format: str = '%.3e', unit: str = '', **kwargs) -> None
2. (gf: ngsolve.comp.GridFunction, sd: int = 2, autoscale: bool = True, min: float = 0.0, max: float = 1.0, **kwargs) -> None
3. (mesh: ngsolve.comp.Mesh, **kwargs) -> None
4. (arg0: object) -> None
Invoked with: <function levelset at 0x7db99cab7d00>, <ngsolve.comp.Mesh object at 0x7db9993f6f70>, 'levelset'
How can I define the levelset function?
I hope you can help me, thank you!!