diff --git a/benchmarks/shell-0001/README.md b/benchmarks/shell-0001/README.md new file mode 100644 index 0000000..bac5d19 --- /dev/null +++ b/benchmarks/shell-0001/README.md @@ -0,0 +1,130 @@ +# shell-0001: Clamped Circular Plate, Center Point Load + +This benchmark compares nodal transverse displacement along a radial profile to the closed-form Reissner–Mindlin solution for a **clamped circular plate** under a **center point load**. Results are reported in **Hughes-normalized** form so that the reference profile is well scaled away from the load singularity. + +The benchmark is intended to catch issues in shell/plate bending response, symmetry and clamp boundary handling, center-load application on a mapped quarter-disk mesh, and the mixed interpolation used by `HeterosisPlate` (where displacement sampling must respect connectivity). + +## Model Geometry + +The finite element model uses **quarter symmetry** of a circular disk of radius `RADIUS`. The patch is built with `xara.Model.surface` on a mapped **Q9** control-point layout (`quarter_disk_points` in `test_shell_0001.py`): one element block with `MESH_SUBDIVISIONS_PER_DIRECTION` by `MESH_SUBDIVISIONS_PER_DIRECTION` subdivisions. + +- Symmetry is enforced on the straight edges meeting at the plate center. +- The **clamped** condition is applied on the **outer circular arc**. +- A **quarter** of the total point load is applied at the center node so that the symmetric model represents a total center load `P_TOTAL`. + +## Files + +- `test_shell_0001.py` + Defines the model, boundary conditions, reference solution, and pytest benchmarks. + +## Elements Tested + +The benchmark is parametrized over the following elements (MITC4/ASDShellQ4 and MITC9 use the same `ELEMENT_ORDER` convention as `shell-0002`): + +```python +ELEMENT_ORDER = { + "ASDShellQ4": 1, + "ShellMITC4": 1, + "HeterosisPlate": 2, + "ShellMITC9": 2, +} +``` + +The interpolation used for the underlying shell formulation follows the same pattern as in `shell-0002`: + +| Element | Transverse displacement interpolation | Rotation interpolation | +|---|---:|---:| +| `ASDShellQ4` | Q4 | Q4 | +| `ShellMITC4` | Q4 | Q4 | +| `ShellMITC9` | Q9 | Q9 | +| `HeterosisPlate` | Q8 | Q9 | + +The `HeterosisPlate` case is exercised in a dedicated test because the ninth connectivity node carries the Q9 rotation field but not the Q8 transverse displacement field; nodal samples at **rotation-only** nodes are excluded from the displacement profile comparison. + +## Material, Thickness, and Load Parameters + +The numerical values in the test module are: + +```python +RADIUS = 5.0 +THICKNESS = 2.0 +E_MOD = 10.92e5 +NU = 0.3 +KAPPA = 5.0 / 6.0 + +P_TOTAL = 1.0 +P_QUARTER = -P_TOTAL / 4.0 +``` + +Flexural rigidity is $D = E t^3 / (12(1-\nu^2))$. Mindlin transverse shear stiffness uses $\kappa G t$ with $G = E / (2(1+\nu))$. + +## Reference Solution + +Transverse deflection (positive downward in the closed form used by the test) combines the **Kirchhoff** clamped-plate term with a **logarithmic shear correction** in the Reissner–Mindlin framework; see `reissner_clamped_center_load_w` in `test_shell_0001.py`. + +Nodal results are compared after **Hughes normalization**: + +$$ +w_{\mathrm{norm}} = w_{\mathrm{down}} \left(\frac{16\pi D\, }{P_{\mathrm{total}}\, R^2}\right) +$$ + +where $w_{\mathrm{down}}$ is the magnitude of downward displacement extracted from the model. + +## What the Test Verifies + +After a linear static analysis (`xara.StaticAnalysis` with `system="Umfpack"`), the benchmark checks that: + +- Nodal transverse displacement along radii in the annulus $0.05 < r/R < 0.95$ matches the analytical profile within tolerances. +- Each sample satisfies either a **relative** or **absolute** error bound (whichever is easier to satisfy where the reference is small). +- `HeterosisPlate` results ignore nodes that are not part of the Q8 transverse displacement connectivity. + +This is an **equilibrium / BVP** check, not a pure kinematic patch test. + +## Error Tolerances and Sampling + +Pointwise checks use: + +```python +MAX_POINTWISE_REL_ERROR = 5.0e-2 +MAX_POINTWISE_ABS_ERROR = 1.0e-2 +REFERENCE_VALUE_FLOOR = 5.0e-2 +``` + +Relative error divides the absolute error by $\max(|y_{\mathrm{ref}}|,\ r_{\mathrm{floor}})$ with $r_{\mathrm{floor}}$ equal to the scalar `REFERENCE_VALUE_FLOOR` in the test module, avoiding blow-ups near zero reference values. + +## Running the Test + +From the `shell-0001` folder: + +```bash +python -m pytest test_shell_0001.py -v -rs +``` + +Expected collection (three parametrized cases plus one Heterosis-specific test): + +```text +test_hughes_clamped_plate_pointwise_error[ASDShellQ4-1] +test_hughes_clamped_plate_pointwise_error[ShellMITC4-1] +test_hughes_clamped_plate_pointwise_error[ShellMITC9-2] +test_hughes_clamped_plate_pointwise_error_heterosis +``` + +The parametrized test **excludes** `HeterosisPlate`; that element is covered only by `test_hughes_clamped_plate_pointwise_error_heterosis`, which applies connectivity-aware displacement sampling. + +If an element type is not compiled into the active OpenSees/Xara build, the corresponding case is **skipped** when an “unknown element type” error is raised during model construction. + +## Notes on Interpretation + +A passing result indicates that the mesh density and element formulation are adequate for this classical plate benchmark under the stated tolerances. + +A failure may indicate: + +- A problem in the potentially new plate/shell element tested (current suite is expected to pass). +- For `HeterosisPlate` or similar elements, incorrect handling of rotation-only nodes when sampling $w$. +- An unavailable or changed element implementation in the active OpenSees/Xara build. + +## References + +[1] T. J. R. Hughes, *The Finite Element Method: Linear Static and Dynamic Finite Element Analysis*. Englewood Cliffs, New Jersey: Prentice-Hall, 1987. + +[2] O. C. Zienkiewicz, R. L. Taylor, and S. Govindjee, *The Finite Element Method: Its Basis and Fundamentals*. Elsevier, 8th ed., Nov. 2024. diff --git a/benchmarks/shell-0001/test_shell_0001.py b/benchmarks/shell-0001/test_shell_0001.py new file mode 100644 index 0000000..b84b4c2 --- /dev/null +++ b/benchmarks/shell-0001/test_shell_0001.py @@ -0,0 +1,343 @@ +"""Self-contained regression test for benchmark `shell-0001`. + +The benchmark compares nodal transverse displacement profiles against the +closed-form Reissner-Mindlin solution for a clamped circular plate under a +center point load. + +The HeterosisPlate case is handled separately because its external OpenSees +connectivity uses 9 shell nodes, while the physical plate interpolation uses +Q8 transverse displacement and Q9 rotations. +""" + +from __future__ import annotations + +import math +from typing import Any + +import numpy as np +import pytest +import xara +from opensees.errors import XaraError +from xara.helpers import find_node + + +RADIUS = 5.0 +THICKNESS = 2.0 +E_MOD = 10.92e5 +NU = 0.3 +KAPPA = 5.0 / 6.0 + +P_TOTAL = 1.0 +P_QUARTER = -P_TOTAL / 4.0 + +SECTION_TAG = 1 +LOAD_PATTERN_TAG = 1 +MESH_SUBDIVISIONS_PER_DIRECTION = 12 + +MAX_POINTWISE_REL_ERROR = 5.0e-2 +MAX_POINTWISE_ABS_ERROR = 1.0e-2 +REFERENCE_VALUE_FLOOR = 5.0e-2 + +ELEMENT_ORDER = { + "ASDShellQ4": 1, + "ShellMITC4": 1, + "HeterosisPlate": 2, + "ShellMITC9": 2, +} + + +def bending_rigidity() -> float: + """Return plate flexural rigidity D = E t^3 / (12 * (1 - nu^2)).""" + return E_MOD * THICKNESS**3 / (12.0 * (1.0 - NU**2)) + + +def shear_rigidity() -> float: + """Return Mindlin shear stiffness kappa * G * t.""" + shear_modulus = E_MOD / (2.0 * (1.0 + NU)) + return KAPPA * shear_modulus * THICKNESS + + +def normalized_displacement(w_down: np.ndarray) -> np.ndarray: + """Return Hughes normalized displacement, with w positive downward.""" + return 16.0 * math.pi * bending_rigidity() * w_down / ( + P_TOTAL * RADIUS**2 + ) + + +def reissner_clamped_center_load_w( + radius_from_center: np.ndarray, +) -> np.ndarray: + """Return Reissner-Mindlin displacement for a clamped circular plate.""" + rho = np.asarray(radius_from_center, dtype=float) / RADIUS + rho = np.clip(rho, 1.0e-14, 1.0) + + flexural = bending_rigidity() + shear = shear_rigidity() + + w_kirchhoff = ( + P_TOTAL + * RADIUS**2 + / (16.0 * math.pi * flexural) + * (1.0 - rho**2 + 2.0 * rho**2 * np.log(rho)) + ) + w_shear = P_TOTAL / (2.0 * math.pi * shear) * np.log(1.0 / rho) + + return w_kirchhoff + w_shear + + +def quarter_disk_points() -> dict[int, list[float]]: + """Return Q9 control points for the mapped quarter disk.""" + r = RADIUS + a = math.pi / 8.0 + + return { + 1: [0.0, 0.0, 0.0], + 2: [r, 0.0, 0.0], + 3: [r / math.sqrt(2.0), r / math.sqrt(2.0), 0.0], + 4: [0.0, r, 0.0], + 5: [0.5 * r, 0.0, 0.0], + 6: [r * math.cos(a), r * math.sin(a), 0.0], + 7: [r * math.cos(3.0 * a), r * math.sin(3.0 * a), 0.0], + 8: [0.0, 0.5 * r, 0.0], + 9: [0.5 * r / math.sqrt(2.0), 0.5 * r / math.sqrt(2.0), 0.0], + } + + +def add_to_fixity(fixity: list[int], dofs_1based: tuple[int, ...]) -> None: + """Set selected 1-based degrees of freedom in a 6-entry fixity vector.""" + for dof in dofs_1based: + fixity[dof - 1] = 1 + + +def nodes_from_edges(edge_segments: list[tuple[int, int]]) -> list[int]: + """Return sorted unique node tags from a sequence of edge node pairs.""" + return sorted({node for edge in edge_segments for node in edge}) + + +def is_heterosis_plate(element_name: str) -> bool: + """Return true for the HeterosisPlate element.""" + return element_name.lower() == "heterosisplate" + + +def heterosis_rotation_only_nodes(model: xara.Model) -> set[int]: + """ + Return nodes that belong to the Q9 rotation field but not to the Q8 w field. + + HeterosisPlate uses the first 8 connectivity nodes for transverse + displacement interpolation and all 9 connectivity nodes for rotation + interpolation. Therefore, the 9th local node is rotation-only unless the + same global node appears in the first 8 positions of another element. + """ + w_nodes: set[int] = set() + rotation_center_nodes: set[int] = set() + + for element_tag in model.getEleTags(): + connectivity = [int(tag) for tag in model.eleNodes(element_tag)] + + if len(connectivity) == 9: + w_nodes.update(connectivity[:8]) + rotation_center_nodes.add(connectivity[8]) + + return rotation_center_nodes - w_nodes + + +def apply_boundary_conditions( + model: xara.Model, + surface: Any, + element_name: str, + n_el: int, +) -> None: + """Apply symmetry, clamp, and Heterosis inactive-degree constraints.""" + edge_node_pairs = list(surface.walk_edge()) + + x_symmetry_nodes = nodes_from_edges(edge_node_pairs[:n_el]) + clamped_outer_arc_nodes = nodes_from_edges(edge_node_pairs[n_el : 3 * n_el]) + y_symmetry_nodes = nodes_from_edges(edge_node_pairs[3 * n_el :]) + + node_fixities = { + node_tag: [0, 0, 0, 0, 0, 0] + for node_tag in model.getNodeTags() + } + + for node_tag in x_symmetry_nodes: + add_to_fixity(node_fixities[node_tag], (2, 4)) + + for node_tag in y_symmetry_nodes: + add_to_fixity(node_fixities[node_tag], (1, 5)) + + for node_tag in clamped_outer_arc_nodes: + add_to_fixity(node_fixities[node_tag], (3, 4, 5)) + + if is_heterosis_plate(element_name): + for node_tag in model.getNodeTags(): + add_to_fixity(node_fixities[node_tag], (1, 2, 6)) + + for node_tag in heterosis_rotation_only_nodes(model): + add_to_fixity(node_fixities[node_tag], (3,)) + + elif "plate" in element_name.lower(): + for node_tag in model.getNodeTags(): + add_to_fixity(node_fixities[node_tag], (6,)) + + for node_tag, fixity in node_fixities.items(): + if any(fixity): + model.fix(node_tag, tuple(fixity)) + + +def build_model(element_name: str, order: int, n_el: int) -> tuple[xara.Model, int]: + """Build the shell-0001 model and return the model and center node tag.""" + model = xara.Model(ndm=3, ndf=6) + + model.section("ElasticShell", SECTION_TAG, E_MOD, NU, THICKNESS) + + surface = model.surface( + (n_el, n_el), + element=element_name, + args={"section": SECTION_TAG}, + order=order, + points=quarter_disk_points(), + ) + + apply_boundary_conditions(model, surface, element_name, n_el) + + center = find_node(model, x=0.0, y=0.0) + + model.pattern("Plain", LOAD_PATTERN_TAG, "Linear") + model.load( + center, + (0.0, 0.0, P_QUARTER, 0.0, 0.0, 0.0), + pattern=LOAD_PATTERN_TAG, + ) + + return model, center + + +def nodal_profile( + model: xara.Model, + excluded_node_tags: set[int] | None = None, +) -> np.ndarray: + """Return sorted nodal samples (rho, radius, w_down).""" + excluded_node_tags = excluded_node_tags or set() + samples = [] + + for tag in model.getNodeTags(): + if tag in excluded_node_tags: + continue + + x, y, _ = model.nodeCoord(tag) + radius = math.hypot(x, y) + rho = radius / RADIUS + + if 0.05 < rho < 0.95: + w_down = -model.nodeDisp(tag)[2] + samples.append((rho, radius, w_down)) + + if not samples: + raise RuntimeError("No nodal profile points found in 0.05 < r/R < 0.95.") + + return np.asarray(sorted(samples, key=lambda row: row[0]), dtype=float) + + +def assert_pointwise_match( + model: xara.Model, + element_name: str, + mesh_order: int, + max_relative_error: float, + max_absolute_error: float, + excluded_node_tags: set[int] | None = None, +) -> None: + """Check the normalized displacement profile against the reference solution.""" + profile = nodal_profile(model, excluded_node_tags=excluded_node_tags) + radius = profile[:, 1] + + y_fem = normalized_displacement(profile[:, 2]) + y_ref = normalized_displacement(reissner_clamped_center_load_w(radius)) + + pointwise_abs_error = np.abs(y_fem - y_ref) + pointwise_rel_error = pointwise_abs_error / np.maximum( + np.abs(y_ref), + REFERENCE_VALUE_FLOOR, + ) + + passing_mask = ( + (pointwise_rel_error <= max_relative_error) + | (pointwise_abs_error <= max_absolute_error) + ) + failing_indices = np.where(~passing_mask)[0] + + assert failing_indices.size == 0, ( + f"{element_name} (order={mesh_order}, " + f"n_el={MESH_SUBDIVISIONS_PER_DIRECTION}) has " + f"{failing_indices.size}/{len(pointwise_abs_error)} nodes violating " + f"rel<={max_relative_error:.3e} OR abs<={max_absolute_error:.3e} " + f"(reference floor={REFERENCE_VALUE_FLOOR:.3e}); " + f"max rel={pointwise_rel_error.max():.6e}, " + f"max abs={pointwise_abs_error.max():.6e}" + ) + + +@pytest.mark.parametrize( + "element_name,mesh_order", + [ + (name, order) + for name, order in sorted(ELEMENT_ORDER.items()) + if name != "HeterosisPlate" + ], +) +def test_hughes_clamped_plate_pointwise_error( + element_name: str, + mesh_order: int, +) -> None: + """Verify shell elements node-by-node against the Reissner reference.""" + try: + model, _ = build_model( + element_name, + mesh_order, + MESH_SUBDIVISIONS_PER_DIRECTION, + ) + except XaraError as exc: + if "unknown element type" in str(exc).lower(): + pytest.skip(f"{element_name} is not available in this OpenSees build.") + raise + + analysis = xara.StaticAnalysis(model, system="Umfpack") + analysis.analyze() + + assert_pointwise_match( + model=model, + element_name=element_name, + mesh_order=mesh_order, + max_relative_error=MAX_POINTWISE_REL_ERROR, + max_absolute_error=MAX_POINTWISE_ABS_ERROR, + ) + + +def test_hughes_clamped_plate_pointwise_error_heterosis() -> None: + """Verify HeterosisPlate using connectivity-aware displacement sampling.""" + element_name = "HeterosisPlate" + mesh_order = ELEMENT_ORDER[element_name] + + try: + model, _ = build_model( + element_name, + mesh_order, + MESH_SUBDIVISIONS_PER_DIRECTION, + ) + except XaraError as exc: + if "unknown element type" in str(exc).lower(): + pytest.skip(f"{element_name} is not available in this OpenSees build.") + raise + + analysis = xara.StaticAnalysis(model, system="Umfpack") + analysis.analyze() + + excluded_nodes = heterosis_rotation_only_nodes(model) + + assert_pointwise_match( + model=model, + element_name=element_name, + mesh_order=mesh_order, + max_relative_error=MAX_POINTWISE_REL_ERROR, + max_absolute_error=MAX_POINTWISE_ABS_ERROR, + excluded_node_tags=excluded_nodes, + ) \ No newline at end of file diff --git a/benchmarks/shell-0002/README.md b/benchmarks/shell-0002/README.md new file mode 100644 index 0000000..dc89d00 --- /dev/null +++ b/benchmarks/shell-0002/README.md @@ -0,0 +1,198 @@ +# shell-0002: Distorted Five-Element Patch Test + +This benchmark verifies generalized strain recovery on a distorted five-element shell/plate patch. The test is a kinematic patch test: a prescribed displacement and rotation field is evaluated through the element interpolation, and the recovered generalized strains are compared against the corresponding analytical strain field. + +The benchmark is intended to detect mistakes in element connectivity, interpolation order, geometry mapping, local node ordering, and the mixed-interpolation treatment used by `HeterosisPlate`. + +## Patch Geometry + +The benchmark uses an enclosing five-element patch. The outer boundary is rectangular, and the center region is a distorted quadrilateral. The five patch regions are: + +- `E_center` +- `E_bottom` +- `E_right` +- `E_top` +- `E_left` + +The same patch geometry is shared by the pytest benchmark and the plotting utility through `_shell_0002_patch.py`. + +![Distorted five-element patch](./_outputs/mesh_plots/HeterosisPlate_distorted_patch.png) + +## Files + +- `_shell_0002_patch.py` + Defines the shared patch geometry and Xara model-building utilities. + +- `test_shell_0002.py` + Defines the pytest benchmark for exact generalized strain recovery. + +- `plot_shell_0002_mesh.py` + Provides a standalone plotting utility for generating VEUx and Matplotlib visualizations of the patch. + +## Elements Tested + +The benchmark is parametrized over the following elements: + +```python +ELEMENT_ORDER = { + "ASDShellQ4": 1, + "ShellMITC4": 1, + "HeterosisPlate": 2, + "ShellMITC9": 2, +} +``` + +The interpolation used for strain recovery depends on the element: + +| Element | Transverse displacement interpolation | Rotation interpolation | +|---|---:|---:| +| `ASDShellQ4` | Q4 | Q4 | +| `ShellMITC4` | Q4 | Q4 | +| `ShellMITC9` | Q9 | Q9 | +| `HeterosisPlate` | Q8 | Q9 | + +The `HeterosisPlate` case is handled inside the recovery logic because its external OpenSees connectivity uses 9 shell nodes, while the physical plate interpolation uses Q8 transverse displacement and Q9 rotations. + +## Prescribed Kinematic Field + +The benchmark prescribes a linear transverse displacement field and linear rotation fields: + +```text +w(x, y) = w0 + wx x + wy y +theta_x(x, y) = theta_x0 + theta_xx x + theta_xy y +theta_y(x, y) = theta_y0 + theta_yx x + theta_yy y +``` + +The current numerical values are: + +```python +w0 = 0.120 +wx = 0.170 +wy = -0.090 + +theta_x0 = 0.030 +theta_xx = 0.041 +theta_xy = -0.022 + +theta_y0 = -0.025 +theta_yx = 0.033 +theta_yy = 0.052 +``` + +This field is exactly representable by Q4, Q8, and Q9 interpolation. Therefore, if the connectivity, geometry mapping, and strain recovery logic are correct, the numerical strain field should match the analytical strain field to machine precision. + +## Strain Components Compared + +At each sampling point, the benchmark recovers and compares the following generalized strain components: + +```text +kappa_xx = theta_x,x +kappa_yy = theta_y,y +kappa_xy = theta_x,y + theta_y,x +gamma_xz = w,x - theta_x +gamma_yz = w,y - theta_y +``` + +The analytical values are obtained directly from the prescribed field: + +```text +kappa_xx = theta_xx +kappa_yy = theta_yy +kappa_xy = theta_xy + theta_yx +gamma_xz = wx - theta_x(x, y) +gamma_yz = wy - theta_y(x, y) +``` + +The bending strain components are constant because the prescribed rotation fields are linear. The transverse shear strain components vary linearly because they include the rotation fields themselves. + +## What the Test Verifies + +The benchmark verifies that: + +- The prescribed kinematic field is reproduced exactly by the selected interpolation basis. +- Shape-function derivatives are correctly transformed from parent coordinates to physical coordinates. +- The physical point associated with the transverse displacement interpolation matches the physical point associated with the rotation interpolation. +- The `HeterosisPlate` Q8/Q9 split is handled correctly. +- Distorted-element geometry does not break exact recovery of a representable field. + +This is not a load-response benchmark. No equilibrium solution is being checked. The test is focused on interpolation and kinematic strain recovery. + +## Sampling Points + +The sampling points depend on the element order: + +- Order 1 elements use 2 by 2 Gauss sampling. +- Order 2 elements use 3 by 3 Gauss sampling. + +For each element, strains are sampled at all corresponding points and compared against the analytical values. + +## Running the Test + +From the `shell-0002` folder: + +```bash +python -m pytest test_shell_0002.py -v -rs +``` + +Expected collection: + +```text +test_distorted_patch_exact_strain_recovery[ASDShellQ4-1] +test_distorted_patch_exact_strain_recovery[HeterosisPlate-2] +test_distorted_patch_exact_strain_recovery[ShellMITC4-1] +test_distorted_patch_exact_strain_recovery[ShellMITC9-2] +test_distorted_patch_output_files_saved[ASDShellQ4-1] +test_distorted_patch_output_files_saved[HeterosisPlate-2] +test_distorted_patch_output_files_saved[ShellMITC4-1] +test_distorted_patch_output_files_saved[ShellMITC9-2] +``` + +The benchmark runs eight cases by default: two tests for each of the four element types. + +## Diagnostic Outputs + +The pytest diagnostics are written under: + +```text +_outputs/patch_diagnostics/ +``` + +The plotting utility writes patch figures under: + +```text +_outputs/mesh_plots/ +``` + +To regenerate the patch visualization: + +```bash +python plot_shell_0002_mesh.py +``` + +By default, the plotting script creates both: + +```text +_outputs/mesh_plots/HeterosisPlate_distorted_patch.html +_outputs/mesh_plots/HeterosisPlate_distorted_patch.png +``` + +The Matplotlib figure is intended for documentation and topology inspection. The VEUx figure is intended for interactive visualization. + +## Notes on Interpretation + +A passing result indicates that the selected interpolation basis can exactly recover the generalized strains of the prescribed representable field on the distorted patch. + +A failure may indicate one of the following issues: + +- Incorrect local node ordering. +- Incorrect Q8/Q9 treatment for `HeterosisPlate`. +- Incorrect parent-to-physical derivative transformation. +- Incompatible displacement and rotation geometry mappings. +- An unintended change in the shared patch geometry. +- An unavailable or changed element implementation in the active OpenSees/Xara build. + +## References + +[1] O. C. Zienkiewicz, R. L. Taylor, and S. Govindjee, *The Finite Element Method: Its Basis and Fundamentals*. Elsevier, 8th ed., Nov. 2024. + +[2] T. J. R. Hughes, *The Finite Element Method: Linear Static and Dynamic Finite Element Analysis*. Englewood Cliffs, New Jersey: Prentice-Hall, 1987. \ No newline at end of file diff --git a/benchmarks/shell-0002/_shell_0002_patch.py b/benchmarks/shell-0002/_shell_0002_patch.py new file mode 100644 index 0000000..e0c995f --- /dev/null +++ b/benchmarks/shell-0002/_shell_0002_patch.py @@ -0,0 +1,118 @@ +# _shell_0002_patch.py +""" +Shared distorted five-element patch definition for benchmark shell-0002. + +This file is imported by both the pytest benchmark and the standalone plotting +script so that the geometry is defined in exactly one place. +""" + +from __future__ import annotations + +import numpy as np +import xara + + +SECTION_TAG = 1 + +E_MOD = 200_000.0 +NU = 0.25 +THICKNESS = 0.2 + +ELEMENT_ORDER = { + "ASDShellQ4": 1, + "ShellMITC4": 1, + "HeterosisPlate": 2, + "ShellMITC9": 2, +} + + +# Enclosing-topology patch: outer boundary plus distorted inner quadrilateral. +OUTER_BL = (0.00, 0.00) +OUTER_BR = (3.40, 0.00) +OUTER_TR = (3.40, 3.00) +OUTER_TL = (0.00, 3.00) + +INNER_BL = (1.12, 1.05) +INNER_BR = (2.24, 1.00) +INNER_TR = (2.15, 2.08) +INNER_TL = (1.02, 2.14) + +PATCH_BLOCKS = { + "center": [INNER_BL, INNER_BR, INNER_TR, INNER_TL], + "bottom": [OUTER_BL, OUTER_BR, INNER_BR, INNER_BL], + "right": [INNER_BR, OUTER_BR, OUTER_TR, INNER_TR], + "top": [INNER_TL, INNER_TR, OUTER_TR, OUTER_TL], + "left": [OUTER_BL, INNER_BL, INNER_TL, OUTER_TL], +} + + +def q4_surface_points(corners: list[tuple[float, float]]) -> dict[int, list[float]]: + """Return Q4 surface control points from four counter-clockwise corners.""" + return { + index + 1: [float(x), float(y), 0.0] + for index, (x, y) in enumerate(corners) + } + + +def q9_surface_points(corners: list[tuple[float, float]]) -> dict[int, list[float]]: + """Return Q9 surface control points from four counter-clockwise corners.""" + p1, p2, p3, p4 = [np.asarray(point, dtype=float) for point in corners] + + points = [ + p1, + p2, + p3, + p4, + 0.5 * (p1 + p2), + 0.5 * (p2 + p3), + 0.5 * (p3 + p4), + 0.5 * (p4 + p1), + 0.25 * (p1 + p2 + p3 + p4), + ] + + return { + index + 1: [float(point[0]), float(point[1]), 0.0] + for index, point in enumerate(points) + } + + +def surface_points( + corners: list[tuple[float, float]], + mesh_order: int, +) -> dict[int, list[float]]: + """Return surface control points for the requested mesh order.""" + if mesh_order == 1: + return q4_surface_points(corners) + + if mesh_order == 2: + return q9_surface_points(corners) + + raise ValueError(f"Unsupported mesh order: {mesh_order}") + + +def build_model(element_name: str, mesh_order: int) -> xara.Model: + """Build the distorted five-element patch.""" + model = xara.Model(ndm=3, ndf=6) + model.section("ElasticShell", SECTION_TAG, E_MOD, NU, THICKNESS) + + for corners in PATCH_BLOCKS.values(): + model.surface( + (1, 1), + element=element_name, + args={"section": SECTION_TAG}, + order=mesh_order, + points=surface_points(corners, mesh_order), + ) + + return model + + +def element_coordinates(model: xara.Model, element_tag: int) -> tuple[list[int], np.ndarray]: + """Return element connectivity and xy coordinates.""" + connectivity = [int(tag) for tag in model.eleNodes(element_tag)] + coordinates = np.array( + [model.nodeCoord(node_tag)[:2] for node_tag in connectivity], + dtype=float, + ) + + return connectivity, coordinates \ No newline at end of file diff --git a/benchmarks/shell-0002/plot_shell_0002_mesh.py b/benchmarks/shell-0002/plot_shell_0002_mesh.py new file mode 100644 index 0000000..7805c87 --- /dev/null +++ b/benchmarks/shell-0002/plot_shell_0002_mesh.py @@ -0,0 +1,451 @@ +""" +Standalone patch-plot utility for benchmark shell-0002. + +This script uses the shared distorted-patch geometry from +``_shell_0002_patch.py``. Therefore, the plotting utility and the pytest +benchmark always use the same patch definition. + +Default behavior +---------------- +Both plotting engines are used by default. Running the script without flags +generates: + + 1. an interactive VEUx/Plotly HTML file + 2. a static Matplotlib PNG file + +For example: + + python plot_shell_0002_mesh.py + +renders the default HeterosisPlate patch with both engines. + +Important plotting convention +----------------------------- +This script plots the patch topology, not the element interpolation topology. +Therefore, it does not distinguish Q4, Q8, Q9, midside nodes, or center nodes. +The same five-region patch is used by all element types in the benchmark. + +To avoid misleading visual overlap: + + - shared patch edges are drawn only once in Matplotlib + - VEUx receives one quadrilateral panel per patch region + - no artificial triangular subdivision is used + - no repeated finite-element nodes are plotted + +Examples +-------- +Render the default HeterosisPlate patch with both engines: + + python plot_shell_0002_mesh.py + +Render a specific element label with both engines: + + python plot_shell_0002_mesh.py --element ShellMITC9 + +Render all element labels with both engines: + + python plot_shell_0002_mesh.py --all + +Render only the VEUx HTML file: + + python plot_shell_0002_mesh.py --renderer veux + +Render only the Matplotlib PNG file: + + python plot_shell_0002_mesh.py --renderer matplotlib +""" + +from __future__ import annotations + +import argparse +from pathlib import Path +from typing import Any + +import matplotlib + +matplotlib.use("Agg") + +import matplotlib.pyplot as plt +import numpy as np +from opensees.errors import XaraError + +from _shell_0002_patch import ( + ELEMENT_ORDER, + PATCH_BLOCKS, + build_model, +) + + +OUTPUT_DIR = Path(__file__).resolve().parent / "_outputs" / "mesh_plots" + + +def default_veux_path(element_name: str, canvas: str) -> Path: + """Return the default VEUx output path.""" + suffix_by_canvas = { + "plotly": ".html", + "gltf": ".glb", + } + + return OUTPUT_DIR / f"{element_name}_distorted_patch{suffix_by_canvas[canvas]}" + + +def default_matplotlib_path(element_name: str) -> Path: + """Return the default Matplotlib output path.""" + return OUTPUT_DIR / f"{element_name}_distorted_patch.png" + + +def coordinate_key(point: tuple[float, float], decimals: int = 12) -> tuple[float, float]: + """Return a stable key for identifying repeated patch coordinates.""" + return (round(float(point[0]), decimals), round(float(point[1]), decimals)) + + +def unique_patch_nodes() -> dict[tuple[float, float], int]: + """Return unique patch coordinates mapped to visualization node tags.""" + coordinate_to_tag: dict[tuple[float, float], int] = {} + + for corners in PATCH_BLOCKS.values(): + for point in corners: + key = coordinate_key(point) + + if key not in coordinate_to_tag: + coordinate_to_tag[key] = len(coordinate_to_tag) + 1 + + return coordinate_to_tag + + +def unique_patch_edges() -> list[tuple[tuple[float, float], tuple[float, float]]]: + """ + Return unique patch edges. + + Adjacent regions share edges. The edge key is orientation-independent so + each shared edge is drawn only once. + """ + edge_by_key: dict[ + tuple[tuple[float, float], tuple[float, float]], + tuple[tuple[float, float], tuple[float, float]], + ] = {} + + for corners in PATCH_BLOCKS.values(): + ordered_points = [coordinate_key(point) for point in corners] + + for start, end in zip( + ordered_points, + ordered_points[1:] + ordered_points[:1], + strict=True, + ): + edge_key = tuple(sorted((start, end))) + + if edge_key not in edge_by_key: + edge_by_key[edge_key] = (start, end) + + return list(edge_by_key.values()) + + +def patch_region_centroid(corners: list[tuple[float, float]]) -> np.ndarray: + """Return the centroid of a patch region.""" + coordinates = np.asarray(corners, dtype=float) + return coordinates.mean(axis=0) + + +def patch_label(region_name: str) -> str: + """Return a formatted label for a patch region.""" + return rf"$E_{{\mathrm{{{region_name}}}}}$" + + +def verify_element_available(element_name: str) -> int: + """ + Build the requested Xara model to verify that the element is available. + + The returned model is not used for plotting. The plot is based on + PATCH_BLOCKS so the displayed geometry stays identical for all elements. + """ + mesh_order = ELEMENT_ORDER[element_name] + + try: + build_model(element_name, mesh_order) + except XaraError as exc: + if "unknown element type" in str(exc).lower(): + raise RuntimeError( + f"{element_name} is not available in this OpenSees build." + ) from exc + raise + + return mesh_order + + +def veux_patch_model() -> dict[str, Any]: + """ + Convert PATCH_BLOCKS into a VEUx-friendly quadrilateral surface model. + + This model is for visualization only. It contains one four-node + quadrilateral panel per patch region. + """ + coordinate_to_tag = unique_patch_nodes() + + nodes = [] + for coordinate, tag in sorted( + coordinate_to_tag.items(), + key=lambda item: item[1], + ): + x, y = coordinate + nodes.append( + { + "name": tag, + "crd": [float(x), float(y), 0.0], + } + ) + + elements = [] + for element_tag, corners in enumerate(PATCH_BLOCKS.values(), start=1): + nodes_for_region = [ + coordinate_to_tag[coordinate_key(point)] + for point in corners + ] + + elements.append( + { + "name": element_tag, + "type": "ShellMITC4", + "nodes": nodes_for_region, + } + ) + + return { + "StructuralAnalysisModel": { + "properties": { + "sections": [], + "nDMaterials": [], + "uniaxialMaterials": [], + "crdTransformations": [], + "patterns": [], + "parameters": [], + }, + "geometry": { + "nodes": nodes, + "elements": elements, + "constraints": [], + }, + } + } + + +def save_veux_artist(artist: Any, output_path: Path) -> None: + """Save a VEUx artist.""" + if not hasattr(artist, "save"): + raise RuntimeError("The VEUx artist does not expose a save() method.") + + artist.save(str(output_path)) + + +def serve_veux_artist(artist: Any) -> None: + """Serve a VEUx artist locally.""" + import veux + + veux.serve(artist) + + +def render_with_veux( + output_path: Path, + canvas: str, + serve: bool, +) -> Path: + """Render the five-region patch with VEUx.""" + import veux + + output_path.parent.mkdir(parents=True, exist_ok=True) + + artist = veux.render( + veux_patch_model(), + canvas=canvas, + vertical=3, + ) + + save_veux_artist(artist, output_path) + + if serve: + serve_veux_artist(artist) + + return output_path + + +def render_with_matplotlib( + element_name: str, + output_path: Path, + show_region_labels: bool, +) -> Path: + """ + Render the five-region patch with Matplotlib. + + The plot uses unique patch edges, so shared internal boundaries are not + drawn repeatedly. + """ + output_path.parent.mkdir(parents=True, exist_ok=True) + + fig, ax = plt.subplots(figsize=(5.2, 4.6)) + + for start, end in unique_patch_edges(): + x_values = [start[0], end[0]] + y_values = [start[1], end[1]] + + ax.plot( + x_values, + y_values, + color="black", + linewidth=1.1, + ) + + if show_region_labels: + for region_name, corners in PATCH_BLOCKS.items(): + centroid = patch_region_centroid(corners) + ax.text( + centroid[0], + centroid[1], + patch_label(region_name), + ha="center", + va="center", + fontsize=8, + ) + + ax.set_aspect("equal", adjustable="box") + ax.set_xlabel("x (unitless)", fontsize=9) + ax.set_ylabel("y (unitless)", fontsize=9) + ax.set_title( + f"{element_name}: distorted 5-element patch", + fontsize=10, + ) + + ax.grid(True, linewidth=0.35, alpha=0.25) + ax.tick_params(labelsize=8) + + fig.tight_layout() + fig.savefig(output_path, dpi=300) + plt.close(fig) + + return output_path + + +def plot_patch( + element_name: str, + renderer: str, + canvas: str, + output_path: Path | None, + serve: bool, + show_region_labels: bool, +) -> list[Path]: + """ + Plot the shared distorted patch. + + When both renderers are requested, the optional output path is ignored so + that each renderer can use its own default file extension. + """ + mesh_order = verify_element_available(element_name) + + print( + f"Rendering shared patch for element label: " + f"{element_name} (mesh order {mesh_order})" + ) + print("Displayed patch geometry is identical for all element types.") + + written_paths: list[Path] = [] + + if renderer in {"both", "veux"}: + veux_path = output_path if renderer == "veux" else None + path = render_with_veux( + output_path=veux_path or default_veux_path(element_name, canvas), + canvas=canvas, + serve=serve, + ) + written_paths.append(path) + print(f"Wrote VEUx patch plot: {path}") + + if renderer in {"both", "matplotlib"}: + matplotlib_path = output_path if renderer == "matplotlib" else None + path = render_with_matplotlib( + element_name=element_name, + output_path=matplotlib_path or default_matplotlib_path(element_name), + show_region_labels=show_region_labels, + ) + written_paths.append(path) + print(f"Wrote Matplotlib patch plot: {path}") + + return written_paths + + +def parse_arguments() -> argparse.Namespace: + """Parse command-line arguments.""" + parser = argparse.ArgumentParser( + description="Generate patch plots for the shell-0002 distorted patch.", + ) + + element_group = parser.add_mutually_exclusive_group() + element_group.add_argument( + "--element", + choices=sorted(ELEMENT_ORDER), + default="HeterosisPlate", + help="Element label to render. Defaults to HeterosisPlate.", + ) + element_group.add_argument( + "--all", + action="store_true", + help="Render plots for all element labels.", + ) + + parser.add_argument( + "--renderer", + choices=("both", "veux", "matplotlib"), + default="both", + help=( + "Renderer to use. Defaults to both, which writes one VEUx HTML " + "file and one Matplotlib PNG file." + ), + ) + parser.add_argument( + "--canvas", + choices=("plotly", "gltf"), + default="plotly", + help="VEUx canvas. Defaults to plotly.", + ) + parser.add_argument( + "--out", + type=Path, + default=None, + help=( + "Output path for single-renderer mode. Ignored when " + "--renderer both or --all is used." + ), + ) + parser.add_argument( + "--serve", + action="store_true", + help="Serve the VEUx rendering locally after writing the output file.", + ) + parser.add_argument( + "--no-region-labels", + action="store_true", + help="Hide patch-region labels in the Matplotlib plot.", + ) + + return parser.parse_args() + + +def main() -> None: + """Run the patch-plot utility.""" + args = parse_arguments() + + element_names = sorted(ELEMENT_ORDER) if args.all else [args.element] + + for element_name in element_names: + output_path = None if args.all else args.out + + plot_patch( + element_name=element_name, + renderer=args.renderer, + canvas=args.canvas, + output_path=output_path, + serve=args.serve, + show_region_labels=not args.no_region_labels, + ) + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/benchmarks/shell-0002/test_shell_0002.py b/benchmarks/shell-0002/test_shell_0002.py new file mode 100644 index 0000000..5373091 --- /dev/null +++ b/benchmarks/shell-0002/test_shell_0002.py @@ -0,0 +1,637 @@ +""" +Self-contained distorted 5-element patch test for shell and plate elements. + +The benchmark builds an enclosing five-element patch with an outer rectangular +boundary and a distorted inner quadrilateral. The patch topology is: + + E_top +E_left E_center E_right + E_bottom + +A representable linear kinematic field is prescribed: + + w(x, y) = w0 + wx x + wy y + theta_x(x, y) = tx0 + txx x + txy y + theta_y(x, y) = ty0 + tyx x + tyy y + +The recovered generalized strains are compared against the analytical values. +The interpolation used for strain recovery depends on the element: + + ASDShellQ4, ShellMITC4 : Q4 w and Q4 rotations + ShellMITC9 : Q9 w and Q9 rotations + HeterosisPlate : Q8 w and Q9 rotations + +The HeterosisPlate case is handled internally because its external OpenSees +connectivity uses 9 shell nodes, while the physical plate interpolation uses +Q8 transverse displacement and Q9 rotations. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path + +import matplotlib + +matplotlib.use("Agg") + +import matplotlib.pyplot as plt +import numpy as np +import pytest +import xara +from opensees.errors import XaraError + +from _shell_0002_patch import ( + ELEMENT_ORDER, + build_model, + element_coordinates, +) + + +ABS_TOL = 1.0e-10 + +OUTPUT_DIR = Path(__file__).resolve().parent / "_outputs" + +STRAIN_NAMES = ( + "kappa_xx", + "kappa_yy", + "kappa_xy", + "gamma_xz", + "gamma_yz", +) + +GAUSS_2 = ( + -1.0 / np.sqrt(3.0), + 1.0 / np.sqrt(3.0), +) + +GAUSS_3 = ( + -np.sqrt(3.0 / 5.0), + 0.0, + np.sqrt(3.0 / 5.0), +) + +GAUSS_POINTS_BY_ORDER = { + 1: tuple((xi, eta) for xi in GAUSS_2 for eta in GAUSS_2), + 2: tuple((xi, eta) for xi in GAUSS_3 for eta in GAUSS_3), +} + + +@dataclass(frozen=True) +class LinearPatchField: + """Linear field exactly representable by all interpolations used here.""" + + w0: float = 0.120 + wx: float = 0.170 + wy: float = -0.090 + + tx0: float = 0.030 + txx: float = 0.041 + txy: float = -0.022 + + ty0: float = -0.025 + tyx: float = 0.033 + tyy: float = 0.052 + + def w(self, x: float, y: float) -> float: + return self.w0 + self.wx * x + self.wy * y + + def theta_x(self, x: float, y: float) -> float: + return self.tx0 + self.txx * x + self.txy * y + + def theta_y(self, x: float, y: float) -> float: + return self.ty0 + self.tyx * x + self.tyy * y + + def expected_at(self, x: float, y: float) -> dict[str, float]: + return { + "kappa_xx": self.txx, + "kappa_yy": self.tyy, + "kappa_xy": self.txy + self.tyx, + "gamma_xz": self.wx - self.theta_x(x, y), + "gamma_yz": self.wy - self.theta_y(x, y), + } + + +def is_heterosis_plate(element_name: str) -> bool: + """Return true for the HeterosisPlate element.""" + return element_name.lower() == "heterosisplate" + + +def build_model_or_skip(element_name: str, mesh_order: int) -> xara.Model: + """Build the patch model, skipping unavailable elements.""" + try: + model = build_model(element_name, mesh_order) + except XaraError as exc: + if "unknown element type" in str(exc).lower(): + pytest.skip(f"{element_name} is not available in this OpenSees build.") + raise + + assert len(model.getEleTags()) == 5 + return model + + +def q4_shape_values_and_derivatives(xi: float, eta: float) -> tuple[np.ndarray, np.ndarray]: + """Return Q4 shape functions and derivatives in local-node order.""" + values = np.array( + [ + 0.25 * (1.0 - xi) * (1.0 - eta), + 0.25 * (1.0 + xi) * (1.0 - eta), + 0.25 * (1.0 + xi) * (1.0 + eta), + 0.25 * (1.0 - xi) * (1.0 + eta), + ], + dtype=float, + ) + + dxi = np.array( + [ + -0.25 * (1.0 - eta), + 0.25 * (1.0 - eta), + 0.25 * (1.0 + eta), + -0.25 * (1.0 + eta), + ], + dtype=float, + ) + + deta = np.array( + [ + -0.25 * (1.0 - xi), + -0.25 * (1.0 + xi), + 0.25 * (1.0 + xi), + 0.25 * (1.0 - xi), + ], + dtype=float, + ) + + return values, np.column_stack((dxi, deta)) + + +def q8_shape_values_and_derivatives(xi: float, eta: float) -> tuple[np.ndarray, np.ndarray]: + """Return Q8 serendipity shape functions and derivatives in local-node order.""" + values = np.array( + [ + 0.25 * (1.0 - xi) * (1.0 - eta) * (-xi - eta - 1.0), + 0.25 * (1.0 + xi) * (1.0 - eta) * (xi - eta - 1.0), + 0.25 * (1.0 + xi) * (1.0 + eta) * (xi + eta - 1.0), + 0.25 * (1.0 - xi) * (1.0 + eta) * (-xi + eta - 1.0), + 0.5 * (1.0 - xi**2) * (1.0 - eta), + 0.5 * (1.0 + xi) * (1.0 - eta**2), + 0.5 * (1.0 - xi**2) * (1.0 + eta), + 0.5 * (1.0 - xi) * (1.0 - eta**2), + ], + dtype=float, + ) + + dxi = np.array( + [ + -0.5 * xi * eta + 0.5 * xi - 0.25 * eta**2 + 0.25 * eta, + -0.5 * xi * eta + 0.5 * xi + 0.25 * eta**2 - 0.25 * eta, + 0.5 * xi * eta + 0.5 * xi + 0.25 * eta**2 + 0.25 * eta, + 0.5 * xi * eta + 0.5 * xi - 0.25 * eta**2 - 0.25 * eta, + xi * eta - xi, + 0.5 - 0.5 * eta**2, + -xi * eta - xi, + 0.5 * eta**2 - 0.5, + ], + dtype=float, + ) + + deta = np.array( + [ + -0.25 * xi**2 - 0.5 * xi * eta + 0.25 * xi + 0.5 * eta, + -0.25 * xi**2 + 0.5 * xi * eta - 0.25 * xi + 0.5 * eta, + 0.25 * xi**2 + 0.5 * xi * eta + 0.25 * xi + 0.5 * eta, + 0.25 * xi**2 - 0.5 * xi * eta - 0.25 * xi + 0.5 * eta, + 0.5 * xi**2 - 0.5, + -xi * eta - eta, + 0.5 - 0.5 * xi**2, + xi * eta - eta, + ], + dtype=float, + ) + + return values, np.column_stack((dxi, deta)) + + +def q9_shape_values_and_derivatives(xi: float, eta: float) -> tuple[np.ndarray, np.ndarray]: + """Return Q9 shape functions and derivatives in local-node order.""" + lx = np.array( + [ + 0.5 * xi * (xi - 1.0), + 1.0 - xi**2, + 0.5 * xi * (xi + 1.0), + ], + dtype=float, + ) + ly = np.array( + [ + 0.5 * eta * (eta - 1.0), + 1.0 - eta**2, + 0.5 * eta * (eta + 1.0), + ], + dtype=float, + ) + dlx = np.array([xi - 0.5, -2.0 * xi, xi + 0.5], dtype=float) + dly = np.array([eta - 0.5, -2.0 * eta, eta + 0.5], dtype=float) + + pairs = ( + (0, 0), + (2, 0), + (2, 2), + (0, 2), + (1, 0), + (2, 1), + (1, 2), + (0, 1), + (1, 1), + ) + + values = np.array([lx[i] * ly[j] for i, j in pairs], dtype=float) + dxi = np.array([dlx[i] * ly[j] for i, j in pairs], dtype=float) + deta = np.array([lx[i] * dly[j] for i, j in pairs], dtype=float) + + return values, np.column_stack((dxi, deta)) + + +def shape_values_and_derivatives( + basis_name: str, + xi: float, + eta: float, +) -> tuple[np.ndarray, np.ndarray]: + """Return shape functions and parent-coordinate derivatives.""" + if basis_name == "q4": + return q4_shape_values_and_derivatives(xi, eta) + + if basis_name == "q8": + return q8_shape_values_and_derivatives(xi, eta) + + if basis_name == "q9": + return q9_shape_values_and_derivatives(xi, eta) + + raise ValueError(f"Unsupported interpolation basis: {basis_name}") + + +def basis_node_count(basis_name: str) -> int: + """Return the number of local nodes used by an interpolation basis.""" + return { + "q4": 4, + "q8": 8, + "q9": 9, + }[basis_name] + + +def interpolation_bases(element_name: str, mesh_order: int) -> tuple[str, str]: + """Return the transverse-displacement and rotation interpolation bases.""" + if is_heterosis_plate(element_name): + return "q8", "q9" + + if mesh_order == 1: + return "q4", "q4" + + if mesh_order == 2: + return "q9", "q9" + + raise ValueError(f"Unsupported mesh order: {mesh_order}") + + +def sample_points(mesh_order: int) -> tuple[tuple[float, float], ...]: + """Return sampling points for the requested element order.""" + return GAUSS_POINTS_BY_ORDER[mesh_order] + + +def shape_gradients_xy( + natural_derivatives: np.ndarray, + coordinates: np.ndarray, +) -> np.ndarray: + """Transform shape-function derivatives from parent to physical coordinates.""" + x = coordinates[:, 0] + y = coordinates[:, 1] + + jacobian = np.array( + [ + [ + natural_derivatives[:, 0] @ x, + natural_derivatives[:, 1] @ x, + ], + [ + natural_derivatives[:, 0] @ y, + natural_derivatives[:, 1] @ y, + ], + ], + dtype=float, + ) + + return natural_derivatives @ np.linalg.inv(jacobian) + + +def field_values( + coordinates: np.ndarray, + field: LinearPatchField, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Return nodal values of w, theta_x, and theta_y.""" + w_values = np.array([field.w(x, y) for x, y in coordinates], dtype=float) + theta_x_values = np.array([field.theta_x(x, y) for x, y in coordinates], dtype=float) + theta_y_values = np.array([field.theta_y(x, y) for x, y in coordinates], dtype=float) + + return w_values, theta_x_values, theta_y_values + + +def recover_strains_at_point( + xi: float, + eta: float, + coordinates_w: np.ndarray, + coordinates_rotation: np.ndarray, + w_values: np.ndarray, + theta_x_values: np.ndarray, + theta_y_values: np.ndarray, + w_basis: str, + rotation_basis: str, +) -> tuple[dict[str, float], tuple[float, float]]: + """Recover generalized strains at one parent-coordinate point.""" + n_w, dn_w_parent = shape_values_and_derivatives(w_basis, xi, eta) + n_rotation, dn_rotation_parent = shape_values_and_derivatives( + rotation_basis, + xi, + eta, + ) + + dn_w_xy = shape_gradients_xy(dn_w_parent, coordinates_w) + dn_rotation_xy = shape_gradients_xy(dn_rotation_parent, coordinates_rotation) + + xy_w = n_w @ coordinates_w + xy_rotation = n_rotation @ coordinates_rotation + + np.testing.assert_allclose( + xy_w, + xy_rotation, + atol=1.0e-12, + err_msg="Displacement and rotation geometries are incompatible.", + ) + + w_x = dn_w_xy[:, 0] @ w_values + w_y = dn_w_xy[:, 1] @ w_values + + theta_x = n_rotation @ theta_x_values + theta_y = n_rotation @ theta_y_values + + theta_x_x = dn_rotation_xy[:, 0] @ theta_x_values + theta_x_y = dn_rotation_xy[:, 1] @ theta_x_values + theta_y_x = dn_rotation_xy[:, 0] @ theta_y_values + theta_y_y = dn_rotation_xy[:, 1] @ theta_y_values + + strains = { + "kappa_xx": theta_x_x, + "kappa_yy": theta_y_y, + "kappa_xy": theta_x_y + theta_y_x, + "gamma_xz": w_x - theta_x, + "gamma_yz": w_y - theta_y, + } + + return strains, (float(xy_rotation[0]), float(xy_rotation[1])) + + +def sample_generalized_strains( + model: xara.Model, + element_name: str, + mesh_order: int, + field: LinearPatchField, +) -> tuple[dict[str, np.ndarray], dict[str, np.ndarray]]: + """Recover generalized strains from the element interpolation.""" + w_basis, rotation_basis = interpolation_bases(element_name, mesh_order) + + sampled = {name: [] for name in STRAIN_NAMES} + expected = {name: [] for name in STRAIN_NAMES} + + for element_tag in model.getEleTags(): + _, coordinates = element_coordinates(model, element_tag) + + w_node_count = basis_node_count(w_basis) + rotation_node_count = basis_node_count(rotation_basis) + + coordinates_w = coordinates[:w_node_count] + coordinates_rotation = coordinates[:rotation_node_count] + + w_values, _, _ = field_values(coordinates_w, field) + _, theta_x_values, theta_y_values = field_values(coordinates_rotation, field) + + for xi, eta in sample_points(mesh_order): + sampled_values, physical_point = recover_strains_at_point( + xi=xi, + eta=eta, + coordinates_w=coordinates_w, + coordinates_rotation=coordinates_rotation, + w_values=w_values, + theta_x_values=theta_x_values, + theta_y_values=theta_y_values, + w_basis=w_basis, + rotation_basis=rotation_basis, + ) + + expected_values = field.expected_at(*physical_point) + + for name in STRAIN_NAMES: + sampled[name].append(sampled_values[name]) + expected[name].append(expected_values[name]) + + return ( + {name: np.asarray(values, dtype=float) for name, values in sampled.items()}, + {name: np.asarray(values, dtype=float) for name, values in expected.items()}, + ) + + +def assert_exact_strain_recovery(element_name: str, mesh_order: int) -> None: + """Verify exact generalized-strain recovery for one element type.""" + model = build_model_or_skip(element_name, mesh_order) + field = LinearPatchField() + + sampled, expected = sample_generalized_strains( + model=model, + element_name=element_name, + mesh_order=mesh_order, + field=field, + ) + + for name in STRAIN_NAMES: + np.testing.assert_allclose( + sampled[name], + expected[name], + atol=ABS_TOL, + err_msg=( + f"{element_name} (order={mesh_order}) failed the distorted " + f"patch strain-recovery test for {name}." + ), + ) + + +def plot_patch_geometry( + model: xara.Model, + element_name: str, + output_directory: Path, +) -> Path: + """Write a geometry plot for manual inspection.""" + output_directory.mkdir(parents=True, exist_ok=True) + path = output_directory / f"{element_name}_distorted_patch_geometry.png" + + fig, ax = plt.subplots(figsize=(7.0, 6.0)) + + for element_tag in model.getEleTags(): + connectivity, coordinates = element_coordinates(model, element_tag) + + corner_coordinates = coordinates[[0, 1, 2, 3, 0], :] + ax.plot(corner_coordinates[:, 0], corner_coordinates[:, 1], linewidth=1.5) + ax.scatter(coordinates[:, 0], coordinates[:, 1], s=16) + + centroid = coordinates[:4].mean(axis=0) + ax.text( + centroid[0], + centroid[1], + str(element_tag), + ha="center", + va="center", + ) + + for node_tag, xy in zip(connectivity, coordinates, strict=True): + ax.text( + xy[0], + xy[1], + str(node_tag), + fontsize=7, + ha="left", + va="bottom", + ) + + ax.set_aspect("equal", adjustable="box") + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_title(f"{element_name} distorted five-element patch") + ax.grid(True, linewidth=0.4, alpha=0.5) + + fig.tight_layout() + fig.savefig(path, dpi=300) + plt.close(fig) + + return path + + +def plot_strain_errors( + sampled: dict[str, np.ndarray], + expected: dict[str, np.ndarray], + element_name: str, + output_directory: Path, +) -> Path: + """Write a diagnostic plot of absolute strain-recovery errors.""" + output_directory.mkdir(parents=True, exist_ok=True) + path = output_directory / f"{element_name}_distorted_patch_strain_errors.png" + + fig, ax = plt.subplots(figsize=(7.0, 4.5)) + + for name in STRAIN_NAMES: + error = np.abs(sampled[name] - expected[name]) + ax.plot( + np.arange(error.size), + error, + marker="o", + linestyle="none", + label=name, + ) + + ax.set_yscale("log") + ax.set_xlabel("sample index") + ax.set_ylabel("absolute error") + ax.set_title(f"{element_name} distorted patch strain-recovery errors") + ax.grid(True, which="both", linewidth=0.4, alpha=0.5) + ax.legend(frameon=False) + + fig.tight_layout() + fig.savefig(path, dpi=300) + plt.close(fig) + + return path + + +def write_strain_report( + sampled: dict[str, np.ndarray], + expected: dict[str, np.ndarray], + element_name: str, + mesh_order: int, + geometry_plot_path: Path, + error_plot_path: Path, + output_directory: Path, +) -> Path: + """Write a compact strain-recovery report.""" + output_directory.mkdir(parents=True, exist_ok=True) + path = output_directory / f"{element_name}_distorted_patch_strain_report.txt" + + lines = [ + f"{element_name} distorted five-element patch test", + "", + f"Element order: {mesh_order}", + f"Geometry plot: {geometry_plot_path}", + f"Error plot: {error_plot_path}", + "", + "Maximum absolute errors:", + ] + + for name in STRAIN_NAMES: + error = np.abs(sampled[name] - expected[name]) + lines.append(f" {name:10s}: {error.max():.16e}") + + path.write_text("\n".join(lines), encoding="utf-8") + return path + + +def save_diagnostics(element_name: str, mesh_order: int) -> None: + """Save geometry, strain-error plot, and text report.""" + model = build_model_or_skip(element_name, mesh_order) + field = LinearPatchField() + + sampled, expected = sample_generalized_strains( + model=model, + element_name=element_name, + mesh_order=mesh_order, + field=field, + ) + + output_directory = OUTPUT_DIR / "patch_diagnostics" + + geometry_plot_path = plot_patch_geometry( + model=model, + element_name=element_name, + output_directory=output_directory, + ) + error_plot_path = plot_strain_errors( + sampled=sampled, + expected=expected, + element_name=element_name, + output_directory=output_directory, + ) + report_path = write_strain_report( + sampled=sampled, + expected=expected, + element_name=element_name, + mesh_order=mesh_order, + geometry_plot_path=geometry_plot_path, + error_plot_path=error_plot_path, + output_directory=output_directory, + ) + + assert geometry_plot_path.exists() + assert error_plot_path.exists() + assert report_path.exists() + + +@pytest.mark.parametrize("element_name,mesh_order", sorted(ELEMENT_ORDER.items())) +def test_distorted_patch_exact_strain_recovery( + element_name: str, + mesh_order: int, +) -> None: + """Verify exact generalized-strain recovery on the distorted patch.""" + assert_exact_strain_recovery(element_name, mesh_order) + + +@pytest.mark.parametrize("element_name,mesh_order", sorted(ELEMENT_ORDER.items())) +def test_distorted_patch_output_files_saved( + element_name: str, + mesh_order: int, +) -> None: + """Save geometry, strain-error plot, and text report for manual inspection.""" + save_diagnostics(element_name, mesh_order) \ No newline at end of file diff --git a/benchmarks/shell-0005/main.py b/benchmarks/shell-0005/main.py index ca87c48..07076c7 100644 --- a/benchmarks/shell-0005/main.py +++ b/benchmarks/shell-0005/main.py @@ -19,6 +19,8 @@ Transverse deflection w at point A = (375, 60), the bottom-right hole corner. """ +from __future__ import annotations + import argparse import xara @@ -26,146 +28,248 @@ from xara.load import Line, SurfaceLoad -# Geometry (mm) -W, H = 500.0, 300.0 -HX0, HX1 = 125.0, 375.0 -HY0, HY1 = 60.0, 240.0 -THICKNESS = 20.0 - -# Material -E_MOD = 200_000.0 # N/mm^2 -NU = 0.25 - -# Load -Q_LINE = -1000.0 # N/mm along -z - +# --- Geometry (mm) ----------------------------------------------------------- +PLATE_WIDTH_MM = 500.0 +PLATE_HEIGHT_MM = 300.0 + +# Rectangular hole: axis-aligned, centred in the plate. +HOLE_X_MIN_MM = 125.0 +HOLE_X_MAX_MM = 375.0 +HOLE_Y_MIN_MM = 60.0 +HOLE_Y_MAX_MM = 240.0 + +THICKNESS_MM = 20.0 + +# Benchmark point A: bottom-right corner of the hole (see module docstring). +POINT_A_X_MM = HOLE_X_MAX_MM +POINT_A_Y_MM = HOLE_Y_MIN_MM + +# --- Material --------------------------------------------------------------- +E_MODULUS_MPA = 200_000.0 # N/mm^2 +POISSON_RATIO = 0.25 + +# --- Loading --------------------------------------------------------------- +# Uniform line load on the inner top edge of the hole, global +z component (N/mm). +LINE_LOAD_Z_N_PER_MM = -1000.0 + +# Nodes closer than this (mm) to a boundary line are treated as on that line. +COORD_MATCH_TOL_MM = 1e-6 + + +# Eight quadrilateral patches cover the plate; the hole is left empty. +# Keys: row (B=bottom, M=middle, T=top) + column (L=left, M=middle, R=right). +# Each value is four (x, y) corners in counter-clockwise order on the mid-surface. +_PATCH_CORNERS_MM: dict[str, list[tuple[float, float]]] = { + "BL": [ + (0.0, 0.0), + (HOLE_X_MIN_MM, 0.0), + (HOLE_X_MIN_MM, HOLE_Y_MIN_MM), + (0.0, HOLE_Y_MIN_MM), + ], + "BM": [ + (HOLE_X_MIN_MM, 0.0), + (HOLE_X_MAX_MM, 0.0), + (HOLE_X_MAX_MM, HOLE_Y_MIN_MM), + (HOLE_X_MIN_MM, HOLE_Y_MIN_MM), + ], + "BR": [ + (HOLE_X_MAX_MM, 0.0), + (PLATE_WIDTH_MM, 0.0), + (PLATE_WIDTH_MM, HOLE_Y_MIN_MM), + (HOLE_X_MAX_MM, HOLE_Y_MIN_MM), + ], + "ML": [ + (0.0, HOLE_Y_MIN_MM), + (HOLE_X_MIN_MM, HOLE_Y_MIN_MM), + (HOLE_X_MIN_MM, HOLE_Y_MAX_MM), + (0.0, HOLE_Y_MAX_MM), + ], + "MR": [ + (HOLE_X_MAX_MM, HOLE_Y_MIN_MM), + (PLATE_WIDTH_MM, HOLE_Y_MIN_MM), + (PLATE_WIDTH_MM, HOLE_Y_MAX_MM), + (HOLE_X_MAX_MM, HOLE_Y_MAX_MM), + ], + "TL": [ + (0.0, HOLE_Y_MAX_MM), + (HOLE_X_MIN_MM, HOLE_Y_MAX_MM), + (HOLE_X_MIN_MM, PLATE_HEIGHT_MM), + (0.0, PLATE_HEIGHT_MM), + ], + "TM": [ + (HOLE_X_MIN_MM, HOLE_Y_MAX_MM), + (HOLE_X_MAX_MM, HOLE_Y_MAX_MM), + (HOLE_X_MAX_MM, PLATE_HEIGHT_MM), + (HOLE_X_MIN_MM, PLATE_HEIGHT_MM), + ], + "TR": [ + (HOLE_X_MAX_MM, HOLE_Y_MAX_MM), + (PLATE_WIDTH_MM, HOLE_Y_MAX_MM), + (PLATE_WIDTH_MM, PLATE_HEIGHT_MM), + (HOLE_X_MAX_MM, PLATE_HEIGHT_MM), + ], +} -# 8-block decomposition of the plate, leaving the centred hole empty. -# Naming: row letter (B=bottom, M=middle, T=top) then column letter (L, M, R). -_BLOCKS = { - "BL": [(0.0, 0.0), (HX0, 0.0), (HX0, HY0), (0.0, HY0)], - "BM": [(HX0, 0.0), (HX1, 0.0), (HX1, HY0), (HX0, HY0)], - "BR": [(HX1, 0.0), (W, 0.0), (W, HY0), (HX1, HY0)], - "ML": [(0.0, HY0), (HX0, HY0), (HX0, HY1), (0.0, HY1)], - "MR": [(HX1, HY0), (W, HY0), (W, HY1), (HX1, HY1)], - "TL": [(0.0, HY1), (HX0, HY1), (HX0, H ), (0.0, H )], - "TM": [(HX0, HY1), (HX1, HY1), (HX1, H ), (HX0, H )], - "TR": [(HX1, HY1), (W, HY1), (W, H ), (HX1, H )], +# Physical width/height of each patch column (L/M/R) and row (B/M/T), in mm. +_PATCH_COLUMN_WIDTH_MM = { + "L": HOLE_X_MIN_MM, + "M": HOLE_X_MAX_MM - HOLE_X_MIN_MM, + "R": PLATE_WIDTH_MM - HOLE_X_MAX_MM, +} +_PATCH_ROW_HEIGHT_MM = { + "B": HOLE_Y_MIN_MM, + "M": HOLE_Y_MAX_MM - HOLE_Y_MIN_MM, + "T": PLATE_HEIGHT_MM - HOLE_Y_MAX_MM, } -_COLUMN_WIDTHS = {"L": HX0, "M": HX1 - HX0, "R": W - HX1} -_ROW_HEIGHTS = {"B": HY0, "M": HY1 - HY0, "T": H - HY1} +def _surface_subdivisions( + patch_key: str, target_edge_length_mm: float +) -> tuple[int, int]: + """ + Number of finite elements along each parent direction for one patch. -def _divs(name, h): + ``patch_key`` is two letters: row (B/M/T) then column (L/M/R). The counts + are chosen so that the shorter patch dimension is near ``target_edge_length_mm``. """ - Return element counts (nx, ny) for a block, targeting edge length h. + row_letter, col_letter = patch_key[0], patch_key[1] + num_elem_x = max( + 1, round(_PATCH_COLUMN_WIDTH_MM[col_letter] / target_edge_length_mm) + ) + num_elem_y = max( + 1, round(_PATCH_ROW_HEIGHT_MM[row_letter] / target_edge_length_mm) + ) + return num_elem_x, num_elem_y + + +def _corner_points_for_xara_surface( + corners_xy: list[tuple[float, float]], +) -> dict[int, list[float]]: + """Map 2-D patch corners to Xara's 1-based point dict with z = 0.""" + return { + index + 1: [float(x), float(y), 0.0] + for index, (x, y) in enumerate(corners_xy) + } + + +def build_model(element: str, order: int, target_edge_length_mm: float) -> xara.Model: + """ + Assemble the mesh, boundary conditions, and line load for this benchmark. + + Parameters + ---------- + element + Shell element type passed to ``Model.surface``. + order + Parent-element order (1 for bilinear, 2 for biquadratic). + target_edge_length_mm + Target physical edge length (mm); actual spacing is rounded per patch. + + Returns + ------- + xara.Model + Fully constrained and loaded 3-D shell model (6 DOFs per node). """ - row, col = name[0], name[1] - nx = max(1, round(_COLUMN_WIDTHS[col] / h)) - ny = max(1, round(_ROW_HEIGHTS[row] / h)) - return nx, ny - - -def _surface_points(corners): - """4 corner tuples -> xara surface points dict, padded to z = 0.""" - return {i + 1: [float(p[0]), float(p[1]), 0.0] for i, p in enumerate(corners)} - - -def build_model(element, order, h): model = xara.Model(ndm=3, ndf=6) + material = xara.TriaxialMaterial( + "ElasticIsotropic", E=E_MODULUS_MPA, nu=POISSON_RATIO + ) + model.material(material) - # - # Material and section - # - mat = xara.TriaxialMaterial("ElasticIsotropic", E=E_MOD, nu=NU) - model.material(mat) + section = xara.ShellSection("Elastic", material, THICKNESS_MM) + model.section(section) - sec = xara.ShellSection("Elastic", mat, THICKNESS) - model.section(sec) - - # - # Mesh each block - # - # xara reuses nodes at coincident block corners. - for name, corners in _BLOCKS.items(): + # Coincident patch corners share nodes automatically. + for patch_key, corners_xy in _PATCH_CORNERS_MM.items(): model.surface( - _divs(name, h), + _surface_subdivisions(patch_key, target_edge_length_mm), element=element, - args={"section": sec}, + args={"section": section}, order=order, - points=_surface_points(corners), + points=_corner_points_for_xara_surface(corners_xy), ) - # - # Boundary - # - - # Clamp the left (x = 0) and top (y = H) outer edges - tol = 1e-6 - for tag in model.getNodeTags(): - x, y, _ = model.nodeCoord(tag) - if abs(x) < tol or abs(y - H) < tol: - model.fix(tag, (1, 1, 1, 1, 1, 1)) + # Clamped outer edges: left (x = 0) and top (y = plate height). + for node_tag in model.getNodeTags(): + x_mm, y_mm, _z_mm = model.nodeCoord(node_tag) + if abs(x_mm) < COORD_MATCH_TOL_MM or abs( + y_mm - PLATE_HEIGHT_MM + ) < COORD_MATCH_TOL_MM: + model.fix(node_tag, (1, 1, 1, 1, 1, 1)) + # Mindlin plate elements carry a drilling DOF; pin it at every node. if "plate" in element.lower(): - for tag in model.getNodeTags(): - # fix drill + for node_tag in model.getNodeTags(): try: - model.fix(tag, (0, 0, 0, 0, 0, 1)) - except: + model.fix(node_tag, (0, 0, 0, 0, 0, 1)) + except Exception: + # Node may already have a conflicting fixity from the edge clamps. pass - # - # Loading - # - - # collect nodes on the inner top edge of the cut-out: - # y = HY1, HX0 <= x <= HX1 - edge_nodes = [ - tag for tag in model.getNodeTags() - if abs(model.nodeCoord(tag)[1] - HY1) < tol - and HX0 - tol <= model.nodeCoord(tag)[0] <= HX1 + tol + # Inner top edge of the hole: y = HOLE_Y_MAX_MM, x between hole x-bounds. + cutout_top_edge_node_tags = [ + tag + for tag in model.getNodeTags() + if abs(model.nodeCoord(tag)[1] - HOLE_Y_MAX_MM) < COORD_MATCH_TOL_MM + and HOLE_X_MIN_MM - COORD_MATCH_TOL_MM + <= model.nodeCoord(tag)[0] + <= HOLE_X_MAX_MM + COORD_MATCH_TOL_MM ] - # Sort by x-coordinate to ensure consistent load application direction - edge_nodes.sort(key=lambda n: model.nodeCoord(n)[0]) + cutout_top_edge_node_tags.sort(key=lambda t: model.nodeCoord(t)[0]) - def q(*_): - return [0.0, 0.0, Q_LINE] + def uniform_line_load_global(*_coords: object) -> list[float]: + """Line-load intensity callback: constant (0, 0, q_z) in global axes.""" + return [0.0, 0.0, LINE_LOAD_Z_N_PER_MM] - load = SurfaceLoad(Line(model, edge_nodes), q) - model.pattern(xara.StaticPattern([load])) + distributed_load = SurfaceLoad( + Line(model, cutout_top_edge_node_tags), uniform_line_load_global + ) + model.pattern(xara.StaticPattern([distributed_load])) return model -def main(): +def main() -> None: parser = argparse.ArgumentParser(description=__doc__.splitlines()[1]) - parser.add_argument("-e", "--element", default="ASDShellQ4", - choices=["ASDShellQ4", "ShellMITC4", "HeterosisPlate", "ShellMITC9"], - help="shell element name (default: ASDShellQ4)") - parser.add_argument("--order", type=int, default=1, choices=(1, 2), - help="mesh order, 1 for Q4 or 2 for Q9 (default: 1)") - parser.add_argument("--h", type=float, default=25.0, - help="target element edge length in mm (default: 25)") + parser.add_argument( + "-e", + "--element", + default="ASDShellQ4", + choices=["ASDShellQ4", "ShellMITC4", "HeterosisPlate", "ShellMITC9"], + help="shell element name (default: ASDShellQ4)", + ) + parser.add_argument( + "--order", + type=int, + default=1, + choices=(1, 2), + help="mesh order, 1 for Q4 or 2 for Q9 (default: 1)", + ) + parser.add_argument( + "--h", + type=float, + default=25.0, + help="target element edge length in mm (default: 25)", + ) args = parser.parse_args() model = build_model(args.element, args.order, args.h) print(f" Element: {args.element}") print(f" {len(model.getNodeTags())} nodes") + analysis = xara.StaticAnalysis(model, system="Umfpack") print(analysis) analysis.analyze() - a_tag = find_node(model, x=HX1, y=HY0) - w_a = model.nodeDisp(a_tag)[2] + node_a = find_node(model, x=POINT_A_X_MM, y=POINT_A_Y_MM) + w_vertical_mm = model.nodeDisp(node_a)[2] - print(f" w at point A = {w_a:.6e} mm") + print(f" w at point A = {w_vertical_mm:.6e} mm") import veux + veux.serve(veux.render(model)) if __name__ == "__main__": main() -