Adaptive refinement of structured sweeps#
Notebook 24 introduced structured sweep bands;
notebooks 30/31
closed the solve → size map → remesh loop for unstructured meshes. This
notebook closes the same loop for models containing structured bands:
remesh_structured adapts every ResolutionSpec — for a band, that means
re-deriving its tangential/normal arrays from the size map (the geometry
.xao never changes), so the band stays an exact tensor grid through
every iteration.
import matplotlib.pyplot as plt
import numpy as np
import shapely
from meshwell.orchestrator import generate_mesh
from meshwell.polysurface import PolySurface
from meshwell.remesh import remesh_structured
from meshwell.resolution import StructuredSweepResolutionSpec
from meshwell.structured.sweep import StructuredSweep
from meshwell.visualization import plot2D
/home/runner/work/meshwell/meshwell/.venv/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdm
Model: two-layer stack with a quantum-well band#
Same fixture as notebook 24: two boxes meeting at y=1, a "qw" band of
thickness 0.4 grown into upper. We checkpoint the CAD to a .xao — the
adaptive loop regenerates from it on every iteration.
sweep = StructuredSweep(name="qw", on="lower___upper", thickness={"upper": 0.4})
specs = {
"qw": [StructuredSweepResolutionSpec(tangential=1.0, normal={"upper": 4})],
}
mesh0 = generate_mesh(
entities=[
PolySurface(
polygons=shapely.box(0, 0, 4, 1), physical_name="lower", mesh_order=2
),
PolySurface(
polygons=shapely.box(0, 1, 4, 2), physical_name="upper", mesh_order=1
),
],
sweeps=[sweep],
dim=2,
output_mesh="adaptive_structured_0.msh",
checkpoint_cad="adaptive_structured.xao",
default_characteristic_length=0.5,
resolution_specs=specs,
)
plot2D(mesh0, title="Iteration 0", wireframe=True)
Info : Clearing all models and views...
Info : Done clearing all models and views
Info : Reading '/tmp/tmpe9x4p3r2/cad.xao'...
Info : Done reading '/tmp/tmpe9x4p3r2/cad.xao'
Info : Writing 'adaptive_structured.xao'...
Info : Done writing 'adaptive_structured.xao'
Signal: an analytic boundary-layer “solution”#
The signal contract is unchanged from remesh.py: the solver supplies
(x, y, z, value) data (or a ready size map). Here we stand in for the
solver with an analytic field u = tanh((y - 1)/0.1) — a sharp transition
at the interface, as a depletion region or thermal boundary layer would
produce — and turn its gradient into a target size.
def size_map_from_solution():
xs, ys = np.meshgrid(
np.linspace(0.0, 4.0, 81), np.linspace(0.0, 2.0, 81), indexing="ij"
)
grad = 10.0 / np.cosh((ys - 1.0) / 0.1) ** 2 # |du/dy|
sizes = np.clip(0.5 / (1.0 + grad), 0.02, 0.5)
return np.column_stack([xs.ravel(), ys.ravel(), np.zeros(xs.size), sizes.ravel()])
The loop: solve → size map → remesh_structured#
Three lines per iteration. The driver damps per-iteration size changes
(change_max) and gradation-limits both the band arrays and the point
cloud (max_ratio), so the loop converges without caller-side hygiene.
mesh_i, specs_i = mesh0, specs
node_counts = [len(mesh0.points)]
for it in range(1, 3):
mesh_i, specs_i = remesh_structured(
input_mesh=mesh_i,
geometry_file="adaptive_structured.xao",
sweeps=[sweep],
resolution_specs=specs_i,
size_map=size_map_from_solution(),
output_mesh=f"adaptive_structured_{it}.msh",
default_characteristic_length=0.5,
)
node_counts.append(len(mesh_i.points))
print("node counts per iteration:", node_counts)
Info : Clearing all models and views...
Info : Done clearing all models and views
Info : Reading 'adaptive_structured.xao'...
Info : Done reading 'adaptive_structured.xao'
Info : Clearing all models and views...
Info : Done clearing all models and views
Info : Reading 'adaptive_structured.xao'...
Info : Done reading 'adaptive_structured.xao'
node counts per iteration: [63, 163, 305]
plot2D(mesh_i, title="Iteration 2: band refined toward the interface", wireframe=True)
The adapted spec is inspectable state — explicit arrays, not hidden mesh:
off = np.asarray(specs_i["qw"][0].normal["upper"])
print("normal offsets:", np.round(off, 4))
plt.figure(figsize=(5, 3))
plt.step(off[:-1], np.diff(off), where="post")
plt.xlabel("η (offset from interface)")
plt.ylabel("normal cell size")
plt.title("Adapted normal grading (no Graded ratio was ever stated)")
plt.show()
normal offsets: [0. 0.043 0.099 0.173 0.27 0.4 ]
Point refinement inside a structured band#
A size map that targets a single point inside a structured band shows the anisotropic character of sweep adaptation: the band’s normal and tangential arrays respond independently (directional min-collapse), so refinement clusters around the point in both directions while the band remains an exact tensor grid — no unstructured island, no admissibility repair.
P = np.array([2.0, 1.2])
def point_size_map():
xs, ys = np.meshgrid(
np.linspace(0.0, 4.0, 121), np.linspace(0.6, 1.8, 61), indexing="ij"
)
dist = np.hypot(xs - P[0], ys - P[1])
sizes = np.clip(0.02 + 0.3 * dist, 0.02, 0.5)
return np.column_stack([xs.ravel(), ys.ravel(), np.zeros(xs.size), sizes.ravel()])
Adapt (three iterations — the per-iteration clamp change_max=2#
bounds how fast sizes may move, so a deep target takes a few steps)
mesh_i, specs_i = mesh0, specs
for it in range(1, 4):
mesh_i, specs_i = remesh_structured(
input_mesh=mesh_i,
geometry_file="point_refine.xao",
sweeps=[sweep],
resolution_specs=specs_i,
size_map=point_size_map(),
output_mesh=f"point_refine_{it}.msh",
default_characteristic_length=0.5,
)
plot2D(mesh_i, title="After: refinement clustered at (2.0, 1.2)", wireframe=True)
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[8], line 3
1 mesh_i, specs_i = mesh0, specs
2 for it in range(1, 4):
----> 3 mesh_i, specs_i = remesh_structured(
4 input_mesh=mesh_i,
5 geometry_file="point_refine.xao",
6 sweeps=[sweep],
File ~/work/meshwell/meshwell/meshwell/remesh.py:1048, in remesh_structured(input_mesh, geometry_file, sweeps, resolution_specs, strategies, size_map, change_max, max_ratio, output_mesh, dim, global_2D_algorithm, global_3D_algorithm, mesh_element_order, optimization_flags, default_characteristic_length, verbosity, n_threads, filename, model, cad_settings)
1043 size_map = gradation_limit_size_map(np.asarray(size_map, float), max_ratio)
1045 remesher = RemeshGMSH(
1046 n_threads=n_threads, filename=filename, model=model, verbosity=verbosity
1047 )
-> 1048 settings = remesher._resolve_intake(geometry_file, cad_settings)
1049 try:
1050 remesher.model_manager.ensure_initialized(str(remesher.model_manager.filename))
File ~/work/meshwell/meshwell/meshwell/remesh.py:175, in Remesher._resolve_intake(self, geometry_file, cad_settings)
167 """Apply the mesh-stage ``.xao`` intake rules (raise on missing/mismatch).
168
169 A caller-supplied model is checked for consistency; an internally
170 created one (``_owns_model``) is not, since its ``point_tolerance``
171 was never chosen by the caller.
172 """
173 from meshwell.mesh import _resolve_mesh_cad_settings
--> 175 settings = _resolve_mesh_cad_settings(
176 input_file=geometry_file,
177 model=None if self._owns_model else self.model_manager,
178 cad_settings=cad_settings,
179 point_tolerance=None,
180 )
181 self.model_manager.cad_settings = settings
182 return settings
File ~/work/meshwell/meshwell/meshwell/mesh.py:718, in _resolve_mesh_cad_settings(input_file, model, cad_settings, point_tolerance)
716 model_settings = model.cad_settings if model is not None else None
717 if input_file is not None:
--> 718 settings = CADSettings.resolve_for_intake(
719 input_file, cad_settings, point_tolerance=point_tolerance
720 )
721 settings.check_matches(
722 model_settings, self_label="xao", other_label="model.cad_settings"
723 )
724 else:
File ~/work/meshwell/meshwell/meshwell/cad_settings.py:418, in CADSettings.resolve_for_intake(cls, xao_path, supplied, point_tolerance)
414 if supplied is not None and not isinstance(supplied, CADSettings):
415 raise TypeError(
416 f"cad_settings must be a CADSettings, got {type(supplied).__name__}"
417 )
--> 418 stored = cls.from_xao(xao_path)
419 if stored is None and supplied is None:
420 raise MissingCADSettingsError(
421 f"{xao_path} carries no meshwell CAD-settings metadata (it was not "
422 "written by meshwell.cad / write_xao(cad_settings=...)). Pass "
423 "cad_settings=CADSettings(...) describing how it was generated."
424 )
File ~/work/meshwell/meshwell/meshwell/cad_settings.py:316, in CADSettings.from_xao(cls, xao_path)
309 """Read embedded settings from a ``.xao``; ``None`` if absent.
310
311 Raises:
312 CADSettingsError: if a meshwell block exists but is malformed or
313 uses an unsupported schema version.
314 """
315 path = Path(xao_path).with_suffix(".xao")
--> 316 with path.open("rb") as fh:
317 fh.seek(0, 2)
318 size = fh.tell()
File ~/work/_temp/uv-python-dir/cpython-3.13.15-linux-x86_64-gnu/lib/python3.13/pathlib/_local.py:537, in Path.open(self, mode, buffering, encoding, errors, newline)
535 if "b" not in mode:
536 encoding = io.text_encoding(encoding)
--> 537 return io.open(self, mode, buffering, encoding, errors, newline)
FileNotFoundError: [Errno 2] No such file or directory: 'point_refine.xao'
The anisotropy, made explicit#
The adapted spec carries the two 1D arrays. The tangential spacing h_t(ξ) dips near ξ = 2.0 and the normal spacing h_n(η) dips near η = 0.2 — each direction saw its own min-collapse of the same isotropic target.
spec = specs_i["qw"][0]
t = np.asarray(spec.tangential)
n = np.asarray(spec.normal["upper"])
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3))
ax1.step(t[:-1], np.diff(t), where="post")
ax1.axvline(P[0], color="k", ls=":", label="target ξ")
ax1.set_xlabel("ξ (tangential)")
ax1.set_ylabel("h_t")
ax1.legend()
ax2.step(n[:-1], np.diff(n), where="post")
ax2.axvline(P[1] - 1.0, color="k", ls=":", label="target η")
ax2.set_xlabel("η (normal, from interface)")
ax2.set_ylabel("h_n")
ax2.legend()
fig.suptitle("Independent tangential / normal response to a point target")
fig.tight_layout()
plt.show()
The band is still an exact tensor grid — every interior node lies on a grid line — which is the admissibility-by-construction witness:
pts = mesh_i.points[:, :2]
band = pts[(pts[:, 1] >= 1.0 - 1e-9) & (pts[:, 1] <= 1.4 + 1e-9)]
xs = np.unique(np.round(band[:, 0], 9))
ys = np.unique(np.round(band[:, 1], 9))
if len(band) != len(xs) * len(ys):
msg = "structured-band nodes do not form a tensor grid"
raise ValueError(msg)
print(f"tensor grid intact: {len(xs)} x {len(ys)} = {len(band)} band nodes")