Download this testcase.
X-Mesh ellipse membrane relaxation
A fixed-connectivity X-Mesh follows an advected level-set interface while enriched pressure and a membrane force relax an ellipse toward a circle.
import argparse
import json
import os
import shutil
import subprocess
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from migflow import fluid, gmsh_io
from _membrane_ale import (interface_edges, mesh_quality,
mesh_shape_quality, ordered_interface,
polygon_area)
HERE = os.path.dirname(os.path.abspath(__file__))
SIGMA = 1.0
RHO = 1.0
ELLIPSE_CENTER = np.array([0.5, 0.5])
ELLIPSE_AXES = np.array([0.32, 0.18])
class XMeshTopology:
"""Fixed connectivity and vertex adjacency required by move_front."""
def __init__(self, elements, n_nodes, boundary_nodes):
self.el = np.asarray(elements, dtype=np.int32)
edge_blocks = [self.el[:, [0, 1]], self.el[:, [1, 2]],
self.el[:, [2, 0]]]
self.edges = np.unique(
np.sort(np.vstack(edge_blocks), axis=1), axis=0).astype(np.int32)
directed = np.vstack((self.edges, self.edges[:, ::-1]))
directed = directed[np.lexsort((directed[:, 1], directed[:, 0]))]
self.neigh_start = np.cumsum(
np.concatenate(([0], np.bincount(directed[:, 0],
minlength=n_nodes))))
self.neigh = directed[:, 1]
self.boundary = np.zeros(n_nodes, dtype=bool)
self.boundary[np.asarray(boundary_nodes, dtype=np.int32)] = True
def relax(front, x, x0, alpha, elements=None, reference_sign=None,
minimum_area=None, minimum_shape=None):
"""Relax off-front nodes toward x0, with optional quality backtracking."""
alpha = float(np.clip(alpha, 0.0, 1.0))
amount = np.where(front, 0.0, alpha)
while True:
relaxed = x * (1.0 - amount[:, None]) + x0 * amount[:, None]
if elements is None:
return relaxed
area_ok = (mesh_quality(relaxed, elements, reference_sign) >=
minimum_area)
shape_ok = mesh_shape_quality(relaxed, elements) >= minimum_shape
if area_ok and shape_ok:
return relaxed
alpha *= 0.5
if alpha < 1e-6:
return x.copy()
amount = np.where(front, 0.0, alpha)
def move_front_python(mesh, front, x, f_old, f, allow_handoff=True):
"""Move front/sign-changing vertices to a zero of the linear level set."""
x_next = x.copy()
front = np.asarray(front, dtype=bool)
candidates = front | ((f_old > 0.0) != (f > 0.0))
candidates[mesh.boundary] = False
front_next = np.zeros_like(front)
proposals = {}
for i in np.flatnonzero(candidates):
neigh = mesh.neigh[mesh.neigh_start[i]:mesh.neigh_start[i + 1]]
opposite = neigh[(~front[neigh]) &
((f[neigh] > 0.0) != (f[i] > 0.0))]
if opposite.size == 0:
continue
denominator = f[i] - f[opposite]
valid = np.abs(denominator) > 1e-14
opposite = opposite[valid]
denominator = denominator[valid]
if opposite.size == 0:
continue
fraction = np.clip(f[i] / denominator, 0.0, 1.0)
roots = ((1.0 - fraction)[:, None] * x[i] +
fraction[:, None] * x[opposite])
target = np.argmin(np.linalg.norm(roots[:, :2] - x[i, :2], axis=1))
neighbour = opposite[target]
root = roots[target]
move_i = np.linalg.norm(root[:2] - x[i, :2])
move_neighbour = np.linalg.norm(root[:2] - x[neighbour, :2])
# Change the owning vertex before the zero becomes too close to the
# opposite end of the edge. Waiting for the midpoint (the purely
# nearest-vertex rule) creates a very thin triangle just before the
# handoff and makes reducing dt ineffective.
owner = (neighbour if allow_handoff and move_neighbour < 2.0 * move_i and
not mesh.boundary[neighbour] else i)
movement = min(move_i, move_neighbour) if owner == neighbour else move_i
if owner not in proposals or movement < proposals[owner][0]:
proposals[owner] = (movement, root)
for node, (_, root) in proposals.items():
front_next[node] = True
x_next[node] = root
return x_next, front_next
def snap_new_front_nodes(mesh, nodes, front, coordinates, levelset):
"""Project newly activated interface nodes onto the P1 zero level set."""
moved = coordinates.copy()
for node in np.flatnonzero(nodes):
if abs(levelset[node]) <= 1e-14:
continue
neighbours = mesh.neigh[
mesh.neigh_start[node]:mesh.neigh_start[node + 1]]
opposite = neighbours[
(~front[neighbours]) &
(levelset[node] * levelset[neighbours] < 0.0)]
if opposite.size == 0:
raise RuntimeError(
f"cannot snap new interface node {node} to the level set")
denominator = levelset[node] - levelset[opposite]
fraction = np.clip(levelset[node] / denominator, 0.0, 1.0)
roots = ((1.0 - fraction)[:, None] * coordinates[node] +
fraction[:, None] * coordinates[opposite])
target = np.argmin(
np.linalg.norm(roots[:, :2] - coordinates[node, :2], axis=1))
moved[node] = roots[target]
return moved
def ellipse_levelset(coordinates):
scaled = (coordinates[:, :2] - ELLIPSE_CENTER) / ELLIPSE_AXES
return np.sum(scaled ** 2, axis=1) - 1.0
def initialize_front(mesh, x0, max_iterations=30):
x = x0.copy()
front = np.zeros(x.shape[0], dtype=bool)
levelset = np.ones(x.shape[0])
for _ in range(max_iterations):
old_x = x.copy()
old_front = front.copy()
old_levelset = levelset.copy()
levelset = ellipse_levelset(x)
x, front = move_front_python(
mesh, front, x, old_levelset, levelset, allow_handoff=False)
x = relax(front, x, x0, 1.0)
levelset = ellipse_levelset(x)
levelset[front] = 0.0
if np.array_equal(front, old_front) and np.max(np.abs(x - old_x)) < 1e-12:
break
else:
raise RuntimeError("X-Mesh front initialization did not converge")
return x, front, levelset
def levelset_gradient(mesh, coordinates, levelset):
triangle = mesh.el
points = coordinates[triangle, :2]
edge1 = points[:, 1] - points[:, 0]
edge2 = points[:, 2] - points[:, 0]
determinant = edge1[:, 0] * edge2[:, 1] - edge1[:, 1] * edge2[:, 0]
rhs1 = levelset[triangle[:, 1]] - levelset[triangle[:, 0]]
rhs2 = levelset[triangle[:, 2]] - levelset[triangle[:, 0]]
gradient = np.column_stack((
(rhs1 * edge2[:, 1] - rhs2 * edge1[:, 1]) / determinant,
(-rhs1 * edge2[:, 0] + rhs2 * edge1[:, 0]) / determinant))
weight = np.abs(determinant) / 2.0
nodal_gradient = np.zeros((coordinates.shape[0], 2))
nodal_weight = np.zeros(coordinates.shape[0])
for local in range(3):
np.add.at(nodal_gradient, triangle[:, local], gradient * weight[:, None])
np.add.at(nodal_weight, triangle[:, local], weight)
nodal_gradient /= nodal_weight[:, None]
return nodal_gradient
def advect_levelset(mesh, coordinates, levelset, velocity, dt):
gradient = levelset_gradient(mesh, coordinates, levelset)
return levelset - dt * np.sum(velocity * gradient, axis=1)
def classify_elements(mesh, levelset, front, previous_phases):
sign = np.where(levelset > 0.0, 1, -1)
sign[front] = 0
element_sign = sign[mesh.el]
phases = (np.sum(element_sign, axis=1) < 0).astype(np.int32)
mixed = np.any(element_sign < 0, axis=1) & np.any(element_sign > 0, axis=1)
phases[mixed] = previous_phases[mixed]
phases[np.all(element_sign == 0, axis=1)] = previous_phases[
np.all(element_sign == 0, axis=1)]
return phases
def _closed_interface_front(mesh, phases):
"""Return the interface-node mask, or None for an invalid topology."""
edges = interface_edges(mesh.el, phases)
if edges.size == 0:
return None
nodes, counts = np.unique(edges, return_counts=True)
if np.any(counts != 2) or nodes.size != edges.shape[0]:
return None
ordered, closed = ordered_interface(edges, np.zeros((mesh.boundary.size, 2)))
if not closed or ordered.size != nodes.size:
return None
result = np.zeros(mesh.boundary.size, dtype=bool)
result[nodes] = True
if np.any(result & mesh.boundary):
return None
return result
def initial_phases(mesh, coordinates):
centroids = coordinates[mesh.el, :2].mean(axis=1)
scaled = (centroids - ELLIPSE_CENTER) / ELLIPSE_AXES
return (np.sum(scaled ** 2, axis=1) < 1.0).astype(np.int32)
def interface_is_conforming(mesh, phases, front):
edges = interface_edges(mesh.el, phases)
return edges.size > 0 and np.all(front[edges])
def interface_phases_match_levelset(mesh, phases, front, levelset):
"""Check the phase on triangles resolved by one off-front vertex."""
front_count = np.count_nonzero(front[mesh.el], axis=1)
resolved = np.flatnonzero(front_count == 2)
for element in resolved:
third = mesh.el[element][~front[mesh.el[element]]][0]
expected = int(levelset[third] < 0.0)
if phases[element] != expected:
return False
return True
def redistance_levelset(mesh, coordinates, phases, front):
"""Build the signed distance to the current piecewise-linear interface."""
edges = interface_edges(mesh.el, phases)
if edges.size == 0 or not np.all(front[edges]):
raise RuntimeError("cannot redistance a non-conforming interface")
point0 = coordinates[edges[:, 0], :2]
tangent = coordinates[edges[:, 1], :2] - point0
length2 = np.sum(tangent ** 2, axis=1)
if np.any(length2 <= 1e-24):
raise RuntimeError("cannot redistance with a zero-length interface edge")
relative = coordinates[:, None, :2] - point0[None, :, :]
coordinate = np.clip(
np.sum(relative * tangent[None, :, :], axis=2) /
length2[None, :], 0.0, 1.0)
projection = point0[None, :, :] + coordinate[:, :, None] * tangent
distance = np.sqrt(np.min(np.sum(
(coordinates[:, None, :2] - projection) ** 2, axis=2), axis=1))
minimum_phase = np.ones(coordinates.shape[0], dtype=np.int32)
maximum_phase = np.zeros(coordinates.shape[0], dtype=np.int32)
for local in range(3):
np.minimum.at(minimum_phase, mesh.el[:, local], phases)
np.maximum.at(maximum_phase, mesh.el[:, local], phases)
ambiguous = (minimum_phase != maximum_phase) & ~front
if np.any(ambiguous):
nodes = np.flatnonzero(ambiguous)
raise RuntimeError(
f"off-front nodes touch both phases: {nodes[:8].tolist()}")
result = np.where(maximum_phase == 1, -distance, distance)
result[front] = 0.0
return result
def capillary_timestep(coordinates, edges, safety=0.8):
length = np.linalg.norm(
coordinates[edges[:, 1], :2] - coordinates[edges[:, 0], :2], axis=1)
return safety * np.sqrt(2.0 * RHO * length.min() ** 3 /
(4.0 * np.pi * SIGMA))
def interface_points(elements, phases, coordinates):
edges = interface_edges(elements, phases)
nodes, closed = ordered_interface(edges, coordinates)
if not closed:
raise RuntimeError("the X-Mesh ellipse interface is not closed")
return coordinates[nodes, :2], edges
def plot_interfaces(initial, final, target_radius, filename):
figure, axis = plt.subplots(figsize=(5.4, 5.0), constrained_layout=True)
initial = np.vstack((initial, initial[0]))
final = np.vstack((final, final[0]))
axis.plot(initial[:, 0], initial[:, 1], "--", color="0.55",
label="initial X-Mesh ellipse")
axis.plot(final[:, 0], final[:, 1], color="#1565c0", lw=2,
label="final X-Mesh interface")
angle = np.linspace(0.0, 2.0 * np.pi, 300)
axis.plot(0.5 + target_radius * np.cos(angle),
0.5 + target_radius * np.sin(angle), ":", color="#d84315",
label="equal-area circle")
axis.set(aspect="equal", xlim=(0.12, 0.88), ylim=(0.12, 0.88),
xlabel="$x$", ylabel="$y$", title="X-Mesh membrane relaxation")
axis.legend(frameon=False)
figure.savefig(filename, dpi=180)
plt.close(figure)
def write_state(problem, output_dir, time, levelset, front):
"""Write the fluid solution and the two X-Mesh diagnostic P1 fields."""
fields = problem.get_default_export()
p1 = problem.get_p1_element()
fields["levelset"] = (levelset.reshape(-1, 1), p1)
fields["front"] = (front.astype(np.float64).reshape(-1, 1), p1)
problem.write_mig(output_dir, time, fields)
def run(output_dir, final_time=0.15, viscosity=0.1, relaxation=0.15):
if final_time <= 0.0:
raise ValueError("final_time must be strictly positive")
shutil.rmtree(output_dir, ignore_errors=True)
os.makedirs(output_dir)
mesh_file = os.path.join(output_dir, "xmesh_ellipse.msh")
subprocess.run(["gmsh", "-2", os.path.join(HERE, "xmesh_ellipse.geo"),
"-o", mesh_file], check=True, stdout=subprocess.DEVNULL)
problem = fluid.FluidProblem2(
np.zeros(2), viscosity, RHO, sigma=SIGMA,
solver=fluid.ILinearSystem.create_pardiso(iterative_refinement=5),
temporal=True, advection=True,
pressure_enrichment=True, surface_tension_model="membrane")
gmsh_io.load_msh(problem, mesh_file)
elements = problem.elements().copy()
boundary_nodes = np.unique(np.concatenate(
[edge.ravel() for edge in problem.mesh_boundaries().values()]))
mesh = XMeshTopology(elements, problem.n_nodes(), boundary_nodes)
x0 = problem.coordinates().copy()
coordinates, front, levelset = initialize_front(mesh, x0)
phases = classify_elements(
mesh, levelset, front, initial_phases(mesh, coordinates))
if not interface_is_conforming(mesh, phases, front):
raise RuntimeError("initial X-Mesh interface is not conforming")
levelset = redistance_levelset(
mesh, coordinates, phases, front)
problem.set_coordinates(coordinates)
problem.mesh_velocity().fill(0.0)
problem.set_materials(phases)
problem.set_mean_pressure(0.0)
problem.set_strong_boundary("Boundary", velocity=[0.0, 0.0])
initial_interface, current_edges = interface_points(
elements, phases, coordinates)
initial_area = polygon_area(initial_interface)
initial_center = initial_interface.mean(axis=0)
initial_radius = np.linalg.norm(
initial_interface - initial_center, axis=1)
initial_radius_variation = initial_radius.std() / initial_radius.mean()
target_radius = np.sqrt(initial_area / np.pi)
triangle = coordinates[elements, :2]
signed_area = ((triangle[:, 1, 0] - triangle[:, 0, 0]) *
(triangle[:, 2, 1] - triangle[:, 0, 1]) -
(triangle[:, 1, 1] - triangle[:, 0, 1]) *
(triangle[:, 2, 0] - triangle[:, 0, 0]))
reference_sign = np.sign(np.median(signed_area))
initial_min_area = np.min(signed_area * reference_sign)
initial_quality = mesh_shape_quality(coordinates, elements)
write_state(problem, output_dir, 0.0, levelset, front)
capillary_dt = capillary_timestep(coordinates, current_edges)
number_of_steps = max(1, int(np.ceil(final_time / capillary_dt)))
dt = final_time / number_of_steps
maximum_vmax = 0.0
minimum_quality = initial_quality
phase_changes = 0
for step in range(number_of_steps):
time = (step + 1) * dt
problem.implicit_euler(dt)
write_state(problem, output_dir, time, levelset, front)
velocity = problem.velocity().copy()
maximum_vmax = max(
maximum_vmax, np.linalg.norm(velocity, axis=1).max())
if step + 1 == number_of_steps:
break
# Prepare the geometry and pressure space for the next fluid solve.
trial_levelset = advect_levelset(
mesh, coordinates, levelset, velocity, dt)
trial_coordinates, moved_front = move_front_python(
mesh, front, coordinates, levelset, trial_levelset)
trial_levelset[moved_front] = 0.0
trial_phases = classify_elements(
mesh, trial_levelset, moved_front, phases)
trial_front = _closed_interface_front(mesh, trial_phases)
if trial_front is None:
raise RuntimeError(
f"invalid phase-interface topology at fixed step {step + 1}")
new_front = trial_front & ~moved_front
trial_coordinates = snap_new_front_nodes(
mesh, new_front, trial_front, trial_coordinates, trial_levelset)
# The final front must be known before relaxation: newly activated
# interface vertices stay fitted, while vertices which just left the
# front are the only former interface vertices allowed to return to x0.
trial_coordinates = relax(
trial_front, trial_coordinates, x0, relaxation, elements,
reference_sign, 0.08 * initial_min_area,
0.12 * initial_quality)
trial_coordinates[mesh.boundary] = x0[mesh.boundary]
# Redistance after the complete ALE update so the next fluid solve
# and its export use the same geometry, front and element phases.
trial_levelset = redistance_levelset(
mesh, trial_coordinates, trial_phases, trial_front)
area_quality = mesh_quality(
trial_coordinates, elements, reference_sign)
shape_quality = mesh_shape_quality(trial_coordinates, elements)
conforming = interface_is_conforming(
mesh, trial_phases, trial_front)
phase_consistent = interface_phases_match_levelset(
mesh, trial_phases, trial_front, trial_levelset)
# A fixed-connectivity X-Mesh transition may briefly create a sliver
# while one interface vertex hands the front to its neighbour. It is
# admissible as long as no element is inverted or numerically singular.
if not (area_quality > 1e-4 * initial_min_area and
shape_quality > 1e-3 * initial_quality and conforming and
phase_consistent):
triangle_points = trial_coordinates[elements, :2]
edge2 = np.sum((triangle_points -
np.roll(triangle_points, -1, axis=1)) ** 2,
axis=2)
twice_area = np.abs(
(triangle_points[:, 1, 0] - triangle_points[:, 0, 0]) *
(triangle_points[:, 2, 1] - triangle_points[:, 0, 1]) -
(triangle_points[:, 1, 1] - triangle_points[:, 0, 1]) *
(triangle_points[:, 2, 0] - triangle_points[:, 0, 0]))
local_shape = (2.0 * np.sqrt(3.0) * twice_area /
np.maximum(np.sum(edge2, axis=1), 1e-30))
worst = int(np.argmin(local_shape))
raise RuntimeError(
"invalid X-Mesh geometry with the fixed capillary timestep: "
f"step={step + 1}, dt={dt:.6g}, "
f"area_ratio={area_quality/initial_min_area:.3e}, "
f"shape_quality={shape_quality:.3e}, conforming={conforming}, "
f"phase_consistent={phase_consistent}, "
f"front_nodes={np.count_nonzero(trial_front)}, "
f"worst_triangle={elements[worst].tolist()}, "
f"front_flags={trial_front[elements[worst]].tolist()}, "
f"coordinates={triangle_points[worst].tolist()}")
old_coordinates = coordinates
coordinates = trial_coordinates
front = trial_front
levelset = trial_levelset
if not np.array_equal(trial_phases, phases):
phase_changes += np.count_nonzero(trial_phases != phases)
phases = trial_phases
problem.set_coordinates(coordinates)
problem.mesh_velocity().set_from_host(np.ascontiguousarray(
(coordinates[:, :2] - old_coordinates[:, :2]) / dt))
# The new interfaces rebuild the pressure DOFs.
problem.set_materials(phases)
else:
problem.set_coordinates(coordinates)
problem.mesh_velocity().set_from_host(np.ascontiguousarray(
(coordinates[:, :2] - old_coordinates[:, :2]) / dt))
minimum_quality = min(minimum_quality, shape_quality)
final_interface, _ = interface_points(elements, phases, coordinates)
center = final_interface.mean(axis=0)
radius = np.linalg.norm(final_interface - center, axis=1)
final_area = polygon_area(final_interface)
result = {
"steps": number_of_steps,
"time": final_time,
"phase_changes": int(phase_changes),
"dt": float(dt),
"capillary_dt_limit": float(capillary_dt),
"initial_radius_variation": float(initial_radius_variation),
"radius_variation": float(radius.std() / radius.mean()),
"relative_area_change": float(abs(final_area - initial_area) /
initial_area),
"minimum_mesh_quality": float(minimum_quality),
"max_velocity": float(maximum_vmax),
}
np.savetxt(os.path.join(output_dir, "interface_initial.txt"),
initial_interface)
np.savetxt(os.path.join(output_dir, "interface_final.txt"), final_interface)
with open(os.path.join(output_dir, "summary_xmesh.json"), "w",
encoding="utf-8") as stream:
json.dump(result, stream, indent=2)
plot_interfaces(initial_interface, final_interface, target_radius,
os.path.join(output_dir, "ellipse_xmesh.png"))
if final_time >= 0.1 and phase_changes == 0:
raise AssertionError("the test did not exercise dynamic phase remapping")
if final_time <= 0.2 + 1e-14 and result["relative_area_change"] > 0.03:
raise AssertionError(result)
if (number_of_steps > 1 and
result["radius_variation"] >= result["initial_radius_variation"]):
raise AssertionError("the X-Mesh interface did not become more circular")
return result
def main():
parser = argparse.ArgumentParser()
parser.add_argument("output", nargs="?", default="output_xmesh_ellipse")
parser.add_argument("--final-time", type=float, default=0.15)
parser.add_argument("--viscosity", type=float, default=0.1)
parser.add_argument("--relaxation", type=float, default=0.15)
args = parser.parse_args()
result = run(os.path.abspath(args.output), args.final_time,
args.viscosity, args.relaxation)
print("X-Mesh ellipse:", result)
if __name__ == "__main__":
main()