checkpoint reached: YOU GOT A MESH, OH YES!

Author: amaurigmartins

.gitignore

@@ -44,3 +44,4 @@ showcase
 !showcase/*.jl
 papers
 sprint
+*/fem_output

examples/fem_output/tutorial/tutorial.geo_unrolled

@@ -36,29 +36,31 @@ addto_linecablesystem!(
 display(cable_system)
 
 # Create a FEMProblemDefinition with custom parameters
+domain_radius = 5.0
 problem = FEMProblemDefinition(
-    domain_radius=5.0,  # 5 m radius
+    domain_radius=domain_radius,
     correct_twisting=true,
-    elements_per_scale_length_conductor=3.0,
-    elements_per_scale_length_insulator=2.0,
-    elements_per_scale_length_semicon=4.0,
-    elements_per_scale_length_interfaces=0.1,
+    elements_per_length_conductor=1,
+    elements_per_length_insulator=2,
+    elements_per_length_semicon=1,
+    elements_per_length_interfaces=5,
+    points_per_circumference=8,
     analysis_type=[0, 1],  # Electrostatic and magnetostatic
-    mesh_size_min=1e-4,
-    mesh_size_max=0.5,
-    mesh_size_default=0.01,
-    mesh_algorithm=6,  # Frontal-Delaunay
+    mesh_size_min=1e-6,
+    mesh_size_max=domain_radius / 5,
+    mesh_size_default=domain_radius / 10,
+    mesh_algorithm=5,
     materials_db=materials_db
 )
 
 # Create a FEMSolver with custom parameters
 solver = FEMSolver(
     force_remesh=true,  # Force remeshing
     run_solver=false,
-    preview_geo=true,  # Preview geometry
-    preview_mesh=false,  # Preview the mesh
+    preview_geo=false,  # Preview geometry
+    preview_mesh=true,  # Preview the mesh
     base_path=joinpath(@__DIR__, "fem_output"),
-    verbosity=2,  # Verbose output
+    verbosity=0,  # Verbose output
     getdp_executable=joinpath("/home/amartins/Applications/onelab-Linux64", "getdp"),  # Path to GetDP executable
 )
 

examples/fem_output/tutorial/tutorial.msh

@@ -203,18 +203,18 @@ struct FEMProblemDefinition <: AbstractProblemDefinition
     domain_radius::Float64
     "Flag to correct for twisting effects \\[dimensionless\\]."
     correct_twisting::Bool
-    "Elements per scale length for conductors \\[dimensionless\\]."
-    elements_per_scale_length_conductor::Float64
-    "Elements per scale length for insulators \\[dimensionless\\]."
-    elements_per_scale_length_insulator::Float64
-    "Elements per scale length for semiconductors \\[dimensionless\\]."
-    elements_per_scale_length_semicon::Float64
-    "Elements per scale length for interfaces \\[dimensionless\\]."
-    elements_per_scale_length_interfaces::Float64
+    "Elements per characteristic length for conductors \\[dimensionless\\]."
+    elements_per_length_conductor::Int
+    "Elements per characteristic length for insulators \\[dimensionless\\]."
+    elements_per_length_insulator::Int
+    "Elements per characteristic length for semiconductors \\[dimensionless\\]."
+    elements_per_length_semicon::Int
+    "Elements per characteristic length for interfaces \\[dimensionless\\]."
+    elements_per_length_interfaces::Int
+    "Points per circumference length (2π radians) \\[dimensionless\\]."
+    points_per_circumference::Int
     "Analysis types to perform \\[dimensionless\\]."
     analysis_type::Vector{Int}
-    "Materials database."
-    materials_db::MaterialsLibrary
 
     "Minimum mesh size \\[m\\]."
     mesh_size_min::Float64
@@ -225,6 +225,9 @@ struct FEMProblemDefinition <: AbstractProblemDefinition
     "Mesh algorithm to use \\[dimensionless\\]."
     mesh_algorithm::Int
 
+    "Materials database."
+    materials_db::MaterialsLibrary
+
     """
     $(TYPEDSIGNATURES)
 
@@ -234,10 +237,10 @@ struct FEMProblemDefinition <: AbstractProblemDefinition
 
     - `domain_radius`: Domain radius for the simulation \\[m\\]. Default: 5.0.
     - `correct_twisting`: Flag to correct for twisting effects \\[dimensionless\\]. Default: true.
-    - `elements_per_scale_length_conductor`: Elements per scale length for conductors \\[dimensionless\\]. Default: 3.0.
-    - `elements_per_scale_length_insulator`: Elements per scale length for insulators \\[dimensionless\\]. Default: 2.0.
-    - `elements_per_scale_length_semicon`: Elements per scale length for semiconductors \\[dimensionless\\]. Default: 4.0.
-    - `elements_per_scale_length_interfaces`: Elements per scale length for interfaces \\[dimensionless\\]. Default: 0.1.
+    - `elements_per_length_conductor`: Elements per scale length for conductors \\[dimensionless\\]. Default: 3.0.
+    - `elements_per_length_insulator`: Elements per scale length for insulators \\[dimensionless\\]. Default: 2.0.
+    - `elements_per_length_semicon`: Elements per scale length for semiconductors \\[dimensionless\\]. Default: 4.0.
+    - `elements_per_length_interfaces`: Elements per scale length for interfaces \\[dimensionless\\]. Default: 0.1.
     - `analysis_type`: Analysis types to perform \\[dimensionless\\]. Default: [0, 1].
     - `mesh_size_min`: Minimum mesh size \\[m\\]. Default: 1e-4.
     - `mesh_size_max`: Maximum mesh size \\[m\\]. Default: 1.0.
@@ -258,18 +261,19 @@ struct FEMProblemDefinition <: AbstractProblemDefinition
     # Create a problem definition with custom parameters
     problem_def = $(FUNCTIONNAME)(
         domain_radius=10.0,
-        elements_per_scale_length_conductor=5.0,
+        elements_per_length_conductor=5.0,
         mesh_algorithm=2
     )
     ```
     """
     function FEMProblemDefinition(;
         domain_radius::Float64=5.0,
         correct_twisting::Bool=true,
-        elements_per_scale_length_conductor::Float64=3.0,
-        elements_per_scale_length_insulator::Float64=2.0,
-        elements_per_scale_length_semicon::Float64=4.0,
-        elements_per_scale_length_interfaces::Float64=0.1,
+        elements_per_length_conductor::Int=3,
+        elements_per_length_insulator::Int=2,
+        elements_per_length_semicon::Int=4,
+        elements_per_length_interfaces::Int=3,
+        points_per_circumference::Int=16,
         analysis_type::Vector{Int}=[0, 1],
         mesh_size_min::Float64=1e-4,
         mesh_size_max::Float64=1.0,
@@ -279,10 +283,8 @@ struct FEMProblemDefinition <: AbstractProblemDefinition
     )
         return new(
             domain_radius, correct_twisting,
-            elements_per_scale_length_conductor, elements_per_scale_length_insulator,
-            elements_per_scale_length_semicon, elements_per_scale_length_interfaces,
-            analysis_type, materials_db,
-            mesh_size_min, mesh_size_max, mesh_size_default, mesh_algorithm
+            elements_per_length_conductor, elements_per_length_insulator,
+            elements_per_length_semicon, elements_per_length_interfaces, points_per_circumference, analysis_type, mesh_size_min, mesh_size_max, mesh_size_default, mesh_algorithm, materials_db
         )
     end
 end
@@ -572,8 +574,8 @@ function run_fem_model(cable_system::LineCableSystem,
             _assign_physical_groups(workspace)
 
             # Mesh sizing
-            # _log(workspace, 1, "Setting up physics-based mesh sizing...")
-            # _config_mesh_sizes(workspace)
+            _log(workspace, 1, "Setting up physics-based mesh sizing...")
+            _config_mesh_sizes(workspace)
 
             # Preview pre-meshing configuration if requested
             if solver.preview_geo

examples/fem_output/tutorial3/tutorial3.geo_unrolled

@@ -130,9 +130,6 @@ function _make_cablepart!(workspace::FEMWorkspace, part::AbstractCablePart,
         material_id     # Material ID from registry
     )
 
-    # Calculate mesh size for this part
-    mesh_size = calc_mesh_size(part, workspace)
-
     # Create physical name
     part_type = lowercase(string(nameof(typeof(part))))
     elementary_name = create_cable_elementary_name(
@@ -148,15 +145,39 @@ function _make_cablepart!(workspace::FEMWorkspace, part::AbstractCablePart,
     radius_in = to_nominal(part.radius_in)
     radius_ext = to_nominal(part.radius_ext)
 
+    # Calculate mesh size for this part
+    if part isa AbstractConductorPart
+        num_elements = workspace.problem_def.elements_per_length_conductor
+    elseif part isa Insulator
+        num_elements = workspace.problem_def.elements_per_length_insulator
+    elseif part isa Semicon
+        num_elements = workspace.problem_def.elements_per_length_semicon
+    end
+
+    mesh_size_current = calc_mesh_size(radius_in, radius_ext, part.material_props, num_elements, workspace)
+
+    # Calculate mesh size for the next part
+    num_layers = length(cabledef.cable.components[comp_idx].conductor_group.layers)
+    next_part = layer_idx < num_layers ? cabledef.cable.components[comp_idx].conductor_group.layers[layer_idx+1] : nothing
+
+    if !isnothing(next_part)
+        next_radius_in = to_nominal(next_part.radius_in)
+        next_radius_ext = to_nominal(next_part.radius_ext)
+        mesh_size_next = calc_mesh_size(next_radius_in, next_radius_ext, next_part.material_props, num_elements, workspace)
+    else
+        mesh_size_next = mesh_size_current
+    end
+
+    mesh_size = min(mesh_size_current, mesh_size_next)
+    num_points_circumference = workspace.problem_def.points_per_circumference
+
     # Create annular shape and assign marker
     if radius_in ≈ 0
         # Solid disk
-        marker = _get_disk_marker(x_center, y_center)
-        entity_tag = _draw_disk(x_center, y_center, radius_ext, mesh_size)
+        _, _, marker = _draw_disk(x_center, y_center, radius_ext, mesh_size, num_points_circumference)
     else
         # Annular shape
-        marker = _get_annular_marker(x_center, y_center, radius_in, radius_ext)
-        entity_tag = _draw_annular(x_center, y_center, radius_in, radius_ext, mesh_size)
+        _, _, marker = _draw_annular(x_center, y_center, radius_in, radius_ext, mesh_size, num_points_circumference)
     end
 
     # Create entity data
@@ -219,34 +240,67 @@ function _make_cablepart!(workspace::FEMWorkspace, part::WireArray,
         material_id     # Material ID from registry
     )
 
-    # Calculate mesh size for this part
-    mesh_size = calc_mesh_size(part, workspace)
-
-    #
-    # First handle the wires
-    #
+    # -------- First handle the wires
 
     # Create physical name
     part_type = lowercase(string(nameof(typeof(part))))
 
     # Extract parameters
     radius_in = to_nominal(part.radius_in)
+    radius_ext = to_nominal(part.radius_ext)
 
-    radius_wire = to_nominal(part.radius_wire)
+    TOL = 5e-6 # Shrink the radius to avoid overlapping boundaries, this must be greater than Gmsh geometry tolerance
+    radius_wire = to_nominal(part.radius_wire) - TOL
     num_wires = part.num_wires
-    lay_radius = num_wires == 1 ? 0 : to_nominal(part.radius_in)
+
+
+    # Calculate mesh size for this part
+    num_elements = workspace.problem_def.elements_per_length_conductor
+    mesh_size_current = calc_mesh_size(radius_in, radius_ext, part.material_props, num_elements, workspace)
+
+    # Calculate mesh size for the next part
+    num_layers = length(cabledef.cable.components[comp_idx].conductor_group.layers)
+    next_part = layer_idx < num_layers ? cabledef.cable.components[comp_idx].conductor_group.layers[layer_idx+1] : nothing
+
+    if !isnothing(next_part)
+        next_radius_in = to_nominal(next_part.radius_in)
+        next_radius_ext = to_nominal(next_part.radius_ext)
+        mesh_size_next = calc_mesh_size(next_radius_in, next_radius_ext, next_part.material_props, num_elements, workspace)
+    else
+        mesh_size_next = mesh_size_current
+    end
+
+    mesh_size = min(mesh_size_current, mesh_size_next)
+    num_points_circumference = workspace.problem_def.points_per_circumference
 
     # Calculate wire positions
-    wire_positions = calc_wirearray_coords(num_wires, radius_wire, lay_radius, C=(x_center, y_center))
+    function _calc_wirearray_coords(
+        num_wires::Number,
+        # radius_wire::Number,
+        radius_in::Number,
+        radius_ext::Number;
+        C=(0.0, 0.0),
+    )
+        wire_coords = []  # Global coordinates of all wires
+        lay_radius = num_wires == 1 ? 0 : (radius_in + radius_ext) / 2
+
+        # Calculate the angle between each wire
+        angle_step = 2 * π / num_wires
+        for i in 0:num_wires-1
+            angle = i * angle_step
+            x = C[1] + lay_radius * cos(angle)
+            y = C[2] + lay_radius * sin(angle)
+            push!(wire_coords, (x, y))  # Add wire center
+        end
+        return wire_coords
+    end
 
+    wire_positions = _calc_wirearray_coords(num_wires, radius_in, radius_ext, C=(x_center, y_center))
 
     # Create wires
     for (wire_idx, (wx, wy)) in enumerate(wire_positions)
-        # Create wire marker
-        marker = _get_disk_marker(wx, wy)
 
-        # Create wire disk
-        entity_tag = _draw_disk(wx, wy, radius_wire, mesh_size)
+        _, _, marker = _draw_disk(wx, wy, radius_wire, mesh_size, num_points_circumference)
 
         # Create wire name
         elementary_name = create_cable_elementary_name(
@@ -267,9 +321,24 @@ function _make_cablepart!(workspace::FEMWorkspace, part::WireArray,
         workspace.unassigned_entities[marker] = entity_data
     end
 
-    #
-    # Then those nasty air gaps, the cause of this entire suffering
-    #
+    # Handle WireArray outermost boundary
+    mesh_size = (radius_ext - radius_in)
+    if !(next_part isa WireArray) && !isnothing(next_part)
+        # step_angle = 2 * pi / num_wires
+        _add_mesh_points(
+            radius_in=radius_ext,
+            radius_ext=radius_ext,
+            theta_0=0,
+            theta_1=2 * pi,
+            mesh_size=mesh_size,
+            num_points_ang=num_points_circumference,
+            num_points_rad=0,
+            C=(x_center, y_center),
+            theta_offset=0 #step_angle / 2
+        )
+    end
+
+    # --------Then those nasty air gaps
 
     # Air gaps will be determined from the boolean fragmentation operation and do not need to be drawn. Only the markers are needed.
     markers_air_gap = _get_air_gap_markers(num_wires, radius_wire, radius_in)
@@ -298,9 +367,6 @@ function _make_cablepart!(workspace::FEMWorkspace, part::WireArray,
         material_id     # Material ID from registry
     )
 
-    # Calculate mesh size for this part
-    mesh_size = calc_mesh_size(part, workspace)
-
     for marker in markers_air_gap
         # elementary names are not assigned to the air gaps because they are not drawn and appear as a result of the boolean operation
         core_data = CoreEntityData(physical_group_tag_air_gap, "", mesh_size)