diff --git a/iga/use_cases/turek_2d_sbm/FluidMaterials.json b/iga/use_cases/turek_2d_sbm/FluidMaterials.json new file mode 100644 index 00000000..8ff36a9c --- /dev/null +++ b/iga/use_cases/turek_2d_sbm/FluidMaterials.json @@ -0,0 +1,34 @@ +{ + "properties": [ + { + "model_part_name": "IgaModelPart", + "properties_id": 1, + "Material": { + "constitutive_law": { + "name": "Newtonian2DLaw" + }, + "Variables": { + "DENSITY": 1000.0, + "DYNAMIC_VISCOSITY": 1.0, + "PENALTY_FACTOR": 0.0 + }, + "Tables": null + } + }, + { + "model_part_name": "IgaModelPart.SBM_Support_inner", + "properties_id": 2, + "Material": { + "constitutive_law": { + "name": "Newtonian2DLaw" + }, + "Variables": { + "DENSITY": 1000.0, + "DYNAMIC_VISCOSITY": 1.0, + "PENALTY_FACTOR": 0.0 + }, + "Tables": null + } + } + ] +} diff --git a/iga/use_cases/turek_2d_sbm/ProjectParameters_2D_fluid.json b/iga/use_cases/turek_2d_sbm/ProjectParameters_2D_fluid.json new file mode 100644 index 00000000..302feb7b --- /dev/null +++ b/iga/use_cases/turek_2d_sbm/ProjectParameters_2D_fluid.json @@ -0,0 +1,203 @@ +{ + "analysis_stage": "KratosMultiphysics.FluidDynamicsApplication.fluid_dynamics_analysis", + "problem_data": { + "problem_name": "turek_sbm_single_patch", + "echo_level": 1, + "parallel_type": "OpenMP", + "start_time": 0.0, + "end_time": 10.0 + }, + "solver_settings": { + "solver_type": "monolithic_iga", + "time_scheme": "bdf2_higher_order_vms", + "analysis_type": "non_linear", + "model_part_name": "IgaModelPart", + "domain_size": 2, + "model_import_settings": { + "input_type": "use_input_model_part" + }, + "material_import_settings": { + "materials_filename": "FluidMaterials.json" + }, + // "linear_solver_settings": { + // "solver_type": "skyline_lu_factorization" + // }, + "linear_solver_settings": { "solver_type": "bicgstab", "tolerance": 1.0e-14, "max_iteration": 1000, "scaling": true, "preconditioner_type": "ilu0"}, + "echo_level": 1, + "enforce_element_and_conditions_replacement": false, + "compute_reactions": false, + "maximum_iterations": 50, + "relative_velocity_tolerance": 1e-11, + "absolute_velocity_tolerance": 1e-11, + "relative_pressure_tolerance": 1e-9, + "absolute_pressure_tolerance": 1e-9, + "time_stepping": { + "time_step": 0.005, + "automatic_time_step": false + }, + "skip_entities_replace_and_check": true + }, + "modelers": [ + { + "modeler_name": "ImportNurbsSbmModeler", + "Parameters": { + "input_filename": "turek_nurbs.json", + "model_part_name": "initial_skin_model_part_in", + "link_layer_to_condition_name": [ + { + "layer_name": "Layer0", + "condition_name": "SbmFluidConditionDirichlet" + } + ] + } + }, + { + "modeler_name": "NurbsGeometryModelerSbm", + "Parameters": { + "model_part_name": "IgaModelPart", + "lower_point_xyz": [0.0, 0.0, 0.0], + "upper_point_xyz": [2.5, 0.41, 0.0], + "lower_point_uvw": [0.0, 0.0, 0.0], + "upper_point_uvw": [2.5, 0.41, 0.0], + "polynomial_order": [1, 1], + // Roughly square background cells for a 2.5 x 0.41 box. + // "number_of_knot_spans": [25, 4], // m0 + // "number_of_knot_spans": [50, 8], // m1 + // "number_of_knot_spans": [100, 16], // m2 + // "number_of_knot_spans": [200, 33], // m3 + "number_of_knot_spans": [400, 66], // m4 + // "number_of_knot_spans": [800, 131], // m5 + // "number_of_knot_spans": [1600, 262], // m6 + "number_initial_points_if_importing_nurbs": 1000, + "lambda_outer": 0.5, + "lambda_inner": 0.0, + "number_of_inner_loops": 1, + "skin_model_part_inner_initial_name": "initial_skin_model_part_in", + "skin_model_part_name": "skin_model_part" + } + }, + { + "modeler_name": "IgaModelerSbm", + "Parameters": { + "echo_level": 0, + "skin_model_part_name": "skin_model_part", + "analysis_model_part_name": "IgaModelPart", + "integrate_on_true_boundary": true, + "element_condition_list": [ + { + "geometry_type": "GeometrySurface", + "iga_model_part": "FluidDomain", + "type": "element", + "name": "NavierStokesElement", + // "name": "StokesElement", + "shape_function_derivatives_order": 3 + }, + { + "geometry_type": "SurfaceEdge", + "iga_model_part": "SBM_Support_outer", + "type": "condition", + "name": "SupportFluidCondition", + "shape_function_derivatives_order": 2, + "brep_ids": [2, 4] + }, + { + "geometry_type": "SurfaceEdge", + "iga_model_part": "Inlet_Velocity", + "type": "condition", + "name": "SupportFluidCondition", + "shape_function_derivatives_order": 2, + "brep_ids": [5] + }, + { + "geometry_type": "SurfaceEdge", + "iga_model_part": "Outlet_Pressure", + "type": "condition", + "name": "SupportPressureCondition", + "shape_function_derivatives_order": 2, + "brep_ids": [3] + }, + { + "boundary_role": "internal", + "geometry_type": "SurfaceEdge", + "iga_model_part": "SBM_Support_inner", + "type": "condition", + "name": "SbmCondition", + "shape_function_derivatives_order": 7, + "sbm_parameters": { + "is_inner": true + } + } + ] + } + } + ], + "processes": { + "additional_processes": [ + { + "python_module": "assign_iga_external_conditions_process", + "kratos_module": "KratosMultiphysics.IgaApplication", + "process_name": "AssignIgaExternalConditionsProcess", + "Parameters": { + "model_part_name": "IgaModelPart", + "echo_level": 0, + "element_condition_list": [ + { + "iga_model_part": "IgaModelPart.FluidDomain", + "variables": [ + { + "variable_name": "BODY_FORCE", + "value": ["0.0", "0.0"] + } + ] + }, + { + "iga_model_part": "IgaModelPart.SBM_Support_outer", + "variables": [ + { + "variable_name": "VELOCITY", + "value": ["0.0", "0.0"] + } + ] + }, + { + "iga_model_part": "IgaModelPart.Inlet_Velocity", + "variables": [ + { + "variable_name": "VELOCITY", + "value": [ + "0.5*(1 - cos(3.141592653589793 * (t + 2.0 - sqrt((t - 2.0)*(t - 2.0))) / 4.0)) * (0.2 * 1.5 * 4 * y * (0.41 - y) / (0.41 * 0.41))", + "0.0" + ] + } + ] + }, + { + "iga_model_part": "IgaModelPart.Outlet_Pressure", + "variables": [ + { + "variable_name": "PRESSURE", + "value": "0.0" + } + ] + } + ] + } + } + ], + "dirichlet_process_list": [ + { + "kratos_module": "KratosMultiphysics", + "python_module": "assign_vector_variable_to_nodes_process", + "Parameters": { + "model_part_name": "skin_model_part.inner", + "variable_name": "VELOCITY", + "value": ["0.0", "0.0", "0.0"] + } + } + ], + "constraints_process_list": [] + }, + "output_processes": { + "output_process_list": [] + } +} diff --git a/iga/use_cases/turek_2d_sbm/README.md b/iga/use_cases/turek_2d_sbm/README.md new file mode 100644 index 00000000..73970866 --- /dev/null +++ b/iga/use_cases/turek_2d_sbm/README.md @@ -0,0 +1,53 @@ +# 2D Turek SBM Example + +Single-patch SBM/IGA version of the classical 2D Turek cylinder benchmark. + +This example uses: +- the unordered Turek NURBS skin for the immersed cylinder +- a single Cartesian background patch generated with `NurbsGeometryModelerSbm` +- Turek-style boundary conditions: ramped parabolic inlet, no-slip walls and cylinder, zero outlet pressure + +## Files + +- `ProjectParameters_2D_fluid.json`: fluid setup and boundary conditions +- `FluidMaterials.json`: 2D Newtonian material +- `turek_nurbs_unordered.json`: immersed cylinder skin +- `run_and_post_nurbs.py`: run the case and create the plots and GIFs +- `plot_geometry.py`: plot the background mesh, skin boundary, and surrogate boundary + +## Python Dependencies + +The post-processing scripts require: +- `numpy` +- `matplotlib` +- `imageio` + +If one of them is missing, the script stops immediately with a clear error message. + +If a `latex` executable is available, matplotlib uses it for text rendering. Otherwise the scripts fall back to Computer Modern mathtext. + +## Run + +```bash +cd /home/nantonelli/Examples/iga/use_cases/turek_2d_sbm +python3 run_and_post_nurbs.py +``` + +To plot only the geometry: + +```bash +cd /home/nantonelli/Examples/iga/use_cases/turek_2d_sbm +python3 plot_geometry.py +``` + +## Outputs + +- `final_step_fields.png` +- `velocity.gif` +- `pressure.gif` +- `geometry_plot.png` + +Internal frame caches are written into: +- `_frame_cache` +- `frames_velocity` +- `frames_pressure` diff --git a/iga/use_cases/turek_2d_sbm/plot_geometry.py b/iga/use_cases/turek_2d_sbm/plot_geometry.py new file mode 100644 index 00000000..2bc3ef7c --- /dev/null +++ b/iga/use_cases/turek_2d_sbm/plot_geometry.py @@ -0,0 +1,323 @@ +#!/usr/bin/env python3 +"""Plot the 2D Turek SBM geometry. + +This script visualizes: +- the background Cartesian mesh +- the skin boundary stored in `skin_model_part` +- the surrogate inner boundary + +It only executes the modelers. It does not run the fluid solve. +""" + +from __future__ import annotations + +import importlib +import os +import shutil +from pathlib import Path + +import KratosMultiphysics +import KratosMultiphysics.IgaApplication # registers IGA modelers + +try: + import numpy as np +except ModuleNotFoundError as exc: + raise SystemExit( + "Missing dependency: numpy\n" + "Install it in the current Python environment before running this script." + ) from exc + +try: + import matplotlib +except ModuleNotFoundError as exc: + raise SystemExit( + "Missing dependency: matplotlib\n" + "Install it in the current Python environment before running this script." + ) from exc + +if not (os.environ.get("DISPLAY") or os.environ.get("WAYLAND_DISPLAY")): + matplotlib.use("Agg") + +try: + import matplotlib.font_manager as font_manager + import matplotlib.pyplot as plt +except ModuleNotFoundError as exc: + raise SystemExit( + "matplotlib is installed but incomplete in this environment.\n" + "Please install a full matplotlib package." + ) from exc + + +PROJECT_PARAMETERS_FILENAME = "ProjectParameters_2D_fluid.json" +OUTPUT_FILENAME = "geometry_plot.png" +SHOW_PLOT = bool(os.environ.get("DISPLAY") or os.environ.get("WAYLAND_DISPLAY")) +MAX_SKIN_SEGMENTS_TO_PLOT = 12000 + + +def _configure_matplotlib() -> None: + def _font_available(font_name: str) -> bool: + try: + font_manager.findfont( + font_manager.FontProperties(family=font_name), + fallback_to_default=False, + ) + return True + except Exception: + return False + + serif_candidates = [ + "Computer Modern Roman", + "CMU Serif", + "DejaVu Serif", + "Times New Roman", + "Times", + ] + serif_font = next((name for name in serif_candidates if _font_available(name)), None) + use_tex = shutil.which("latex") is not None + + rc_params = { + "text.usetex": use_tex, + "font.family": "serif", + } + if serif_font: + rc_params["font.serif"] = [serif_font] + if not use_tex and serif_font in {"Computer Modern Roman", "CMU Serif"}: + rc_params["mathtext.fontset"] = "cm" + plt.rcParams.update(rc_params) + + +_configure_matplotlib() + + +def _require_file(path: Path, description: str) -> None: + if not path.is_file(): + raise FileNotFoundError(f"Required {description} was not found: {path}") + + +def _load_parameters(path: Path) -> KratosMultiphysics.Parameters: + _require_file(path, "project parameter file") + with path.open("r", encoding="utf-8") as parameter_file: + return KratosMultiphysics.Parameters(parameter_file.read()) + + +def _normalize_solver_settings(parameters: KratosMultiphysics.Parameters) -> None: + if not parameters.Has("solver_settings"): + raise RuntimeError('Project parameters must contain "solver_settings".') + + solver_settings = parameters["solver_settings"] + if solver_settings.Has("solver_type"): + solver_settings["solver_type"].SetString("monolithic_iga") + else: + solver_settings.AddEmptyValue("solver_type").SetString("monolithic_iga") + + if solver_settings.Has("skip_entities_replace_and_check"): + solver_settings["skip_entities_replace_and_check"].SetBool(True) + else: + solver_settings.AddEmptyValue("skip_entities_replace_and_check").SetBool(True) + + if solver_settings.Has("volume_model_part_name"): + solver_settings["volume_model_part_name"].SetString("IgaModelPart") + else: + solver_settings.AddEmptyValue("volume_model_part_name").SetString("IgaModelPart") + + if solver_settings.Has("skin_parts"): + solver_settings.RemoveValue("skin_parts") + if solver_settings.Has("no_skin_parts"): + solver_settings.RemoveValue("no_skin_parts") + + +def _execute_modelers(parameters: KratosMultiphysics.Parameters, model: KratosMultiphysics.Model) -> None: + if not model.HasModelPart("IgaModelPart"): + model.CreateModelPart("IgaModelPart") + + for modeler_data in parameters["modelers"].values(): + if not modeler_data.Has("modeler_name"): + continue + + modeler_name = modeler_data["modeler_name"].GetString() + modeler_parameters = modeler_data["Parameters"] + + if KratosMultiphysics.HasModeler(modeler_name): + modeler = KratosMultiphysics.CreateModeler(modeler_name, model, modeler_parameters) + else: + if not modeler_data.Has("kratos_module"): + raise RuntimeError(f'Python modeler "{modeler_name}" is missing "kratos_module".') + + kratos_module_name = modeler_data["kratos_module"].GetString() + if not kratos_module_name.startswith("KratosMultiphysics"): + kratos_module_name = "KratosMultiphysics." + kratos_module_name + python_modeler_name = kratos_module_name + ".modelers." + modeler_name + python_module = importlib.import_module(python_modeler_name) + modeler = python_module.Factory(model, modeler_parameters) + + setup_geometry = getattr(modeler, "SetupGeometryModel", None) + if callable(setup_geometry): + setup_geometry() + + prepare_geometry = getattr(modeler, "PrepareGeometryModel", None) + if callable(prepare_geometry): + prepare_geometry() + + setup_model_part = getattr(modeler, "SetupModelPart", None) + if callable(setup_model_part): + setup_model_part() + + +def _read_background_bounds(parameters: KratosMultiphysics.Parameters) -> tuple[float, float, float, float]: + for modeler_data in parameters["modelers"].values(): + if not modeler_data.Has("modeler_name"): + continue + if modeler_data["modeler_name"].GetString() != "NurbsGeometryModelerSbm": + continue + + modeler_parameters = modeler_data["Parameters"] + lower = modeler_parameters["lower_point_xyz"] + upper = modeler_parameters["upper_point_xyz"] + return ( + lower[0].GetDouble(), + upper[0].GetDouble(), + lower[1].GetDouble(), + upper[1].GetDouble(), + ) + + raise RuntimeError('Could not find "NurbsGeometryModelerSbm" bounds in project parameters.') + + +def _knot_lines(model_part) -> tuple[np.ndarray, np.ndarray]: + if not model_part.Has(KratosMultiphysics.KratosGlobals.GetVariable("KNOT_VECTOR_U")): + raise RuntimeError("IgaModelPart is missing KNOT_VECTOR_U.") + if not model_part.Has(KratosMultiphysics.KratosGlobals.GetVariable("KNOT_VECTOR_V")): + raise RuntimeError("IgaModelPart is missing KNOT_VECTOR_V.") + + knot_vector_u = np.array(model_part.GetValue(KratosMultiphysics.KratosGlobals.GetVariable("KNOT_VECTOR_U")), dtype=float) + knot_vector_v = np.array(model_part.GetValue(KratosMultiphysics.KratosGlobals.GetVariable("KNOT_VECTOR_V")), dtype=float) + return np.unique(np.round(knot_vector_u, decimals=12)), np.unique(np.round(knot_vector_v, decimals=12)) + + +def _collect_surrogate_segments(model) -> list[np.ndarray]: + candidate_names = ( + "IgaModelPart.surrogate_inner", + "IgaModelPart.SBM_Support_inner", + ) + + segments = [] + seen = set() + for name in candidate_names: + if not model.HasModelPart(name): + continue + + for condition in model[name].Conditions: + geometry = condition.GetGeometry() + if geometry.PointsNumber() != 2: + continue + + segment = np.array( + [[geometry[0].X, geometry[0].Y], [geometry[1].X, geometry[1].Y]], + dtype=float, + ) + key = tuple(np.round(segment.reshape(-1), 12)) + if key in seen: + continue + seen.add(key) + segments.append(segment) + + return segments + + +def _collect_skin_boundary_segments(model) -> list[np.ndarray]: + candidate_names = ( + "skin_model_part.inner", + "skin_model_part.inner.Layer0", + ) + + segments = [] + seen = set() + for name in candidate_names: + if not model.HasModelPart(name): + continue + + model_part = model[name] + for condition in model_part.Conditions: + geometry = condition.GetGeometry() + if geometry.PointsNumber() != 2: + continue + + segment = np.array( + [[geometry[0].X, geometry[0].Y], [geometry[1].X, geometry[1].Y]], + dtype=float, + ) + p0 = tuple(np.round(segment[0], 12)) + p1 = tuple(np.round(segment[1], 12)) + key = tuple(sorted((p0, p1))) + if key in seen: + continue + seen.add(key) + segments.append(segment) + + if len(segments) > MAX_SKIN_SEGMENTS_TO_PLOT: + stride = int(np.ceil(len(segments) / MAX_SKIN_SEGMENTS_TO_PLOT)) + segments = segments[::stride] + + return segments + + +def main() -> int: + script_dir = Path(__file__).resolve().parent + project_parameters_path = script_dir / PROJECT_PARAMETERS_FILENAME + output_path = script_dir / OUTPUT_FILENAME + + parameters = _load_parameters(project_parameters_path) + _normalize_solver_settings(parameters) + + model = KratosMultiphysics.Model() + _execute_modelers(parameters, model) + + if not model.HasModelPart("IgaModelPart"): + raise RuntimeError('Modeler execution did not create "IgaModelPart".') + + iga_model_part = model["IgaModelPart"] + x_min, x_max, y_min, y_max = _read_background_bounds(parameters) + knot_lines_u, knot_lines_v = _knot_lines(iga_model_part) + skin_boundary_segments = _collect_skin_boundary_segments(model) + surrogate_segments = _collect_surrogate_segments(model) + + fig, ax = plt.subplots(figsize=(9, 4.8)) + + for x_value in knot_lines_u: + ax.plot([x_value, x_value], [y_min, y_max], color="#d9d9d9", linewidth=0.55, zorder=1) + for y_value in knot_lines_v: + ax.plot([x_min, x_max], [y_value, y_value], color="#d9d9d9", linewidth=0.55, zorder=1) + + for segment in skin_boundary_segments: + ax.plot(segment[:, 0], segment[:, 1], color="#1f77b4", linewidth=0.7, alpha=0.7, zorder=2) + + for segment in surrogate_segments: + ax.plot(segment[:, 0], segment[:, 1], color="#d62728", linewidth=1.2, zorder=3) + + ax.set_title(r"2D Turek geometry: mesh, skin boundary, surrogate boundary") + ax.set_xlabel(r"$x$") + ax.set_ylabel(r"$y$") + ax.set_aspect("equal") + ax.set_xlim(x_min, x_max) + ax.set_ylim(y_min, y_max) + ax.grid(False) + + legend_handles = [ + plt.Line2D([0], [0], color="#d9d9d9", linewidth=1.0, label="Background mesh"), + plt.Line2D([0], [0], color="#1f77b4", linewidth=0.7, alpha=0.7, label="Skin boundary"), + plt.Line2D([0], [0], color="#d62728", linewidth=1.2, label="Surrogate boundary"), + ] + ax.legend(handles=legend_handles, loc="upper right", frameon=True) + + fig.tight_layout() + fig.savefig(output_path, dpi=180) + if SHOW_PLOT: + plt.show() + plt.close(fig) + + print(f"Saved {output_path}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/iga/use_cases/turek_2d_sbm/run_and_post_nurbs.py b/iga/use_cases/turek_2d_sbm/run_and_post_nurbs.py new file mode 100644 index 00000000..955b7dc3 --- /dev/null +++ b/iga/use_cases/turek_2d_sbm/run_and_post_nurbs.py @@ -0,0 +1,636 @@ +#!/usr/bin/env python3 +"""Run the 2D Turek SBM example and post-process the solution. + +This example produces: +- one static image of the final step +- one GIF for velocity magnitude +- one GIF for pressure + +Required third-party Python packages: +- numpy +- matplotlib +- imageio + +If a LaTeX installation is available, matplotlib uses it for text rendering. +Otherwise the script falls back to Computer Modern mathtext. +""" + +from __future__ import annotations + +import importlib +import os +import shutil +import sys +import time +from pathlib import Path + +import KratosMultiphysics +import KratosMultiphysics.IgaApplication # registers IGA modelers and conditions + +try: + import numpy as np +except ModuleNotFoundError as exc: + raise SystemExit( + "Missing dependency: numpy\n" + "Install it in the current Python environment before running this example." + ) from exc + +try: + import matplotlib +except ModuleNotFoundError as exc: + raise SystemExit( + "Missing dependency: matplotlib\n" + "Install it in the current Python environment before running this example." + ) from exc + +if not (os.environ.get("DISPLAY") or os.environ.get("WAYLAND_DISPLAY")): + matplotlib.use("Agg") + +try: + import matplotlib.font_manager as font_manager + import matplotlib.pyplot as plt + from matplotlib.colors import Normalize +except ModuleNotFoundError as exc: + raise SystemExit( + "matplotlib is installed but incomplete in this environment.\n" + "Please install a full matplotlib package." + ) from exc + +try: + import imageio.v2 as imageio +except ModuleNotFoundError as exc: + raise SystemExit( + "Missing dependency: imageio\n" + "Install it in the current Python environment before running this example." + ) from exc + +PROJECT_PARAMETERS_FILENAME = "ProjectParameters_2D_fluid.json" +FINAL_SHOT_FILENAME = "final_step_fields.png" +VELOCITY_GIF_FILENAME = "velocity.gif" +PRESSURE_GIF_FILENAME = "pressure.gif" +FRAME_DATA_DIRNAME = "_frame_cache" +VELOCITY_FRAME_DIRNAME = "frames_velocity" +PRESSURE_FRAME_DIRNAME = "frames_pressure" + +FRAME_EVERY = 1 +GIF_DURATION = 0.05 +GIF_PALETTE_SIZE = 64 +MIN_VELOCITY_COLOR_MAX = 0.3 +FRAME_DPI = 120 +def _configure_matplotlib() -> None: + def _font_available(font_name: str) -> bool: + try: + font_manager.findfont( + font_manager.FontProperties(family=font_name), + fallback_to_default=False, + ) + return True + except Exception: + return False + + serif_candidates = [ + "Computer Modern Roman", + "CMU Serif", + "DejaVu Serif", + "Times New Roman", + "Times", + ] + serif_font = next((name for name in serif_candidates if _font_available(name)), None) + use_tex = shutil.which("latex") is not None + + rc_params = { + "text.usetex": use_tex, + "font.family": "serif", + } + if serif_font: + rc_params["font.serif"] = [serif_font] + if not use_tex and serif_font in {"Computer Modern Roman", "CMU Serif"}: + rc_params["mathtext.fontset"] = "cm" + + plt.rcParams.update(rc_params) + + +_configure_matplotlib() + + +def _require_file(path: Path, description: str) -> None: + if not path.is_file(): + raise FileNotFoundError(f"Required {description} was not found: {path}") + + +def _load_parameters(path: Path) -> KratosMultiphysics.Parameters: + _require_file(path, "project parameter file") + with path.open("r", encoding="utf-8") as parameter_file: + return KratosMultiphysics.Parameters(parameter_file.read()) + + +def _load_analysis_stage_class(parameters: KratosMultiphysics.Parameters): + module_name = parameters["analysis_stage"].GetString() + class_name = "".join(part.title() for part in module_name.split(".")[-1].split("_")) + module = importlib.import_module(module_name) + return getattr(module, class_name) + + +def _normalize_solver_settings(parameters: KratosMultiphysics.Parameters) -> None: + if not parameters.Has("solver_settings"): + raise RuntimeError('Project parameters must contain "solver_settings".') + + solver_settings = parameters["solver_settings"] + + if solver_settings.Has("solver_type"): + solver_settings["solver_type"].SetString("monolithic_iga") + else: + solver_settings.AddEmptyValue("solver_type").SetString("monolithic_iga") + + if solver_settings.Has("time_scheme"): + if solver_settings["time_scheme"].GetString() == "bdf2": + solver_settings["time_scheme"].SetString("bdf2_higher_order_vms") + else: + solver_settings.AddEmptyValue("time_scheme").SetString("bdf2_higher_order_vms") + + if solver_settings.Has("skip_entities_replace_and_check"): + solver_settings["skip_entities_replace_and_check"].SetBool(True) + else: + solver_settings.AddEmptyValue("skip_entities_replace_and_check").SetBool(True) + + if solver_settings.Has("volume_model_part_name"): + solver_settings["volume_model_part_name"].SetString("IgaModelPart") + else: + solver_settings.AddEmptyValue("volume_model_part_name").SetString("IgaModelPart") + + if solver_settings.Has("skin_parts"): + solver_settings.RemoveValue("skin_parts") + if solver_settings.Has("no_skin_parts"): + solver_settings.RemoveValue("no_skin_parts") + + +def CreateAnalysisStageWithFlushInstance(cls, global_model, parameters): + class AnalysisStageWithFlush(cls): + def __init__(self, model, project_parameters, flush_frequency=10.0): + super().__init__(model, project_parameters) + self.flush_frequency = flush_frequency + self.last_flush = time.time() + sys.stdout.flush() + + def Initialize(self): + super().Initialize() + sys.stdout.flush() + + def InitializeSolutionStep(self): + current_time = self._GetSolver().GetComputingModelPart().ProcessInfo[ + KratosMultiphysics.TIME + ] + if global_model.HasModelPart("skin_model_part"): + global_model["skin_model_part"].ProcessInfo[KratosMultiphysics.TIME] = current_time + super().InitializeSolutionStep() + + def FinalizeSolutionStep(self): + super().FinalizeSolutionStep() + if self.parallel_type == "OpenMP": + now = time.time() + if now - self.last_flush > self.flush_frequency: + sys.stdout.flush() + self.last_flush = now + + return AnalysisStageWithFlush(global_model, parameters) + + +def _sample_element_fields(model_part): + x_coord = [] + y_coord = [] + velx = [] + vely = [] + pressure = [] + + for element in model_part.Elements: + geometry = element.GetGeometry() + center = geometry.Center() + x_coord.append(center.X) + y_coord.append(center.Y) + + n_nodes = len(geometry) + if n_nodes == 0: + velx.append(0.0) + vely.append(0.0) + pressure.append(0.0) + continue + + vx = 0.0 + vy = 0.0 + p = 0.0 + for node in geometry: + vx += node.GetSolutionStepValue(KratosMultiphysics.VELOCITY_X, 0) + vy += node.GetSolutionStepValue(KratosMultiphysics.VELOCITY_Y, 0) + p += node.GetSolutionStepValue(KratosMultiphysics.PRESSURE, 0) + + inv_n = 1.0 / n_nodes + velx.append(vx * inv_n) + vely.append(vy * inv_n) + pressure.append(p * inv_n) + + if not x_coord: + return ( + np.empty((0, 2), dtype=float), + np.empty((0, 2), dtype=float), + np.array([], dtype=float), + ) + + points = np.column_stack((x_coord, y_coord)) + velocity = np.column_stack((velx, vely)) + return points, velocity, np.array(pressure, dtype=float) + + +def _prepare_output_directory(path: Path) -> None: + if path.is_dir(): + shutil.rmtree(path) + path.mkdir(parents=True, exist_ok=True) + + +def _cached_frame_paths(frame_data_dir: Path) -> list[Path]: + return sorted(frame_data_dir.glob("frame_*.npz")) + + +def _write_gif(frame_dir: Path, gif_path: Path) -> None: + frames = sorted(frame_dir.glob("frame_*.png")) + if not frames: + raise RuntimeError(f"No PNG frames were written in {frame_dir}. Cannot create {gif_path.name}.") + + with imageio.get_writer( + gif_path, + mode="I", + duration=GIF_DURATION, + palettesize=GIF_PALETTE_SIZE, + subrectangles=True, + ) as writer: + for frame_path in frames: + writer.append_data(imageio.imread(frame_path)) + + +def _extract_inner_boundary_outline(model) -> np.ndarray | None: + def _is_reasonable_outline(points: np.ndarray) -> bool: + if points.shape[0] < 16: + return False + + center = points.mean(axis=0) + radii = np.linalg.norm(points - center, axis=1) + radii = radii[radii > 1e-12] + if radii.size < 16: + return False + + radius_p10 = float(np.percentile(radii, 10.0)) + radius_p90 = float(np.percentile(radii, 90.0)) + if radius_p10 <= 1e-12: + return False + + # Reject point clouds that contain long connector segments or large radial outliers. + return (radius_p90 / radius_p10) < 1.35 + + candidate_names = ( + "skin_model_part.inner", + "initial_skin_model_part_in", + "IgaModelPart.SBM_Support_inner", + ) + + for name in candidate_names: + if not model.HasModelPart(name): + continue + + model_part = model[name] + node_map = {} + for condition in model_part.Conditions: + for node in condition.GetGeometry(): + node_map[node.Id] = (node.X, node.Y) + if not node_map and model_part.NumberOfNodes() > 0: + for node in model_part.Nodes: + node_map[node.Id] = (node.X, node.Y) + if len(node_map) < 3: + continue + + points = np.array(list(node_map.values()), dtype=float) + if not _is_reasonable_outline(points): + continue + + center = points.mean(axis=0) + angles = np.arctan2(points[:, 1] - center[1], points[:, 0] - center[0]) + return points[np.argsort(angles)] + + return None + + +def _draw_inner_boundary(ax, boundary_outline: np.ndarray | None) -> None: + if boundary_outline is None or boundary_outline.shape[0] < 3: + return + + ax.fill( + boundary_outline[:, 0], + boundary_outline[:, 1], + facecolor="white", + edgecolor="black", + linewidth=1.1, + zorder=4, + ) + + +def _square_marker_size(ax, x_values: np.ndarray, y_values: np.ndarray) -> float: + unique_x = np.unique(np.round(x_values, decimals=12)) + unique_y = np.unique(np.round(y_values, decimals=12)) + + def _spacing(values: np.ndarray) -> float: + if values.size <= 1: + return 1.0 + diffs = np.diff(values) + diffs = diffs[diffs > 1e-12] + return float(np.min(diffs)) if diffs.size > 0 else 1.0 + + dx = _spacing(unique_x) + dy = _spacing(unique_y) + + ax.figure.canvas.draw() + origin = ax.transData.transform((0.0, 0.0)) + dx_px = abs(ax.transData.transform((dx, 0.0))[0] - origin[0]) if dx > 0.0 else 6.0 + dy_px = abs(ax.transData.transform((0.0, dy))[1] - origin[1]) if dy > 0.0 else 6.0 + marker_size_pts = max(dx_px, dy_px) * 72.0 / ax.figure.dpi * 1.08 + return marker_size_pts ** 2 + + +def _plot_field(ax, points, values, title, color_norm, colorbar_label, boundary_outline) -> None: + x_values = points[:, 0] + y_values = points[:, 1] + marker_size = _square_marker_size(ax, x_values, y_values) + + artist = ax.scatter( + x_values, + y_values, + c=values, + cmap="jet", + norm=color_norm, + s=marker_size, + marker="s", + linewidths=0.0, + zorder=3, + ) + ax.set_title(title) + ax.set_xlabel(r"$x$") + ax.set_ylabel(r"$y$") + ax.set_facecolor("white") + ax.set_aspect("equal") + ax.set_xlim(float(x_values.min()), float(x_values.max())) + ax.set_ylim(float(y_values.min()), float(y_values.max())) + _draw_inner_boundary(ax, boundary_outline) + plt.colorbar(artist, ax=ax, shrink=0.92, pad=0.02, label=colorbar_label) + + +def _save_single_field_frame( + frame_path: Path, + points, + values, + title, + color_norm, + colorbar_label, + boundary_outline, +) -> None: + fig, ax = plt.subplots(figsize=(9, 4.8)) + _plot_field(ax, points, values, title, color_norm, colorbar_label, boundary_outline) + fig.tight_layout() + fig.savefig(frame_path, dpi=FRAME_DPI) + plt.close(fig) + + +def _save_final_shot( + figure_path: Path, + points, + velocity_magnitude, + pressure, + velocity_norm, + pressure_norm, + boundary_outline, +) -> None: + fig, axes = plt.subplots(1, 2, figsize=(14, 5)) + _plot_field( + axes[0], + points, + velocity_magnitude, + r"Final velocity magnitude $|\mathbf{v}|$", + velocity_norm, + r"$|\mathbf{v}|$", + boundary_outline, + ) + _plot_field( + axes[1], + points, + pressure, + r"Final pressure $p$", + pressure_norm, + r"$p$", + boundary_outline, + ) + fig.tight_layout() + fig.savefig(figure_path, dpi=160) + plt.close(fig) + + +def _cache_frame_data( + frame_data_dir: Path, + frame_index: int, + time_value: float, + points, + velocity_magnitude, + pressure, +) -> None: + cache_path = frame_data_dir / f"frame_{frame_index:06d}.npz" + np.savez_compressed( + cache_path, + time=np.array([time_value], dtype=float), + points=points, + velocity_magnitude=velocity_magnitude, + pressure=pressure, + ) + + +def main() -> int: + script_dir = Path(__file__).resolve().parent + project_parameters_path = script_dir / PROJECT_PARAMETERS_FILENAME + final_shot_path = script_dir / FINAL_SHOT_FILENAME + velocity_frame_dir = script_dir / VELOCITY_FRAME_DIRNAME + pressure_frame_dir = script_dir / PRESSURE_FRAME_DIRNAME + frame_data_dir = script_dir / FRAME_DATA_DIRNAME + velocity_gif_path = script_dir / VELOCITY_GIF_FILENAME + pressure_gif_path = script_dir / PRESSURE_GIF_FILENAME + + parameters = _load_parameters(project_parameters_path) + _normalize_solver_settings(parameters) + + _prepare_output_directory(frame_data_dir) + _prepare_output_directory(velocity_frame_dir) + _prepare_output_directory(pressure_frame_dir) + + analysis_stage_class = _load_analysis_stage_class(parameters) + global_model = KratosMultiphysics.Model() + simulation = CreateAnalysisStageWithFlushInstance( + analysis_stage_class, + global_model, + parameters, + ) + + simulation.Initialize() + boundary_outline = _extract_inner_boundary_outline(global_model) + if boundary_outline is None: + print( + "Warning: could not reconstruct a clean immersed-boundary outline for plotting. " + "The obstacle fill will be skipped." + ) + + pressure_min = None + pressure_max = None + velocity_max = MIN_VELOCITY_COLOR_MAX + final_frame_path = None + frame_index = 0 + + while simulation.KeepAdvancingSolutionLoop(): + simulation.time = simulation._AdvanceTime() + simulation.InitializeSolutionStep() + simulation._GetSolver().Predict() + simulation._GetSolver().SolveSolutionStep() + simulation.FinalizeSolutionStep() + simulation.OutputSolutionStep() + + if not global_model.HasModelPart("IgaModelPart"): + raise RuntimeError('Expected model part "IgaModelPart" was not created.') + + model_part = global_model["IgaModelPart"] + points, velocity, pressure = _sample_element_fields(model_part) + if points.size == 0: + raise RuntimeError("No element data were found for plotting.") + + velocity_magnitude = np.linalg.norm(velocity, axis=1) + velocity_max = max(velocity_max, float(np.max(velocity_magnitude))) + + current_pressure_min = float(np.min(pressure)) + current_pressure_max = float(np.max(pressure)) + pressure_min = current_pressure_min if pressure_min is None else min(pressure_min, current_pressure_min) + pressure_max = current_pressure_max if pressure_max is None else max(pressure_max, current_pressure_max) + + if frame_index % FRAME_EVERY == 0: + time_value = float(model_part.ProcessInfo[KratosMultiphysics.TIME]) + _cache_frame_data( + frame_data_dir, + frame_index, + time_value, + points, + velocity_magnitude, + pressure, + ) + final_frame_path = frame_data_dir / f"frame_{frame_index:06d}.npz" + + # Write preview frames immediately so users can inspect progress while the solve runs. + preview_velocity_norm = Normalize(vmin=0.0, vmax=velocity_max) + preview_pressure_min = pressure_min if pressure_min is not None else current_pressure_min + preview_pressure_max = pressure_max if pressure_max is not None else current_pressure_max + if abs(preview_pressure_max - preview_pressure_min) < 1e-14: + preview_pressure_min -= 1.0 + preview_pressure_max += 1.0 + preview_pressure_norm = Normalize(vmin=preview_pressure_min, vmax=preview_pressure_max) + + frame_suffix = f"frame_{frame_index:06d}.png" + time_label = rf"$t = {time_value:.3f}$" + _save_single_field_frame( + velocity_frame_dir / frame_suffix, + points, + velocity_magnitude, + rf"Velocity magnitude $|\mathbf{{v}}|$, {time_label}", + preview_velocity_norm, + r"$|\mathbf{v}|$", + boundary_outline, + ) + _save_single_field_frame( + pressure_frame_dir / frame_suffix, + points, + pressure, + rf"Pressure $p$, {time_label}", + preview_pressure_norm, + r"$p$", + boundary_outline, + ) + + frame_index += 1 + + simulation.Finalize() + + if final_frame_path is None or not final_frame_path.is_file(): + raise RuntimeError("No frame data were cached during the simulation.") + + if pressure_min is None or pressure_max is None: + raise RuntimeError("Pressure bounds could not be computed from the simulation output.") + if abs(pressure_max - pressure_min) < 1e-14: + pressure_min -= 1.0 + pressure_max += 1.0 + + velocity_norm = Normalize(vmin=0.0, vmax=velocity_max) + pressure_norm = Normalize(vmin=pressure_min, vmax=pressure_max) + + cached_paths = _cached_frame_paths(frame_data_dir) + if not cached_paths: + raise RuntimeError(f"No cached frame data were found in {frame_data_dir}.") + + print( + "Simulation finished. Rendering " + f"{len(cached_paths)} cached frames for velocity and pressure GIFs..." + ) + + for frame_counter, cache_path in enumerate(cached_paths, start=1): + with np.load(cache_path) as data: + points = data["points"] + velocity_magnitude = data["velocity_magnitude"] + pressure = data["pressure"] + time_value = float(data["time"][0]) + + frame_suffix = cache_path.stem + ".png" + time_label = rf"$t = {time_value:.3f}$" + + _save_single_field_frame( + velocity_frame_dir / frame_suffix, + points, + velocity_magnitude, + rf"Velocity magnitude $|\mathbf{{v}}|$, {time_label}", + velocity_norm, + r"$|\mathbf{v}|$", + boundary_outline, + ) + _save_single_field_frame( + pressure_frame_dir / frame_suffix, + points, + pressure, + rf"Pressure $p$, {time_label}", + pressure_norm, + r"$p$", + boundary_outline, + ) + + if frame_counter % 25 == 0 or frame_counter == len(cached_paths): + print(f"Rendered {frame_counter}/{len(cached_paths)} frames") + + with np.load(final_frame_path) as final_data: + final_points = final_data["points"] + final_velocity_magnitude = final_data["velocity_magnitude"] + final_pressure = final_data["pressure"] + + _save_final_shot( + final_shot_path, + final_points, + final_velocity_magnitude, + final_pressure, + velocity_norm, + pressure_norm, + boundary_outline, + ) + _write_gif(velocity_frame_dir, velocity_gif_path) + _write_gif(pressure_frame_dir, pressure_gif_path) + + print(f"Saved {final_shot_path}") + print(f"Saved {velocity_gif_path}") + print(f"Saved {pressure_gif_path}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/iga/use_cases/turek_2d_sbm/turek_nurbs.json b/iga/use_cases/turek_2d_sbm/turek_nurbs.json new file mode 100644 index 00000000..8df2b198 --- /dev/null +++ b/iga/use_cases/turek_2d_sbm/turek_nurbs.json @@ -0,0 +1,259 @@ +{ + "Surfaces": [ + { + "Weights": [ + 0, + 0, + 0, + 0 + ], + "CPCoordinates": [ + [ + 0.12395495292560721, + 0.144582386160613, + 0.0 + ], + [ + 0.6983853929574487, + 0.144582386160613, + 0.0 + ], + [ + 0.12395495292560721, + 0.2641872773841938, + 0.0 + ], + [ + 0.6983853929574487, + 0.2641872773841938, + 0.0 + ] + ], + "knotVectors": { + "Xi": [ + 0.0, + 0.0, + 1.0, + 1.0 + ], + "Eta": [ + 0.0, + 0.0, + 1.0, + 1.0 + ] + }, + "id": 1, + "pDegreeU": 1, + "pDegreeV": 1 + } + ], + "Lines": [ + { + "Layer": "Layer0", + "Weights": [ + 1, + 1 + ], + "CPCoordinates": [ + [ + 0.6, + 0.2, + 0.0 + ], + [ + 0.6, + 0.21, + 0.0 + ] + ], + "knotVector": [ + 0.0, + 0.0, + 1.0, + 1.0 + ], + "id": 3, + "pDegree": 1 + }, + { + "Layer": "Layer0", + "Weights": [ + 1, + 1 + ], + "CPCoordinates": [ + [ + 0.6, + 0.19, + 0.0 + ], + [ + 0.6, + 0.2, + 0.0 + ] + ], + "knotVector": [ + 0.0, + 0.0, + 1.0, + 1.0 + ], + "id": 7, + "pDegree": 1 + }, + { + "Layer": "Layer0", + "Weights": [ + 1.0, + 0.7745966692414529, + 1.0 + ], + "CPCoordinates": [ + [ + 0.2489897948556646, + 0.20999999999999525, + 0.0 + ], + [ + 0.24082482904639038, + 0.24999999999999997, + 0.0 + ], + [ + 0.2, + 0.25, + 0.0 + ] + ], + "knotVector": [ + 0.0, + 0.0, + 0.0, + 1.0, + 1.0, + 1.0 + ], + "id": 9, + "pDegree": 2 + }, + { + "Layer": "Layer0", + "Weights": [ + 1, + 1 + ], + "CPCoordinates": [ + [ + 0.6, + 0.21, + 0.0 + ], + [ + 0.2489897948556646, + 0.20999999999999525, + 0.0 + ] + ], + "knotVector": [ + 0.0, + 0.0, + 1.0, + 1.0 + ], + "id": 10, + "pDegree": 1 + }, + { + "Layer": "Layer0", + "Weights": [ + 1, + 1 + ], + "CPCoordinates": [ + [ + 0.24898979485848902, + 0.19, + 0.0 + ], + [ + 0.6, + 0.19, + 0.0 + ] + ], + "knotVector": [ + 0.0, + 0.0, + 1.0, + 1.0 + ], + "id": 13, + "pDegree": 1 + }, + { + "Layer": "Layer0", + "Weights": [ + 1.0, + 0.5000000000000001, + 1.0, + 0.5000000000000001, + 1.0, + 0.9870481592612333, + 1.0 + ], + "CPCoordinates": [ + [ + 0.2, + 0.25, + 0.0 + ], + [ + 0.1133974596215562, + 0.25, + 0.0 + ], + [ + 0.1566987298107781, + 0.17500000000000002, + 0.0 + ], + [ + 0.2, + 0.10000000000000006, + 0.0 + ], + [ + 0.24330127018922193, + 0.175, + 0.0 + ], + [ + 0.24736450209470942, + 0.18203772410323926, + 0.0 + ], + [ + 0.2489897948563544, + 0.1900000000033843, + 0.0 + ] + ], + "knotVector": [ + 0.0, + 0.0, + 0.0, + 0.4629844867437035, + 0.4629844867437035, + 0.925968973487407, + 0.925968973487407, + 1.0, + 1.0, + 1.0 + ], + "id": 14, + "pDegree": 2 + } + ] +} \ No newline at end of file