A cantilever under a tip load in 2D and 3D, with four stress invariants of the 2D solve.
uv run python examples/cli.py run linear_elastic
2D triangles 9452 3D tetrahedra 4860 3D degrees of freedom 4116 3D peak deflection 0.6144
The functions that pose and solve the problem. The figures are below the fold.
"""A cantilever under a tip load, in 2D and 3D. `bend_2d` and `bend_3d` each state and solve one problem; `run` calls them and returns a `CantileverStudy` of plain results. Nothing here draws: `figures.py` does that from the `CantileverStudy`, and this file is what the gallery shows. """ from dataclasses import dataclass import numpy as np from fem.algebra.backends import IterativeBackend from fem.boundary import Dirichlet, Neumann from fem.conditions import Conditions from fem.mesh.mesh import Mesh from fem.mesh.structured import box_mesh from fem.physics.equations import LinearElastic from fem.post.solution import ElasticSolution from fem.regions import in_box, intersect, on_plane E, NU = 200.0, 0.4 def clamp_and_tip_load(width) -> Conditions: """Clamped on the left, pulled down over the middle of the right edge.""" return Conditions( Dirichlet(on_plane(0, 0.0), [0, 0]), # Transverse, so the beam bends. Sized for a tip deflection near 9% of the span, # inside the small-strain regime. Neumann(intersect(on_plane(0, width), in_box([None, 0.2], [None, 0.8])), [0, -0.5]), ) def bend_2d(mesh: Mesh, bc: Conditions) -> ElasticSolution: """The 2D cantilever solve.""" return LinearElastic(E, NU).problem(mesh, bc).solve() def bend_3d(n_3d) -> tuple[Mesh, ElasticSolution]: """The same clamp-and-load, one dimension up. The same assembly, the equation reading the tetrahedron off the connectivity. AMG-CG rather than a direct factorization, whose fill-in hurts in 3D. """ box = box_mesh(corners=[[0, 0, 0], [4, 1, 1]], resolution=(4 * n_3d // 2, n_3d // 2, n_3d // 2)) bc_3d = Conditions( Dirichlet(on_plane(0, 0.0), [0, 0, 0]), Neumann(on_plane(0, 4.0), [0, 0, -0.5]), ) solution = LinearElastic(E, NU).problem(box, bc_3d).with_backend(IterativeBackend()).solve() return box, solution @dataclass class CantileverStudy: """Everything `run` computed, for the figures and the summary to read.""" mesh: Mesh bc: Conditions solution: ElasticSolution box: Mesh solution_3d: ElasticSolution @property def tip_3d(self) -> float: """The 3D solve's largest vertical deflection.""" return float(np.abs(self.solution_3d.component(2)).max()) @property def invariants(self) -> list[tuple[str, np.ndarray]]: """Rotation-invariant reductions of the 2D stress tensor: von Mises, mean normal stress, the Tresca measure, and the largest tensile principal value.""" s = self.solution return [ ('Von Mises', s.von_mises), ('Pressure', s.pressure), ('Max shear', s.max_shear), ('Max principal', s.principal_stress[:, -1]), ] def run(mesh: Mesh, n_3d=14) -> CantileverStudy: """Solve the cantilever in 2D on `mesh` and in 3D on a box `n_3d` deep.""" bc = clamp_and_tip_load(np.max(mesh.vertices[:, 0])) solution = bend_2d(mesh, bc) box, solution_3d = bend_3d(n_3d) return CantileverStudy(mesh, bc, solution, box, solution_3d)
"""The figures and summary of the cantilever demo, drawn from a `CantileverStudy`.""" from functools import partial from demo_registry import Demo, DemoResult, Figure from demos._charts import conditions_figure from demos.linear_elastic import physics from demos.linear_elastic.physics import CantileverStudy, run from fem.mesh.structured import box_mesh from fem.plot.plotter import Plotter def _fields_figure(s: CantileverStudy) -> Figure: fields = Plotter(1, 2, figsize=(10.5, 4.2), title='Linear elasticity in 2D and 3D') fields.plot(s.solution.deformed_mesh(), s.solution.von_mises, mode='colored', idx=(0, 0), label='von Mises stress', title=f'2D: {len(s.mesh.elements)} triangles') # Only the boundary surface is drawn. fields.plot(s.solution_3d.deformed_mesh(), s.solution_3d.von_mises, mode='solid', idx=(0, 1), label='von Mises stress', title=f'3D: {len(s.box.elements)} tetrahedra') return Figure( fields, 'The same clamp and load solved in 2D and 3D. The bending stress is largest ' 'at the clamp, with tension over the neutral axis and compression under it. ' 'The 3D solve carries a three-component displacement and recovers stress the ' 'same way, drawn on its boundary surface.', 'fields') def _invariants_figure(s: CantileverStudy) -> Figure: deformed = s.solution.deformed_mesh() invariants = Plotter(2, 2, title='Stress invariants of the same solve', panel_aspect=4.0) for i, (name, values) in enumerate(s.invariants): invariants.plot(deformed, values, mode='colored', idx=divmod(i, 2), title=name) return Figure( invariants, 'Four rotation-invariant reductions of the same 2D stress tensor: von Mises, ' 'mean normal stress, the Tresca measure, and the largest tensile principal ' 'value.', 'invariants') def _conditions_figure(s: CantileverStudy) -> Figure: return conditions_figure( s.mesh, s.bc, 'Clamped along the left edge, pulled down over the middle of the right one; ' 'everything between is traction-free. The 3D solve imposes the same clamp ' 'and tip load, one dimension up.', panel_aspect=4.0) def _summary(s: CantileverStudy) -> str: return (f'2D triangles {len(s.mesh.elements)}\n' f'3D tetrahedra {len(s.box.elements)}\n' f'3D degrees of freedom {3 * len(s.box.vertices)}\n' f'3D peak deflection {s.tip_3d:.4f}') def demo(mesh, **kwargs) -> DemoResult: """A cantilever under a tip load in 2D and 3D, with four stress invariants of the 2D solve.""" s = run(mesh, **kwargs) return DemoResult([ _fields_figure(s), _invariants_figure(s), _conditions_figure(s), ], text=_summary(s)) # The 2D cantilever whose domain this is, plus a 3D box the demo builds for itself. DEMO = Demo('linear_elastic', demo, section='Solids & structures', domain=partial(box_mesh, [[0.0, 0.0], [4.0, 1.0]], (140, 35)), smoke_kwargs={'n_3d': 6}, show_source=physics)