API Reference - Mesh
MeshInputError
Bases: ValueError
Raised for mesh inputs gmsh cannot handle — e.g. a reinforcement/constraint line lying on the domain boundary — so callers can report a clear message instead of a raw gmsh 'Unable to recover edge' crash.
Source code in xslope/mesh.py
class MeshInputError(ValueError):
"""Raised for mesh inputs gmsh cannot handle — e.g. a reinforcement/constraint
line lying on the domain boundary — so callers can report a clear message
instead of a raw gmsh 'Unable to recover edge' crash."""
pass
add_dload_points_to_polygons(polygons, slope_data)
Add distributed load points to polygon edges if they are coincident with edges but not existing vertices.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def add_dload_points_to_polygons(polygons, slope_data):
"""
Add distributed load points to polygon edges if they are coincident with edges
but not existing vertices.
Parameters:
polygons: List of polygons (list of (x,y) tuples) or dicts with "coords"
slope_data: Dictionary containing slope data
Returns:
Updated list of polygons with added points
"""
import numpy as np
tol = 1e-8
# Collect distributed load points to check
points_to_check = []
# Add distributed load points from 'dloads' and 'dloads2' keys
# Each is a list of load lines; each load line is a list of dicts with X, Y, Normal
for key in ('dloads', 'dloads2'):
for load_line in slope_data.get(key, []):
if isinstance(load_line, list):
for pt in load_line:
if isinstance(pt, dict) and 'X' in pt and 'Y' in pt:
points_to_check.append((pt['X'], pt['Y']))
elif isinstance(load_line, dict) and 'xy' in load_line:
for point in load_line['xy']:
points_to_check.append(point)
return _add_points_to_polygons(polygons, points_to_check, tol)
add_intersection_points_to_polygons(polygons, lines, debug=False)
Add intersection points between reinforcement lines and polygon edges to the polygon vertex lists. This ensures that polygons have vertices at all intersection points with reinforcement lines.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def add_intersection_points_to_polygons(polygons, lines, debug=False):
"""
Add intersection points between reinforcement lines and polygon edges to the polygon vertex lists.
This ensures that polygons have vertices at all intersection points with reinforcement lines.
Parameters:
polygons: List of polygons (lists of (x,y) tuples) or dicts with "coords"
lines: List of reinforcement lines (lists of (x,y) tuples)
debug: Enable debug output
Returns:
Updated list of polygons with intersection points added
"""
if not lines:
return polygons
if debug:
print("Adding intersection points to polygons...")
# Make a copy of polygons to modify
updated_polygons = []
for poly in polygons:
if isinstance(poly, dict):
updated_entry = dict(poly)
updated_entry["coords"] = list(poly.get("coords", []))
updated_polygons.append(updated_entry)
else:
updated_polygons.append(list(poly)) # Convert to list for modification
# Find all intersections
for line_idx, line_pts in enumerate(lines):
line_pts_clean = remove_duplicate_endpoint(list(line_pts))
if debug:
print(f"Processing line {line_idx}: {line_pts_clean}")
# Check each segment of the reinforcement line
for i in range(len(line_pts_clean) - 1):
line_seg_start = line_pts_clean[i]
line_seg_end = line_pts_clean[i + 1]
# Check intersection with each polygon
for poly_idx, poly in enumerate(updated_polygons):
poly_coords = poly.get("coords", []) if isinstance(poly, dict) else poly
# Check each edge of this polygon
for j in range(len(poly_coords)):
poly_edge_start = poly_coords[j]
poly_edge_end = poly_coords[(j + 1) % len(poly_coords)]
# Find intersection point if it exists
intersection = line_segment_intersection(
line_seg_start, line_seg_end,
poly_edge_start, poly_edge_end
)
if intersection:
if debug:
print(f"Found intersection {intersection} between line {line_idx} segment {i} and polygon {poly_idx} edge {j}")
# Check if intersection point is already a vertex of this polygon
is_vertex = False
for vertex in poly_coords:
if abs(vertex[0] - intersection[0]) < 1e-8 and abs(vertex[1] - intersection[1]) < 1e-8:
is_vertex = True
break
if not is_vertex:
# Insert intersection point into polygon at the correct position
# Insert after vertex j (which is the start of the edge)
insert_idx = j + 1
if isinstance(updated_polygons[poly_idx], dict):
updated_polygons[poly_idx]["coords"].insert(insert_idx, intersection)
else:
updated_polygons[poly_idx].insert(insert_idx, intersection)
if debug:
print(f"Added intersection point {intersection} to polygon {poly_idx} at position {insert_idx}")
return updated_polygons
add_seep_bc_points_to_polygons(polygons, slope_data)
Insert the vertices of every seepage BC polyline as polygon vertices.
Exactly the same need the distributed loads have. A specified-flux BC is applied EDGE-wise — a boundary edge carries load only when BOTH its corner nodes lie on the polyline — so an endpoint landing part-way along an element edge drops that whole edge and silently applies less water than the user asked for, the same way a dload end that straddles an edge under-applies its traction. Pinning a node at each vertex makes the loaded length equal the specified length exactly.
Head and exit-face BCs are matched node-wise rather than edge-wise, so a stray endpoint only rounds their extent to the nearest node. Pinning their vertices makes those extents exact too, which matters most for an exit face: its endpoints decide where the phreatic surface is allowed to emerge.
Interior vertices are pinned as well as endpoints — a polyline can kink at a point that is not a corner of the material polygons, and an edge spanning that kink would otherwise cut across it.
Source code in xslope/mesh.py
def add_seep_bc_points_to_polygons(polygons, slope_data):
"""Insert the vertices of every seepage BC polyline as polygon vertices.
Exactly the same need the distributed loads have. A specified-flux BC is applied
EDGE-wise — a boundary edge carries load only when BOTH its corner nodes lie on
the polyline — so an endpoint landing part-way along an element edge drops that
whole edge and silently applies less water than the user asked for, the same way a
dload end that straddles an edge under-applies its traction. Pinning a node at
each vertex makes the loaded length equal the specified length exactly.
Head and exit-face BCs are matched node-wise rather than edge-wise, so a stray
endpoint only rounds their extent to the nearest node. Pinning their vertices
makes those extents exact too, which matters most for an exit face: its endpoints
decide where the phreatic surface is allowed to emerge.
Interior vertices are pinned as well as endpoints — a polyline can kink at a point
that is not a corner of the material polygons, and an edge spanning that kink would
otherwise cut across it.
"""
pts = []
for key in ('seepage_bc', 'seepage_bc2'):
bc = slope_data.get(key) or {}
for block in ((bc.get('specified_heads') or [])
+ (bc.get('specified_fluxes') or [])):
pts.extend(tuple(c) for c in (block.get('coords') or []))
pts.extend(tuple(c) for c in (bc.get('exit_face') or []))
return _add_points_to_polygons(polygons, pts)
build_mesh_from_polygons(polygons, target_size, element_type='tri3', lines=None, debug=False, mesh_params=None, target_size_1d=None, profile_lines=None, point_constraints=None, refine_factor=None, refine_features=None, material_k=None)
Build a finite element mesh with material regions using Gmsh. Fixed version that properly handles shared boundaries between polygons.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def build_mesh_from_polygons(polygons, target_size, element_type='tri3', lines=None, debug=False, mesh_params=None, target_size_1d=None, profile_lines=None, point_constraints=None, refine_factor=None, refine_features=None, material_k=None):
"""
Build a finite element mesh with material regions using Gmsh.
Fixed version that properly handles shared boundaries between polygons.
Parameters:
polygons : List of polygon coordinate lists or dicts with "coords"/"mat_id"
target_size : Desired element size
element_type : 'tri3' (3-node triangles), 'tri6' (6-node triangles),
'quad4' (4-node quadrilaterals), 'quad8' (8-node quadrilaterals),
'quad9' (9-node quadrilaterals)
lines : Optional list of lines, each defined by list of (x, y) tuples for 1D elements
debug : Enable debug output
mesh_params : Optional dictionary of GMSH meshing parameters to override defaults
target_size_1d : Optional target size for 1D elements (default None, which is set to target_size if None)
profile_lines: Optional list of profile line dicts with 'mat_id' keys for material assignment
refine_factor: Optional feature-aware auto-refinement. None (default) = OFF; the
mesh is byte-identical to the historical output. A value > 1 drives the
local element size down to target_size/refine_factor near model features
(reinforcement/pile lines, crack/notch tips, thin material zones) using
gmsh native size fields, growing smoothly back to target_size away from
them. Crack tips refine twice as strongly.
refine_features: Optional list selecting which feature classes to refine near, from
{'reinforcement','piles','cracks','thin_zones','interfaces'}. None =
the default four (NOT 'interfaces', which is seepage-specific and opt-in).
Ignored when refine_factor is None.
material_k : Optional mapping region_id -> major hydraulic conductivity k1 (a dict, or
a sequence indexed by region_id). Only consulted when 'interfaces' is in
refine_features: material boundaries whose two sides differ in k1 by >= 100x
get the same Distance+Threshold band as the line features. Ignored otherwise.
Returns:
mesh dict containing:
nodes : np.ndarray of node coordinates (n_nodes, 2)
elements : np.ndarray of 2D element vertex indices (n_elements, 9) - unused nodes set to 0
element_types: np.ndarray indicating number of nodes per 2D element (3, 4, 6, 8, or 9)
element_materials: np.ndarray of material ID for each 2D element
If lines is provided, also includes:
elements_1d : np.ndarray of 1D element vertex indices (n_elements_1d, 3) - unused nodes set to 0
element_types_1d: np.ndarray indicating element type (2 for linear, 3 for quadratic)
element_materials_1d: np.ndarray of material ID for each 1D element (line index)
"""
gmsh = _get_gmsh()
from collections import defaultdict
# Set default target_size_1d if None
if target_size_1d is None:
target_size_1d = target_size
if debug:
print(f"Using default target_size_1d = target_size = {target_size_1d}")
# Validate / normalize feature-aware refinement. refine_factor is None => OFF
# and NOTHING below changes (byte-identical to the historical mesh).
if refine_factor is not None:
if not (refine_factor > 1.0):
raise ValueError(
"refine_factor must be > 1.0 (local size = target_size/refine_factor); "
f"got {refine_factor!r}. Use None to disable refinement.")
if refine_features is None:
# Default set excludes the opt-in 'interfaces' (seepage-specific).
refine_set = set(_REFINE_FEATURES)
else:
refine_set = set(refine_features)
unknown = refine_set - set(_REFINE_FEATURES_ALL)
if unknown:
raise ValueError(
f"refine_features contains unknown entries {sorted(unknown)}; "
f"valid options are {sorted(_REFINE_FEATURES_ALL)}")
else:
refine_set = set()
# Normalize polygons to coordinate lists and optional mat_id
polygon_coords = []
polygon_mat_ids = []
for i, polygon in enumerate(polygons):
if isinstance(polygon, dict):
polygon_coords.append(polygon.get("coords", []))
polygon_mat_ids.append(polygon.get("mat_id"))
else:
polygon_coords.append(polygon)
polygon_mat_ids.append(None)
# Point constraints (e.g. line-load application points): insert each point as
# a vertex into every polygon edge that contains it, so gmsh places a node
# exactly there. Done BEFORE the conforming pass so shared edges stay welded.
if point_constraints:
_insert_point_constraints(polygon_coords, point_constraints,
tol=1e-6 * max(1.0, target_size), debug=debug)
# Make adjacent zones conforming: split any edge at a neighbour's vertex that
# lies in its interior (a T-junction), so shared interfaces mesh without slits.
make_polygons_conforming(polygon_coords, debug=debug)
# Reject constraint/reinforcement lines gmsh can't embed (on the boundary /
# outside the domain) up front, with a clear message instead of a gmsh crash.
if lines:
_validate_constraint_lines(lines, polygon_coords)
# Build a list of region ids (list of material IDs - one per polygon)
if any(mat_id is not None for mat_id in polygon_mat_ids):
region_ids = [
mat_id if mat_id is not None else i
for i, mat_id in enumerate(polygon_mat_ids)
]
elif profile_lines and len(profile_lines) >= len(polygon_coords):
region_ids = []
for i in range(len(polygon_coords)):
mat_id = profile_lines[i].get('mat_id')
if mat_id is not None:
region_ids.append(mat_id)
else:
# Fallback to polygon index if no mat_id
region_ids.append(i)
else:
# Fallback to sequential IDs if no profile_lines provided
region_ids = [i for i in range(len(polygon_coords))]
if element_type not in ['tri3', 'tri6', 'quad4', 'quad8', 'quad9']:
raise ValueError("element_type must be 'tri3', 'tri6', 'quad4', 'quad8', or 'quad9'")
# Determine if we need quadratic elements - but always generate linear first
quadratic = element_type in ['tri6', 'quad8', 'quad9']
# For quadratic elements, always start with linear base element
if quadratic:
if element_type == 'tri6':
base_element_type = 'tri3'
elif element_type in ['quad8', 'quad9']:
base_element_type = 'quad4'
if debug:
print(f"Quadratic element '{element_type}' requested: generating '{base_element_type}' first, then post-processing")
else:
base_element_type = element_type
# Adjust target_size for quads to compensate for recombination creating finer meshes
if element_type.startswith('quad'):
# Different adjustment factors based on meshing parameters
if mesh_params and 'size_factor' in mesh_params:
size_factor = mesh_params['size_factor']
else:
# Default size factors for different approaches
if mesh_params and mesh_params.get("Mesh.RecombinationAlgorithm") == 0:
size_factor = 1.2 # Fast algorithm needs less adjustment
elif mesh_params and mesh_params.get("Mesh.RecombineOptimizeTopology", 0) > 50:
size_factor = 1.8 # High optimization creates more elements
else:
size_factor = 1.4 # Default
adjusted_target_size = target_size * size_factor
if debug:
print(f"Adjusted target size for quads: {target_size} -> {adjusted_target_size} (factor: {size_factor})")
else:
adjusted_target_size = target_size
# gmsh installs a SIGINT handler when interruptible; signal handlers can only be
# set from the main thread, so disable it when meshing off-thread (e.g. a GUI
# worker). On the main thread keep it so Ctrl+C still interrupts long meshes.
import threading
gmsh.initialize(interruptible=(threading.current_thread() is threading.main_thread()))
gmsh.option.setNumber("General.Verbosity", 4) # Reduce verbosity
gmsh.model.add("multi_region_mesh")
# Global point map to ensure shared boundaries use the same points
point_map = {} # maps (x, y) to Gmsh point tag
# Track all unique edges and their usage
edge_map = {} # maps (pt1, pt2) tuple to line tag
edge_usage = defaultdict(list) # maps edge to list of (region_id, orientation)
def add_point(x, y, size_override=None):
key = (x, y)
if key not in point_map:
point_size = size_override if size_override is not None else adjusted_target_size
tag = gmsh.model.geo.addPoint(x, y, 0, point_size)
point_map[key] = tag
return point_map[key]
def get_edge_key(pt1, pt2):
"""Get canonical edge key (always smaller point first)"""
return (min(pt1, pt2), max(pt1, pt2))
# Feature-aware refinement, thin-zone boundary sizing (opt-in). A whole-thin
# material polygon (e.g. a soft band) refines to a finer element size, but the
# long-edge transfinite constraints below would otherwise pin its boundary at
# the coarse target size and block the size field from resolving across it. So
# pre-compute a finer transfinite element size for the boundary edges of each
# whole-thin polygon, keyed by the (rounded) coordinate pair, and honor it when
# the long-edge transfinite constraint is set. Empty (no override) when OFF.
thin_edge_size = {}
if refine_factor is not None and 'thin_zones' in refine_set:
_floor = target_size / (refine_factor * _REFINE_CRACK_TIP_MULT)
for _zone in detect_thin_zones(polygon_coords, target_size):
if _zone['kind'] != 'whole':
continue
_i = _zone['poly_index']
if not (0 <= _i < len(polygon_coords)):
continue
_size = min(max(_zone['size'], _floor), target_size)
_ring = remove_duplicate_endpoint(list(polygon_coords[_i]))
for _j in range(len(_ring)):
_a = _ring[_j]
_b = _ring[(_j + 1) % len(_ring)]
_key = frozenset([(round(_a[0], 9), round(_a[1], 9)),
(round(_b[0], 9), round(_b[1], 9))])
thin_edge_size[_key] = min(thin_edge_size.get(_key, _size), _size)
# First pass: Create all points and identify short edges
polygon_data = []
short_edge_points = set() # Points that are endpoints of short edges
# Pre-pass to identify short edges - improved logic
for idx, (poly_pts, region_id) in enumerate(zip(polygon_coords, region_ids)):
poly_pts_clean = remove_duplicate_endpoint(list(poly_pts))
for i in range(len(poly_pts_clean)):
p1 = poly_pts_clean[i]
p2 = poly_pts_clean[(i + 1) % len(poly_pts_clean)]
edge_length = ((p2[0] - p1[0])**2 + (p2[1] - p1[1])**2)**0.5
# Only mark as short edge if it's genuinely short AND not a major boundary
# Major boundaries should maintain consistent mesh sizing
is_major_boundary = False
# Check if this edge is part of a major boundary (long horizontal or vertical edge)
if abs(p2[0] - p1[0]) > adjusted_target_size * 5: # Long horizontal edge
is_major_boundary = True
elif abs(p2[1] - p1[1]) > adjusted_target_size * 5: # Long vertical edge
is_major_boundary = True
# Only apply short edge sizing if edge is genuinely short AND not a major boundary
if edge_length < adjusted_target_size and not is_major_boundary:
short_edge_points.add(p1)
short_edge_points.add(p2)
if debug:
print(f"Short edge found: {p1} to {p2}, length={edge_length:.2f}")
elif debug and edge_length < adjusted_target_size:
print(f"Short edge ignored (major boundary): {p1} to {p2}, length={edge_length:.2f}")
# Main pass: Create points with appropriate sizes
for idx, (poly_pts, region_id) in enumerate(zip(polygon_coords, region_ids)):
poly_pts_clean = remove_duplicate_endpoint(list(poly_pts)) # make a copy
pt_tags = []
for x, y in poly_pts_clean:
# Use larger size for points on short edges to discourage subdivision
# But be more conservative about when to apply this
if (x, y) in short_edge_points:
point_size = adjusted_target_size * 2.0 # Reduced from 3.0 to 2.0
pt_tags.append(add_point(x, y, point_size))
else:
pt_tags.append(add_point(x, y))
# Track edges for this polygon
edges = []
for i in range(len(pt_tags)):
pt1 = pt_tags[i]
pt2 = pt_tags[(i + 1) % len(pt_tags)]
edge_key = get_edge_key(pt1, pt2)
# Determine orientation: True if pt1 < pt2, False otherwise
forward = (pt1 < pt2)
# Store edge usage
edge_usage[edge_key].append((region_id, forward))
edges.append((pt1, pt2, edge_key, forward))
polygon_data.append({
'region_id': region_id,
'pt_tags': pt_tags,
'edges': edges
})
# Second pass: Create all unique lines and track short edges
short_edges = [] # Track short edges for later processing
for edge_key in edge_usage.keys():
pt1, pt2 = edge_key
line_tag = gmsh.model.geo.addLine(pt1, pt2)
edge_map[edge_key] = line_tag
# Calculate edge length from point coordinates
pt1_coords = None
pt2_coords = None
for (x, y), tag in point_map.items():
if tag == pt1:
pt1_coords = (x, y)
if tag == pt2:
pt2_coords = (x, y)
if pt1_coords and pt2_coords:
edge_length = ((pt2_coords[0] - pt1_coords[0])**2 + (pt2_coords[1] - pt1_coords[1])**2)**0.5
# Add transfinite constraints for long boundary edges to ensure consistent mesh sizing
# This prevents the creation of overly coarse elements along major boundaries
if edge_length > adjusted_target_size * 3: # Long edge
# Calculate how many elements should be along this edge
num_elements = max(3, int(edge_length / adjusted_target_size))
# Feature-aware refinement: a whole-thin polygon's boundary edge is
# divided at its finer element size so the size field can resolve
# across the band instead of being pinned to the coarse target.
if thin_edge_size:
_ek = frozenset([(round(pt1_coords[0], 9), round(pt1_coords[1], 9)),
(round(pt2_coords[0], 9), round(pt2_coords[1], 9))])
_sz = thin_edge_size.get(_ek)
if _sz:
num_elements = max(num_elements, int(round(edge_length / _sz)))
try:
gmsh.model.geo.mesh.setTransfiniteCurve(line_tag, num_elements)
if debug:
print(f"Set transfinite constraint on long edge: {pt1_coords} to {pt2_coords}, length={edge_length:.2f}, num_elements={num_elements}")
except Exception as e:
if debug:
print(f"Warning: Could not set transfinite constraint on edge {pt1_coords} to {pt2_coords}: {e}")
# Add transfinite constraints for short edges to prevent subdivision
# This forces GMSH to use exactly 2 nodes (start and end) for short edges
elif edge_length < adjusted_target_size:
try:
gmsh.model.geo.mesh.setTransfiniteCurve(line_tag, 2) # Exactly 2 nodes
if debug:
print(f"Set transfinite constraint on short edge: {pt1_coords} to {pt2_coords}, length={edge_length:.2f}, exactly 2 nodes")
except Exception as e:
if debug:
print(f"Warning: Could not set transfinite constraint on short edge {pt1_coords} to {pt2_coords}: {e}")
# Short edges are now handled by point sizing, no need for transfinite curves
# Ensure all polygon points (including intersection points) are created as GMSH points
# The intersection points were already added to polygons in build_polygons(),
# so we just need to ensure they exist as GMSH geometric entities
if lines is not None:
if debug:
print("Ensuring all polygon points (including intersections) are created as GMSH points...")
# Collect all points from all polygons to ensure they exist in GMSH
all_polygon_points = set()
for poly_data in polygon_data:
pt_tags = poly_data['pt_tags']
for tag in pt_tags:
# Find the coordinates for this point tag
for (x, y), point_tag in point_map.items():
if point_tag == tag:
all_polygon_points.add((x, y))
break
# Create any missing GMSH points
for x, y in all_polygon_points:
key = (x, y)
if key not in point_map:
pt_tag = gmsh.model.geo.addPoint(x, y, 0.0, adjusted_target_size * 0.5)
point_map[key] = pt_tag
if debug:
print(f"Created GMSH point for polygon vertex {key}: tag {pt_tag}")
if debug:
print(f"Ensured {len(all_polygon_points)} polygon points exist as GMSH entities")
# Create enhanced reinforcement lines that include intersection points from polygons
# This is essential for proper mesh generation with embedded 1D elements
# Snap constraint line endpoints to nearby polygon boundary points.
# This prevents near-zero-length elements when a line endpoint is close
# to (but not exactly on) a polygon boundary (e.g., pile top near ground surface).
snap_tol = adjusted_target_size * 0.05
poly_pts_list = list(all_polygon_points)
if lines is not None and poly_pts_list:
poly_pts_arr = np.array(poly_pts_list)
for line_idx in range(len(lines)):
snapped = list(lines[line_idx])
for i in [0, len(snapped) - 1]: # snap endpoints only
px, py = snapped[i]
dists = np.sqrt((poly_pts_arr[:, 0] - px)**2 + (poly_pts_arr[:, 1] - py)**2)
j = np.argmin(dists)
if dists[j] < snap_tol and dists[j] > 1e-12:
if debug:
print(f"Snapped line {line_idx} endpoint ({px:.4f},{py:.4f}) "
f"-> ({poly_pts_list[j][0]:.4f},{poly_pts_list[j][1]:.4f}) "
f"dist={dists[j]:.6f}")
snapped[i] = poly_pts_list[j]
lines[line_idx] = snapped
enhanced_lines = []
for line_idx, line_pts in enumerate(lines):
line_pts_clean = remove_duplicate_endpoint(list(line_pts))
# Collect all points for this line: original + intersection points from polygons
all_line_points = []
# Add original line points
for x, y in line_pts_clean:
all_line_points.append((x, y, 'original'))
# Add intersection points that are on this line (from polygon data)
for poly_data in polygon_data:
pt_tags = poly_data['pt_tags']
for tag in pt_tags:
# Find the coordinates for this point tag
for (x, y), point_tag in point_map.items():
if point_tag == tag:
# Check if this point is on the reinforcement line
if is_point_on_line_segments((x, y), line_pts_clean, tolerance=1e-6):
all_line_points.append((x, y, 'intersection'))
break
# Sort all points along the line to maintain proper order
if len(all_line_points) > 1:
all_line_points.sort(key=lambda p: line_segment_parameter((p[0], p[1]), line_pts_clean[0], line_pts_clean[-1]))
# Remove duplicates (keep first occurrence)
unique_points = []
seen = set()
for x, y, point_type in all_line_points:
point_key = (round(x, 8), round(y, 8)) # Round to avoid floating point issues
if point_key not in seen:
seen.add(point_key)
unique_points.append((x, y, point_type))
# Create the enhanced line
enhanced_line = [(x, y) for x, y, _ in unique_points]
enhanced_lines.append(enhanced_line)
if debug:
print(f"Enhanced line {line_idx}: {len(line_pts_clean)} original points -> {len(enhanced_line)} total points")
# Replace original lines with enhanced lines
lines = enhanced_lines
# Create reinforcement lines as geometric constraints to force 2D mesh edges
line_data = []
if lines is not None:
for line_idx, line_pts in enumerate(lines):
# Use the enhanced line coordinates (which include intersection points)
line_pts_clean = remove_duplicate_endpoint(list(line_pts))
# Create points for this reinforcement line
line_point_tags = []
# Create all points for this line (original + intersection points)
for x, y in line_pts_clean:
key = (x, y)
if key in point_map:
line_point_tags.append((x, y, point_map[key]))
else:
# Create new point with small mesh size to ensure it's preserved
pt_tag = gmsh.model.geo.addPoint(x, y, 0.0, adjusted_target_size * 0.5)
point_map[key] = pt_tag
line_point_tags.append((x, y, pt_tag))
# Sort points along the line to maintain proper order
line_point_tags.sort(key=lambda p: line_segment_parameter((p[0], p[1]), line_pts_clean[0], line_pts_clean[-1]))
# Extract just the point tags in order
pt_tags = [tag for _, _, tag in line_point_tags]
if debug:
print(f" Line {line_idx} points: {[(x, y) for x, y, _ in line_point_tags]}")
# Create line segments as geometric constraints with controlled meshing
line_tags = []
for i in range(len(pt_tags) - 1):
pt1, pt2 = pt_tags[i], pt_tags[i + 1]
# Calculate segment length to determine number of subdivisions
coord1 = None
coord2 = None
for (x, y), tag in point_map.items():
if tag == pt1:
coord1 = (x, y)
if tag == pt2:
coord2 = (x, y)
if coord1 and coord2:
segment_length = ((coord2[0] - coord1[0])**2 + (coord2[1] - coord1[1])**2)**0.5
# Calculate number of elements needed to achieve target_size_1d
# For segments longer than target_size_1d, we want multiple elements
# For segments shorter than target_size_1d, we still want at least 2 elements
if segment_length > target_size_1d:
num_elements = max(3, int(round(segment_length / target_size_1d)))
else:
num_elements = 2
if debug:
print(f" Segment {i}: length {segment_length:.2f}, creating {num_elements} elements")
line_tag = gmsh.model.geo.addLine(pt1, pt2)
line_tags.append(line_tag)
# Set transfinite constraint to create appropriate number of nodes
try:
gmsh.model.geo.mesh.setTransfiniteCurve(line_tag, num_elements)
if debug:
print(f" Set transfinite constraint on line segment {i}: {num_elements} nodes")
except Exception as e:
if debug:
print(f" Warning: Could not set transfinite constraint on segment {i}: {e}")
else:
# Fallback: create line with default 2 nodes
line_tag = gmsh.model.geo.addLine(pt1, pt2)
line_tags.append(line_tag)
try:
gmsh.model.geo.mesh.setTransfiniteCurve(line_tag, 2)
if debug:
print(f" Set transfinite constraint on line segment {i}: 2 nodes (fallback)")
except Exception as e:
if debug:
print(f" Warning: Could not set transfinite constraint on segment {i}: {e}")
# Store line data for later 1D element extraction
# Use the enhanced line coordinates (which include intersection points)
line_data.append({
'line_idx': line_idx,
'line_tags': line_tags,
'point_coords': line_pts_clean # This now contains the enhanced coordinates
})
if debug:
print(f"Created reinforcement constraint line {line_idx} with {len(line_tags)} segments: {line_pts_clean}")
# Third pass: Create surfaces using the shared lines
surface_to_region = {}
# Surface tag per polygon (in polygon_coords order); None where creation failed.
# Used by feature-aware thin-zone refinement to Restrict a size field to a
# whole-thin material polygon's own surface (follows an inclined band exactly).
surface_tags_by_polygon = []
for poly_data in polygon_data:
region_id = poly_data['region_id']
edges = poly_data['edges']
line_tags = []
for pt1, pt2, edge_key, forward in edges:
line_tag = edge_map[edge_key]
# Use positive or negative line tag based on orientation
if forward:
line_tags.append(line_tag)
else:
line_tags.append(-line_tag)
# Create curve loop and surface
try:
loop = gmsh.model.geo.addCurveLoop(line_tags)
surface = gmsh.model.geo.addPlaneSurface([loop])
surface_to_region[surface] = region_id
surface_tags_by_polygon.append(surface)
except Exception as e:
print(f"Warning: Could not create surface for region {region_id}: {e}")
surface_tags_by_polygon.append(None)
continue
# Synchronize geometry
gmsh.model.geo.synchronize()
# Force mesh edges along reinforcement lines by creating additional geometric constraints
if lines is not None:
for line_info in line_data:
line_idx = line_info['line_idx']
line_tags = line_info['line_tags']
line_pts = line_info['point_coords']
# Set transfinite constraints to force mesh edges along each line segment
# REMOVED: This was conflicting with the target_size_1d calculations above
# for i, line_tag in enumerate(line_tags):
# try:
# # Force exactly 2 nodes (start and end) to prevent subdivision
# gmsh.model.geo.mesh.setTransfiniteCurve(line_tag, 2)
# if debug:
# print(f"Set transfinite constraint on line {line_idx} segment {i}: exactly 2 nodes")
# except Exception as e:
# if debug:
# print(f"Warning: Could not set transfinite constraint on line {line_idx} segment {i}: {e}")
# Embed reinforcement lines in all surfaces to ensure they're part of the mesh
for surface in surface_to_region.keys():
try:
# Embed all line segments of this reinforcement line
gmsh.model.mesh.embed(1, line_tags, 2, surface)
if debug:
print(f"Embedded reinforcement line {line_idx} in surface {surface}")
except Exception as e:
if debug:
print(f"Could not embed line {line_idx} in surface {surface}: {e}")
# CRITICAL: Set mesh coherence to ensure shared nodes along boundaries
# This forces Gmsh to use the same nodes for shared geometric entities
gmsh.model.mesh.removeDuplicateNodes()
# Create physical groups for material regions (this helps with mesh consistency)
physical_surfaces = []
for surface, region_id in surface_to_region.items():
physical_tag = gmsh.model.addPhysicalGroup(2, [surface])
physical_surfaces.append((physical_tag, region_id))
# Create physical groups for embedded reinforcement lines
physical_lines = []
if lines is not None:
for line_info in line_data:
line_idx = line_info['line_idx']
line_tags = line_info['line_tags']
physical_tag = gmsh.model.addPhysicalGroup(1, line_tags)
physical_lines.append((physical_tag, line_idx))
# Check for potential quad4 + reinforcement line conflicts
has_reinforcement_lines = lines is not None and len(lines) > 0
wants_quads = base_element_type.startswith('quad')
# Set mesh algorithm and recombination options BEFORE generating mesh
if base_element_type.startswith('quad'):
# Check if we need to use a more robust algorithm for reinforcement lines
if has_reinforcement_lines:
if debug:
print(f"Detected quad elements with reinforcement lines.")
print(f"Using robust recombination algorithm to handle embedded line constraints.")
# Use 'fast' algorithm which is more robust with embedded constraints
default_params = {
"Mesh.Algorithm": 8, # Frontal-Delaunay for quads
"Mesh.RecombineAll": 1, # Recombine triangles into quads
"Mesh.RecombinationAlgorithm": 0, # Standard (more robust than simple)
"Mesh.SubdivisionAlgorithm": 0, # Mixed tri/quad where needed
"Mesh.RecombineOptimizeTopology": 0, # Minimal optimization
"Mesh.RecombineNodeRepositioning": 1, # Still reposition nodes
"Mesh.RecombineMinimumQuality": 0.01, # Keep quality threshold
"Mesh.Smoothing": 5, # Reduced smoothing
"Mesh.SmoothNormals": 1, # Keep smooth normals
"Mesh.SmoothRatio": 1.8, # Keep smoothing ratio
}
else:
# Standard quad meshing parameters for cases without reinforcement lines
default_params = {
"Mesh.Algorithm": 8, # Frontal-Delaunay for quads (try 5, 6, 8)
"Mesh.RecombineAll": 1, # Recombine triangles into quads
"Mesh.RecombinationAlgorithm": 1, # Simple recombination (try 0, 1, 2, 3)
"Mesh.SubdivisionAlgorithm": 1, # All quads (try 0, 1, 2)
"Mesh.RecombineOptimizeTopology": 5, # Optimize topology (0-100)
"Mesh.RecombineNodeRepositioning": 1, # Reposition nodes (0 or 1)
"Mesh.RecombineMinimumQuality": 0.01, # Minimum quality threshold
"Mesh.Smoothing": 10, # Number of smoothing steps (try 0-100)
"Mesh.SmoothNormals": 1, # Smooth normals
"Mesh.SmoothRatio": 1.8, # Smoothing ratio (1.0-3.0)
}
# Override with user-provided parameters
if mesh_params:
default_params.update(mesh_params)
# Apply all parameters (except our custom ones)
for param, value in default_params.items():
if param not in ['size_factor']: # Skip our custom parameters
gmsh.option.setNumber(param, value)
# Set recombination for each surface
for surface in surface_to_region.keys():
gmsh.model.mesh.setRecombine(2, surface)
else:
gmsh.option.setNumber("Mesh.Algorithm", 6) # Frontal-Delaunay for triangles
# Always generate linear elements first - quadratic conversion is done in post-processing
# This avoids gmsh issues with quadratic elements and embedded 1D lines
gmsh.option.setNumber("Mesh.ElementOrder", 1)
# Force mesh coherence before generation
gmsh.option.setNumber("Mesh.ToleranceInitialDelaunay", 1e-12)
# Short edge control is now handled by point sizing during geometry creation
# Feature-aware auto refinement (opt-in). With refine_factor is None this block
# is skipped entirely and no size field is installed, so the generated mesh is
# byte-identical to the historical output. When set, native gmsh size fields
# drive the element size down near features; gmsh composes them with the existing
# point/boundary sizing by taking the minimum, so far-field sizing is unchanged.
if refine_factor is not None:
all_line_curve_tags = [t for info in line_data for t in info['line_tags']]
_apply_feature_refinement(gmsh, target_size, refine_factor, refine_set,
polygon_coords, all_line_curve_tags, point_map,
surface_tags_by_polygon, debug=debug,
region_ids=region_ids, material_k=material_k,
edge_map=edge_map)
# Generate mesh
gmsh.model.mesh.generate(2)
# Remove duplicate nodes again after mesh generation (belt and suspenders)
gmsh.model.mesh.removeDuplicateNodes()
# Get nodes
node_tags, coords, _ = gmsh.model.mesh.getNodes()
nodes = np.array(coords).reshape(-1, 3)[:, :2]
# Create node tag to index mapping
node_tag_to_index = {tag: i for i, tag in enumerate(node_tags)}
elements = []
mat_ids = []
element_node_counts = []
# For quad8: track center nodes to delete later
center_nodes_to_delete = set() if element_type == 'quad8' else None
# Extract elements using physical groups for better region identification
for physical_tag, region_id in physical_surfaces:
try:
# Get entities in this physical group
entities = gmsh.model.getEntitiesForPhysicalGroup(2, physical_tag)
for entity in entities:
# Get all elements for this entity
elem_types, elem_tags_list, node_tags_list = gmsh.model.mesh.getElements(2, entity)
for elem_type, elem_tags, node_tags in zip(elem_types, elem_tags_list, node_tags_list):
# Gmsh element type mapping:
# 2: 3-node triangle, 9: 6-node triangle
# 3: 4-node quadrilateral, 10: 8-node quadrilateral
if elem_type == 2: # 3-node triangle
elements_array = np.array(node_tags).reshape(-1, 3)
for element in elements_array:
idxs = [node_tag_to_index[tag] for tag in element]
# GMSH returns clockwise triangles - reorder to counter-clockwise
idxs[1], idxs[2] = idxs[2], idxs[1]
# Pad to 9 columns with zeros
padded_idxs = idxs + [0] * (9 - len(idxs))
elements.append(padded_idxs)
mat_ids.append(region_id)
element_node_counts.append(3)
elif elem_type == 9: # 6-node triangle
elements_array = np.array(node_tags).reshape(-1, 6)
for element in elements_array:
idxs = [node_tag_to_index[tag] for tag in element]
# GMSH returns clockwise tri6 elements - reorder to counter-clockwise
# Swap corner nodes 1 and 2
idxs[1], idxs[2] = idxs[2], idxs[1]
# Fix midpoint assignments after corner swap 1<->2:
# GMSH gives: n3=edge(0-1), n4=edge(1-2), n5=edge(2-0)
# After swap: n3=edge(0-2), n4=edge(2-1), n5=edge(1-0)
# Standard requires: n3=edge(0-1), n4=edge(1-2), n5=edge(2-0)
# So remap: new_n3=old_n5, new_n4=old_n4, new_n5=old_n3
old_3, old_4, old_5 = idxs[3], idxs[4], idxs[5]
idxs[3] = old_5 # standard edge(0-1) gets GMSH edge(2-0) midpoint
idxs[4] = old_4 # standard edge(1-2) gets GMSH edge(1-2) midpoint
idxs[5] = old_3 # standard edge(2-0) gets GMSH edge(0-1) midpoint
# Pad to 9 columns with zeros
padded_idxs = idxs + [0] * (9 - len(idxs))
elements.append(padded_idxs)
mat_ids.append(region_id)
element_node_counts.append(6)
elif elem_type == 3: # 4-node quadrilateral
elements_array = np.array(node_tags).reshape(-1, 4)
for element in elements_array:
idxs = [node_tag_to_index[tag] for tag in element]
# Fix node ordering for quadrilateral elements
if element_type.startswith('quad'):
idxs = idxs[::-1] # Simple reversal of node order
# Pad to 9 columns with zeros
padded_idxs = idxs + [0] * (9 - len(idxs))
elements.append(padded_idxs)
mat_ids.append(region_id)
element_node_counts.append(4)
elif elem_type == 10: # Quadratic quadrilateral (gmsh generates 9-node Lagrange)
# Gmsh always generates 9-node Lagrange quads for order 2
elements_array = np.array(node_tags).reshape(-1, 9)
for element in elements_array:
idxs = [node_tag_to_index[tag] for tag in element]
if element_type in ['quad8', 'quad9']:
# Both quad8 and quad9 need CW to CCW conversion for first 8 nodes
# Convert from Gmsh CW to CCW ordering for quadrilateral
# Corner nodes: reverse order (0,1,2,3) -> (0,3,2,1)
# Midpoint nodes need to be reordered accordingly:
# GMSH: n4=edge(0-1), n5=edge(1-2), n6=edge(2-3), n7=edge(3-0)
# After corner reversal: need n4=edge(0-3), n5=edge(3-2), n6=edge(2-1), n7=edge(1-0)
# So: new_n4=old_n7, new_n5=old_n6, new_n6=old_n5, new_n7=old_n4
reordered_first8 = [
idxs[0], # corner 0 stays
idxs[3], # corner 1 -> corner 3
idxs[2], # corner 2 stays
idxs[1], # corner 3 -> corner 1
idxs[7], # edge(0-1) -> edge(0-3) = old edge(3-0)
idxs[6], # edge(1-2) -> edge(3-2) = old edge(2-3)
idxs[5], # edge(2-3) -> edge(2-1) = old edge(1-2)
idxs[4] # edge(3-0) -> edge(1-0) = old edge(0-1)
]
if element_type == 'quad8':
# For quad8, skip center node and mark for deletion
center_node_idx = idxs[8] # Mark center node for deletion
center_nodes_to_delete.add(center_node_idx)
padded_idxs = reordered_first8 + [0] # Skip center node, pad to 9
elements.append(padded_idxs)
mat_ids.append(region_id)
element_node_counts.append(8)
else: # quad9
# For quad9, keep center node (9th node unchanged)
full_idxs = reordered_first8 + [idxs[8]] # Add center node
elements.append(full_idxs)
mat_ids.append(region_id)
element_node_counts.append(9)
else:
# This should never happen since element_type is validated earlier
raise ValueError(f"Unexpected element_type '{element_type}' for Gmsh elem_type {elem_type}")
except Exception as e:
print(f"Warning: Could not extract elements for physical group {physical_tag} (region {region_id}): {e}")
continue
# Convert to numpy arrays
elements_array = np.array(elements, dtype=int)
element_types = np.array(element_node_counts, dtype=int)
element_materials = np.array(mat_ids, dtype=int)
# Extract 1D elements from Gmsh-generated 1D mesh along reinforcement lines
elements_1d = []
mat_ids_1d = []
element_node_counts_1d = []
if lines is not None:
# Extract 1D elements from physical groups for each reinforcement line
for physical_tag, line_idx in physical_lines:
try:
# Get entities in this physical group
entities = gmsh.model.getEntitiesForPhysicalGroup(1, physical_tag)
if debug:
print(f" Physical group {physical_tag} (line {line_idx}): found {len(entities)} entities")
for entity in entities:
# Get all 1D elements for this entity
elem_types, elem_tags_list, node_tags_list = gmsh.model.mesh.getElements(1, entity)
for elem_type, elem_tags, node_tags in zip(elem_types, elem_tags_list, node_tags_list):
# Gmsh 1D element type mapping:
# 1: 2-node line (linear), 8: 3-node line (quadratic)
if elem_type == 1: # Linear 1D elements (2 nodes)
elements_array = np.array(node_tags).reshape(-1, 2)
for element in elements_array:
try:
# Convert numpy arrays to regular Python scalars
element_list = element.tolist() # Convert to Python list
if len(element_list) >= 2:
tag1 = int(element_list[0])
tag2 = int(element_list[1])
# Get node indices
idx1 = node_tag_to_index[tag1]
idx2 = node_tag_to_index[tag2]
# Skip zero-length elements (degenerate)
if idx1 == idx2:
continue
coord1 = nodes[idx1]
coord2 = nodes[idx2]
seg_len = ((coord2[0]-coord1[0])**2 + (coord2[1]-coord1[1])**2)**0.5
if seg_len < 1e-6:
continue
# Create 1D element
padded_idxs = [idx1, idx2, 0]
elements_1d.append(padded_idxs)
mat_ids_1d.append(line_idx)
element_node_counts_1d.append(2)
if debug:
coord1 = nodes[idx1]
coord2 = nodes[idx2]
print(f" Created 1D element: {coord1} -> {coord2}")
except (KeyError, TypeError, ValueError, IndexError) as e:
if debug:
print(f" Skipping 1D element due to error: {e}")
continue
elif elem_type == 8: # Quadratic 1D elements (3 nodes)
elements_array = np.array(node_tags).reshape(-1, 3)
for element in elements_array:
try:
# Convert numpy arrays to regular Python scalars
element_list = element.tolist() # Convert to Python list
if len(element_list) >= 3:
tag1 = int(element_list[0])
tag2 = int(element_list[1])
tag3 = int(element_list[2])
# Get node indices
idx1 = node_tag_to_index[tag1]
idx2 = node_tag_to_index[tag2]
idx3 = node_tag_to_index[tag3]
# Skip zero-length elements (degenerate)
coord1 = nodes[idx1]
coord2 = nodes[idx2]
seg_len = ((coord2[0]-coord1[0])**2 + (coord2[1]-coord1[1])**2)**0.5
if seg_len < 1e-6:
continue
# Create 1D element
padded_idxs = [idx1, idx2, idx3]
elements_1d.append(padded_idxs)
mat_ids_1d.append(line_idx)
element_node_counts_1d.append(3)
except (KeyError, TypeError, ValueError, IndexError) as e:
if debug:
print(f" Skipping quadratic 1D element due to error: {e}")
continue
except Exception as e:
if debug:
print(f" Error extracting 1D elements for line {line_idx}: {e}")
continue
gmsh.finalize()
# Clean up center nodes for quad8 elements
if element_type == 'quad8' and center_nodes_to_delete:
print(f"Quad8 cleanup: removing {len(center_nodes_to_delete)} center nodes from {len(nodes)} total nodes")
# c) Create array tracking original node numbering
original_node_count = len(nodes)
nodes_to_keep = [i for i in range(original_node_count) if i not in center_nodes_to_delete]
# d) Delete center nodes - create new nodes array
new_nodes = nodes[nodes_to_keep]
# e) Create mapping from old node indices to new node indices
old_to_new_mapping = {old_idx: new_idx for new_idx, old_idx in enumerate(nodes_to_keep)}
# f) Update element topology to use new node numbering
new_elements = []
for element in elements_array:
new_element = []
for node_idx in element:
if node_idx == 0: # Keep padding zeros
new_element.append(0)
elif node_idx in center_nodes_to_delete:
# This should not happen since we set center nodes to 0
new_element.append(0)
else:
# Map to new node index
new_element.append(old_to_new_mapping[node_idx])
new_elements.append(new_element)
# g) Replace arrays with consolidated versions
elements_array = np.array(new_elements, dtype=int)
nodes = new_nodes
print(f"Quad8 cleanup complete: {len(nodes)} nodes, {len(elements_array)} elements")
# Convert lists to arrays
elements_array = np.array(elements, dtype=int)
element_types = np.array(element_node_counts, dtype=int)
element_materials = np.array(mat_ids, dtype=int) + 1 # Make 1-based
# Normalize winding: downstream assembly routines (e.g. the linear-triangle
# seepage stiffness) assume CCW elements with positive signed area. Gmsh
# inherits the orientation of the input polygon ring, so polygon coordinates
# entered clockwise would otherwise propagate CW elements and silently drop
# out of the assembly.
n_flipped = ensure_ccw_elements(nodes, elements_array, element_types)
if debug and n_flipped:
print(f"Flipped {n_flipped} clockwise elements to CCW winding")
mesh = {
"nodes": nodes,
"elements": elements_array,
"element_types": element_types,
"element_materials": element_materials,
}
# Add 1D element data if lines were provided
if lines is not None and len(elements_1d) > 0:
elements_1d_array = np.array(elements_1d, dtype=int)
element_types_1d = np.array(element_node_counts_1d, dtype=int)
element_materials_1d = np.array(mat_ids_1d, dtype=int) + 1 # Make 1-based
mesh["elements_1d"] = elements_1d_array
mesh["element_types_1d"] = element_types_1d
mesh["element_materials_1d"] = element_materials_1d
# Conforming fallback for constraint lines the geo-kernel embed could not recover.
# Long inclined reinforcement rooted on / crossing a boundary or material edge —
# VP60's soil nails root on the vertical wall face — leaves 1D nodes floating
# inside 2D triangles (orphans -> singular stiffness) or produces no elements at
# all. Detect that and rebuild via the OCC kernel's boolean fragment, which splits
# the surfaces along the lines so they become real, conforming mesh edges. Only
# triangular meshes are rebuilt (the reinforcement corpus is tri3/tri6); this is
# gated on FAILURE, so a mesh that already conformed is byte-identical (it never
# enters this branch).
if (lines is not None and base_element_type == 'tri3'
and (len(mesh.get("elements", [])) == 0 or _has_orphan_1d_nodes(mesh))):
if debug:
print("Primary embed left orphan 1D nodes / no elements — rebuilding "
"with OCC boolean fragment (conforming fallback)")
mesh = _remesh_with_occ_fragment(polygon_coords, region_ids, lines,
target_size, base_element_type, debug=debug)
# Post-process to convert linear elements to quadratic if requested
if quadratic:
if debug:
print(f"Converting linear {base_element_type} mesh to quadratic {element_type}")
mesh = convert_linear_to_quadratic_mesh(mesh, element_type, debug=debug)
# An empty mesh means gmsh failed to recover the geometry (e.g. a constraint
# line without matching polygon vertices). Fail loudly here — downstream the
# symptom is a cryptic singular-matrix error in the solver.
if len(mesh.get("elements", [])) == 0:
raise ValueError(
"Mesh generation produced 0 elements — gmsh could not recover the "
"geometry. Check that material zones tile the section without gaps "
"and that reinforcement/pile lines have their intersection points "
"inserted into the zone boundaries (get_material_polygons(slope_data, "
"reinf_lines=...) does this automatically).")
return mesh
build_polygons(slope_data, reinf_lines=None, tol=1e-06, debug=False)
Build material zone polygons from slope_data.
Extracts profile lines and max depth, then creates polygons for each material zone. Also integrates distributed load points and reinforcement line endpoints that are coincident with polygon edges.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def build_polygons(slope_data, reinf_lines=None, tol = 0.000001, debug=False):
"""
Build material zone polygons from slope_data.
Extracts profile lines and max depth, then creates polygons for each material zone.
Also integrates distributed load points and reinforcement line endpoints that are
coincident with polygon edges.
Parameters:
slope_data: Dictionary containing slope geometry data
Returns:
List of polygons as dicts with keys:
"coords": list of (x, y) coordinate tuples
"mat_id": optional material ID (0-based) or None
"""
import numpy as np
import copy
# Extract profile lines and max depth from slope_data
profile_lines = slope_data.get('profile_lines', [])
max_depth = slope_data.get('max_depth', None)
if not profile_lines:
raise ValueError("Need at least 1 profile line to create material zones")
# For single profile line, max_depth serves as the bottom boundary
if len(profile_lines) == 1:
if max_depth is None:
raise ValueError("When using only 1 profile line, max_depth must be specified")
n = len(profile_lines)
lines = [list(line['coords']) for line in copy.deepcopy(profile_lines)]
for i in range(n - 1):
top = lines[i]
for endpoint in [0, -1]: # left and right
x_top, y_top = top[endpoint]
# Find the highest lower profile at this x
best_j = None
best_y = -np.inf
for j in range(i + 1, n):
lower = lines[j]
xs_lower = np.array([x for x, y in lower])
ys_lower = np.array([y for x, y in lower])
if xs_lower[0] - tol <= x_top <= xs_lower[-1] + tol:
y_proj = np.interp(x_top, xs_lower, ys_lower)
if y_proj > best_y:
best_y = y_proj
best_j = j
if best_j is not None:
lower = lines[best_j]
xs_lower = np.array([x for x, y in lower])
ys_lower = np.array([y for x, y in lower])
y_proj = np.interp(x_top, xs_lower, ys_lower)
# Check if lower profile already has a point at this x (within tol)
found = False
for (x_l, y_l) in lower:
if abs(x_l - x_top) < tol:
found = True
break
if abs(y_proj - y_top) < tol:
# Coincident: insert (x_top, y_top) if not present
if not found:
insert_idx = np.searchsorted(xs_lower, x_top)
lower.insert(insert_idx, (round(x_top, 6), round(y_top, 6)))
else:
# Not coincident: insert (x_top, y_proj) if not present
if not found:
insert_idx = np.searchsorted(xs_lower, x_top)
lower.insert(insert_idx, (round(x_top, 6), round(y_proj, 6)))
def clean_polygon(poly, tol=1e-8):
# Remove consecutive duplicate points (except for closing point)
if not poly:
return poly
cleaned = [poly[0]]
for pt in poly[1:]:
if abs(pt[0] - cleaned[-1][0]) > tol or abs(pt[1] - cleaned[-1][1]) > tol:
cleaned.append(pt)
# Ensure closed
if abs(cleaned[0][0] - cleaned[-1][0]) > tol or abs(cleaned[0][1] - cleaned[-1][1]) > tol:
cleaned.append(cleaned[0])
return cleaned
# Now build polygons as before
polygons = []
for i, top_line in enumerate(lines):
xs_top, ys_top = zip(*top_line)
xs_top = np.array(xs_top)
ys_top = np.array(ys_top)
left_x, left_y = xs_top[0], ys_top[0]
right_x, right_y = xs_top[-1], ys_top[-1]
# Initialize variables for debug output
lower_left_x = None
lower_right_x = None
proj_left_x = None
proj_right_x = None
bottom_cleaned = []
# Initialize vertical edge points (used for intermediate points on vertical edges)
left_vertical_points = [] # Intermediate points on left vertical edge (bottom to top)
right_vertical_points = [] # Intermediate points on right vertical edge (top to bottom)
left_y_bot = -np.inf
right_y_bot = -np.inf
if i < n - 1:
# Use the immediate next line as the lower boundary
lower_line = lines[i + 1]
xs_bot, ys_bot = zip(*lower_line)
xs_bot = np.array(xs_bot)
ys_bot = np.array(ys_bot)
lower_left_x = xs_bot[0]
lower_right_x = xs_bot[-1]
# Collect actual points from all lower lines within the top line's x-range
# But only include a point if it's actually on the highest lower profile at that x
bottom_points = [] # List of (x, y, line_idx) tuples
for j in range(i + 1, n):
lower_candidate = lines[j]
xs_cand = np.array([x for x, y in lower_candidate])
ys_cand = np.array([y for x, y in lower_candidate])
# Only include points that are within the top line's x-range
mask = (xs_cand >= left_x - tol) & (xs_cand <= right_x + tol)
for x, y in zip(xs_cand[mask], ys_cand[mask]):
# Check if this point is actually on the highest lower profile at this x
# Compare with all other lower lines at this x-coordinate
is_highest = True
for k in range(i + 1, n):
if k == j:
continue
other_line = lines[k]
xs_other = np.array([x_o for x_o, y_o in other_line])
ys_other = np.array([y_o for x_o, y_o in other_line])
if xs_other[0] - tol <= x <= xs_other[-1] + tol:
y_other = np.interp(x, xs_other, ys_other)
if y_other > y + tol: # Other line is higher
is_highest = False
break
if is_highest:
bottom_points.append((x, y, j))
# Build bottom_cleaned: ordered path along the bottom boundary
# Check if all bottom points come from a single lower profile
line_indices_set = set(line_idx for _, _, line_idx in bottom_points)
if len(line_indices_set) == 1:
# Single lower profile: extract ordered sub-path directly from the profile.
# This preserves vertical segments (multiple y-values at same x) that would
# be lost by the x-based grouping approach.
j = next(iter(line_indices_set))
lower_path = lines[j]
bottom_cleaned = [
(x, y) for x, y in lower_path
if left_x - tol <= x <= right_x + tol
]
else:
# Multiple lower profiles: group by x-coordinate and keep only the highest y
# at each x. This handles cases where multiple lines have points at the same x.
bottom_dict = {} # x_key -> (y, line_idx, orig_x, orig_y)
for x, y, line_idx in bottom_points:
x_key = round(x / tol) * tol # Round to tolerance to group nearby points
if x_key not in bottom_dict or y > bottom_dict[x_key][0]:
bottom_dict[x_key] = (y, line_idx, x, y)
# Convert to sorted list
bottom_cleaned = sorted([(orig_x, orig_y) for _, _, orig_x, orig_y in bottom_dict.values()])
# Helper function to check if a point already exists in a list
def point_exists(point_list, x, y, tol=1e-8):
"""Check if a point (x, y) already exists in the point list within tolerance."""
for px, py in point_list:
if abs(px - x) < tol and abs(py - y) < tol:
return True
return False
# Helper function to find the lowest y value at a given x by checking all segments
def find_lowest_y_at_x(line_points, x_query, tol=1e-8):
"""
Find the lowest y value at x_query by checking all segments of the line.
Handles vertical segments properly by finding all y values at that x and returning the minimum.
Returns:
tuple: (y_value, is_at_endpoint) where is_at_endpoint indicates if x_query is at an endpoint
"""
if not line_points:
return None, False
xs = np.array([x for x, y in line_points])
ys = np.array([y for x, y in line_points])
# Check if x_query is within the line's x-range
if xs[0] - tol > x_query or xs[-1] + tol < x_query:
return None, False
# Check if x_query is at an endpoint
is_at_left_endpoint = abs(x_query - xs[0]) < tol
is_at_right_endpoint = abs(x_query - xs[-1]) < tol
is_at_endpoint = is_at_left_endpoint or is_at_right_endpoint
# Find all y values at x_query by checking all segments
y_values = []
# Check all points that are exactly at x_query
for k in range(len(line_points)):
if abs(xs[k] - x_query) < tol:
y_values.append(ys[k])
# Check all segments that contain x_query
for k in range(len(line_points) - 1):
x1, y1 = line_points[k]
x2, y2 = line_points[k + 1]
# Check if segment is vertical and contains x_query
if abs(x1 - x_query) < tol and abs(x2 - x_query) < tol:
# Vertical segment - include both y values
y_values.append(y1)
y_values.append(y2)
# Check if segment is horizontal or sloped and contains x_query
elif min(x1, x2) - tol <= x_query <= max(x1, x2) + tol:
# Interpolate y value
if abs(x2 - x1) < tol:
# Segment is vertical (should have been caught above, but just in case)
y_values.append(y1)
y_values.append(y2)
else:
# Linear interpolation
t = (x_query - x1) / (x2 - x1)
if 0 <= t <= 1:
y_interp = y1 + t * (y2 - y1)
y_values.append(y_interp)
if not y_values:
return None, False
# Return the lowest y value
y_min = min(y_values)
return y_min, is_at_endpoint
def find_projected_y_at_x(line_points, x_query, y_ref, side, tol=1e-8):
"""
For vertical endpoint projections: choose the intersection y at x_query that is
closest *below* the point we're projecting from.
This fixes the case where a candidate profile has a vertical segment at x_query
(e.g., (260,229) then (260,202)). In that situation, using the "lowest y" (202)
is wrong; we want the first hit when projecting downward (229).
Behavior is intentionally conservative:
- If there is at least one intersection strictly below y_ref, return the highest of those.
- Otherwise fall back to the original behavior (lowest y), preserving legacy behavior
in edge cases (e.g., coincident/above intersections).
"""
# Reuse the exact same intersection enumeration logic as find_lowest_y_at_x,
# but keep the full set of y-values.
if not line_points:
return None, False
xs = np.array([x for x, y in line_points])
ys = np.array([y for x, y in line_points])
if xs[0] - tol > x_query or xs[-1] + tol < x_query:
return None, False
is_at_left_endpoint = abs(x_query - xs[0]) < tol
is_at_right_endpoint = abs(x_query - xs[-1]) < tol
is_at_endpoint = is_at_left_endpoint or is_at_right_endpoint
y_values = []
for k in range(len(line_points)):
if abs(xs[k] - x_query) < tol:
y_values.append(float(ys[k]))
for k in range(len(line_points) - 1):
x1, y1 = line_points[k]
x2, y2 = line_points[k + 1]
if abs(x1 - x_query) < tol and abs(x2 - x_query) < tol:
y_values.append(float(y1))
y_values.append(float(y2))
elif min(x1, x2) - tol <= x_query <= max(x1, x2) + tol:
if abs(x2 - x1) < tol:
y_values.append(float(y1))
y_values.append(float(y2))
else:
t = (x_query - x1) / (x2 - x1)
if 0 <= t <= 1:
y_values.append(float(y1 + t * (y2 - y1)))
if not y_values:
return None, False
# If the polyline has multiple *vertices* exactly at this x (vertical segment / duplicate-x),
# use a deterministic selection based on which side we are projecting from:
# - projecting from LEFT endpoint of the upper line: keep the LAST y encountered
# - projecting from RIGHT endpoint of the upper line: keep the FIRST y encountered
#
# This matches the intended "walk along the lower boundary" behavior and fixes cases like:
# - right projection at x=260 with vertices (260,229) then (260,202): choose 229 (first)
# - left projection at x=240 with vertices (240,140) then (240,190): choose 190 (last)
vertex_y_at_x = [float(y) for (x, y) in line_points if abs(x - x_query) < tol]
if len(vertex_y_at_x) >= 2:
if side == "right":
# first encountered vertex at this x
y_pick = vertex_y_at_x[0]
# If we are exactly on a vertex at y_ref, that is the first hit.
if abs(y_pick - y_ref) < tol:
return float(y_ref), is_at_endpoint
if y_pick < (y_ref - tol):
return y_pick, is_at_endpoint
elif side == "left":
# last encountered vertex at this x
y_pick = vertex_y_at_x[-1]
# If we are exactly on a vertex at y_ref, that is the first hit.
if abs(y_pick - y_ref) < tol:
return float(y_ref), is_at_endpoint
if y_pick < (y_ref - tol):
return y_pick, is_at_endpoint
y_below = [y for y in y_values if y < (y_ref - tol)]
if y_below:
return max(y_below), is_at_endpoint
# Fall back to legacy behavior
return min(y_values), is_at_endpoint
# Project endpoints - find highest lower profile or use max_depth
# When projecting right side: if intersection is at left end of lower line,
# add that point but continue projecting down
# When projecting left side: if intersection is at right end of lower line,
# add that point but continue projecting down
for j in range(i + 1, n):
lower_candidate = lines[j]
xs_cand = np.array([x for x, y in lower_candidate])
ys_cand = np.array([y for x, y in lower_candidate])
# Check left endpoint projection
if xs_cand[0] - tol <= left_x <= xs_cand[-1] + tol:
y_cand, is_at_endpoint = find_projected_y_at_x(lower_candidate, left_x, left_y, side="left", tol=tol)
if y_cand is not None:
# If intersection is at the right end of the lower line, add point but continue
if is_at_endpoint and abs(left_x - xs_cand[-1]) < tol: # At right endpoint
# Only add if not duplicate of the endpoint being projected and not already in list
if abs(y_cand - left_y) > tol and not point_exists(left_vertical_points, left_x, y_cand, tol):
left_vertical_points.append((left_x, y_cand))
else: # Not at endpoint, use as stopping point
if y_cand > left_y_bot:
left_y_bot = y_cand
# Check right endpoint projection
if xs_cand[0] - tol <= right_x <= xs_cand[-1] + tol:
y_cand, is_at_endpoint = find_projected_y_at_x(lower_candidate, right_x, right_y, side="right", tol=tol)
if y_cand is not None:
# If intersection is at the left end of the lower line, add point but continue
if is_at_endpoint and abs(right_x - xs_cand[0]) < tol: # At left endpoint
# Only add if not duplicate of the endpoint being projected and not already in list
if abs(y_cand - right_y) > tol and not point_exists(right_vertical_points, right_x, y_cand, tol):
right_vertical_points.append((right_x, y_cand))
else: # Not at endpoint, use as stopping point
if y_cand > right_y_bot:
right_y_bot = y_cand
# If no lower profile at endpoints, use max_depth
if left_y_bot == -np.inf:
left_y_bot = max_depth if max_depth is not None else -np.inf
if right_y_bot == -np.inf:
right_y_bot = max_depth if max_depth is not None else -np.inf
# Filter vertical-edge "continue projecting" points so we only keep points that
# actually lie on the final vertical edge between the top and bottom of this zone.
#
# Without this, a deeper left-endpoint intersection (e.g., (240,190) at the left
# endpoint of some deeper line) can be appended to right_vertical_points even after
# we've already found the correct bottom (e.g., right_y_bot=229). That creates the
# dangling vertical segment you observed.
if right_y_bot != -np.inf:
right_vertical_points = [
(x, y) for (x, y) in right_vertical_points
if (y < right_y - tol) and (y > right_y_bot + tol)
]
if left_y_bot != -np.inf:
# Left edge runs from bottom up to top; keep points strictly between bottom and top.
left_vertical_points = [
(x, y) for (x, y) in left_vertical_points
if (y > left_y_bot + tol) and (y < left_y - tol)
]
# Deduplicate vertical points (remove points that are too close to each other)
def deduplicate_points(points, tol=1e-8):
"""Remove duplicate points within tolerance."""
if not points:
return []
unique_points = [points[0]]
for p in points[1:]:
# Check if this point is too close to any existing unique point
is_duplicate = False
for up in unique_points:
if abs(p[0] - up[0]) < tol and abs(p[1] - up[1]) < tol:
is_duplicate = True
break
if not is_duplicate:
unique_points.append(p)
return unique_points
right_vertical_points = deduplicate_points(right_vertical_points, tol)
left_vertical_points = deduplicate_points(left_vertical_points, tol)
# Sort vertical points: right edge top to bottom, left edge bottom to top
right_vertical_points.sort(key=lambda p: -p[1]) # Sort by y descending (top to bottom)
left_vertical_points.sort(key=lambda p: p[1]) # Sort by y ascending (bottom to top)
# Build bottom boundary: right projection, intermediate points (right to left), left projection
# The bottom should go from right to left to close the polygon
bottom = []
# Start with right endpoint
if right_y_bot != -np.inf:
bottom.append((right_x, right_y_bot))
# Add intermediate points in reverse order (right to left)
# Filter out points too close to endpoints
for x, y in reversed(bottom_cleaned):
if abs(x - left_x) > tol and abs(x - right_x) > tol:
bottom.append((x, y))
# End with left endpoint
if left_y_bot != -np.inf:
bottom.append((left_x, left_y_bot))
# Store for debug output
proj_left_x = left_x
proj_right_x = right_x
else:
# For the lowest polygon, bottom is at max_depth
# Only need endpoints - no intermediate points
left_y_bot = max_depth if max_depth is not None else -np.inf
right_y_bot = max_depth if max_depth is not None else -np.inf
bottom = []
bottom.append((right_x, max_depth))
bottom.append((left_x, max_depth))
# Build polygon: top left-to-right, right vertical edge (with intermediate points),
# bottom right-to-left, left vertical edge (with intermediate points)
poly = []
# Top edge: left to right along profile line
for x, y in zip(xs_top, ys_top):
poly.append((round(x, 6), round(y, 6)))
# Right vertical edge: from (right_x, right_y) down to (right_x, right_y_bot)
# Include intermediate points where we intersect left endpoints of lower lines
# Note: (right_x, right_y_bot) will be added as part of the bottom edge, so don't add it here
if i < n - 1:
for x, y in right_vertical_points:
# Only add if it's between top and bottom (not duplicate of endpoints)
if abs(y - right_y) > tol and abs(y - right_y_bot) > tol:
poly.append((round(x, 6), round(y, 6)))
# Bottom edge: right to left (already includes (right_x, right_y_bot) and (left_x, left_y_bot))
for x, y in bottom:
poly.append((round(x, 6), round(y, 6)))
# Left vertical edge: from (left_x, left_y_bot) up to (left_x, left_y)
# Include intermediate points where we intersect right endpoints of lower lines
# Note: (left_x, left_y_bot) was already added as part of the bottom edge
if i < n - 1:
for x, y in reversed(left_vertical_points): # Reverse to go bottom to top
# Only add if it's between bottom and top (not duplicate of endpoints)
if abs(y - left_y_bot) > tol and abs(y - left_y) > tol:
poly.append((round(x, 6), round(y, 6)))
# Clean up polygon (should rarely do anything)
poly = clean_polygon(poly)
mat_id = profile_lines[i].get("mat_id") if i < len(profile_lines) else None
polygons.append({
"coords": poly,
"mat_id": mat_id
})
# Distributed-load endpoints are inserted as polygon vertices by
# get_material_polygons() (via add_dload_points_to_polygons) — not here, since
# build_polygons() also produces LEM slicing geometry, which doesn't need them.
# That keeps element edges aligned to the loaded span for the FEM load assembly.
# Add intersection points with reinforcement lines if provided
if reinf_lines is not None:
polygons = add_intersection_points_to_polygons(polygons, reinf_lines, debug=debug)
return polygons
convert_linear_to_quadratic_mesh(mesh, target_element_type, debug=False)
Convert a linear mesh (tri3/quad4) to quadratic (tri6/quad8/quad9) by adding midside nodes.
This is much more robust than gmsh's built-in quadratic generation, especially when dealing with embedded 1D elements (reinforcement lines).
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def convert_linear_to_quadratic_mesh(mesh, target_element_type, debug=False):
"""
Convert a linear mesh (tri3/quad4) to quadratic (tri6/quad8/quad9) by adding midside nodes.
This is much more robust than gmsh's built-in quadratic generation, especially
when dealing with embedded 1D elements (reinforcement lines).
Parameters:
mesh: Dictionary containing linear mesh data
target_element_type: 'tri6', 'quad8', or 'quad9'
debug: Enable debug output
Returns:
Updated mesh dictionary with quadratic elements
"""
if debug:
print(f"Converting to {target_element_type} elements...")
nodes = mesh["nodes"].copy()
elements = mesh["elements"].copy()
element_types = mesh["element_types"].copy()
element_materials = mesh["element_materials"].copy()
# Handle 1D elements if present
elements_1d = mesh.get("elements_1d")
element_types_1d = mesh.get("element_types_1d")
element_materials_1d = mesh.get("element_materials_1d")
has_1d_elements = elements_1d is not None
if has_1d_elements:
elements_1d = elements_1d.copy()
element_types_1d = element_types_1d.copy()
element_materials_1d = element_materials_1d.copy()
# Dictionary to store midside nodes: (node1_idx, node2_idx) -> midside_node_idx
# Always store with node1_idx < node2_idx for consistency
midside_nodes = {}
next_node_idx = len(nodes)
new_node_coords = [] # collect new coords, append to nodes at the end
def get_or_create_midside_node(n1_idx, n2_idx):
"""Get existing midside node or create new one between n1 and n2"""
nonlocal next_node_idx
# Ensure consistent ordering
if n1_idx > n2_idx:
n1_idx, n2_idx = n2_idx, n1_idx
edge_key = (n1_idx, n2_idx)
if edge_key in midside_nodes:
return midside_nodes[edge_key]
# Create new midside node at edge center
midside_coord = (nodes[n1_idx] + nodes[n2_idx]) / 2.0
new_node_coords.append(midside_coord)
midside_idx = next_node_idx
midside_nodes[edge_key] = midside_idx
next_node_idx += 1
if debug and len(midside_nodes) <= 10: # Only print first few
print(f" Created midside node {midside_idx} between {n1_idx}-{n2_idx} at {midside_coord}")
return midside_idx
# Convert 2D elements
new_elements = []
new_element_types = []
for elem_idx, element in enumerate(elements):
elem_type = element_types[elem_idx]
if target_element_type == 'tri6' and elem_type == 3:
# Convert tri3 to tri6
n0, n1, n2 = element[0], element[1], element[2]
# Get/create midside nodes
n3 = get_or_create_midside_node(n0, n1) # edge 0-1
n4 = get_or_create_midside_node(n1, n2) # edge 1-2
n5 = get_or_create_midside_node(n2, n0) # edge 2-0
# tri6 node ordering: [corner_nodes, midside_nodes]
new_element = [n0, n1, n2, n3, n4, n5, 0, 0, 0]
new_elements.append(new_element)
new_element_types.append(6)
elif target_element_type == 'quad8' and elem_type == 4:
# Convert quad4 to quad8
n0, n1, n2, n3 = element[0], element[1], element[2], element[3]
# Get/create midside nodes on edges
n4 = get_or_create_midside_node(n0, n1) # edge 0-1
n5 = get_or_create_midside_node(n1, n2) # edge 1-2
n6 = get_or_create_midside_node(n2, n3) # edge 2-3
n7 = get_or_create_midside_node(n3, n0) # edge 3-0
# quad8 node ordering: [corner_nodes, midside_nodes]
new_element = [n0, n1, n2, n3, n4, n5, n6, n7, 0]
new_elements.append(new_element)
new_element_types.append(8)
elif target_element_type == 'quad9' and elem_type == 4:
# Convert quad4 to quad9
n0, n1, n2, n3 = element[0], element[1], element[2], element[3]
# Get/create midside nodes on edges
n4 = get_or_create_midside_node(n0, n1) # edge 0-1
n5 = get_or_create_midside_node(n1, n2) # edge 1-2
n6 = get_or_create_midside_node(n2, n3) # edge 2-3
n7 = get_or_create_midside_node(n3, n0) # edge 3-0
# Create center node (append to new_node_coords like midside nodes)
center_coord = (nodes[n0] + nodes[n1] + nodes[n2] + nodes[n3]) / 4.0
new_node_coords.append(center_coord)
n8 = next_node_idx
next_node_idx += 1
# quad9 node ordering: [corner_nodes, midside_nodes, center_node]
new_element = [n0, n1, n2, n3, n4, n5, n6, n7, n8]
new_elements.append(new_element)
new_element_types.append(9)
elif elem_type == 3 and target_element_type in ['quad8', 'quad9']:
# Convert leftover tri3 to tri6 in quad-dominant meshes
n0, n1, n2 = element[0], element[1], element[2]
n3 = get_or_create_midside_node(n0, n1)
n4 = get_or_create_midside_node(n1, n2)
n5 = get_or_create_midside_node(n2, n0)
new_element = [n0, n1, n2, n3, n4, n5, 0, 0, 0]
new_elements.append(new_element)
new_element_types.append(6)
else:
# Keep original element unchanged
new_elements.append(element.tolist())
new_element_types.append(elem_type)
# Keep 1D elements as linear (2-node). Truss stiffness uses only end nodes,
# so midside nodes add no physical fidelity and can cause singular K if they
# are not shared with a 2D element edge.
new_elements_1d = []
new_element_types_1d = []
if has_1d_elements:
for elem_idx, element in enumerate(elements_1d):
new_elements_1d.append(element.tolist())
new_element_types_1d.append(element_types_1d[elem_idx])
# Append all new midside node coordinates at once
if new_node_coords:
nodes = np.vstack([nodes, np.array(new_node_coords)])
if debug:
print(f" Added {len(midside_nodes)} midside nodes")
print(f" Total nodes: {len(nodes)} (was {len(mesh['nodes'])})")
# Create updated mesh
updated_mesh = {
"nodes": nodes,
"elements": np.array(new_elements, dtype=int),
"element_types": np.array(new_element_types, dtype=int),
"element_materials": element_materials
}
if has_1d_elements:
updated_mesh["elements_1d"] = np.array(new_elements_1d, dtype=int)
updated_mesh["element_types_1d"] = np.array(new_element_types_1d, dtype=int)
updated_mesh["element_materials_1d"] = element_materials_1d
return updated_mesh
detect_crack_tips(polygon_coords, angle_thresh_deg=_REFINE_CRACK_ANGLE_DEG)
Detect V-notch / crack tips in the material polygons.
A boundary vertex whose two incident edges meet at an interior angle sharper
than angle_thresh_deg forms a thin re-entrant spike — exactly the wall-crack
idiom the seepage sheet-pile / clay-blanket samples build (a ground surface that
dips straight to the wall tip and back up, leaving a narrow slit; the tip is the
deepest vertex of that notch). Returns a de-duplicated, sorted list of (x, y)
tip coordinates so downstream field construction is deterministic.
Source code in xslope/mesh.py
def detect_crack_tips(polygon_coords, angle_thresh_deg=_REFINE_CRACK_ANGLE_DEG):
"""Detect V-notch / crack tips in the material polygons.
A boundary vertex whose two incident edges meet at an interior angle sharper
than ``angle_thresh_deg`` forms a thin re-entrant spike — exactly the wall-crack
idiom the seepage sheet-pile / clay-blanket samples build (a ground surface that
dips straight to the wall tip and back up, leaving a narrow slit; the tip is the
deepest vertex of that notch). Returns a de-duplicated, sorted list of ``(x, y)``
tip coordinates so downstream field construction is deterministic.
"""
import math
cos_thresh = math.cos(math.radians(angle_thresh_deg))
tips = {}
for poly in polygon_coords:
pts = remove_duplicate_endpoint(list(poly))
n = len(pts)
if n < 3:
continue
for i in range(n):
a = pts[(i - 1) % n]
v = pts[i]
b = pts[(i + 1) % n]
ux, uy = a[0] - v[0], a[1] - v[1]
wx, wy = b[0] - v[0], b[1] - v[1]
lu = math.hypot(ux, uy)
lw = math.hypot(wx, wy)
if lu < 1e-12 or lw < 1e-12:
continue
# cos(angle) near 1 => the two edges point nearly the same way from v,
# i.e. a narrow spike (a crack/notch tip), not an ordinary corner.
cos_ang = (ux * wx + uy * wy) / (lu * lw)
if cos_ang > cos_thresh:
tips[(round(v[0], 9), round(v[1], 9))] = (v[0], v[1])
return [tips[k] for k in sorted(tips.keys())]
detect_interface_edges(polygon_coords, region_ids, material_k, contrast_thresh=_REFINE_INTERFACE_CONTRAST)
Detect material-boundary edges across which the hydraulic conductivity contrast
is at least contrast_thresh (default 100x) — the seepage-critical interfaces
where a sharp permeability jump (a clay core in a pervious shell, an impermeable
cutoff) forces a steep head gradient the mesh must resolve to converge.
Pure geometry + adjacency: a segment shared by two material polygons is a boundary;
its k1 ratio (max/min of the two sides' material_k) decides whether it is
high-contrast. material_k maps region_id (the value in region_ids, one
per polygon) to k1 (a dict) or is a sequence indexed by region_id. Edges whose
either side has no k1 entry, or a non-positive k1, are skipped.
Returns a de-duplicated, sorted list of ((x1,y1),(x2,y2)) coordinate-pair edges
with each endpoint rounded to 9 decimals and the pair in canonical (sorted) order,
so downstream field construction is deterministic.
Source code in xslope/mesh.py
def detect_interface_edges(polygon_coords, region_ids, material_k,
contrast_thresh=_REFINE_INTERFACE_CONTRAST):
"""Detect material-boundary edges across which the hydraulic conductivity contrast
is at least ``contrast_thresh`` (default 100x) — the seepage-critical interfaces
where a sharp permeability jump (a clay core in a pervious shell, an impermeable
cutoff) forces a steep head gradient the mesh must resolve to converge.
Pure geometry + adjacency: a segment shared by two material polygons is a boundary;
its k1 ratio (``max/min`` of the two sides' ``material_k``) decides whether it is
high-contrast. ``material_k`` maps ``region_id`` (the value in ``region_ids``, one
per polygon) to k1 (a dict) or is a sequence indexed by region_id. Edges whose
either side has no k1 entry, or a non-positive k1, are skipped.
Returns a de-duplicated, sorted list of ``((x1,y1),(x2,y2))`` coordinate-pair edges
with each endpoint rounded to 9 decimals and the pair in canonical (sorted) order,
so downstream field construction is deterministic.
"""
if material_k is None or not region_ids:
return []
edge_regions = {}
for idx, poly in enumerate(polygon_coords):
if idx >= len(region_ids):
break
rid = region_ids[idx]
pts = remove_duplicate_endpoint(list(poly))
n = len(pts)
if n < 3:
continue
for i in range(n):
a = pts[i]
b = pts[(i + 1) % n]
ka = (round(a[0], 9), round(a[1], 9))
kb = (round(b[0], 9), round(b[1], 9))
if ka == kb:
continue
key = (ka, kb) if ka <= kb else (kb, ka)
edge_regions.setdefault(key, set()).add(rid)
out = []
for key, regions in edge_regions.items():
if len(regions) != 2:
continue # unshared boundary or (degenerate) 3+ regions on an edge
ks = [_k_of(material_k, r) for r in regions]
if any(k is None or k <= 0 for k in ks):
continue
if max(ks) / min(ks) >= contrast_thresh:
out.append(key)
out.sort()
return out
detect_thin_zones(polygon_coords, target_size, min_elems=_REFINE_THIN_MIN_ELEMS)
Detect thin material zones — pinches and slender (possibly inclined) bands whose
local width is too small to fit min_elems elements at target_size.
Erosion test: a material polygon that vanishes when eroded by min_elems*target_size/2
is thin EVERYWHERE (its local width < min_elems*target_size) — a soft band is the
canonical case. Because such a band is often inclined, its axis-aligned bounding box is
a poor proxy; it is flagged as a 'whole' zone so the caller can restrict a size
field to the polygon's own surface and follow its shape exactly.
For polygons that survive erosion, a morphological opening (erode then dilate back)
isolates any LOCAL pinch — the residue polygon - opening keeps only the parts
thinner than min_elems*target_size. A pinch is emitted as a 'box' zone ONLY
when it is compact (both bbox dimensions within a few element sizes); spread-out
residues of otherwise-thick polygons are corner artifacts and are skipped, since an
axis-aligned box over them would explode the node count without resolving a real
feature.
Returns a deterministic list of dicts. 'whole' zones:
{kind:'whole', poly_index:i, width, size}. 'box' zones:
{kind:'box', bbox:(xmin,ymin,xmax,ymax), width, size}. size = width/min_elems.
Source code in xslope/mesh.py
def detect_thin_zones(polygon_coords, target_size, min_elems=_REFINE_THIN_MIN_ELEMS):
"""Detect thin material zones — pinches and slender (possibly inclined) bands whose
local width is too small to fit ``min_elems`` elements at ``target_size``.
Erosion test: a material polygon that vanishes when eroded by ``min_elems*target_size/2``
is thin EVERYWHERE (its local width < ``min_elems*target_size``) — a soft band is the
canonical case. Because such a band is often inclined, its axis-aligned bounding box is
a poor proxy; it is flagged as a ``'whole'`` zone so the caller can restrict a size
field to the polygon's own surface and follow its shape exactly.
For polygons that survive erosion, a morphological opening (erode then dilate back)
isolates any LOCAL pinch — the residue ``polygon - opening`` keeps only the parts
thinner than ``min_elems*target_size``. A pinch is emitted as a ``'box'`` zone ONLY
when it is compact (both bbox dimensions within a few element sizes); spread-out
residues of otherwise-thick polygons are corner artifacts and are skipped, since an
axis-aligned box over them would explode the node count without resolving a real
feature.
Returns a deterministic list of dicts. ``'whole'`` zones:
``{kind:'whole', poly_index:i, width, size}``. ``'box'`` zones:
``{kind:'box', bbox:(xmin,ymin,xmax,ymax), width, size}``. ``size = width/min_elems``.
"""
from shapely.geometry import Polygon, MultiPolygon, GeometryCollection
w_half = min_elems * target_size / 2.0
min_area = target_size ** 2 # ignore residues smaller than one element
max_box = _REFINE_BAND_ELEMS * min_elems * target_size # compactness limit for a pinch box
whole, boxes = [], []
for idx, poly in enumerate(polygon_coords):
pts = remove_duplicate_endpoint(list(poly))
if len(pts) < 3:
continue
try:
P = Polygon(pts)
if not P.is_valid:
P = P.buffer(0)
if P.is_empty or P.area <= 0 or P.length <= 0:
continue
eroded = P.buffer(-w_half)
except Exception:
continue
if eroded.is_empty:
# Thin everywhere -> whole-polygon zone (handles inclined bands exactly).
width = 2.0 * P.area / P.length
whole.append({'kind': 'whole', 'poly_index': idx,
'width': width, 'size': width / min_elems})
continue
# Survived erosion: isolate any compact local pinch.
try:
residue = P.difference(eroded.buffer(+w_half))
except Exception:
continue
if residue.is_empty:
continue
geoms = residue.geoms if isinstance(residue, (MultiPolygon, GeometryCollection)) else [residue]
for g in geoms:
area = getattr(g, 'area', 0.0)
length = getattr(g, 'length', 0.0)
if area < min_area or length < 1e-9:
continue
xmin, ymin, xmax, ymax = g.bounds
if (xmax - xmin) > max_box or (ymax - ymin) > max_box:
continue # not a localized pinch — skip (node-safety)
width = 2.0 * area / length
boxes.append({'kind': 'box',
'bbox': (round(xmin, 9), round(ymin, 9),
round(xmax, 9), round(ymax, 9)),
'width': width, 'size': width / min_elems})
whole.sort(key=lambda d: d['poly_index'])
boxes.sort(key=lambda d: d['bbox'])
return whole + boxes
ensure_ccw_elements(nodes, elements, element_types)
Flip any clockwise 2D element to counter-clockwise winding, in place.
Downstream assembly routines assume positive (CCW) signed areas — the linear-triangle seepage stiffness, for example, skips elements with area <= 0, so a fully CW mesh assembles a singular system. Corner nodes are reversed and midside nodes remapped per element type.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def ensure_ccw_elements(nodes, elements, element_types):
"""
Flip any clockwise 2D element to counter-clockwise winding, in place.
Downstream assembly routines assume positive (CCW) signed areas — the
linear-triangle seepage stiffness, for example, skips elements with
area <= 0, so a fully CW mesh assembles a singular system. Corner nodes
are reversed and midside nodes remapped per element type.
Parameters:
nodes : (n_nodes, 2+) array of node coordinates
elements : (n_elements, 9) padded element node indices (modified in place)
element_types : (n_elements,) number of nodes per element (3/4/6/8/9)
Returns:
int: number of elements flipped
"""
# node-order permutations that reverse winding per element type
flip_orders = {
3: [0, 2, 1],
4: [0, 3, 2, 1],
6: [0, 2, 1, 5, 4, 3],
8: [0, 3, 2, 1, 7, 6, 5, 4],
9: [0, 3, 2, 1, 7, 6, 5, 4, 8],
}
nodes = np.asarray(nodes)
n_flipped = 0
for i in range(len(elements)):
etype = int(element_types[i])
n_corner = 3 if etype in (3, 6) else 4
corners = elements[i, :n_corner]
p = nodes[corners, :2]
# shoelace signed area (x2)
signed = 0.0
for a in range(n_corner):
b = (a + 1) % n_corner
signed += p[a, 0] * p[b, 1] - p[b, 0] * p[a, 1]
if signed < 0:
elem = elements[i, :etype].copy()
elements[i, :etype] = elem[flip_orders[etype]]
n_flipped += 1
return n_flipped
export_mesh_to_json(mesh, filename)
Save mesh dictionary to JSON file.
Source code in xslope/mesh.py
def export_mesh_to_json(mesh, filename):
"""Save mesh dictionary to JSON file."""
import json
import numpy as np
# Convert numpy arrays to lists for JSON serialization
mesh_json = {}
for key, value in mesh.items():
if isinstance(value, np.ndarray):
mesh_json[key] = value.tolist()
else:
mesh_json[key] = value
with open(filename, 'w') as f:
json.dump(mesh_json, f, indent=2)
print(f"Mesh saved to {filename}")
extract_1d_elements_from_2d_edges(nodes, elements_2d, element_types_2d, lines, debug=False)
Extract 1D elements from 2D element edges that lie along reinforcement lines. This ensures proper finite element integration where 1D elements are shared edges of 2D elements.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def extract_1d_elements_from_2d_edges(nodes, elements_2d, element_types_2d, lines, debug=False):
"""
Extract 1D elements from 2D element edges that lie along reinforcement lines.
This ensures proper finite element integration where 1D elements are shared edges of 2D elements.
Parameters:
nodes: np.ndarray of node coordinates (n_nodes, 2)
elements_2d: np.ndarray of 2D element vertex indices (n_elements, 9)
element_types_2d: np.ndarray indicating 2D element type (3, 4, 6, 8, or 9 nodes)
lines: List of reinforcement lines, each defined by list of (x, y) tuples
debug: Enable debug output
Returns:
tuple: (elements_1d, mat_ids_1d, element_node_counts_1d)
"""
import numpy as np
from collections import defaultdict
elements_1d = []
mat_ids_1d = []
element_node_counts_1d = []
# Build edge-to-element mapping from 2D elements
edge_to_element = defaultdict(list) # edge (n1, n2) -> list of element indices
element_edges = {} # element_idx -> list of edges
for elem_idx, (element, elem_type) in enumerate(zip(elements_2d, element_types_2d)):
edges = []
if elem_type in [3, 6]: # Triangle
# Triangle edges: (0,1), (1,2), (2,0)
corner_nodes = [element[0], element[1], element[2]]
edge_pairs = [(0, 1), (1, 2), (2, 0)]
for i, j in edge_pairs:
n1, n2 = corner_nodes[i], corner_nodes[j]
edge_key = (min(n1, n2), max(n1, n2)) # Canonical edge representation
edges.append(edge_key)
edge_to_element[edge_key].append(elem_idx)
elif elem_type in [4, 8, 9]: # Quadrilateral
# Quadrilateral edges: (0,1), (1,2), (2,3), (3,0)
corner_nodes = [element[0], element[1], element[2], element[3]]
edge_pairs = [(0, 1), (1, 2), (2, 3), (3, 0)]
for i, j in edge_pairs:
n1, n2 = corner_nodes[i], corner_nodes[j]
edge_key = (min(n1, n2), max(n1, n2)) # Canonical edge representation
edges.append(edge_key)
edge_to_element[edge_key].append(elem_idx)
element_edges[elem_idx] = edges
if debug:
print(f"Built edge map with {len(edge_to_element)} unique edges from {len(elements_2d)} 2D elements")
# For each reinforcement line, find 2D element edges that lie along it
for line_idx, line_pts in enumerate(lines):
line_pts_clean = remove_duplicate_endpoint(list(line_pts))
if len(line_pts_clean) < 2:
continue
if debug:
print(f"Processing reinforcement line {line_idx}: {line_pts_clean}")
# Find all 2D element edges that lie along this reinforcement line
line_edges = []
for edge_key, elem_indices in edge_to_element.items():
n1, n2 = edge_key
# Get coordinates of edge endpoints
coord1 = nodes[n1]
coord2 = nodes[n2]
# Check if this edge lies along the reinforcement line
if is_edge_on_reinforcement_line(coord1, coord2, line_pts_clean, tolerance=1e-6):
line_edges.append((n1, n2))
if debug:
print(f" Found edge ({n1}, {n2}) at coords {coord1} -> {coord2}")
# Sort edges to form continuous 1D elements along the line
if line_edges:
sorted_edges = sort_edges_along_line(line_edges, nodes, line_pts_clean, debug)
# Create 1D elements from sorted edges
for n1, n2 in sorted_edges:
# For linear elements, just use the two nodes
elements_1d.append([n1, n2, 0]) # Pad to 3 columns
mat_ids_1d.append(line_idx)
element_node_counts_1d.append(2)
if debug:
print(f" Created {len(sorted_edges)} 1D elements for line {line_idx}")
if debug:
print(f"Total 1D elements extracted: {len(elements_1d)}")
return elements_1d, mat_ids_1d, element_node_counts_1d
extract_constraint_line_geometry(slope_data)
Extract all constraint line geometry (reinforcement + piles) for mesh generation.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def extract_constraint_line_geometry(slope_data):
"""
Extract all constraint line geometry (reinforcement + piles) for mesh generation.
Parameters:
slope_data: Dictionary containing slope data
Returns:
lines: Combined list of constraint lines (reinforcement first, then piles)
n_reinf: Number of reinforcement lines
n_pile: Number of pile lines
"""
reinf_lines = extract_reinforcement_line_geometry(slope_data)
pile_lines = extract_pile_line_geometry(slope_data)
return reinf_lines + pile_lines, len(reinf_lines), len(pile_lines)
extract_pile_line_geometry(slope_data)
Extract pile line geometry from slope_data in the format needed for mesh generation.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def extract_pile_line_geometry(slope_data):
"""
Extract pile line geometry from slope_data in the format needed for mesh generation.
Parameters:
slope_data: Dictionary containing slope data with 'pile_lines' key
Returns:
List of pile lines, where each line is a list of (x, y) coordinate tuples
"""
lines = []
if 'pile_lines' in slope_data and slope_data['pile_lines']:
for pile in slope_data['pile_lines']:
lines.append([(pile['x1'], pile['y1']), (pile['x2'], pile['y2'])])
return lines
extract_point_constraints(slope_data)
Mesh point constraints from slope_data: currently the line-load application points (a node must land exactly at each so the load can be applied as a nodal force). Returns a list of (x, y) tuples for build_mesh_from_polygons(point_constraints=...).
Seepage BC vertices are NOT here: they are inserted as polygon vertices inside get_material_polygons(), the same way distributed-load points are, so that every caller gets them without having to opt in.
Source code in xslope/mesh.py
def extract_point_constraints(slope_data):
"""Mesh point constraints from slope_data: currently the line-load
application points (a node must land exactly at each so the load can be
applied as a nodal force). Returns a list of (x, y) tuples for
build_mesh_from_polygons(point_constraints=...).
Seepage BC vertices are NOT here: they are inserted as polygon vertices inside
get_material_polygons(), the same way distributed-load points are, so that every
caller gets them without having to opt in."""
return [(ll['x'], ll['y']) for ll in (slope_data.get('line_loads') or [])]
extract_reinforcement_line_geometry(slope_data)
Extract reinforcement line geometry from slope_data in the format needed for mesh generation.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def extract_reinforcement_line_geometry(slope_data):
"""
Extract reinforcement line geometry from slope_data in the format needed for mesh generation.
Parameters:
slope_data: Dictionary containing slope data with 'reinforce_lines' key
Returns:
List of reinforcement lines, where each line is a list of (x, y) coordinate tuples
"""
lines = []
if 'reinforce_lines' in slope_data and slope_data['reinforce_lines']:
for line in slope_data['reinforce_lines']:
# Convert from dict format to tuple format
line_coords = [(point['X'], point['Y']) for point in line]
lines.append(line_coords)
return lines
find_element_containing_point(nodes, elements, element_types, point)
Find which element contains the given point using spatial indexing for efficiency.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def find_element_containing_point(nodes, elements, element_types, point):
"""
Find which element contains the given point using spatial indexing for efficiency.
Parameters:
nodes: np.ndarray of node coordinates (n_nodes, 2)
elements: np.ndarray of element vertex indices (n_elements, 9) - unused nodes set to 0
element_types: np.ndarray indicating element type (3, 4, 6, 8, or 9 nodes)
point: tuple (x, y) coordinates of the point to find
Returns:
int: Index of the element containing the point, or -1 if not found
"""
x, y = point
# Use spatial indexing to find candidate elements quickly
# Build spatial hash grid if not already built, or rebuild if mesh changed
mesh_id = id(nodes)
cache = getattr(find_element_containing_point, '_cache', None)
if cache is None or cache['mesh_id'] != mesh_id:
find_element_containing_point._cache = {
'mesh_id': mesh_id,
'spatial_grid': _build_spatial_grid(nodes, elements, element_types)
}
spatial_grid = find_element_containing_point._cache['spatial_grid']
# Find grid cell containing the point
grid_x = int((x - spatial_grid['x_min']) / spatial_grid['cell_size'])
grid_y = int((y - spatial_grid['y_min']) / spatial_grid['cell_size'])
def _point_in_element(elem_idx):
element = elements[elem_idx]
elem_type = element_types[elem_idx]
if elem_type in [3, 6]: # Triangle (linear or quadratic)
# For point-in-element testing, use only corner nodes
x1, y1 = nodes[element[0]]
x2, y2 = nodes[element[1]]
x3, y3 = nodes[element[2]]
# Calculate barycentric coordinates
det = (y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3)
if abs(det) < 1e-12: # Degenerate triangle
return False
lambda1 = ((y2 - y3) * (x - x3) + (x3 - x2) * (y - y3)) / det
lambda2 = ((y3 - y1) * (x - x3) + (x1 - x3) * (y - y3)) / det
lambda3 = 1.0 - lambda1 - lambda2
# Check if point is inside triangle (all barycentric coordinates >= 0)
return lambda1 >= -1e-12 and lambda2 >= -1e-12 and lambda3 >= -1e-12
elif elem_type in [4, 8, 9]: # Quadrilateral (linear or quadratic)
# For point-in-element testing, use only corner nodes
x1, y1 = nodes[element[0]]
x2, y2 = nodes[element[1]]
x3, y3 = nodes[element[2]]
x4, y4 = nodes[element[3]]
# Use point-in-polygon test for quadrilaterals
# Check if point is inside by counting crossings
vertices = [(x1, y1), (x2, y2), (x3, y3), (x4, y4)]
inside = False
for j in range(len(vertices)):
xi, yi = vertices[j]
xj, yj = vertices[(j + 1) % len(vertices)]
if ((yi > y) != (yj > y)) and (x < (xj - xi) * (y - yi) / (yj - yi) + xi):
inside = not inside
return inside
return False
cells = spatial_grid['cells']
# Elements are registered into every cell their bounding box overlaps, so the
# point's own cell already holds the containing element in the common case.
# Check it first and return immediately — this avoids the wide neighbor sweep.
for elem_idx in cells.get((grid_x, grid_y), ()):
if _point_in_element(elem_idx):
return elem_idx
# Fallback: scan the surrounding ring. Rarely needed (e.g. a point landing
# exactly on a cell boundary with floating-point round-off), but keeps the
# result identical to a full search.
for dx in (-1, 0, 1):
for dy in (-1, 0, 1):
if dx == 0 and dy == 0:
continue
for elem_idx in cells.get((grid_x + dx, grid_y + dy), ()):
if _point_in_element(elem_idx):
return elem_idx
return -1 # Point not found in any element
get_material_polygons(slope_data, reinf_lines=None)
Return material-zone polygons for mesh generation, in the format the mesher expects: a list of dicts with 'coords' (closed (x, y) list) and 'mat_id'.
Unified entry point that works for both geometry inputs (plan_polygons.md §6): - Polygon input: use slope_data['polygons'] directly (the polygons ARE the geometry — no profile lines to convert). - Profile input: convert profile lines via build_polygons(), which also integrates reinforcement-line vertices when reinf_lines is provided.
Callers (main_seep/main_fem/main_mesh) should use this instead of calling build_polygons() directly, so they work on polygon-sheet inputs.
Distributed-load endpoints (and intermediate vertices) are inserted as polygon vertices when they fall on a boundary edge, so the mesh places nodes there and element edges align to the loaded span. Without this, the consistent edge-load integrator in fem.py (Pass 2a) drops the element edges that straddle a load end and under-applies the load (~one element edge of traction at each end).
Seepage BC polyline vertices are inserted the same way and for the same reason: a specified-flux BC loads a boundary edge only when both its corner nodes lie on the polyline, so an endpoint straddling an edge would silently drop it and apply less water than specified. See add_seep_bc_points_to_polygons.
Source code in xslope/mesh.py
def get_material_polygons(slope_data, reinf_lines=None):
"""
Return material-zone polygons for mesh generation, in the format the mesher
expects: a list of dicts with 'coords' (closed (x, y) list) and 'mat_id'.
Unified entry point that works for both geometry inputs (plan_polygons.md §6):
- Polygon input: use slope_data['polygons'] directly (the polygons ARE the
geometry — no profile lines to convert).
- Profile input: convert profile lines via build_polygons(), which also
integrates reinforcement-line vertices when reinf_lines is provided.
Callers (main_seep/main_fem/main_mesh) should use this instead of calling
build_polygons() directly, so they work on polygon-sheet inputs.
Distributed-load endpoints (and intermediate vertices) are inserted as polygon
vertices when they fall on a boundary edge, so the mesh places nodes there and
element edges align to the loaded span. Without this, the consistent edge-load
integrator in fem.py (Pass 2a) drops the element edges that straddle a load end
and under-applies the load (~one element edge of traction at each end).
Seepage BC polyline vertices are inserted the same way and for the same reason:
a specified-flux BC loads a boundary edge only when both its corner nodes lie on
the polyline, so an endpoint straddling an edge would silently drop it and apply
less water than specified. See add_seep_bc_points_to_polygons.
"""
if slope_data.get('profile_lines'):
polygons = build_polygons(slope_data, reinf_lines=reinf_lines)
else:
polygons = [
{'coords': list(p['polygon'].exterior.coords), 'mat_id': p['mat_id']}
for p in (slope_data.get('polygons') or [])
]
# Polygon-sheet input skips build_polygons, so the reinforcement/pile
# vertex integration must happen here too — otherwise gmsh cannot recover
# the constraint-line edges and silently produces an empty mesh.
if reinf_lines:
polygons = add_intersection_points_to_polygons(polygons, reinf_lines)
polygons = add_dload_points_to_polygons(_clean_pinchouts(polygons), slope_data)
return add_seep_bc_points_to_polygons(polygons, slope_data)
get_quad_mesh_presets()
Returns dictionary of preset quad meshing parameter combinations to try.
Source code in xslope/mesh.py
def get_quad_mesh_presets():
"""
Returns dictionary of preset quad meshing parameter combinations to try.
"""
presets = {
'default': {
"Mesh.Algorithm": 8,
"Mesh.RecombinationAlgorithm": 1,
"Mesh.SubdivisionAlgorithm": 1,
"Mesh.RecombineOptimizeTopology": 5,
"Mesh.Smoothing": 10,
"size_factor": 1.4, # Target size adjustment
},
'blossom': {
"Mesh.Algorithm": 6,
"Mesh.RecombinationAlgorithm": 2, # Blossom
"Mesh.SubdivisionAlgorithm": 1,
"Mesh.RecombineOptimizeTopology": 20,
"Mesh.Smoothing": 20,
"size_factor": 1.6, # Slightly larger for better recombination
},
'blossom_full': {
"Mesh.Algorithm": 5,
"Mesh.RecombinationAlgorithm": 3, # Blossom full-quad
"Mesh.SubdivisionAlgorithm": 1,
"Mesh.RecombineOptimizeTopology": 50,
"Mesh.Smoothing": 30,
"size_factor": 1.7, # Larger for complex recombination
},
'high_quality': {
"Mesh.Algorithm": 6,
"Mesh.RecombinationAlgorithm": 1,
"Mesh.SubdivisionAlgorithm": 1,
"Mesh.RecombineOptimizeTopology": 100,
"Mesh.RecombineNodeRepositioning": 1,
"Mesh.RecombineMinimumQuality": 0.1,
"Mesh.Smoothing": 50,
"Mesh.SmoothRatio": 2.0,
"size_factor": 2.0, # Much larger due to heavy optimization
},
'fast': {
"Mesh.Algorithm": 8,
"Mesh.RecombinationAlgorithm": 0, # Standard (fastest)
"Mesh.SubdivisionAlgorithm": 0,
"Mesh.RecombineOptimizeTopology": 0,
"Mesh.Smoothing": 5,
"size_factor": 0.7, # Smaller adjustment = more elements
}
}
return presets
import_mesh_from_json(filename)
Load mesh dictionary from JSON file.
Source code in xslope/mesh.py
def import_mesh_from_json(filename):
"""Load mesh dictionary from JSON file."""
import json
import numpy as np
with open(filename, 'r') as f:
mesh_json = json.load(f)
# Convert lists back to numpy arrays
mesh = {}
for key, value in mesh_json.items():
if isinstance(value, list):
mesh[key] = np.array(value)
else:
mesh[key] = value
# Defensive: normalize winding for meshes exported before CCW normalization
# was added to build_mesh_from_polygons (or edited externally).
if 'elements' in mesh and 'element_types' in mesh and 'nodes' in mesh:
ensure_ccw_elements(mesh['nodes'], mesh['elements'], mesh['element_types'])
return mesh
insert_point_into_polygon_edge(intersection, edge_start, edge_end, poly_data, point_map, target_size)
Insert an intersection point into a polygon edge, updating the polygon's coordinate list.
Source code in xslope/mesh.py
def insert_point_into_polygon_edge(intersection, edge_start, edge_end, poly_data, point_map, target_size):
"""Insert an intersection point into a polygon edge, updating the polygon's coordinate list."""
x, y = intersection
# Ensure the point exists in the point_map (for Gmsh)
if (x, y) not in point_map:
tag = len(point_map) + 1 # Simple tag assignment
point_map[(x, y)] = tag
# Insert the intersection point into the polygon's coordinate list at the correct edge
# poly_data['pt_tags'] is a list of Gmsh point tags, but we need to update the coordinate list used to build the polygon
# We'll reconstruct the coordinate list from the tags and point_map
pt_tags = poly_data['pt_tags']
# Build coordinate list for the polygon
coords = []
tag_to_coord = {v: k for k, v in point_map.items()}
for tag in pt_tags:
if tag in tag_to_coord:
coords.append(tag_to_coord[tag])
else:
# Fallback: try to find the coordinate in point_map
found = False
for (cx, cy), t in point_map.items():
if t == tag:
coords.append((cx, cy))
found = True
break
if not found:
coords.append((None, None)) # Should not happen
# Find the edge to insert after
insert_idx = None
for i in range(len(coords)):
a = coords[i]
b = coords[(i + 1) % len(coords)]
if (abs(a[0] - edge_start[0]) < 1e-8 and abs(a[1] - edge_start[1]) < 1e-8 and
abs(b[0] - edge_end[0]) < 1e-8 and abs(b[1] - edge_end[1]) < 1e-8):
insert_idx = i + 1
break
# Also check reversed edge
if (abs(a[0] - edge_end[0]) < 1e-8 and abs(a[1] - edge_end[1]) < 1e-8 and
abs(b[0] - edge_start[0]) < 1e-8 and abs(b[1] - edge_start[1]) < 1e-8):
insert_idx = i + 1
break
if insert_idx is not None:
# Insert the intersection point into the coordinate list
coords.insert(insert_idx, (x, y))
# Now update pt_tags to match
tag = point_map[(x, y)]
pt_tags.insert(insert_idx, tag)
# Update poly_data
poly_data['pt_tags'] = pt_tags
interpolate_at_point(nodes, elements, element_types, values, point, return_found=False, signed=False)
Interpolate values at a given point using the mesh.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def interpolate_at_point(nodes, elements, element_types, values, point, return_found=False, signed=False):
"""
Interpolate values at a given point using the mesh.
Parameters:
nodes: np.ndarray of node coordinates (n_nodes, 2)
elements: np.ndarray of element vertex indices (n_elements, 8)
element_types: np.ndarray indicating element type (3, 4, 6, or 8 nodes)
values: np.ndarray of values at nodes (n_nodes,)
point: tuple (x, y) coordinates of the point to interpolate at
return_found: if True, return (value, found) so callers can distinguish a
genuine interpolated zero from a point that fell outside the mesh
signed: if False (default) the interpolated value is clamped at 0.0 on the
historical assumption that pore pressure cannot be negative -- every
existing caller relies on this. Pass signed=True to receive the raw
interpolated value INCLUDING negative pressures; an unsaturated seepage
field carries genuine matric suction (u < 0) above the water table, and
the caller is then responsible for its own clamping where an
effective-normal pore pressure is required.
Returns:
float: Interpolated value at the point, or 0.0 if point not found
(or (value, found) when return_found=True)
"""
# Find the element containing the point
element_idx = find_element_containing_point(nodes, elements, element_types, point)
if element_idx == -1:
return (0.0, False) if return_found else 0.0 # Point not found in any element
element = elements[element_idx]
elem_type = element_types[element_idx]
x, y = point
if elem_type == 3: # Linear triangle
# Get triangle vertices and values
x1, y1 = nodes[element[0]]
x2, y2 = nodes[element[1]]
x3, y3 = nodes[element[2]]
v1 = values[element[0]]
v2 = values[element[1]]
v3 = values[element[2]]
# Calculate barycentric coordinates
det = (y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3)
lambda1 = ((y2 - y3) * (x - x3) + (x3 - x2) * (y - y3)) / det
lambda2 = ((y3 - y1) * (x - x3) + (x1 - x3) * (y - y3)) / det
lambda3 = 1.0 - lambda1 - lambda2
# Interpolate using barycentric coordinates
interpolated_value = lambda1 * v1 + lambda2 * v2 + lambda3 * v3
elif elem_type == 6: # Quadratic triangle
# Get all 6 nodes: corners (0,1,2) and midpoints (3,4,5)
# Node ordering: 0-1-2 corners, 3 midpoint of 0-1, 4 midpoint of 1-2, 5 midpoint of 2-0
corner_nodes = [element[0], element[1], element[2]]
midpoint_nodes = [element[3], element[4], element[5]]
# Get coordinates
x1, y1 = nodes[corner_nodes[0]] # Node 0
x2, y2 = nodes[corner_nodes[1]] # Node 1
x3, y3 = nodes[corner_nodes[2]] # Node 2
# Calculate barycentric coordinates (L1, L2, L3)
det = (y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3)
L1 = ((y2 - y3) * (x - x3) + (x3 - x2) * (y - y3)) / det
L2 = ((y3 - y1) * (x - x3) + (x1 - x3) * (y - y3)) / det
L3 = 1.0 - L1 - L2
# Quadratic shape functions for 6-node triangle
N = np.zeros(6)
N[0] = L1 * (2*L1 - 1) # Corner node 0
N[1] = L2 * (2*L2 - 1) # Corner node 1
N[2] = L3 * (2*L3 - 1) # Corner node 2
N[3] = 4 * L1 * L2 # Midpoint node 0-1
N[4] = 4 * L2 * L3 # Midpoint node 1-2
N[5] = 4 * L3 * L1 # Midpoint node 2-0
# Interpolate using quadratic shape functions
interpolated_value = 0.0
for i in range(6):
interpolated_value += N[i] * values[element[i]]
elif elem_type == 4: # Linear quadrilateral
# Get quadrilateral vertices and values
x1, y1 = nodes[element[0]]
x2, y2 = nodes[element[1]]
x3, y3 = nodes[element[2]]
x4, y4 = nodes[element[3]]
v1 = values[element[0]]
v2 = values[element[1]]
v3 = values[element[2]]
v4 = values[element[3]]
# Use proper bilinear shape functions for quadrilaterals
# Map to natural coordinates (xi, eta) in [-1, 1] x [-1, 1]
# For bilinear quad4, use iterative Newton-Raphson to find natural coordinates
# Initial guess at element center
xi, eta = 0.0, 0.0
# Newton-Raphson iteration to find (xi, eta) such that physical coordinates match
for _ in range(10): # Max 10 iterations
# Bilinear shape functions
N = np.array([
0.25 * (1-xi) * (1-eta), # Node 0
0.25 * (1+xi) * (1-eta), # Node 1
0.25 * (1+xi) * (1+eta), # Node 2
0.25 * (1-xi) * (1+eta) # Node 3
])
# Shape function derivatives
dN_dxi = np.array([
-0.25 * (1-eta), # Node 0
0.25 * (1-eta), # Node 1
0.25 * (1+eta), # Node 2
-0.25 * (1+eta) # Node 3
])
dN_deta = np.array([
-0.25 * (1-xi), # Node 0
-0.25 * (1+xi), # Node 1
0.25 * (1+xi), # Node 2
0.25 * (1-xi) # Node 3
])
# Current physical coordinates
x_curr = N[0]*x1 + N[1]*x2 + N[2]*x3 + N[3]*x4
y_curr = N[0]*y1 + N[1]*y2 + N[2]*y3 + N[3]*y4
# Residual
fx = x_curr - x
fy = y_curr - y
if abs(fx) < 1e-10 and abs(fy) < 1e-10:
break
# Jacobian
dx_dxi = dN_dxi[0]*x1 + dN_dxi[1]*x2 + dN_dxi[2]*x3 + dN_dxi[3]*x4
dx_deta = dN_deta[0]*x1 + dN_deta[1]*x2 + dN_deta[2]*x3 + dN_deta[3]*x4
dy_dxi = dN_dxi[0]*y1 + dN_dxi[1]*y2 + dN_dxi[2]*y3 + dN_dxi[3]*y4
dy_deta = dN_deta[0]*y1 + dN_deta[1]*y2 + dN_deta[2]*y3 + dN_deta[3]*y4
det_J = dx_dxi * dy_deta - dx_deta * dy_dxi
if abs(det_J) < 1e-12:
break
# Newton-Raphson update
dxi = (dy_deta * fx - dx_deta * fy) / det_J
deta = (-dy_dxi * fx + dx_dxi * fy) / det_J
xi -= dxi
eta -= deta
# Clamp to [-1,1]
xi = max(-1, min(1, xi))
eta = max(-1, min(1, eta))
# Final bilinear shape functions
N = np.array([
0.25 * (1-xi) * (1-eta), # Node 0
0.25 * (1+xi) * (1-eta), # Node 1
0.25 * (1+xi) * (1+eta), # Node 2
0.25 * (1-xi) * (1+eta) # Node 3
])
# Interpolate using bilinear shape functions
interpolated_value = N[0]*v1 + N[1]*v2 + N[2]*v3 + N[3]*v4
elif elem_type == 8: # Quadratic quadrilateral
# Get all 8 nodes: corners (0,1,2,3) and midpoints (4,5,6,7)
# Node ordering: 0-1-2-3 corners, 4 midpoint of 0-1, 5 midpoint of 1-2,
# 6 midpoint of 2-3, 7 midpoint of 3-0
# Get corner coordinates for mapping to natural coordinates
x1, y1 = nodes[element[0]] # Node 0
x2, y2 = nodes[element[1]] # Node 1
x3, y3 = nodes[element[2]] # Node 2
x4, y4 = nodes[element[3]] # Node 3
# For quadratic quads, we need to map from physical (x,y) to natural coordinates (xi,eta)
# This is complex for general quadrilaterals, so use simplified approach:
# Map to unit square [-1,1] x [-1,1] using bilinear mapping of corners
# Bilinear inverse mapping (approximate for general quads)
# Solve for natural coordinates xi, eta in [-1,1] x [-1,1]
# For simplicity, use area coordinate method similar to linear quad
# but with quadratic shape functions
# Calculate area coordinates (this is an approximation)
A_total = 0.5 * abs((x3-x1)*(y4-y2) - (x4-x2)*(y3-y1))
if A_total < 1e-12:
# Degenerate element, fall back to linear
A1 = abs((x - x1) * (y2 - y1) - (x2 - x1) * (y - y1)) / 2
A2 = abs((x - x2) * (y3 - y2) - (x3 - x2) * (y - y2)) / 2
A3 = abs((x - x3) * (y4 - y3) - (x4 - x3) * (y - y3)) / 2
A4 = abs((x - x4) * (y1 - y4) - (x1 - x4) * (y - y4)) / 2
A_sum = A1 + A2 + A3 + A4
if A_sum > 1e-12:
w1, w2, w3, w4 = A1/A_sum, A2/A_sum, A3/A_sum, A4/A_sum
else:
w1 = w2 = w3 = w4 = 0.25
# Linear interpolation as fallback
interpolated_value = (w1 * values[element[0]] + w2 * values[element[1]] +
w3 * values[element[2]] + w4 * values[element[3]])
else:
# For proper quadratic interpolation, we need natural coordinates
# This is a simplified implementation - full implementation would solve
# the nonlinear system for xi,eta
# Use parametric coordinates estimation
# Map point to approximate natural coordinates
xi_approx = 2 * (x - 0.5*(x1+x3)) / (x2+x3-x1-x4) if abs(x2+x3-x1-x4) > 1e-12 else 0
eta_approx = 2 * (y - 0.5*(y1+y3)) / (y2+y4-y1-y3) if abs(y2+y4-y1-y3) > 1e-12 else 0
# Clamp to [-1,1]
xi = max(-1, min(1, xi_approx))
eta = max(-1, min(1, eta_approx))
# Quadratic shape functions for 8-node quad in natural coordinates
N = np.zeros(8)
# Corner nodes
N[0] = 0.25 * (1-xi) * (1-eta) * (-xi-eta-1) # Node 0
N[1] = 0.25 * (1+xi) * (1-eta) * (xi-eta-1) # Node 1
N[2] = 0.25 * (1+xi) * (1+eta) * (xi+eta-1) # Node 2
N[3] = 0.25 * (1-xi) * (1+eta) * (-xi+eta-1) # Node 3
# Midpoint nodes
N[4] = 0.5 * (1-xi*xi) * (1-eta) # Node 4 (midpoint 0-1)
N[5] = 0.5 * (1+xi) * (1-eta*eta) # Node 5 (midpoint 1-2)
N[6] = 0.5 * (1-xi*xi) * (1+eta) # Node 6 (midpoint 2-3)
N[7] = 0.5 * (1-xi) * (1-eta*eta) # Node 7 (midpoint 3-0)
# Interpolate using quadratic shape functions
interpolated_value = 0.0
for i in range(8):
interpolated_value += N[i] * values[element[i]]
elif elem_type == 9: # Biquadratic quadrilateral (9-node Lagrange)
# Get all 9 nodes: corners (0,1,2,3), edges (4,5,6,7), and center (8)
# Node ordering: 0-1-2-3 corners, 4 midpoint of 0-1, 5 midpoint of 1-2,
# 6 midpoint of 2-3, 7 midpoint of 3-0, 8 center
# Get corner coordinates for mapping to natural coordinates
x1, y1 = nodes[element[0]] # Node 0
x2, y2 = nodes[element[1]] # Node 1
x3, y3 = nodes[element[2]] # Node 2
x4, y4 = nodes[element[3]] # Node 3
# Newton-Raphson iteration to find natural coordinates (xi, eta)
xi, eta = 0.0, 0.0 # Initial guess at element center
for _ in range(10): # Max 10 iterations
# Biquadratic Lagrange shape functions for all 9 nodes
N = np.zeros(9)
# Corner nodes
N[0] = 0.25 * xi * (xi-1) * eta * (eta-1) # Node 0: (-1,-1)
N[1] = 0.25 * xi * (xi+1) * eta * (eta-1) # Node 1: (1,-1)
N[2] = 0.25 * xi * (xi+1) * eta * (eta+1) # Node 2: (1,1)
N[3] = 0.25 * xi * (xi-1) * eta * (eta+1) # Node 3: (-1,1)
# Edge nodes
N[4] = 0.5 * (1-xi*xi) * eta * (eta-1) # Node 4: (0,-1)
N[5] = 0.5 * xi * (xi+1) * (1-eta*eta) # Node 5: (1,0)
N[6] = 0.5 * (1-xi*xi) * eta * (eta+1) # Node 6: (0,1)
N[7] = 0.5 * xi * (xi-1) * (1-eta*eta) # Node 7: (-1,0)
# Center node
N[8] = (1-xi*xi) * (1-eta*eta) # Node 8: (0,0)
# Shape function derivatives w.r.t. xi
dN_dxi = np.zeros(9)
dN_dxi[0] = 0.25 * (2*xi-1) * eta * (eta-1)
dN_dxi[1] = 0.25 * (2*xi+1) * eta * (eta-1)
dN_dxi[2] = 0.25 * (2*xi+1) * eta * (eta+1)
dN_dxi[3] = 0.25 * (2*xi-1) * eta * (eta+1)
dN_dxi[4] = -xi * eta * (eta-1)
dN_dxi[5] = 0.5 * (2*xi+1) * (1-eta*eta)
dN_dxi[6] = -xi * eta * (eta+1)
dN_dxi[7] = 0.5 * (2*xi-1) * (1-eta*eta)
dN_dxi[8] = -2*xi * (1-eta*eta)
# Shape function derivatives w.r.t. eta
dN_deta = np.zeros(9)
dN_deta[0] = 0.25 * xi * (xi-1) * (2*eta-1)
dN_deta[1] = 0.25 * xi * (xi+1) * (2*eta-1)
dN_deta[2] = 0.25 * xi * (xi+1) * (2*eta+1)
dN_deta[3] = 0.25 * xi * (xi-1) * (2*eta+1)
dN_deta[4] = 0.5 * (1-xi*xi) * (2*eta-1)
dN_deta[5] = -eta * xi * (xi+1)
dN_deta[6] = 0.5 * (1-xi*xi) * (2*eta+1)
dN_deta[7] = -eta * xi * (xi-1)
dN_deta[8] = -2*eta * (1-xi*xi)
# Current physical coordinates using all 9 nodes
node_coords = nodes[element[:9]]
x_curr = np.sum(N * node_coords[:, 0])
y_curr = np.sum(N * node_coords[:, 1])
# Residual
fx = x_curr - x
fy = y_curr - y
if abs(fx) < 1e-10 and abs(fy) < 1e-10:
break
# Jacobian
dx_dxi = np.sum(dN_dxi * node_coords[:, 0])
dx_deta = np.sum(dN_deta * node_coords[:, 0])
dy_dxi = np.sum(dN_dxi * node_coords[:, 1])
dy_deta = np.sum(dN_deta * node_coords[:, 1])
det_J = dx_dxi * dy_deta - dx_deta * dy_dxi
if abs(det_J) < 1e-12:
break
# Newton-Raphson update
dxi = (dy_deta * fx - dx_deta * fy) / det_J
deta = (-dy_dxi * fx + dx_dxi * fy) / det_J
xi -= dxi
eta -= deta
# Clamp to [-1,1]
xi = max(-1, min(1, xi))
eta = max(-1, min(1, eta))
# Final biquadratic shape functions
N = np.zeros(9)
N[0] = 0.25 * xi * (xi-1) * eta * (eta-1) # Node 0
N[1] = 0.25 * xi * (xi+1) * eta * (eta-1) # Node 1
N[2] = 0.25 * xi * (xi+1) * eta * (eta+1) # Node 2
N[3] = 0.25 * xi * (xi-1) * eta * (eta+1) # Node 3
N[4] = 0.5 * (1-xi*xi) * eta * (eta-1) # Node 4
N[5] = 0.5 * xi * (xi+1) * (1-eta*eta) # Node 5
N[6] = 0.5 * (1-xi*xi) * eta * (eta+1) # Node 6
N[7] = 0.5 * xi * (xi-1) * (1-eta*eta) # Node 7
N[8] = (1-xi*xi) * (1-eta*eta) # Node 8
# Interpolate using biquadratic shape functions
interpolated_value = 0.0
for i in range(9):
interpolated_value += N[i] * values[element[i]]
else:
return (0.0, False) if return_found else 0.0 # Unknown element type
# Clamp negative interpolated values to zero (pore pressure cannot be negative)
# unless the caller explicitly opts into the signed field via signed=True.
result = interpolated_value if signed else max(0.0, interpolated_value)
return (result, True) if return_found else result
is_edge_on_reinforcement_line(coord1, coord2, line_pts, tolerance=1e-06)
Check if an edge lies along a reinforcement line.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def is_edge_on_reinforcement_line(coord1, coord2, line_pts, tolerance=1e-6):
"""
Check if an edge lies along a reinforcement line.
Parameters:
coord1, coord2: Edge endpoint coordinates (x, y)
line_pts: List of (x, y) points defining the reinforcement line
tolerance: Tolerance for coincidence checking
Returns:
bool: True if edge lies along the reinforcement line
"""
x1, y1 = coord1
x2, y2 = coord2
# Check if both endpoints lie on the reinforcement line
point1_on_line = is_point_on_line_segments(coord1, line_pts, tolerance)
point2_on_line = is_point_on_line_segments(coord2, line_pts, tolerance)
if not (point1_on_line and point2_on_line):
return False
# Additional check: ensure edge direction is consistent with line direction
# This prevents selecting edges that cross the reinforcement line
edge_vector = np.array([x2 - x1, y2 - y1])
edge_length = np.linalg.norm(edge_vector)
if edge_length < tolerance:
return False
edge_unit = edge_vector / edge_length
# Check alignment with any segment of the reinforcement line
# This allows edges to span multiple segments after intersection preprocessing
for i in range(len(line_pts) - 1):
seg_start = np.array(line_pts[i])
seg_end = np.array(line_pts[i + 1])
seg_vector = seg_end - seg_start
seg_length = np.linalg.norm(seg_vector)
if seg_length < tolerance:
continue
seg_unit = seg_vector / seg_length
# Check if edge is aligned with this segment (or opposite direction)
dot_product = abs(np.dot(edge_unit, seg_unit))
if dot_product > 0.95: # Nearly parallel (cos(18°) ≈ 0.95)
# More flexible check: edge should be collinear with the reinforcement line
# and both endpoints should lie on the line (but not necessarily on the same segment)
return True
return False
is_point_on_edge(point, edge_start, edge_end, tol=1e-08)
Check if a point lies on a line segment (edge).
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def is_point_on_edge(point, edge_start, edge_end, tol=1e-8):
"""
Check if a point lies on a line segment (edge).
Parameters:
point: (x, y) tuple of point to check
edge_start: (x, y) tuple of edge start
edge_end: (x, y) tuple of edge end
tol: Tolerance for coincidence
Returns:
bool: True if point lies on edge segment
"""
px, py = point
x1, y1 = edge_start
x2, y2 = edge_end
# Check if point is within bounding box of edge
if not (min(x1, x2) - tol <= px <= max(x1, x2) + tol and
min(y1, y2) - tol <= py <= max(y1, y2) + tol):
return False
# Check if point is collinear with edge
# Use cross product to check collinearity
cross_product = abs((py - y1) * (x2 - x1) - (px - x1) * (y2 - y1))
# If cross product is close to zero, point is on the line
# Also check that it's within the segment bounds
if cross_product < tol:
# Check if point is between edge endpoints
dot_product = (px - x1) * (x2 - x1) + (py - y1) * (y2 - y1)
edge_length_sq = (x2 - x1) ** 2 + (y2 - y1) ** 2
if edge_length_sq < tol: # Edge is essentially a point
return abs(px - x1) < tol and abs(py - y1) < tol
# Parameter t should be between 0 and 1 for point to be on segment
t = dot_product / edge_length_sq
return -tol <= t <= 1 + tol
return False
is_point_on_line_segment(point, seg_start, seg_end, tolerance=1e-06)
Check if a point lies on a line segment.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def is_point_on_line_segment(point, seg_start, seg_end, tolerance=1e-6):
"""
Check if a point lies on a line segment.
Parameters:
point: (x, y) coordinates of point to check
seg_start: (x, y) coordinates of segment start
seg_end: (x, y) coordinates of segment end
tolerance: Tolerance for coincidence checking
Returns:
bool: True if point lies on the line segment
"""
px, py = point
x1, y1 = seg_start
x2, y2 = seg_end
# Check if point is within bounding box of segment
if not (min(x1, x2) - tolerance <= px <= max(x1, x2) + tolerance and
min(y1, y2) - tolerance <= py <= max(y1, y2) + tolerance):
return False
# Check collinearity using cross product
cross_product = abs((py - y1) * (x2 - x1) - (px - x1) * (y2 - y1))
# Check if cross product is close to zero (collinear)
if cross_product < tolerance:
# Verify point is between segment endpoints using dot product
dot_product = (px - x1) * (x2 - x1) + (py - y1) * (y2 - y1)
segment_length_sq = (x2 - x1) ** 2 + (y2 - y1) ** 2
if segment_length_sq < tolerance: # Degenerate segment
return abs(px - x1) < tolerance and abs(py - y1) < tolerance
# Parameter t should be between 0 and 1 for point to be on segment
t = dot_product / segment_length_sq
return -tolerance <= t <= 1 + tolerance
return False
is_point_on_line_segments(point, line_pts, tolerance=1e-06)
Check if a point lies on any segment of a multi-segment line.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def is_point_on_line_segments(point, line_pts, tolerance=1e-6):
"""
Check if a point lies on any segment of a multi-segment line.
Parameters:
point: (x, y) coordinates of point to check
line_pts: List of (x, y) points defining the line segments
tolerance: Tolerance for coincidence checking
Returns:
bool: True if point lies on any line segment
"""
for i in range(len(line_pts) - 1):
if is_point_on_line_segment(point, line_pts[i], line_pts[i + 1], tolerance):
return True
return False
line_segment_intersection(p1, p2, p3, p4, tol=1e-08)
Find intersection point between two line segments. Returns intersection point (x, y) if it exists, None otherwise.
Source code in xslope/mesh.py
def line_segment_intersection(p1, p2, p3, p4, tol=1e-8):
"""
Find intersection point between two line segments.
Returns intersection point (x, y) if it exists, None otherwise.
"""
x1, y1 = p1
x2, y2 = p2
x3, y3 = p3
x4, y4 = p4
# Calculate direction vectors
d1x, d1y = x2 - x1, y2 - y1
d2x, d2y = x4 - x3, y4 - y3
# Calculate determinant
det = d1x * d2y - d1y * d2x
if abs(det) < tol: # Lines are parallel
return None
# Calculate parameters
t1 = ((x3 - x1) * d2y - (y3 - y1) * d2x) / det
t2 = ((x3 - x1) * d1y - (y3 - y1) * d1x) / det
# Check if intersection is within both segments
if 0 <= t1 <= 1 and 0 <= t2 <= 1:
# Calculate intersection point
ix = x1 + t1 * d1x
iy = y1 + t1 * d1y
return (round(ix, 6), round(iy, 6))
return None
line_segment_parameter(point, line_start, line_end)
Calculate the parameter t (0 to 1) of a point along a line segment. Returns t where point = line_start + t * (line_end - line_start)
Source code in xslope/mesh.py
def line_segment_parameter(point, line_start, line_end):
"""
Calculate the parameter t (0 to 1) of a point along a line segment.
Returns t where point = line_start + t * (line_end - line_start)
"""
px, py = point
x1, y1 = line_start
x2, y2 = line_end
# Calculate parameter t
dx = x2 - x1
dy = y2 - y1
if abs(dx) > abs(dy):
t = (px - x1) / dx
else:
t = (py - y1) / dy
return t
make_polygons_conforming(polygon_coords, tol=1e-08, debug=False)
Make adjacent material-zone polygons share matching vertices along every common boundary, so the mesh is conforming across interfaces.
The mesh builder shares an interface between two zones only when both define it with identical endpoints (the gmsh edge is keyed on its endpoint points). If one zone places a vertex partway along an edge that its neighbour treats as a single straight segment — a T-junction / hanging node — the two sides become independent geometric paths and the mesh comes out non-conforming across that interface (a slit that blocks flow/stress transfer). This function detects any polygon vertex lying in the interior of another polygon's edge and splits that edge at the vertex, so both sides emit the same sub-edges and gmsh deduplicates them into one conforming interface.
Operates on a list of per-polygon (x, y) coordinate lists (the closing
endpoint, if duplicated, is dropped). Mutates polygon_coords in place and
returns it. Inserted coordinates are the exact foreign-vertex tuples, so they
match the neighbour's point under the downstream exact-coordinate de-duplication.
Source code in xslope/mesh.py
def make_polygons_conforming(polygon_coords, tol=1e-8, debug=False):
"""Make adjacent material-zone polygons share matching vertices along every
common boundary, so the mesh is conforming across interfaces.
The mesh builder shares an interface between two zones only when both define it
with identical endpoints (the gmsh edge is keyed on its endpoint points). If one
zone places a vertex partway along an edge that its neighbour treats as a single
straight segment — a **T-junction / hanging node** — the two sides become
independent geometric paths and the mesh comes out non-conforming across that
interface (a slit that blocks flow/stress transfer). This function detects any
polygon vertex lying in the *interior* of another polygon's edge and splits that
edge at the vertex, so both sides emit the same sub-edges and gmsh deduplicates
them into one conforming interface.
Operates on a list of per-polygon ``(x, y)`` coordinate lists (the closing
endpoint, if duplicated, is dropped). Mutates ``polygon_coords`` in place and
returns it. Inserted coordinates are the exact foreign-vertex tuples, so they
match the neighbour's point under the downstream exact-coordinate de-duplication.
"""
# Normalise to tuples and drop any duplicate closing endpoint.
polys = []
for pts in polygon_coords:
pts = [tuple(p) for p in pts]
if len(pts) > 1 and abs(pts[0][0] - pts[-1][0]) < tol and abs(pts[0][1] - pts[-1][1]) < tol:
pts = pts[:-1]
polys.append(pts)
# All distinct vertices across every polygon = candidate split points.
all_verts, seen = [], set()
for pts in polys:
for p in pts:
key = (round(p[0], 9), round(p[1], 9))
if key not in seen:
seen.add(key)
all_verts.append(p)
def _same(u, v):
return abs(u[0] - v[0]) < tol and abs(u[1] - v[1]) < tol
def _t(v, a, b): # projection parameter of v along a->b
dx, dy = b[0] - a[0], b[1] - a[1]
denom = dx * dx + dy * dy
return ((v[0] - a[0]) * dx + (v[1] - a[1]) * dy) / denom if denom else 0.0
n_inserted = 0
# Each pass only inserts already-existing vertices, so it converges quickly.
for _ in range(len(all_verts) + 1):
changed = False
for pi, pts in enumerate(polys):
n = len(pts)
if n < 2:
continue
new_pts = []
for i in range(n):
a, b = pts[i], pts[(i + 1) % n]
new_pts.append(a)
on_edge = [v for v in all_verts
if not _same(v, a) and not _same(v, b)
and is_point_on_edge(v, a, b, tol)]
for v in sorted(on_edge, key=lambda v: _t(v, a, b)):
if not _same(v, new_pts[-1]):
new_pts.append(v)
changed = True
n_inserted += 1
polys[pi] = new_pts
if not changed:
break
if debug and n_inserted:
print(f"make_polygons_conforming: inserted {n_inserted} vertex(es) to "
f"resolve T-junctions on shared edges")
polygon_coords[:] = polys
return polygon_coords
point_near_existing(point, existing_points, tol=1e-08)
Check if a point is near any existing points.
Source code in xslope/mesh.py
def point_near_existing(point, existing_points, tol=1e-8):
"""Check if a point is near any existing points."""
px, py = point
for ex, ey in existing_points:
if abs(px - ex) < tol and abs(py - ey) < tol:
return True
return False
print_mesh_connectivity_report(mesh, tolerance=1e-08)
Print a detailed report about mesh connectivity.
| Parameters: |
|
|---|
Source code in xslope/mesh.py
def print_mesh_connectivity_report(mesh, tolerance=1e-8):
"""
Print a detailed report about mesh connectivity.
Parameters:
mesh: Mesh dictionary
tolerance: Tolerance for considering nodes as duplicates
"""
results = verify_mesh_connectivity(mesh, tolerance)
print("=== MESH CONNECTIVITY REPORT ===")
print(f"Total nodes: {results['total_nodes']}")
print(f"Total elements: {results['total_elements']}")
print(f"Mesh is properly connected: {results['is_connected']}")
print()
if results['duplicate_node_groups']:
print(f"WARNING: Found {len(results['duplicate_node_groups'])} groups of duplicate nodes:")
for i, group in enumerate(results['duplicate_node_groups']):
print(f" Group {i+1}: Nodes {group} at position {mesh['nodes'][group[0]]}")
print()
if results['isolated_nodes']:
print(f"WARNING: Found {len(results['isolated_nodes'])} isolated nodes:")
for node_idx in results['isolated_nodes']:
print(f" Node {node_idx} at position {mesh['nodes'][node_idx]}")
print()
if results['elements_with_duplicates']:
print(f"WARNING: Found {len(results['elements_with_duplicates'])} elements with duplicate nodes:")
for elem_idx in results['elements_with_duplicates']:
print(f" Element {elem_idx}: {mesh['elements'][elem_idx]}")
print()
if results['is_connected']:
print("✓ Mesh connectivity is good - no duplicate nodes or isolated nodes found.")
else:
print("✗ Mesh connectivity issues detected. Consider regenerating the mesh.")
print_polygon_summary(polygons)
Prints a summary of the generated polygons for diagnostic purposes.
| Parameters: |
|
|---|
Source code in xslope/mesh.py
def print_polygon_summary(polygons):
"""
Prints a summary of the generated polygons for diagnostic purposes.
Parameters:
polygons: List of polygon coordinate lists or dicts with "coords"
"""
print("=== POLYGON SUMMARY ===")
print(f"Number of material zones: {len(polygons)}")
print()
for i, polygon in enumerate(polygons):
coords = polygon.get("coords") if isinstance(polygon, dict) else polygon
mat_id = polygon.get("mat_id") if isinstance(polygon, dict) else i
if mat_id is None:
mat_id = i
print(f"Material Zone {i+1} (Material ID: {mat_id}):")
print(f" Number of vertices: {len(coords)}")
# Calculate area (simple shoelace formula)
area = 0
for j in range(len(coords) - 1):
x1, y1 = coords[j]
x2, y2 = coords[j + 1]
area += (x2 - x1) * (y2 + y1) / 2
area = abs(area)
print(f" Approximate area: {area:.2f} square units")
# Print bounding box
xs = [x for x, y in coords]
ys = [y for x, y in coords]
print(f" Bounding box: x=[{min(xs):.2f}, {max(xs):.2f}], y=[{min(ys):.2f}, {max(ys):.2f}]")
print()
sort_edges_along_line(edges, nodes, line_pts, debug=False)
Sort edges to form a continuous sequence along a reinforcement line.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def sort_edges_along_line(edges, nodes, line_pts, debug=False):
"""
Sort edges to form a continuous sequence along a reinforcement line.
Parameters:
edges: List of (n1, n2) edge tuples
nodes: Node coordinates array
line_pts: Reinforcement line points
debug: Enable debug output
Returns:
list: Sorted list of (n1, n2) edge tuples
"""
if not edges:
return []
if len(edges) == 1:
return edges
# Build connectivity graph
node_connections = defaultdict(list)
for n1, n2 in edges:
node_connections[n1].append(n2)
node_connections[n2].append(n1)
# Find start node (should have only one connection, or be closest to line start)
line_start = np.array(line_pts[0])
line_end = np.array(line_pts[-1])
start_candidates = []
for node in node_connections:
if len(node_connections[node]) == 1: # End node
start_candidates.append(node)
if not start_candidates:
# No clear end nodes, use node closest to line start
min_dist = float('inf')
start_node = list(node_connections.keys())[0]
for node in node_connections:
dist = np.linalg.norm(nodes[node] - line_start)
if dist < min_dist:
min_dist = dist
start_node = node
else:
# Choose end node closest to line start
min_dist = float('inf')
start_node = start_candidates[0]
for node in start_candidates:
dist = np.linalg.norm(nodes[node] - line_start)
if dist < min_dist:
min_dist = dist
start_node = node
# Trace path from start node
sorted_edges = []
used_edges = set()
current_node = start_node
while True:
# Find next unused edge from current node
next_node = None
for neighbor in node_connections[current_node]:
edge_key = (min(current_node, neighbor), max(current_node, neighbor))
if edge_key not in used_edges:
next_node = neighbor
used_edges.add(edge_key)
sorted_edges.append((current_node, next_node))
break
if next_node is None:
break
current_node = next_node
if debug:
print(f" Sorted {len(sorted_edges)} edges along line")
return sorted_edges
test_1d_element_alignment(mesh, reinforcement_lines, tolerance=1e-06, debug=True)
Test that 1D elements correctly align with reinforcement lines.
This function verifies that: 1. Each reinforcement line is represented by a sequence of 1D elements 2. The 1D elements form continuous paths along each reinforcement line 3. The element endpoints match the expected line segment endpoints
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def test_1d_element_alignment(mesh, reinforcement_lines, tolerance=1e-6, debug=True):
"""
Test that 1D elements correctly align with reinforcement lines.
This function verifies that:
1. Each reinforcement line is represented by a sequence of 1D elements
2. The 1D elements form continuous paths along each reinforcement line
3. The element endpoints match the expected line segment endpoints
Parameters:
mesh: Dictionary containing nodes and 1D element data
reinforcement_lines: List of reinforcement lines, each containing coordinate tuples
tolerance: Tolerance for coordinate comparison (default 1e-6)
debug: Enable detailed debug output
Returns:
bool: True if all tests pass, False otherwise
"""
if debug:
print("\n=== Testing 1D Element Alignment ===")
if 'elements_1d' not in mesh:
print("ERROR: No 1D elements found in mesh")
return False
elements_1d = mesh['elements_1d']
if elements_1d is None or len(elements_1d) == 0:
print("ERROR: No 1D elements found in mesh")
return False
nodes = np.array(mesh['nodes'])
elements_1d = mesh['elements_1d']
if debug:
print(f"Testing {len(reinforcement_lines)} reinforcement lines")
print(f"Found {len(elements_1d)} 1D elements")
success = True
for line_idx, line_pts in enumerate(reinforcement_lines):
if debug:
print(f"\nTesting line {line_idx}: {line_pts}")
# Remove duplicate endpoints and get expected segments
line_pts_clean = remove_duplicate_endpoint(list(line_pts))
if len(line_pts_clean) < 2:
if debug:
print(f" Skipping line {line_idx}: insufficient points")
continue
# Expected segments for this line
expected_segments = []
for i in range(len(line_pts_clean) - 1):
expected_segments.append((line_pts_clean[i], line_pts_clean[i + 1]))
if debug:
print(f" Expected {len(expected_segments)} segments:")
for i, (start, end) in enumerate(expected_segments):
print(f" Segment {i}: {start} -> {end}")
# Find 1D elements that belong to this reinforcement line using material IDs
line_elements = []
if 'element_materials_1d' in mesh:
element_materials_1d = mesh['element_materials_1d']
for elem_idx, (element, material_id) in enumerate(zip(elements_1d, element_materials_1d)):
# Skip zero-padded elements
if len(element) < 2 or element[1] == 0:
continue
# Check if this element belongs to the current line
if material_id == line_idx + 1: # Material IDs are 1-based
# Get element coordinates
try:
coord1 = nodes[element[0]]
coord2 = nodes[element[1]]
except IndexError:
if debug:
print(f" WARNING: Element {elem_idx} has invalid node indices {element[0]}, {element[1]}")
continue
line_elements.append((elem_idx, coord1, coord2))
else:
# Fallback: use the old method if material IDs are not available
for elem_idx, element in enumerate(elements_1d):
# Skip zero-padded elements
if len(element) < 2 or element[1] == 0:
continue
# Get element coordinates
try:
coord1 = nodes[element[0]]
coord2 = nodes[element[1]]
except IndexError:
if debug:
print(f" WARNING: Element {elem_idx} has invalid node indices {element[0]}, {element[1]}")
continue
# Check if this element lies on the current reinforcement line
if is_edge_on_reinforcement_line(coord1, coord2, line_pts_clean, tolerance):
line_elements.append((elem_idx, coord1, coord2))
if debug:
print(f" Found {len(line_elements)} 1D elements on this line:")
for elem_idx, coord1, coord2 in line_elements:
print(f" Element {elem_idx}: {coord1} -> {coord2}")
# Test 1: Check that we have at least some 1D elements for this line
if len(line_elements) == 0:
print(f"ERROR: Line {line_idx} has no 1D elements")
success = False
continue
# Test 2: Check that we have reasonable number of elements
# After intersection preprocessing, we may have more elements than original segments
# But we should have at least some elements for each line
if len(line_elements) == 0:
print(f"ERROR: Line {line_idx} has no 1D elements")
success = False
continue
# Test 2: Check if elements form continuous path
if len(line_elements) > 1:
# Sort elements to form continuous sequence
sorted_elements = []
remaining_elements = line_elements.copy()
# Start with first element
current_elem = remaining_elements.pop(0)
sorted_elements.append(current_elem)
# Build chain by finding connecting elements
while remaining_elements:
last_coord = sorted_elements[-1][2] # End coordinate of last element
# Find next element that starts where last one ended
found_next = False
for i, (elem_idx, coord1, coord2) in enumerate(remaining_elements):
if np.linalg.norm(np.array(coord1) - np.array(last_coord)) < tolerance:
sorted_elements.append((elem_idx, coord1, coord2))
remaining_elements.pop(i)
found_next = True
break
elif np.linalg.norm(np.array(coord2) - np.array(last_coord)) < tolerance:
# Element is reversed, flip it
sorted_elements.append((elem_idx, coord2, coord1))
remaining_elements.pop(i)
found_next = True
break
if not found_next:
print(f"ERROR: Line {line_idx} elements do not form continuous path")
print(f" Cannot connect from {last_coord}")
print(f" Remaining elements: {remaining_elements}")
success = False
break
line_elements = sorted_elements
# Test 3: Check that the 1D elements cover the reinforcement line from start to end
if len(line_elements) > 0:
# Get the start and end points of the reinforcement line
line_start = line_pts_clean[0]
line_end = line_pts_clean[-1]
# Find the first and last 1D elements
first_elem = line_elements[0]
last_elem = line_elements[-1]
# Check if the first element starts near the line start
first_start_dist = np.linalg.norm(np.array(first_elem[1]) - np.array(line_start))
first_end_dist = np.linalg.norm(np.array(first_elem[2]) - np.array(line_start))
# Check if the last element ends near the line end
last_start_dist = np.linalg.norm(np.array(last_elem[1]) - np.array(line_end))
last_end_dist = np.linalg.norm(np.array(last_elem[2]) - np.array(line_end))
# The first element should start near the line start (either direction)
# Be more flexible due to intersection preprocessing
if first_start_dist > tolerance * 10 and first_end_dist > tolerance * 10:
print(f"WARNING: Line {line_idx} first element does not start at line start")
print(f" Line start: {line_start}")
print(f" First element: {first_elem[1]} -> {first_elem[2]}")
print(f" Start distances: {first_start_dist:.2e}, {first_end_dist:.2e}")
# Don't fail the test for this - just warn
# The last element should end near the line end (either direction)
# Be more flexible due to intersection preprocessing
if last_start_dist > tolerance * 10 and last_end_dist > tolerance * 10:
print(f"WARNING: Line {line_idx} last element does not end at line end")
print(f" Line end: {line_end}")
print(f" Last element: {last_elem[1]} -> {last_elem[2]}")
print(f" End distances: {last_start_dist:.2e}, {last_end_dist:.2e}")
# Don't fail the test for this - just warn
# Test 4: Check that line path is continuous
if len(line_elements) > 1:
for i in range(len(line_elements) - 1):
end_coord = line_elements[i][2] # End of current element
start_coord = line_elements[i + 1][1] # Start of next element
gap = np.linalg.norm(np.array(end_coord) - np.array(start_coord))
if gap > tolerance:
print(f"ERROR: Line {line_idx} has gap between elements {i} and {i+1}")
print(f" Gap size: {gap:.2e}")
print(f" Element {i} end: {end_coord}")
print(f" Element {i+1} start: {start_coord}")
success = False
if debug and success:
print(f" ✓ Line {line_idx} passes all alignment tests")
if debug:
if success:
print("\n=== All 1D Element Alignment Tests PASSED ===")
else:
print("\n=== 1D Element Alignment Tests FAILED ===")
return success
verify_mesh_connectivity(mesh, tolerance=1e-08)
Verify that the mesh is properly connected by checking for duplicate nodes at shared boundaries.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in xslope/mesh.py
def verify_mesh_connectivity(mesh, tolerance=1e-8):
"""
Verify that the mesh is properly connected by checking for duplicate nodes at shared boundaries.
Parameters:
mesh: Mesh dictionary with 'nodes' and 'elements' keys
tolerance: Tolerance for considering nodes as duplicates
Returns:
dict: Connectivity verification results
"""
import numpy as np
from collections import defaultdict
nodes = mesh["nodes"]
elements = mesh["elements"]
# Find duplicate nodes (nodes at same location)
duplicate_groups = []
used_indices = set()
for i in range(len(nodes)):
if i in used_indices:
continue
duplicates = [i]
for j in range(i + 1, len(nodes)):
if j in used_indices:
continue
if np.linalg.norm(nodes[i] - nodes[j]) < tolerance:
duplicates.append(j)
used_indices.add(j)
if len(duplicates) > 1:
duplicate_groups.append(duplicates)
used_indices.add(i)
# Check element connectivity
element_connectivity = defaultdict(set)
for elem_idx, element in enumerate(elements):
for node_idx in element:
element_connectivity[node_idx].add(elem_idx)
# Find isolated nodes (nodes not used by any element)
isolated_nodes = []
for i in range(len(nodes)):
if i not in element_connectivity:
isolated_nodes.append(i)
# Find elements with duplicate nodes
elements_with_duplicates = []
for elem_idx, element in enumerate(elements):
unique_nodes = set(element)
if len(unique_nodes) != len(element):
elements_with_duplicates.append(elem_idx)
results = {
"total_nodes": len(nodes),
"total_elements": len(elements),
"duplicate_node_groups": duplicate_groups,
"isolated_nodes": isolated_nodes,
"elements_with_duplicates": elements_with_duplicates,
"is_connected": len(duplicate_groups) == 0 and len(isolated_nodes) == 0
}
return results