From 3cb332af03364386513f038f861ccf4c7fb1e8ee Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Fri, 13 Sep 2024 18:35:16 +0200
Subject: [PATCH 01/18] rm deps

---
 Project.toml | 13 +++++--------
 1 file changed, 5 insertions(+), 8 deletions(-)

diff --git a/Project.toml b/Project.toml
index e2530df..2571368 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,13 +6,12 @@ version = "0.1.2"
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Chain = "8be319e6-bccf-4806-a6f7-6fae938471bc"
+CubicSplines = "9c784101-8907-5a6d-9be6-98f00873c89b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
-GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationBBO = "3e6eede4-6085-4f62-9a71-46d9bc1eb92b"
@@ -29,16 +28,14 @@ Documenter = "0.27.24"
 DocumenterTools = "0.1.16"
 Downloads = "1.6.0"
 FileIO = "1.16.0"
-Gridap = "0.17.17"
-GridapGmsh = "0.6.1"
-Plots = "1.38.16"
-Revise = "3.5.1"
-XLSX = "0.8.4, 0.9"
-julia = "1.8.2"
 Optim = "1.7.6"
 Optimization = "3.15.2"
 OptimizationBBO = "0.1.4"
 Parameters = "0.12.3"
+Plots = "1.38.16"
+Revise = "3.5.1"
+XLSX = "0.8.4, 0.9"
+julia = "1.8.2"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

From 0ead29f43d803261c1cf445a3728d55f73b9125a Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Fri, 13 Sep 2024 18:37:01 +0200
Subject: [PATCH 02/18] updates

---
 src/AirfoilGmsh.jl    |  1 +
 src/CreateGeoFile.jl  | 16 ++++++++----
 src/GmshUtils.jl      | 17 ++++++------
 src/MapLines.jl       | 61 +++++++++++++++++++++++++++++++------------
 src/WriteFileUtils.jl |  2 +-
 5 files changed, 66 insertions(+), 31 deletions(-)

diff --git a/src/AirfoilGmsh.jl b/src/AirfoilGmsh.jl
index 21c0612..61863a2 100644
--- a/src/AirfoilGmsh.jl
+++ b/src/AirfoilGmsh.jl
@@ -10,6 +10,7 @@ using Plots
 using Optim
 using Optimization, OptimizationBBO
 using Parameters
+using CubicSplines
 
 export from_url_to_csv
 include("ReadWeb.jl")
diff --git a/src/CreateGeoFile.jl b/src/CreateGeoFile.jl
index 226994f..8dc25c8 100644
--- a/src/CreateGeoFile.jl
+++ b/src/CreateGeoFile.jl
@@ -22,7 +22,8 @@ It is possible to create a mesh with the following options:
 
 """
 function create_geofile(filename::String; Reynolds = -1, h0 = -1, leading_edge_points = Int64[], trailing_edge_points = Int64[], chord=1.0, dimension=2, elements = :QUAD, open_geo = false)
-    
+
+
 refinement_params = refinement_parameters(Reynolds, h0, chord)
     
 if elements == :QUAD || elements == :HEX
@@ -49,9 +50,14 @@ Airfoil.points.trailing_edge[Airfoil.sharp_idx]
 N_edge = 5
 # Create airfoil lines
     
-spline_airfoil_top = addSpline(Airfoil.points.trailing_edge[1] : Airfoil.points.leading_edge[1], Lines, io)[end][1]
-spline_airfoil_le = addSpline(Airfoil.points.leading_edge[1] : Airfoil.points.leading_edge[2], Lines, io)[end][1]
-    
+spline_airfoil_top = addSpline(collect(Airfoil.points.trailing_edge[1] : Airfoil.points.leading_edge[1]), Lines, io;all=true)[end][1]
+
+spline_airfoil_le = addSpline(collect(Airfoil.points.leading_edge[1] : Airfoil.points.leading_edge[2]), Lines, io; all=true)[end][1]
+
+Airfoil.points.leading_edge[1] 
+Airfoil.points.leading_edge[2]
+Lines
+
 if ! is_sharp(Airfoil)
     spline_airfoil_te = addLine(Airfoil.points.trailing_edge[1], Airfoil.points.trailing_edge[2], Lines, io)[end][1]
     spline_airfoil_bottom = addSpline(Airfoil.points.leading_edge[2] : Airfoil.points.trailing_edge[Airfoil.sharp_idx] , Lines, io)[end][1]
@@ -154,7 +160,7 @@ point7 = addPoint("L", "-L* " * string(x_tmp) * "*Sin(AoA) + " * string(y_tmp) *
     l33r = addLine(point3, point3r, Lines, io)[end][1]
     l44r = addLine(point4, point4r, Lines, io)[end][1]
     
-
+    
     loop1 = LoopfromPoints([point1, point1r, point2r, point2], Lines)
     loop1r = LoopfromPoints([point1r, Airfoil.points.leading_edge[1], Airfoil.points.leading_edge[2], point2r], Lines)
     
diff --git a/src/GmshUtils.jl b/src/GmshUtils.jl
index 8d5edc5..2c2d81b 100644
--- a/src/GmshUtils.jl
+++ b/src/GmshUtils.jl
@@ -51,13 +51,13 @@ function addLine(a1, a2, Lines::Vector{Vector}, io::IOStream; tag="")
     return Lines
 end
 
-function addSpline(a, Lines::Vector{Vector}, io::IOStream; tag="")
+function addSpline(a, Lines::Vector{Vector}, io::IOStream; tag="", all=false)
     nn = length(Lines) + 1
-    if typeof(a) <:Vector
+    if all==false
         str_tmp = "Spline($nn) = {$(a[1]), $(a[2])};\n"
         push!(Lines, [nn, a[1][1], a[2], tag])
     else
-        str_tmp = "Spline($nn) = {$a};\n"
+        str_tmp = "Spline($nn) = {$(string(a)[2:end-1])};\n"
         push!(Lines, [nn, a[1], a[end], tag])
     end
     write(io, str_tmp)
@@ -83,11 +83,10 @@ function LoopfromPoints(a::Vector{Int64}, Lines::Vector{Vector})
     push!(lines_id, LinefromPoints(a[1], a[2], Lines))
     loop = Any[]
     push!(loop, lines_id[1])
-    #the second node
 
+    #the second node
     loop[end] > 0 ? ctrl_sing = 2 : ctrl_sing = 1
     count = 2
-
     while count <= length(a)
 
         loop[end] > 0 ? ctrl_sing = 2 : ctrl_sing = 1
@@ -131,7 +130,7 @@ function LinefromPoints(p1::Int64, p2::Int64, Lines::Vector{Vector})
     if line_found
         return line
     else
-        return "Line not found"
+        return @error("Line $(p1) $p2 not found")
     end
 end
 
@@ -153,13 +152,13 @@ function addPlaneSurface(a, Surfaces::Vector{Vector}, io::IOStream)
     push!(Surfaces, [nn, a])
 end
 
-function TransfiniteCurve(curves::Vector, nodes::String, progression::Union{Float64,String}, io::IOStream)
+function TransfiniteCurve(curves::Vector, nodes::String, progression::Union{Float64,String}, io::IOStream; ptype="Progression")
     str_curves = "$(curves[1])"
     for i = 2:1:length(curves)
         str_curves = str_curves * ", $(curves[i])"
 
     end
-    str_tmp = "Transfinite Curve {$str_curves} = $nodes Using Progression $progression; \n"
+    str_tmp = "Transfinite Curve {$str_curves} = $nodes Using $ptype $progression; \n"
     write(io, str_tmp)
 
 end
@@ -194,7 +193,7 @@ end
 
 function addPhysicalGroup(name::String, entities::Vector, type::String, PhysicalGroups::DataFrame, io::IOStream; add=false)
 
-    if type != "Point" && type != "Curve" && type != "Surface"
+    if type != "Point" && type != "Curve" && type != "Surface" && type != "Volume"
         error("Type of Physical Group not recognized, available:Point, Curve, Surface")
     end
 
diff --git a/src/MapLines.jl b/src/MapLines.jl
index 38b5234..043152f 100644
--- a/src/MapLines.jl
+++ b/src/MapLines.jl
@@ -1,4 +1,4 @@
-function map_entities(airfoil::AirfoilParams, PhysicalGroups::DataFrame, io::IOStream)
+function map_entities(airfoil::AirfoilParams, PhysicalGroups::DataFrame, io::IOStream; pfc=false)
     
     N_airfoil_points = airfoil.points.num
     sharp_end = airfoil.sharp_end
@@ -10,11 +10,11 @@ function map_entities(airfoil::AirfoilParams, PhysicalGroups::DataFrame, io::IOS
         sheet = "NonSharp"
     end
     
-    periodicmap = get_map_periodic_surfaces(sharp_end)
+    periodicmap = get_map_periodic_surfaces(sharp_end,pfc)
     
-    points_physical_map = get_map_points(sharp_end) 
-    lines_physical_map = get_map_lines(sharp_end) 
-    surfaces_physical_map = get_map_surfaces(sharp_end) 
+    points_physical_map = get_map_points(sharp_end,pfc) 
+    lines_physical_map = get_map_lines(sharp_end,pfc) 
+    surfaces_physical_map = get_map_surfaces(sharp_end,pfc) 
     
     
     points_physical = Vector[]
@@ -80,7 +80,8 @@ end
 
 
 
-function get_map_points(sharp_end::Bool)
+function get_map_points(sharp_end::Bool,pfc::Bool)
+    if pfc==false
 if sharp_end
     v1 = 7, 13, 12,  26, 30, 24
     v2 = 2, 6, 4, 18, 21, 25, 1, 5, 3, 15, 19, 23
@@ -92,11 +93,21 @@ else
     v3 = 28, 29, 30, 31
 end
 
-    return DataFrame(Points =[v1, v2, v3], Physical=["outlet", "limits", "airfoil"])
+else
+    v1 = 22,    16,    21,    35,    39,    32    
+    v2 = 11,    15,    27,    28,    13,    34,    14,    30,    10,    24,    7,    8,    9,    44,    42,    40,    12,    31
+    
+    v3 = 36,    37,    46,    47,    48,    38
+    
+
+end
+
+return DataFrame(Points =[v1, v2, v3], Physical=["outlet", "limits", "airfoil"])
 end
 
 
-function get_map_lines(sharp_end::Bool)
+function get_map_lines(sharp_end::Bool,pfc::Bool)
+    if pfc==false
     if sharp_end
         v1 = 1, 2, 3, 58, 62, 64, 57, 55, 60
         v2 = 4, 34
@@ -112,12 +123,19 @@ function get_map_lines(sharp_end::Bool)
         v4= 16, 17, 18, 19, 20, 55, 72, 69, 50, 56, 73, 75, 70, 51
 
     end
+else
+    v1 =1,2,3,4,5,6,98,104,101,68,70,107,100,71,75,106,103,73
+    v2= 7,49
+    v3= 12,16,51,64,47,50,63,8,11,13,82,96,57,42,55,56,87,92,81,86,91,9,10
+    v4= 24,25,26,27,59,77,80,66,65,76,58
+end
     
         return DataFrame(Lines =[v1, v2, v3, v4], Physical=["airfoil", "inlet", "limits", "outlet"])
 end
 
 
-function get_map_surfaces(sharp_end::Bool)
+function get_map_surfaces(sharp_end::Bool,pfc::Bool)
+    if pfc==false
     if sharp_end
         v1 =14
         v2 =29, 45, 42, 25
@@ -135,12 +153,21 @@ function get_map_surfaces(sharp_end::Bool)
 
 
     end
-    
+else
+    v1 =20
+    v2= 27,42,46,32
+    v3= 31,22,26,60,56,52,48
+    v4= 72,64,40,36,67,70
+
+end
+
         return DataFrame(Surfaces =[v1, v2, v3, v4], Physical=["inlet", "outlet", "limits", "airfoil"])
 end
 
 
-function get_map_periodic_surfaces(sharp_end::Bool)
+function get_map_periodic_surfaces(sharp_end::Bool,pfc::Bool)
+    if pfc==false
+
 if sharp_end
     from = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
     to = 15,    19,    23,27,31,35,38,41,44,46
@@ -150,13 +177,15 @@ else
     to = 16,    20,    24,    28,    32,    36,    39,    42,    45,    48,    51
     
 end
+
+else
+    from = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12,13,14,15,16
+    to=21,25,30,34,38,41,45,47,51,55,59,62,65,68,71,73
+    
+end
+
 return DataFrame(From =from, To=to)
 end
 
-sa = get_map_periodic_surfaces(false)
-sa
-for i = 1:1:length(sa.To[1])
 
-    println(sa.To[1][i]) 
-end
 
diff --git a/src/WriteFileUtils.jl b/src/WriteFileUtils.jl
index 025857a..753ea02 100644
--- a/src/WriteFileUtils.jl
+++ b/src/WriteFileUtils.jl
@@ -47,7 +47,7 @@ function start_writing(Airfoil::AirfoilParams, dimension::Int64, chord::Float64,
 
     write(io, "AoA_deg = DefineNumber[ 0, Name \"Parameters/AoA\" ];\n")
     write(io, "AoA = AoA_deg*3.14159/180;\n")
-    write(io, "a_dim = 0.2;\n")
+    write(io, "a_dim = 2;\n")
 
 
         

From 3b3ec5ccd38b63b3c489bdede461b51652be0977 Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Mon, 16 Sep 2024 14:34:48 +0200
Subject: [PATCH 03/18] interface BL

---
 src/AirfoilGmsh.jl    |  8 +++++---
 src/BLanalysis.jl     | 46 ++++++++++++++++++++++++++++++++-----------
 src/CreateGeoFile.jl  |  7 ++++---
 src/WriteFileUtils.jl |  9 ++++++---
 4 files changed, 50 insertions(+), 20 deletions(-)

diff --git a/src/AirfoilGmsh.jl b/src/AirfoilGmsh.jl
index 61863a2..4bb305c 100644
--- a/src/AirfoilGmsh.jl
+++ b/src/AirfoilGmsh.jl
@@ -24,8 +24,6 @@ export get_airfoil_name
 include("AirfoilUtils.jl")
 
 
-export start_writing
-include("WriteFileUtils.jl")
 
 export addAirfoilPoints
 export addShearPoint
@@ -44,9 +42,13 @@ export RecombineSurfaces
 export addPhysicalGroup
 include("GmshUtils.jl")
 
-export refinement_parameters
+export BoundaryLayer
 include("BLanalysis.jl")
 
+
+export start_writing
+include("WriteFileUtils.jl")
+
 export map_entities
 include("MapLines.jl")
 
diff --git a/src/BLanalysis.jl b/src/BLanalysis.jl
index 0eec8d6..a9012f0 100644
--- a/src/BLanalysis.jl
+++ b/src/BLanalysis.jl
@@ -91,17 +91,41 @@ function boundary_layer_characteristics(Re::Real, H::Real, h0::Real, chord::Floa
 end
 
 
-function refinement_parameters(Reynolds::Real, h0::Real, chord::Real)
-    if Reynolds < 0 && h0 < 0 #If no reynolds or height specified
-        return 0.35, 100, 1.12, h0
-    else
-        H = 0.35 * chord
-        if h0 < 0
-            h0 = chord * sqrt(74) * Reynolds^(-13 / 14)
-            println("Extimated h0 = $h0 m")
-        end
-        H_levels, N_levels, G, h0 = boundary_layer_characteristics(Reynolds, H, h0, chord)
+mutable struct BoundaryLayer
+    H_levels::Real 
+    N_levels::Real
+    G::Real
+    h0::Real
+end
+
+function BoundaryLayer(H_levels::Float64,N_levels::Int64,G::Float64)
+    BoundaryLayer(H_levels,N_levels,G,0.0)
+end
+
+### Providing the Reynolds Number or Height.
+## If H0<0.0, then is not provided
+function BoundaryLayer(Reynolds::Real,h0::Float64, chord::Float64)
+    H = 0.35 * chord
+    
+    if h0<0.0 ### I am assuming it is a Reynolds number
+        h0 = chord * sqrt(74) * Reynolds^(-13 / 14)
+        println("Extimated h0 = $h0 m")
+    end
+
+    H_levels, N_levels, G, h0 = boundary_layer_characteristics(Reynolds, H, h0, chord)
+    BoundaryLayer(H_levels, N_levels, G)
+end
 
-        return H_levels, N_levels, G, h0
+function BoundaryLayer(H_levels,N_levels,G, Reynolds,h0, chord)
+if H_levels<0.0 && N_levels<0.0 && G<0.0
+    @assert chord>0.0
+    if Reynolds<0.0 && h0<0.0
+        return BoundaryLayer(0.35*chord, 85, 1.12, 0.0) ## Default values
     end
+    return BoundaryLayer(Reynolds,h0, chord)
+else
+    return BoundaryLayer(H_levels,N_levels,G)
 end
+
+end
+
diff --git a/src/CreateGeoFile.jl b/src/CreateGeoFile.jl
index 8dc25c8..ab0fa6c 100644
--- a/src/CreateGeoFile.jl
+++ b/src/CreateGeoFile.jl
@@ -21,10 +21,11 @@ It is possible to create a mesh with the following options:
 | Thetraedreal   | 3D        | :TETRA    |
 
 """
-function create_geofile(filename::String; Reynolds = -1, h0 = -1, leading_edge_points = Int64[], trailing_edge_points = Int64[], chord=1.0, dimension=2, elements = :QUAD, open_geo = false)
+function create_geofile(filename::String; Reynolds = -1, h0 = -1, leading_edge_points = Int64[], trailing_edge_points = Int64[], chord=1.0, dimension=2, elements = :QUAD,  H_levels=-1, N_levels=-1, G=-1)
 
+    open_geo = false
 
-refinement_params = refinement_parameters(Reynolds, h0, chord)
+boundary_layer = BoundaryLayer(H_levels,N_levels,G, Reynolds,h0, chord)
     
 if elements == :QUAD || elements == :HEX
     recombine = true
@@ -41,7 +42,7 @@ PhysicalGroups = DataFrame(number=Int64[], name=String[], entities=Vector[], typ
 
 Airfoil = AirfoilParams(filename, chord, trailing_edge_points, leading_edge_points)
 
-io = start_writing(Airfoil, dimension, chord, refinement_params)
+io = start_writing(Airfoil, dimension, chord, boundary_layer)
 
 addAirfoilPoints(Airfoil, Points, io)
 Airfoil.points.leading_edge
diff --git a/src/WriteFileUtils.jl b/src/WriteFileUtils.jl
index 753ea02..06720d9 100644
--- a/src/WriteFileUtils.jl
+++ b/src/WriteFileUtils.jl
@@ -17,12 +17,15 @@ const N_edge = 7 #minimum value, it will be overwritten
 
 
 """
-    start_writing(Airfoil::AirfoilParams, dimension::Int64, chord::Float64, refinement_params::Tuple)
+    start_writing(Airfoil::AirfoilParams, dimension::Int64, chord::Float64, boundary_layer::BoundaryLayer)
 
 It starts writing a new .geo file. It writes all the custom parameters that can be later modified when the file is opened in Gmsh.
 """
-function start_writing(Airfoil::AirfoilParams, dimension::Int64, chord::Float64, refinement_params::Tuple)
-    Refinement_offset, N_refinement, P_refinement, h0 = refinement_params
+function start_writing(Airfoil::AirfoilParams, dimension::Int64, chord::Float64, boundary_layer::BoundaryLayer)
+
+    Refinement_offset = boundary_layer.H_levels
+    N_refinement = boundary_layer.N_levels
+    P_refinement = boundary_layer.G
 
     io = open("$(Airfoil.name)_$(dimension)D.geo", "w")
     write(io, "SetFactory(\"OpenCASCADE\");\n")

From efe6cef25aa7c6ab9791b924dcc4fd4b0c54a866 Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Mon, 16 Sep 2024 16:08:21 +0200
Subject: [PATCH 04/18] re-parametrizing

---
 src/BLanalysis.jl     | 55 +++++++++-----------------
 src/CreateGeoFile.jl  | 14 ++++---
 src/WriteFileUtils.jl | 92 ++++++++++++++++++++++++++++---------------
 3 files changed, 88 insertions(+), 73 deletions(-)

diff --git a/src/BLanalysis.jl b/src/BLanalysis.jl
index a9012f0..1dba186 100644
--- a/src/BLanalysis.jl
+++ b/src/BLanalysis.jl
@@ -1,3 +1,13 @@
+@with_kw mutable struct BoundaryLayer
+    H_levels::Float64 = 0.35
+    N_levels::Float64 = 80
+    G::Real = 1.12
+    Reynolds::Int64=1e6
+    h0::Float64=1e-4
+end
+
+
+
 
 """
     yt(yh::Float64, G::Float64, N::Int)
@@ -91,41 +101,14 @@ function boundary_layer_characteristics(Re::Real, H::Real, h0::Real, chord::Floa
 end
 
 
-mutable struct BoundaryLayer
-    H_levels::Real 
-    N_levels::Real
-    G::Real
-    h0::Real
-end
-
-function BoundaryLayer(H_levels::Float64,N_levels::Int64,G::Float64)
-    BoundaryLayer(H_levels,N_levels,G,0.0)
-end
-
-### Providing the Reynolds Number or Height.
-## If H0<0.0, then is not provided
-function BoundaryLayer(Reynolds::Real,h0::Float64, chord::Float64)
-    H = 0.35 * chord
-    
-    if h0<0.0 ### I am assuming it is a Reynolds number
-        h0 = chord * sqrt(74) * Reynolds^(-13 / 14)
-        println("Extimated h0 = $h0 m")
-    end
-
-    H_levels, N_levels, G, h0 = boundary_layer_characteristics(Reynolds, H, h0, chord)
-    BoundaryLayer(H_levels, N_levels, G)
-end
-
-function BoundaryLayer(H_levels,N_levels,G, Reynolds,h0, chord)
-if H_levels<0.0 && N_levels<0.0 && G<0.0
-    @assert chord>0.0
-    if Reynolds<0.0 && h0<0.0
-        return BoundaryLayer(0.35*chord, 85, 1.12, 0.0) ## Default values
-    end
-    return BoundaryLayer(Reynolds,h0, chord)
-else
-    return BoundaryLayer(H_levels,N_levels,G)
-end
-
+function BoundaryLayer(Reynolds::Int64,ds::DomainSize )
+    @unpack chord = ds
+    h0 = chord * sqrt(74) * Reynolds^(-13 / 14)
+    BoundaryLayer(h0, ds)
 end
 
+function BoundaryLayer(h0::Float64, ds::DomainSize)
+    @unpack chord, H_offset = ds
+    H_levels, N_levels, G, h0 = boundary_layer_characteristics(Reynolds, H_offset, h0, chord)
+    BoundaryLayer(H_levels, N_levels, G,Reynolds, h0)
+end
\ No newline at end of file
diff --git a/src/CreateGeoFile.jl b/src/CreateGeoFile.jl
index ab0fa6c..9ea5a37 100644
--- a/src/CreateGeoFile.jl
+++ b/src/CreateGeoFile.jl
@@ -23,7 +23,7 @@ It is possible to create a mesh with the following options:
 """
 function create_geofile(filename::String; Reynolds = -1, h0 = -1, leading_edge_points = Int64[], trailing_edge_points = Int64[], chord=1.0, dimension=2, elements = :QUAD,  H_levels=-1, N_levels=-1, G=-1)
 
-    open_geo = false
+open_geo = false
 
 boundary_layer = BoundaryLayer(H_levels,N_levels,G, Reynolds,h0, chord)
     
@@ -264,14 +264,18 @@ point7 = addPoint("L", "-L* " * string(x_tmp) * "*Sin(AoA) + " * string(y_tmp) *
     
     
     
-    airfoil_lines = [LinefromPoints(point1, point3, Lines),
+    airfoil_lines_top = [LinefromPoints(point1, point3, Lines),
         LinefromPoints(point1r, point3r, Lines),
-        LinefromPoints(Airfoil.points.trailing_edge[1], Airfoil.points.leading_edge[1], Lines),
+        LinefromPoints(Airfoil.points.trailing_edge[1], Airfoil.points.leading_edge[1], Lines)]
+
+    airfoil_lines_bottom = [
         LinefromPoints(Airfoil.points.trailing_edge[Airfoil.sharp_idx], Airfoil.points.leading_edge[2], Lines),
         LinefromPoints(point2, point4, Lines),
         LinefromPoints(point2r, point4r, Lines)]
-    TransfiniteCurve(airfoil_lines, "N_airfoil", 1.0, io)
-    
+
+    TransfiniteCurve(airfoil_lines_top, "N_airfoil_top", 1.0, io)
+    TransfiniteCurve(airfoil_lines_bottom, "N_airfoil_bottom", 1.0, io)
+
     if is_sharp(Airfoil)
         shear_lines = [LinefromPoints(point3, point5, Lines),
             LinefromPoints(point3r, point7r, Lines),
diff --git a/src/WriteFileUtils.jl b/src/WriteFileUtils.jl
index 06720d9..ed0fe91 100644
--- a/src/WriteFileUtils.jl
+++ b/src/WriteFileUtils.jl
@@ -1,18 +1,46 @@
-#Default Values
-const N_inlet = 30
-const N_vertical = 30
-const P_vertical = 1.1
+using Parameters
+mutable struct DomainDivision
+    numdiv::Int64
+    progression::Float64
+end
+
+@with_kw mutable struct DomainMeshDivisions
+    Inlet::DomainDivision
+    Vertical::DomainDivision
+    Shear::DomainDivision
+    Airfoil_Top::DomainDivision
+    Airfoil_Bottom::DomainDivision
+    Edge::DomainDivision
+    Nz::DomainDivision
+end
+
+@with_kw mutable struct DomainSize
+    AoA::Float64 = 0.0
+    chord::Float64 = 1.0
+    C::Float64 = 6.0*chord
+    L::Float64 = 6.0*chord
+    H_offset::Float64=0.35*chord
+    Hz::Float64=0.2
+    dimension::Int64=3
+end
+
+### DEFAULT VALUES FOR DIVISION AND SIZE
+
+domain_size_default = DomainSize(chord=3.0)
+
+inlet = DomainDivision(30,1.0)
+vertical = DomainDivision(30,1.1)
+shear = DomainDivision(30,1.2)
+airfoil_top = DomainDivision(200,1.0)
+airfoil_bottom = DomainDivision(50,1.0)
+edge = DomainDivision(7,1.0)
+nz = DomainDivision(16,1.0)
+
+domain_mesh_divisions_default = DomainMeshDivisions(inlet, vertical, shear, airfoil_top, airfoil_bottom, edge, nz)
+
 
-const N_airfoil = 50
-const N_shear = 30
-const P_shear = 1.2
-const L_domain = 6
-const C_domain = 6
 
-const Hz_coeff = 0.2
-const Nz_default = 22
 
-const N_edge = 7 #minimum value, it will be overwritten
 
 
 
@@ -21,41 +49,41 @@ const N_edge = 7 #minimum value, it will be overwritten
 
 It starts writing a new .geo file. It writes all the custom parameters that can be later modified when the file is opened in Gmsh.
 """
-function start_writing(Airfoil::AirfoilParams, dimension::Int64, chord::Float64, boundary_layer::BoundaryLayer)
+function start_writing(Airfoil::AirfoilParams, domain_size::DomainSize, domain_mesh_divisions::DomainMeshDivisions, boundary_layer::BoundaryLayer)
+    @unpack H_levels,N_levels,G = boundary_layer
+    @unpack AoA, chord, dimension,C,L, Hz  = domain_size
 
-    Refinement_offset = boundary_layer.H_levels
-    N_refinement = boundary_layer.N_levels
-    P_refinement = boundary_layer.G
 
     io = open("$(Airfoil.name)_$(dimension)D.geo", "w")
     write(io, "SetFactory(\"OpenCASCADE\");\n")
     
-    write(io, "N_inlet = DefineNumber[ $(N_inlet), Name \"Parameters/N_inlet\" ];\n")
-    write(io, "N_vertical = DefineNumber[ $(N_vertical), Name \"Parameters/N_vertical\" ];\n")
-    write(io, "P_vertical = DefineNumber[ $(P_vertical), Name \"Parameters/P_vertical\" ];\n")
-    
-    write(io, "N_airfoil = DefineNumber[ $(N_airfoil), Name \"Parameters/N_airfoil\" ];\n")
+    write(io, "N_inlet = DefineNumber[ $(domain_mesh_divisions.Inlet.numdiv), Name \"Parameters/N_inlet\" ];\n")
+    write(io, "N_vertical = DefineNumber[ $(domain_mesh_divisions.Vertical.numdiv), Name \"Parameters/N_vertical\" ];\n")
+    write(io, "P_vertical = DefineNumber[ $(domain_mesh_divisions.Vertical.progression), Name \"Parameters/P_vertical\" ];\n")
     
-    write(io, "N_shear = DefineNumber[ $(N_shear), Name \"Parameters/N_shear\" ];\n")
-    write(io, "P_shear = DefineNumber[ $(P_shear), Name \"Parameters/P_shear\" ];\n")
-    write(io, "L = DefineNumber[ $(L_domain), Name \"Parameters/L\" ];\n")
-    write(io, "C = DefineNumber[ $(C_domain), Name \"Parameters/C\" ];\n")
+    write(io, "N_airfoil = DefineNumber[ $(domain_mesh_divisions.Airfoil_top.numdiv), Name \"Parameters/N_airfoil_top\" ];\n")
+    write(io, "N_airfoil = DefineNumber[ $(domain_mesh_divisions.Airfoil_bottom.numdiv), Name \"Parameters/N_airfoil_bottom\" ];\n")
+
+    write(io, "N_shear = DefineNumber[ $(domain_mesh_divisions.Shear.numdiv), Name \"Parameters/N_shear\" ];\n")
+    write(io, "P_shear = DefineNumber[ $(domain_mesh_divisions.Shear.progression), Name \"Parameters/P_shear\" ];\n")
+    write(io, "L = DefineNumber[ $(domain_mesh_divisions.L.numdiv), Name \"Parameters/L\" ];\n")
+    write(io, "C = DefineNumber[ $(domain_mesh_divisions.C.numdiv), Name \"Parameters/C\" ];\n")
     
-    write(io, "Hz = DefineNumber[ $(chord*Hz_coeff), Name \"Parameters/Hz\" ];\n")
-    write(io, "Nz = DefineNumber[ $(Nz_default), Name \"Parameters/Nz\" ];\n")
+    write(io, "Hz = DefineNumber[ $(chord*Hz), Name \"Parameters/Hz\" ];\n")
+    write(io, "Nz = DefineNumber[ $(domain_mesh_divisions.Nz.numdiv), Name \"Parameters/Nz\" ];\n")
     
-    write(io, "Refinement_offset = DefineNumber[ $Refinement_offset, Name \"Parameters/Refinement_offset\" ];\n")
-    write(io, "N_refinement = DefineNumber[ $N_refinement, Name \"Parameters/N_refinement\" ];\n")
-    write(io, "P_refinement = DefineNumber[ $P_refinement, Name \"Parameters/P_refinement\" ];\n")
+    write(io, "Refinement_offset = DefineNumber[ $H_levels, Name \"Parameters/Refinement_offset\" ];\n")
+    write(io, "N_refinement = DefineNumber[ $N_levels, Name \"Parameters/N_refinement\" ];\n")
+    write(io, "P_refinement = DefineNumber[ $G, Name \"Parameters/P_refinement\" ];\n")
 
-    write(io, "AoA_deg = DefineNumber[ 0, Name \"Parameters/AoA\" ];\n")
+    write(io, "AoA_deg = DefineNumber[ $(AoA), Name \"Parameters/AoA\" ];\n")
     write(io, "AoA = AoA_deg*3.14159/180;\n")
     write(io, "a_dim = 2;\n")
 
 
         
     if !is_sharp(Airfoil) 
-        write(io, "N_edge = DefineNumber[ $N_edge, Name \"Parameters/N_edge\" ];\n")
+        write(io, "N_edge = DefineNumber[ $(domain_mesh_divisions.Edge.numdiv), Name \"Parameters/N_edge\" ];\n")
     end
 
 return io

From 8edd01b62b00225e080a7b15f08588acab1dfb57 Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Tue, 17 Sep 2024 10:24:49 +0200
Subject: [PATCH 05/18] add parametrization

---
 src/AirfoilGmsh.jl    | 18 ++++++++-
 src/BLanalysis.jl     | 18 ++-------
 src/CreateGeoFile.jl  | 90 +++++++++++++++++++++++++++----------------
 src/DefaultValues.jl  | 46 ++++++++++++++++++++++
 src/WriteFileUtils.jl | 60 ++++-------------------------
 test/test_driver.jl   | 10 ++---
 6 files changed, 134 insertions(+), 108 deletions(-)
 create mode 100644 src/DefaultValues.jl

diff --git a/src/AirfoilGmsh.jl b/src/AirfoilGmsh.jl
index 4bb305c..9708e60 100644
--- a/src/AirfoilGmsh.jl
+++ b/src/AirfoilGmsh.jl
@@ -12,9 +12,23 @@ using Optimization, OptimizationBBO
 using Parameters
 using CubicSplines
 
+
+export DomainDivision
+export DomainMeshDivisions
+export DomainInfo
+export BoundaryLayer
+export ElementShape
+export QUAD
+export TRI
+export HEX
+export TETRA
+include("DefaultValues.jl")
+
+
 export from_url_to_csv
 include("ReadWeb.jl")
 
+
 export AirfoilPoints
 export AirfoilParams
 export get_airfoil_features
@@ -42,10 +56,10 @@ export RecombineSurfaces
 export addPhysicalGroup
 include("GmshUtils.jl")
 
-export BoundaryLayer
-include("BLanalysis.jl")
 
 
+include("BLanalysis.jl")
+
 export start_writing
 include("WriteFileUtils.jl")
 
diff --git a/src/BLanalysis.jl b/src/BLanalysis.jl
index 1dba186..7bba485 100644
--- a/src/BLanalysis.jl
+++ b/src/BLanalysis.jl
@@ -1,14 +1,3 @@
-@with_kw mutable struct BoundaryLayer
-    H_levels::Float64 = 0.35
-    N_levels::Float64 = 80
-    G::Real = 1.12
-    Reynolds::Int64=1e6
-    h0::Float64=1e-4
-end
-
-
-
-
 """
     yt(yh::Float64, G::Float64, N::Int)
 
@@ -101,14 +90,15 @@ function boundary_layer_characteristics(Re::Real, H::Real, h0::Real, chord::Floa
 end
 
 
-function BoundaryLayer(Reynolds::Int64,ds::DomainSize )
+function BoundaryLayer(Reynolds::Int64,ds::DomainInfo )
     @unpack chord = ds
     h0 = chord * sqrt(74) * Reynolds^(-13 / 14)
     BoundaryLayer(h0, ds)
 end
 
-function BoundaryLayer(h0::Float64, ds::DomainSize)
+function BoundaryLayer(h0::Float64, ds::DomainInfo)
     @unpack chord, H_offset = ds
     H_levels, N_levels, G, h0 = boundary_layer_characteristics(Reynolds, H_offset, h0, chord)
     BoundaryLayer(H_levels, N_levels, G,Reynolds, h0)
-end
\ No newline at end of file
+end
+
diff --git a/src/CreateGeoFile.jl b/src/CreateGeoFile.jl
index 9ea5a37..2354426 100644
--- a/src/CreateGeoFile.jl
+++ b/src/CreateGeoFile.jl
@@ -1,5 +1,6 @@
 """
-    create_geofile(filename::String; Reynolds = -1, h0 = -1, leading_edge_points = Int64[], trailing_edge_points = Int64[], chord=1.0, dimension=2, elements = :QUAD, open_geo = true)
+  create_geofile(filename::String, domain_mesh_divisions::DomainMeshDivisions, domain_info::DomainInfo, boundary_layer::BoundaryLayer)
+
 
 It is the main function of the package. From a csv file containing the airfoil points it creates a .geo file.
 The .geo file can be created using the function [`from_url_to_csv`](@ref).
@@ -7,27 +8,34 @@ The user can specify just the file name.
 ```julia
     create_geofile("naca0012.csv")
 ```
-It is also possibile to provide extra arguments such as the Reynolds number and/or the first layer height for a better extimation of the boundary-cell properties.
-It is possible to overwrite the extimation of the trailing edge and leading edge made by the code providing the relative points numbers.
+It is possible to provide just the string of the filename.csv where the points are stored.
+It is also possible to customize the Mesh Division, Domain Size and provide information on the boundary layer.
+See the default parameters in the "DefaultValues.jl" file.
+
+
 The mesh can be created in 2D or 3D. In 3D case by default are created periodic boundary conditions in the `z` direction.
+
 It is possible to create a mesh with the following options:
 
 
-| Type of element | Dimension | Symbol   |
+| Type of element | Dimension | struct   |
 | ---------------|-----------|-----------|
-| Quadrilateral  | 2D        | :QUAD     |
-| Hexaedral      | 3D        | :HEX      |
-| Triangular     | 2D        | :TRI      |
-| Thetraedreal   | 3D        | :TETRA    |
+| Quadrilateral  | 2D        | QUAD()     |
+| Hexaedral      | 3D        | HEX()      |
+| Triangular     | 2D        | TRI()      |
+| Thetraedreal   | 3D        | TETRA()    |
 
 """
-function create_geofile(filename::String; Reynolds = -1, h0 = -1, leading_edge_points = Int64[], trailing_edge_points = Int64[], chord=1.0, dimension=2, elements = :QUAD,  H_levels=-1, N_levels=-1, G=-1)
 
-open_geo = false
+function create_geofile(filename::String, domain_mesh_divisions::DomainMeshDivisions, domain_info::DomainInfo, boundary_layer::BoundaryLayer)
 
-boundary_layer = BoundaryLayer(H_levels,N_levels,G, Reynolds,h0, chord)
+
+leading_edge_points = Int64[]
+trailing_edge_points = Int64[]
+
+@unpack chord,dimension = domain_info
     
-if elements == :QUAD || elements == :HEX
+if domain_info.elements == QUAD() || domain_info.elements == HEX()
     recombine = true
     else
     recombine = false
@@ -40,35 +48,30 @@ Loops = Vector[]
 PhysicalGroups = DataFrame(number=Int64[], name=String[], entities=Vector[], type=String[])
 
 
-Airfoil = AirfoilParams(filename, chord, trailing_edge_points, leading_edge_points)
+Airfoil = AirfoilParams( filename, chord, trailing_edge_points, leading_edge_points)
 
-io = start_writing(Airfoil, dimension, chord, boundary_layer)
+io = start_writing(Airfoil, domain_info, domain_mesh_divisions, boundary_layer)
 
 addAirfoilPoints(Airfoil, Points, io)
 Airfoil.points.leading_edge
 Airfoil.points.trailing_edge[Airfoil.sharp_idx]
 
-N_edge = 5
+N_edge = domain_mesh_divisions.Edge.numdiv
 # Create airfoil lines
-    
 spline_airfoil_top = addSpline(collect(Airfoil.points.trailing_edge[1] : Airfoil.points.leading_edge[1]), Lines, io;all=true)[end][1]
 
 spline_airfoil_le = addSpline(collect(Airfoil.points.leading_edge[1] : Airfoil.points.leading_edge[2]), Lines, io; all=true)[end][1]
 
-Airfoil.points.leading_edge[1] 
-Airfoil.points.leading_edge[2]
-Lines
-
 if ! is_sharp(Airfoil)
     spline_airfoil_te = addLine(Airfoil.points.trailing_edge[1], Airfoil.points.trailing_edge[2], Lines, io)[end][1]
-    spline_airfoil_bottom = addSpline(Airfoil.points.leading_edge[2] : Airfoil.points.trailing_edge[Airfoil.sharp_idx] , Lines, io)[end][1]
+    spline_airfoil_bottom = addSpline(collect(Airfoil.points.leading_edge[2] : Airfoil.points.trailing_edge[Airfoil.sharp_idx] ), Lines, io; all=true)[end][1]
 else
     b_points = [Airfoil.points.leading_edge[2]: Airfoil.points.num, Airfoil.points.trailing_edge[Airfoil.sharp_idx]]
+    
     spline_airfoil_bottom = addSpline(b_points, Lines, io)[end][1]
 end
 
 
-
 #External Domain points
 point1 = addPoint(0, "C", 0, Points, io; tag ="external")[end][1]
 point2 = addPoint(0, "-C", 0, Points, io; tag ="external")[end][1]
@@ -171,7 +174,8 @@ point7 = addPoint("L", "-L* " * string(x_tmp) * "*Sin(AoA) + " * string(y_tmp) *
     loop3 = LoopfromPoints([point2, point4, point4r, point2r], Lines)
     loop3r = LoopfromPoints([point2r, point4r, Airfoil.points.trailing_edge[Airfoil.sharp_idx], Airfoil.points.leading_edge[2]], Lines)
     
-    LinefromPoints(point8r, point6, Lines)
+
+    LinefromPoints(Airfoil.points.trailing_edge[Airfoil.sharp_idx], Airfoil.points.leading_edge[2], Lines)
     
     loop4 = LoopfromPoints([point3, point5, point7r, point3r], Lines)
     loop4r = LoopfromPoints([point3r, point7r, point7, Airfoil.points.trailing_edge[1]], Lines)
@@ -372,20 +376,38 @@ point7 = addPoint("L", "-L* " * string(x_tmp) * "*Sin(AoA) + " * string(y_tmp) *
     map_entities(Airfoil, PhysicalGroups, io)
     
     end
-    
-    
-    
-    
-    
-    close(io)
-    
-    ## Open GMSH for visualization of the .geo file and changing the parameters
-    if open_geo
-        open_gmsh(Airfoil, dimension)
-    end
+        
+   close(io)
+
 
     return io    
 
 end
 
 
+### DEFAULT
+function create_geofile(filename::String )
+    domain_mesh_divisions=DomainMeshDivisions()
+    domain_info=DomainInfo()
+    boundary_layer=BoundaryLayer()
+    return create_geofile(filename, domain_mesh_divisions, domain_info, boundary_layer)
+end
+
+function create_geofile(filename::String, domain_mesh_divisions::DomainMeshDivisions)
+    domain_info=DomainInfo()
+    boundary_layer=BoundaryLayer()
+    return create_geofile(filename, domain_mesh_divisions, domain_info, boundary_layer)
+end
+
+function create_geofile(filename::String, domain_info::DomainInfo)
+    domain_mesh_divisions=DomainMeshDivisions()
+    boundary_layer=BoundaryLayer()
+    return create_geofile(filename, domain_mesh_divisions, domain_info, boundary_layer)
+end
+
+
+function create_geofile(filename::String,      boundary_layer::BoundaryLayer)
+    domain_mesh_divisions=DomainMeshDivisions()
+    domain_info=DomainInfo()
+    return create_geofile(filename, domain_mesh_divisions, domain_info, boundary_layer)
+end
\ No newline at end of file
diff --git a/src/DefaultValues.jl b/src/DefaultValues.jl
new file mode 100644
index 0000000..840833a
--- /dev/null
+++ b/src/DefaultValues.jl
@@ -0,0 +1,46 @@
+abstract type ElementShape end
+
+# Define concrete types for specific shapes
+struct QUAD <: ElementShape end
+struct TRI <: ElementShape end
+struct HEX <: ElementShape end
+struct TETRA <: ElementShape end
+
+
+mutable struct DomainDivision
+    numdiv::Int64
+    progression::Float64
+end
+
+### Domain Mesh Divisions Default Parameters
+@with_kw mutable struct DomainMeshDivisions
+    Inlet::DomainDivision= DomainDivision(30,1.0)
+    Vertical::DomainDivision= DomainDivision(30,1.1)
+    Shear::DomainDivision= DomainDivision(30,1.2)
+    Airfoil_Top::DomainDivision= DomainDivision(200,1.0)
+    Airfoil_Bottom::DomainDivision = DomainDivision(50,1.0)
+    Edge::DomainDivision= DomainDivision(7,1.0)
+    Nz::DomainDivision= DomainDivision(16,1.0)
+end
+
+### Domain Size Default Parameters
+@with_kw mutable struct DomainInfo
+    AoA::Float64 = 0.0
+    chord::Float64 = 1.0
+    C::Float64 = 6.0*chord
+    L::Float64 = 6.0*chord
+    H_offset::Float64=0.35*chord
+    Hz::Float64=0.2
+    dimension::Int64=3
+    elements::ElementShape=QUAD()
+end
+
+### Boundary Layer Default Parameters
+@with_kw mutable struct BoundaryLayer
+    H_levels::Float64 = 0.35
+    N_levels::Float64 = 80
+    G::Real = 1.12
+    Reynolds::Int64=1e6
+    h0::Float64=1e-4
+end
+
diff --git a/src/WriteFileUtils.jl b/src/WriteFileUtils.jl
index ed0fe91..a8a0ca1 100644
--- a/src/WriteFileUtils.jl
+++ b/src/WriteFileUtils.jl
@@ -1,57 +1,11 @@
-using Parameters
-mutable struct DomainDivision
-    numdiv::Int64
-    progression::Float64
-end
-
-@with_kw mutable struct DomainMeshDivisions
-    Inlet::DomainDivision
-    Vertical::DomainDivision
-    Shear::DomainDivision
-    Airfoil_Top::DomainDivision
-    Airfoil_Bottom::DomainDivision
-    Edge::DomainDivision
-    Nz::DomainDivision
-end
-
-@with_kw mutable struct DomainSize
-    AoA::Float64 = 0.0
-    chord::Float64 = 1.0
-    C::Float64 = 6.0*chord
-    L::Float64 = 6.0*chord
-    H_offset::Float64=0.35*chord
-    Hz::Float64=0.2
-    dimension::Int64=3
-end
-
-### DEFAULT VALUES FOR DIVISION AND SIZE
-
-domain_size_default = DomainSize(chord=3.0)
-
-inlet = DomainDivision(30,1.0)
-vertical = DomainDivision(30,1.1)
-shear = DomainDivision(30,1.2)
-airfoil_top = DomainDivision(200,1.0)
-airfoil_bottom = DomainDivision(50,1.0)
-edge = DomainDivision(7,1.0)
-nz = DomainDivision(16,1.0)
-
-domain_mesh_divisions_default = DomainMeshDivisions(inlet, vertical, shear, airfoil_top, airfoil_bottom, edge, nz)
-
-
-
-
-
-
-
 """
-    start_writing(Airfoil::AirfoilParams, dimension::Int64, chord::Float64, boundary_layer::BoundaryLayer)
+    start_writing(Airfoil::AirfoilParams, domain_size::DomainInfo, domain_mesh_divisions::DomainMeshDivisions, boundary_layer::BoundaryLayer)
 
 It starts writing a new .geo file. It writes all the custom parameters that can be later modified when the file is opened in Gmsh.
 """
-function start_writing(Airfoil::AirfoilParams, domain_size::DomainSize, domain_mesh_divisions::DomainMeshDivisions, boundary_layer::BoundaryLayer)
+function start_writing(Airfoil::AirfoilParams, domain_info::DomainInfo, domain_mesh_divisions::DomainMeshDivisions, boundary_layer::BoundaryLayer)
     @unpack H_levels,N_levels,G = boundary_layer
-    @unpack AoA, chord, dimension,C,L, Hz  = domain_size
+    @unpack AoA, chord, dimension,C,L, Hz  = domain_info
 
 
     io = open("$(Airfoil.name)_$(dimension)D.geo", "w")
@@ -61,13 +15,13 @@ function start_writing(Airfoil::AirfoilParams, domain_size::DomainSize, domain_m
     write(io, "N_vertical = DefineNumber[ $(domain_mesh_divisions.Vertical.numdiv), Name \"Parameters/N_vertical\" ];\n")
     write(io, "P_vertical = DefineNumber[ $(domain_mesh_divisions.Vertical.progression), Name \"Parameters/P_vertical\" ];\n")
     
-    write(io, "N_airfoil = DefineNumber[ $(domain_mesh_divisions.Airfoil_top.numdiv), Name \"Parameters/N_airfoil_top\" ];\n")
-    write(io, "N_airfoil = DefineNumber[ $(domain_mesh_divisions.Airfoil_bottom.numdiv), Name \"Parameters/N_airfoil_bottom\" ];\n")
+    write(io, "N_airfoil_top = DefineNumber[ $(domain_mesh_divisions.Airfoil_Top.numdiv), Name \"Parameters/N_airfoil_top\" ];\n")
+    write(io, "N_airfoil_bottom = DefineNumber[ $(domain_mesh_divisions.Airfoil_Bottom.numdiv), Name \"Parameters/N_airfoil_bottom\" ];\n")
 
     write(io, "N_shear = DefineNumber[ $(domain_mesh_divisions.Shear.numdiv), Name \"Parameters/N_shear\" ];\n")
     write(io, "P_shear = DefineNumber[ $(domain_mesh_divisions.Shear.progression), Name \"Parameters/P_shear\" ];\n")
-    write(io, "L = DefineNumber[ $(domain_mesh_divisions.L.numdiv), Name \"Parameters/L\" ];\n")
-    write(io, "C = DefineNumber[ $(domain_mesh_divisions.C.numdiv), Name \"Parameters/C\" ];\n")
+    write(io, "L = DefineNumber[ $(L), Name \"Parameters/L\" ];\n")
+    write(io, "C = DefineNumber[ $(C), Name \"Parameters/C\" ];\n")
     
     write(io, "Hz = DefineNumber[ $(chord*Hz), Name \"Parameters/Hz\" ];\n")
     write(io, "Nz = DefineNumber[ $(domain_mesh_divisions.Nz.numdiv), Name \"Parameters/Nz\" ];\n")
diff --git a/test/test_driver.jl b/test/test_driver.jl
index 2f4df0a..62c56ee 100644
--- a/test/test_driver.jl
+++ b/test/test_driver.jl
@@ -3,16 +3,16 @@ function test_airfoil(url::String)
     csvname = fname*".csv"
     @test from_url_to_csv(url) == csvname
     @test typeof(create_geofile(csvname)) == IOStream
-    @test typeof(create_geofile(csvname; Reynolds=200e3)) == IOStream
-    @test typeof(create_geofile(csvname; dimension = 3 )) == IOStream
-    x,y,wl,wu = increase_resolution_airfoil(csvname,500; maxiters = 10, maxtime = 10)
-    @test typeof(x) <: Vector{Float64}
+
+    @test typeof(create_geofile(csvname, BoundaryLayer(Reynolds=200e3))) == IOStream
+    @test typeof(create_geofile(csvname, DomainInfo(dimension = 2) )) == IOStream
+    # x,y,wl,wu = increase_resolution_airfoil(csvname,500; maxiters = 10, maxtime = 10)
+    # @test typeof(x) <: Vector{Float64}
     # rm(csvname)
     # rm(fname*"_2D.geo")
     # rm(fname*"_3D.geo")
 end
 
-
 function get_airfoil_name_test(url::String)
     s = findlast("/", url)[1]
     return url[s+1:end-4]

From baaffe4d216baa64d43d972925efb627ef8bfddd Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Tue, 17 Sep 2024 11:14:49 +0200
Subject: [PATCH 06/18] rm Hoffset

---
 src/DefaultValues.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/DefaultValues.jl b/src/DefaultValues.jl
index 840833a..94b7c1d 100644
--- a/src/DefaultValues.jl
+++ b/src/DefaultValues.jl
@@ -29,7 +29,6 @@ end
     chord::Float64 = 1.0
     C::Float64 = 6.0*chord
     L::Float64 = 6.0*chord
-    H_offset::Float64=0.35*chord
     Hz::Float64=0.2
     dimension::Int64=3
     elements::ElementShape=QUAD()

From caae7fe55340995f4c6a44dafea32294edfcbe5a Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Tue, 17 Sep 2024 15:18:57 +0200
Subject: [PATCH 07/18] new API for BLanalysis

---
 src/BLanalysis.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/BLanalysis.jl b/src/BLanalysis.jl
index 7bba485..0c795c0 100644
--- a/src/BLanalysis.jl
+++ b/src/BLanalysis.jl
@@ -99,6 +99,9 @@ end
 function BoundaryLayer(h0::Float64, ds::DomainInfo)
     @unpack chord, H_offset = ds
     H_levels, N_levels, G, h0 = boundary_layer_characteristics(Reynolds, H_offset, h0, chord)
-    BoundaryLayer(H_levels, N_levels, G,Reynolds, h0)
+    BoundaryLayer(H_levels=H_levels, N_levels=N_levels, G=G, Reynolds=Reynolds, h0=h0)
 end
 
+function BoundaryLayer(H_levels::Float64,N_levels::Int64,G::Float64)
+    BoundaryLayer(H_levels=H_levels, N_levels=N_levels, G=G)
+end
\ No newline at end of file

From 4508f7979bd53a5baa494fbcad400990adba91f0 Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Wed, 18 Sep 2024 09:28:32 +0200
Subject: [PATCH 08/18] export geo_filename

---
 src/CreateGeoFile.jl  | 4 ++--
 src/WriteFileUtils.jl | 6 +++---
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/CreateGeoFile.jl b/src/CreateGeoFile.jl
index 2354426..3c0d589 100644
--- a/src/CreateGeoFile.jl
+++ b/src/CreateGeoFile.jl
@@ -50,7 +50,7 @@ PhysicalGroups = DataFrame(number=Int64[], name=String[], entities=Vector[], typ
 
 Airfoil = AirfoilParams( filename, chord, trailing_edge_points, leading_edge_points)
 
-io = start_writing(Airfoil, domain_info, domain_mesh_divisions, boundary_layer)
+io, geo_filename = start_writing(Airfoil, domain_info, domain_mesh_divisions, boundary_layer)
 
 addAirfoilPoints(Airfoil, Points, io)
 Airfoil.points.leading_edge
@@ -380,7 +380,7 @@ point7 = addPoint("L", "-L* " * string(x_tmp) * "*Sin(AoA) + " * string(y_tmp) *
    close(io)
 
 
-    return io    
+    return io,geo_filename
 
 end
 
diff --git a/src/WriteFileUtils.jl b/src/WriteFileUtils.jl
index a8a0ca1..b703cdc 100644
--- a/src/WriteFileUtils.jl
+++ b/src/WriteFileUtils.jl
@@ -7,8 +7,8 @@ function start_writing(Airfoil::AirfoilParams, domain_info::DomainInfo, domain_m
     @unpack H_levels,N_levels,G = boundary_layer
     @unpack AoA, chord, dimension,C,L, Hz  = domain_info
 
-
-    io = open("$(Airfoil.name)_$(dimension)D.geo", "w")
+    geo_filename = "$(Airfoil.name)_$(dimension)D.geo"
+    io = open(geo_filename, "w")
     write(io, "SetFactory(\"OpenCASCADE\");\n")
     
     write(io, "N_inlet = DefineNumber[ $(domain_mesh_divisions.Inlet.numdiv), Name \"Parameters/N_inlet\" ];\n")
@@ -40,7 +40,7 @@ function start_writing(Airfoil::AirfoilParams, domain_info::DomainInfo, domain_m
         write(io, "N_edge = DefineNumber[ $(domain_mesh_divisions.Edge.numdiv), Name \"Parameters/N_edge\" ];\n")
     end
 
-return io
+return io,geo_filename
 end
 
 

From 52a2a995d2abd8120b2747f1f7b1e895e2062b17 Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Wed, 18 Sep 2024 10:46:16 +0200
Subject: [PATCH 09/18] convert to sharp

---
 src/AirfoilUtils.jl | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/src/AirfoilUtils.jl b/src/AirfoilUtils.jl
index b6de379..637f728 100644
--- a/src/AirfoilUtils.jl
+++ b/src/AirfoilUtils.jl
@@ -78,7 +78,8 @@ function get_airfoil_features(filename::String, c::Float64, trailing_edge_points
     airfoil_points_list = CSV.File(filename, header=true) |> Tables.matrix
     formatting_airfoil_points!(airfoil_points_list,c)
     
-    trailing_edge_points = findTE(trailing_edge_points, c, airfoil_points_list)
+    trailing_edge_points,airfoil_points_list = findTE(trailing_edge_points, c, airfoil_points_list)
+
     sharp_end, sharp_idx = detect_end(trailing_edge_points)
     
     leading_edge_points = findLE(leading_edge_points, c, airfoil_points_list)
@@ -131,6 +132,8 @@ function verify_trailing_edge(airfoil_points_list::Matrix{Float64}, c::Float64)
 return airfoil_points_list
 end
 
+
+
 """
     findTE(trailing_edge_points, c::Float64, airfoil_points_list::Matrix{Float64})
 
@@ -138,12 +141,17 @@ Automatically detects the trailing edge points indexes
 """
 function findTE(trailing_edge_points, c::Float64, airfoil_points_list::Matrix{Float64})
     atol = 1e-6
+    if norm(airfoil_points_list[1,:] - airfoil_points_list[end,:])<c .* 0.5e-2
+        airfoil_points_list = airfoil_points_list[1:end-1,:]
+    end
+
     while isempty(trailing_edge_points)
         trailing_edge_points = findall(x -> isapprox(x, c; atol=atol), airfoil_points_list[:, 1])
         atol = atol * 2
     end
     sort!(trailing_edge_points)
-    trailing_edge_points
+    
+    return trailing_edge_points,airfoil_points_list
 end
 
 """

From 8eec0926b013621b9d97a23aa526421c7cac54c4 Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Wed, 18 Sep 2024 11:05:20 +0200
Subject: [PATCH 10/18] add LinearAlgebra

---
 Project.toml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 2571368..f2e1c9d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "AirfoilGmsh"
 uuid = "e68ec710-8e48-43c5-af1a-c38c4248fd7b"
 authors = ["Carlo Brunelli"]
-version = "0.1.2"
+version = "0.1.3"
 
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
@@ -12,6 +12,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationBBO = "3e6eede4-6085-4f62-9a71-46d9bc1eb92b"

From e0404422835631ed2a9e9b2d19da7c0f8f26a75c Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Wed, 18 Sep 2024 11:10:21 +0200
Subject: [PATCH 11/18] LinearAlgebra

---
 src/AirfoilGmsh.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/AirfoilGmsh.jl b/src/AirfoilGmsh.jl
index 9708e60..86b6994 100644
--- a/src/AirfoilGmsh.jl
+++ b/src/AirfoilGmsh.jl
@@ -11,7 +11,7 @@ using Optim
 using Optimization, OptimizationBBO
 using Parameters
 using CubicSplines
-
+using LinearAlgebra
 
 export DomainDivision
 export DomainMeshDivisions

From 980f47f240f552250061dd542336cd7e8755a794 Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Mon, 23 Sep 2024 15:20:20 +0200
Subject: [PATCH 12/18] fix tests

---
 test/runtests.jl    | 1 +
 test/test_driver.jl | 6 +++---
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index a47d03c..2c1c1cc 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,7 @@
 using AirfoilGmsh
 using Test
 
+
 url_test = ["https://m-selig.ae.illinois.edu/ads/coord/c141a.dat", 
 "https://m-selig.ae.illinois.edu/ads/coord/e1098.dat",
 "https://m-selig.ae.illinois.edu/ads/coord/n0012.dat",
diff --git a/test/test_driver.jl b/test/test_driver.jl
index 62c56ee..4238837 100644
--- a/test/test_driver.jl
+++ b/test/test_driver.jl
@@ -2,10 +2,10 @@ function test_airfoil(url::String)
     fname= get_airfoil_name_test(url)
     csvname = fname*".csv"
     @test from_url_to_csv(url) == csvname
-    @test typeof(create_geofile(csvname)) == IOStream
+    @test typeof(create_geofile(csvname)) == Tuple{IOStream, String}
 
-    @test typeof(create_geofile(csvname, BoundaryLayer(Reynolds=200e3))) == IOStream
-    @test typeof(create_geofile(csvname, DomainInfo(dimension = 2) )) == IOStream
+    @test typeof(create_geofile(csvname, BoundaryLayer(Reynolds=200e3))) == Tuple{IOStream, String}
+    @test typeof(create_geofile(csvname, DomainInfo(dimension = 2) )) == Tuple{IOStream, String}
     # x,y,wl,wu = increase_resolution_airfoil(csvname,500; maxiters = 10, maxtime = 10)
     # @test typeof(x) <: Vector{Float64}
     # rm(csvname)

From 565e1b092ec6a389e2e4736e2515200f63f82164 Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Mon, 23 Sep 2024 15:21:09 +0200
Subject: [PATCH 13/18] doc create_geofile

---
 src/CreateGeoFile.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/CreateGeoFile.jl b/src/CreateGeoFile.jl
index 3c0d589..ed15312 100644
--- a/src/CreateGeoFile.jl
+++ b/src/CreateGeoFile.jl
@@ -26,7 +26,6 @@ It is possible to create a mesh with the following options:
 | Thetraedreal   | 3D        | TETRA()    |
 
 """
-
 function create_geofile(filename::String, domain_mesh_divisions::DomainMeshDivisions, domain_info::DomainInfo, boundary_layer::BoundaryLayer)
 
 

From dbd035c702ff2860256fa37370819bab67ccfcb2 Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Fri, 27 Sep 2024 09:53:55 +0200
Subject: [PATCH 14/18] bugfix CST

---
 src/OptimizationCST.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/OptimizationCST.jl b/src/OptimizationCST.jl
index cba00dc..a3dfdc6 100644
--- a/src/OptimizationCST.jl
+++ b/src/OptimizationCST.jl
@@ -131,8 +131,10 @@ function increase_resolution_airfoil(filename::String, N::Int64; dz = 0.0, w0 =
 
     split_idx = 4
     
+    ub = vcat(zeros(Int64, count(w0.<0)),    ones(Int64,count(w0.>0)))
+    lb = ub .- 1
     params = ( split_idx, xl,xu,dz,y0)
-    prob = Optimization.OptimizationProblem(error_function, w0, params, lb = [-1,-1,-1,0,0,0,0,0,0], ub = [0,0,0,1,1,1,1,1,1])
+    prob = Optimization.OptimizationProblem(error_function, w0, params, lb =lb, ub = ub)
     sol = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited(), maxiters = maxiters,  maxtime = maxtime)
     sol = collect(sol)
     wl, wu = compute_wlwu(sol, split_idx)
@@ -142,4 +144,4 @@ function increase_resolution_airfoil(filename::String, N::Int64; dz = 0.0, w0 =
     write_csv_cst(x,y,filename)
    
     return x,y,wl,wu
-end
\ No newline at end of file
+end

From 3c2f40e7655aad005b7cb95bd8e2c5f10a1f3ac1 Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Fri, 27 Sep 2024 10:39:38 +0200
Subject: [PATCH 15/18] add save_cst

---
 src/OptimizationCST.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/OptimizationCST.jl b/src/OptimizationCST.jl
index a3dfdc6..649de41 100644
--- a/src/OptimizationCST.jl
+++ b/src/OptimizationCST.jl
@@ -122,7 +122,7 @@ plot!(xlims =(0.0,1.0), ylims =(-0.2,0.65))
 plot!(xlabel = "x", ylabel = "y")
 ```
 """
-function increase_resolution_airfoil(filename::String, N::Int64; dz = 0.0, w0 = [-0.1294, -0.0036, -0.0666, -0.01, 0.206, 0.2728, 0.2292, 0.1, 0.1,0.1], maxiters = 100.0, maxtime=100.0)
+function increase_resolution_airfoil(filename::String, N::Int64; dz = 0.0, w0 = [-0.1294, -0.0036, -0.0666, -0.01, 0.206, 0.2728, 0.2292, 0.1, 0.1,0.1], maxiters = 100.0, maxtime=100.0, save_cst=false)
     xu,xl,yu,yl = get_airfoil_coordinates(filename)
     y0 = [yu;yl]
 
@@ -141,7 +141,7 @@ function increase_resolution_airfoil(filename::String, N::Int64; dz = 0.0, w0 =
    
     #Using wl,wu the new coordinates x,y are computed
     x,y = CST_airfoil(wl,wu,dz,N)
-    write_csv_cst(x,y,filename)
+    save_cst && write_csv_cst(x,y,filename)
    
     return x,y,wl,wu
-end
+end
\ No newline at end of file

From ddd7329c011c658347a4e148033d00b74852595d Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Fri, 27 Sep 2024 10:39:58 +0200
Subject: [PATCH 16/18] cn

---
 src/OptimizationCST.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/OptimizationCST.jl b/src/OptimizationCST.jl
index 649de41..e54fa66 100644
--- a/src/OptimizationCST.jl
+++ b/src/OptimizationCST.jl
@@ -122,7 +122,7 @@ plot!(xlims =(0.0,1.0), ylims =(-0.2,0.65))
 plot!(xlabel = "x", ylabel = "y")
 ```
 """
-function increase_resolution_airfoil(filename::String, N::Int64; dz = 0.0, w0 = [-0.1294, -0.0036, -0.0666, -0.01, 0.206, 0.2728, 0.2292, 0.1, 0.1,0.1], maxiters = 100.0, maxtime=100.0, save_cst=false)
+function increase_resolution_airfoil(filename::String, N::Int64; dz = 0.0, w0 = [-0.1294, -0.0036, -0.0666, -0.01, 0.206, 0.2728, 0.2292, 0.1, 0.1,0.1], maxiters = 100.0, maxtime=100.0, write_cst=false)
     xu,xl,yu,yl = get_airfoil_coordinates(filename)
     y0 = [yu;yl]
 
@@ -141,7 +141,7 @@ function increase_resolution_airfoil(filename::String, N::Int64; dz = 0.0, w0 =
    
     #Using wl,wu the new coordinates x,y are computed
     x,y = CST_airfoil(wl,wu,dz,N)
-    save_cst && write_csv_cst(x,y,filename)
+    write_cst && write_csv_cst(x,y,filename)
    
     return x,y,wl,wu
 end
\ No newline at end of file

From d98c7fd869102bc9630d75db2bd7edf95524c8e3 Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Fri, 27 Sep 2024 10:44:13 +0200
Subject: [PATCH 17/18] robust origin find

---
 src/OptimizationCST.jl | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/src/OptimizationCST.jl b/src/OptimizationCST.jl
index e54fa66..8e112ec 100644
--- a/src/OptimizationCST.jl
+++ b/src/OptimizationCST.jl
@@ -14,7 +14,9 @@ end
 Distinguish the upper and lower coordinates of the airfoil
 """
 function find_lower_upper(x::Vector{Float64},y::Vector{Float64})
-    origin_idx = findall(isapprox.(x,0.0))[1]
+    _,origin_idx = findmin(abs.(x) )
+
+    # origin_idx = findall(isapprox.(x,0.0))[1]
     n = length(x)
     if y[origin_idx+1]<0
         idx_upper = 1:origin_idx
@@ -26,6 +28,10 @@ function find_lower_upper(x::Vector{Float64},y::Vector{Float64})
     return x[idx_upper],x[idx_lower],y[idx_upper],y[idx_lower]
 end
 
+
+z = rand(10) .- 0.5
+
+
 """
     get_airfoil_coordinates(filename::String)
 Get the airfoil points coordinates form .csv file. It distinguish the top and bottom side

From ae6822db61c41a356b935462995dce0c6104cf58 Mon Sep 17 00:00:00 2001
From: carlodev <carlobrunellibs@gmail.com>
Date: Fri, 25 Oct 2024 10:50:33 +0200
Subject: [PATCH 18/18] fix

---
 src/OptimizationCST.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/OptimizationCST.jl b/src/OptimizationCST.jl
index 8e112ec..66fae69 100644
--- a/src/OptimizationCST.jl
+++ b/src/OptimizationCST.jl
@@ -18,7 +18,7 @@ function find_lower_upper(x::Vector{Float64},y::Vector{Float64})
 
     # origin_idx = findall(isapprox.(x,0.0))[1]
     n = length(x)
-    if y[origin_idx+1]<0
+    if y[origin_idx+1]<y[origin_idx-1]
         idx_upper = 1:origin_idx
         idx_lower = origin_idx+1:n
     else