{ "cells": [ { "cell_type": "markdown", "id": "eb3f78b0", "metadata": {}, "source": [ "# Multiscale elasticity: resolve material structure through local equations\n", "\n", "\n", "The main workflow is **meshes → spaces → local equations → global balance → assemble → solve → fields and errors**. `MeshHierarchy` associates macro and local meshes; `bind_interface` and `LocalContext` own supported numbering, geometric orientation and native coordinate conversions. The mathematical forms, physical trace meaning, local modes and gauge remain explicit in the cells below. Fully manual/custom spaces use the same numerical owners; see the [custom-interface notebook](https://github.com/ipes-lncc/pymhm/blob/main/notebooks/foundations/operators/custom_interface.ipynb).\n", "\n", "Build a plane-strain primal MHM model explicitly with UFL energy forms, a `LocalContext`, a global `Equation`, and `bind_problem`. A heterogeneous solid under horizontal extension creates nonuniform strains that a small classical mesh misses. MHM keeps a small global mesh while resolving the material inside each macroelement.\n", "\n", "The baseline is separately assembled classical conforming displacement Galerkin, on three successively refined meshes. There is no analytical solution for this physical case. We measure the baseline's own refinement in displacement and physical Cauchy stress before comparing methods.\n", "\n", "Install `pymhm[notebooks,visualization]` and the compatible native DOLFINx/UFL backend, then open this notebook with `jupyter lab` in a writable directory. Physical coefficients, weak forms, local rigid-motion modes and conforming reference forms are declared below. Importable vector helpers handle field evaluation, norms, figures and archives. Mesh-associated bindings provide the supported coordinate maps. Optional DOLFINx/UFL imports belong to this notebook environment; the portable PyMHM core does not import them.\n", "\n", "The local rigid-motion complement and global traction coupling follow [Harder, Madureira and Valentin (2016)](https://doi.org/10.1051/m2an/2015046). The oscillatory material and extension loading define an original application here.\n", "\n", "The PyMHM distribution contains only the library. The first cell explicitly\n", "downloads a checksum-verified companion archive and prepares the declared data\n", "using its local Python helpers. You can inspect these support files in `ROOT`.\n", "Complete studies and historical replay retain their existing opt-in flags.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "fab23fb8", "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "import os\n", "import sys\n", "from pymhm.io.workspace import workspace_from_archive\n", "\n", "# Download verified support files; this operation does not execute them.\n", "COMPANION_URL = \"https://ipes-lncc.github.io/pymhm/downloads/2b46d35d598a626ec0cf24aae3095fc315ba3bb6687dd33814116e45d3e19a1a/multiscale_elasticity-companion.zip\"\n", "COMPANION_SHA256 = \"2b46d35d598a626ec0cf24aae3095fc315ba3bb6687dd33814116e45d3e19a1a\"\n", "WORKSPACE = Path(\n", " os.environ.get(\"PYMHM_WORKSPACE\", Path.cwd() / \".pymhm-companions\" / COMPANION_SHA256)\n", ")\n", "ROOT = workspace_from_archive(COMPANION_URL, sha256=COMPANION_SHA256, directory=WORKSPACE)\n", "os.environ[\"PYMHM_WORKSPACE\"] = str(ROOT)\n", "if str(ROOT) not in sys.path:\n", " sys.path.insert(0, str(ROOT))\n", "\n", "# Explicitly prepare declared data with the downloaded Python helpers.\n", "from scripts.notebook_reproduction import notebook_workspace\n", "\n", "ROOT = notebook_workspace(\"introduction/multiscale_elasticity.ipynb\", directory=ROOT)\n", "root = ROOT\n", "print(\"Workspace:\", ROOT)\n", "\n", "from functools import partial\n", "import numpy as np\n", "from numpy.typing import NDArray\n", "import matplotlib.pyplot as plt\n", "import ufl\n", "\n", "\n", "from pymhm import (\n", " Equation,\n", " LocalEquations,\n", " MeshHierarchy,\n", " LocalContext,\n", " assemble,\n", " columns,\n", " bind_interface,\n", " bind_problem,\n", ")\n", "from pymhm.execution.cpu import ExecutionConfig\n", "from pymhm.meshes.triangle import TriangleMesh\n", "from pymhm.fem.scalar.triangle import nodal_space\n", "from pymhm.fem.scalar.operators import boundary_data\n", "from pymhm.fem.traces.interval import FaceSpace, SkeletonSpace\n", "from examples.introduction.vector import (\n", " BrokenVectorEvaluator,\n", " evaluate_named_displacement,\n", " elasticity_reference,\n", " elasticity_errors,\n", " triangle_grid_quadrature,\n", " plot_convergence,\n", " plot_elasticity_fields,\n", " inspect_elasticity_solution,\n", ")\n", "\n", "Array = NDArray[np.float64]\n" ] }, { "cell_type": "markdown", "id": "aabc0ce0", "metadata": {}, "source": [ "## 1. Declare the material and physical loading\n", "\n", "On the unit square, impose a one-percent horizontal extension on the whole exterior, with no body force:\n", "\n", "$$\n", "-\\nabla\\cdot\\sigma=0,\\qquad\n", "\\sigma=2\\mu\\varepsilon(\\boldsymbol u)+\\lambda\\operatorname{div}(\\boldsymbol u)I,\n", "\\qquad\\boldsymbol u_D=(0.01x,0).\n", "$$\n", "\n", "Here $\\varepsilon(\\boldsymbol u)=(\\nabla\\boldsymbol u+\\nabla\\boldsymbol u^T)/2$. Use a smooth material with short wavelength and Lamé-modulus contrast $e^3\\simeq20.1$:\n", "\n", "$$\n", "\\lambda(x,y)=\\mu(x,y)=\\exp\\!\\left(1.5\\sin(16\\pi x)\\sin(16\\pi y)\\right).\n", "$$\n", "\n", "The material wavelength is $1/8$. Taking $\\lambda=\\mu$ gives Poisson ratio $1/4$ in plane strain, so this case isolates multiscale heterogeneity rather than near-incompressibility. The stiffness is uniformly positive on symmetric strains. Both methods solve the same tensor, domain and boundary values." ] }, { "cell_type": "code", "execution_count": null, "id": "8c0a8b6b", "metadata": {}, "outputs": [], "source": [ "def micro_modulus(points: Array) -> Array:\n", " \"\"\"Return smooth plane-strain Lamé values in [exp(-1.5),exp(1.5)], wavelength 1/8.\n", "\n", " Both lambda and mu equal this value, hence Poisson ratio is 1/4 and\n", " heterogeneity is separated from near-incompressibility and locking.\n", " \"\"\"\n", " return np.exp(1.5 * np.sin(16 * np.pi * points[:, 0]) * np.sin(16 * np.pi * points[:, 1]))\n", "\n", "\n", "def extension_boundary(points: Array) -> Array:\n", " \"\"\"Prescribe a one-percent horizontal extension, with zero vertical displacement.\"\"\"\n", " return np.column_stack((0.01 * points[:, 0], np.zeros(len(points))))\n" ] }, { "cell_type": "code", "execution_count": null, "id": "87ca3b96", "metadata": {}, "outputs": [], "source": [ "macro = TriangleMesh.unit_square(4)\n", "local_degree, local_refinement, trace_segments = 3, 16, 4\n", "skeleton = SkeletonSpace(\n", " macro, tuple(FaceSpace.uniform(1, trace_segments) for _ in macro.faces), components=2\n", ")\n", "print(\n", " {\n", " \"macro_triangles\": len(macro.cells),\n", " \"material_wavelength\": 1 / 8,\n", " \"material_contrast\": float(np.exp(3.0)),\n", " \"local_degree\": local_degree,\n", " \"local_refinement\": local_refinement,\n", " \"trace_degree\": 1,\n", " \"trace_segments\": trace_segments,\n", " }\n", ")\n" ] }, { "cell_type": "markdown", "id": "52ffb39f", "metadata": {}, "source": [ "## 2. Declare the finite element and bind it to the local mesh\n", "\n", "Choose a vector Basix/UFL element for displacement. `local.native_space` binds it to the current local mesh and provides its native function space. PyMHM owns the topology and coefficient conversions, so the local UFL forms do not contain numbering or orientation arrays. The later physical-norm and archive sections state the executed nodal basis explicitly.\n", "\n", "The mesh is the actual independent local refinement of a macro triangle. No global refined mesh is passed to MHM. Native forms use `MPI.COMM_SELF`; the generic form compiler copies their assembled arrays and releases native assembly resources." ] }, { "cell_type": "code", "execution_count": null, "id": "daf30ba6", "metadata": {}, "outputs": [], "source": [ "def rigid_modes(points: Array, center: Array) -> Array:\n", " \"\"\"Return two translations and centered counterclockwise rotation, shape (n,2,3).\"\"\"\n", " modes = np.zeros((len(points), 2, 3))\n", " modes[:, :, :2] = np.eye(2)\n", " modes[:, 0, 2] = -(points[:, 1] - center[1])\n", " modes[:, 1, 2] = points[:, 0] - center[0]\n", " return modes\n" ] }, { "cell_type": "markdown", "id": "b4a0b22e", "metadata": {}, "source": [ "## 3. Derive each local equation and identify its kernel\n", "\n", "Integration by parts yields\n", "\n", "$$\n", "\\int_T \\left[2\\mu\\varepsilon(\\boldsymbol u):\\varepsilon(\\boldsymbol v)\n", "+\\lambda\\operatorname{div}\\boldsymbol u\\operatorname{div}\\boldsymbol v\\right]\n", "-\\langle\\sigma\\boldsymbol n_T,\\boldsymbol v\\rangle_{\\partial T}=0.\n", "$$\n", "\n", "Declare the skeletal multiplier $\\boldsymbol\\lambda=-\\sigma\\boldsymbol n_F$: it is **negative physical Cauchy traction** in the unique macroface orientation. Thus the local equation is $A_T\\boldsymbol u+B_T\\boldsymbol\\lambda=0$, with the signed trace map $B_T$.\n", "\n", "The local energy has three exact rigid modes: translations $(1,0)$ and $(0,1)$, and centered rotation $(-(y-y_T),x-x_T)$. Their physical $L^2$ moments select a unique response in the complement, while three retained coordinates represent the rigid component. We declare these modes, rather than infer them from a numerically rotated nullspace.\n", "\n", "The trace contains the restrictions of rigid motions: degree one. The local degree is three. [Harder, Madureira and Valentin (2016), Lemma 6.6](https://doi.org/10.1051/m2an/2015046) give the sufficient rule $k\\ge\\ell+2$ for odd-degree traces on a single triangular local element. Our refined and segmented spaces differ from that single-element setting, so the rule alone is not their stability theorem. Segments align with the fine boundary and each covers four fine edges. We additionally check the rank of the actual trace pairing; an algebraic residual alone would not establish compatibility.\n", "\n", "| Weak-form term | Visible code |\n", "| --- | --- |\n", "| $2\\mu\\varepsilon(u):\\varepsilon(v)+\\lambda\\operatorname{div}u\\operatorname{div}v$ | UFL expression `a` |\n", "| $\\langle s_{TF}\\lambda,v\\rangle$ | `local.trace_pairings(lambda phi, ds: inner(phi,v)*ds)` |\n", "| Rigid kernel | `kernel` |\n", "| Physical rigid moments | `columns` of UFL linear forms |\n", "| Global displacement continuity | independent UFL pairing `-inner(phi,u)*ds`, with `axis=\"rows\"` |" ] }, { "cell_type": "code", "execution_count": null, "id": "a2ec35ce", "metadata": {}, "outputs": [], "source": [ "def local_elasticity(local: LocalContext) -> LocalEquations:\n", " \"\"\"Declare plane-strain UFL energy, negative Cauchy traction and rigid moments.\"\"\"\n", " cell, fine = local.cell, local.mesh\n", " import basix\n", " import basix.ufl\n", "\n", " element = basix.ufl.element(\n", " \"Lagrange\",\n", " \"triangle\",\n", " local_degree,\n", " lagrange_variant=basix.LagrangeVariant.equispaced,\n", " shape=(2,),\n", " )\n", " binding = local.native_space(element)\n", " domain, space, mapping = binding.mesh, binding.space, binding.mapping\n", " u, v = ufl.TrialFunction(space), ufl.TestFunction(space)\n", " x = ufl.SpatialCoordinate(domain)\n", " mu = ufl.exp(1.5 * ufl.sin(16 * np.pi * x[0]) * ufl.sin(16 * np.pi * x[1]))\n", " lam = mu\n", " dx = ufl.Measure(\"dx\", domain=domain, metadata={\"quadrature_degree\": 16})\n", " # This is exactly the physical energy written above.\n", " a = (\n", " 2 * mu * ufl.inner(ufl.sym(ufl.grad(u)), ufl.sym(ufl.grad(v)))\n", " + lam * ufl.div(u) * ufl.div(v)\n", " ) * dx\n", " _, nodes = nodal_space(fine, local_degree)\n", " center = macro.points[macro.cells[cell]].mean(axis=0)\n", " portable_kernel = rigid_modes(nodes, center).reshape(-1, 3)\n", " kernel = binding.to_native(portable_kernel)\n", " rigid = (\n", " ufl.as_vector((1.0, 0.0)),\n", " ufl.as_vector((0.0, 1.0)),\n", " ufl.as_vector((-(x[1] - center[1]), x[0] - center[0])),\n", " )\n", " moments = columns(*(ufl.inner(v, mode) * dx for mode in rigid))\n", " local.field(\"displacement\", binding)\n", " b = local.trace_pairings(lambda phi, ds: ufl.inner(phi, v) * ds)\n", " c = local.trace_pairings(lambda phi, ds: -ufl.inner(phi, u) * ds, axis=\"rows\")\n", " return local.equations(\n", " a=a,\n", " L=np.zeros(len(mapping)),\n", " b=b,\n", " c=c,\n", " kernel=kernel,\n", " moments=moments,\n", " metadata=(fine, mapping),\n", " )\n" ] }, { "cell_type": "markdown", "id": "ad58cac9", "metadata": {}, "source": [ "## 4. Declare the global displacement equation and solve\n", "\n", "The global weak equation imposes displacement-jump moments on interior macrofaces and prescribed displacement moments on the exterior. With $C=-B^T$, its additional load is $-\\langle\\boldsymbol\\mu,\\boldsymbol u_D\\rangle_{\\partial\\Omega}$. Full displacement data remove the three global rigid ambiguities, so no additional physical gauge is required.\n", "\n", "Only the macroface and three retained rigid coordinates per macroelement enter the global matrix. The local fields are recovered in their executed native basis, then mapped to the portable nodal representation for physical evaluation." ] }, { "cell_type": "code", "execution_count": null, "id": "1b8236a0", "metadata": {}, "outputs": [], "source": [ "boundary, fixed = boundary_data(skeleton, extension_boundary, order=8)\n", "hierarchy = MeshHierarchy(\n", " macro, tuple(macro.submesh(cell, local_refinement) for cell in range(len(macro.cells)))\n", ")\n", "problem = bind_problem(\n", " hierarchy,\n", " bind_interface(skeleton, convention=\"normal\"),\n", " local_elasticity,\n", " global_equation=lambda global_problem: Equation(0, global_problem.trace_load(-boundary)),\n", " retained=3,\n", " fixed=fixed,\n", ")\n", "system = assemble(problem, execution=ExecutionConfig(\"serial\", native_threads=1))\n", "solution = system.solve()\n", "displacement_fields = solution.field(\"displacement\")\n", "inspect_elasticity_solution(system, solution, macro)\n" ] }, { "cell_type": "markdown", "id": "479467b5", "metadata": {}, "source": [ "## 5. Physical evaluation and norms in the executed bases\n", "\n", "The evaluators below locate points in the original triangular connectivity and use the package's Basix basis tables. Separate local fields retain one-sided macroface values. The physical stress is recomputed from the **same** material and symmetric strain, without smoothing.\n", "\n", "Measure displacement in $L^2$, Cauchy stress in Frobenius $L^2$, and the strain-energy difference:\n", "\n", "$$\n", "\\left(\\int_\\Omega\n", "2\\mu\\varepsilon(\\boldsymbol e_u):\\varepsilon(\\boldsymbol e_u)\n", "+\\lambda(\\operatorname{div}\\boldsymbol e_u)^2\\right)^{1/2}.\n", "$$\n", "\n", "Raw primal stress is symmetric but is not claimed to belong to $H(\\mathrm{div})$ or to equilibrate every fine cell. Skeletal rigid equations impose macro force and moment compatibility." ] }, { "cell_type": "code", "execution_count": null, "id": "3e0a236d", "metadata": {}, "outputs": [], "source": [ "# Norms and one-sided field sampling use the importable helpers above.\n", "# Their implementation is in examples/introduction/vector.py.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c6d67988", "metadata": {}, "outputs": [], "source": [ "mhm_evaluator = BrokenVectorEvaluator(\n", " macro, tuple(partial(evaluate_named_displacement, field) for field in displacement_fields)\n", ")\n", "# Resolve every reference/local interface on a common 256×256 square grid.\n", "error_points, error_weights = triangle_grid_quadrature(256, order=5)\n" ] }, { "cell_type": "markdown", "id": "3cf83169", "metadata": {}, "source": [ "## 6. Independent classical assembly, including a coarse comparison\n", "\n", "Assemble the UFL energy globally on one continuous $P_2$ displacement mesh. Strong Dirichlet elimination uses the prescribed extension. There is no MHM skeleton or local condensation. A coarse solve on the macro mesh illustrates the cost of unresolved material; it is not used as the reference.\n", "\n", "Three fine meshes, $64\\times64$, $128\\times128$, and $256\\times256$ squares split into triangles, verify the baseline's own refinement in physical fields. The reference shares DOLFINx/Basix with local integration but has an independent global conforming assembly; no independent external code is executed." ] }, { "cell_type": "code", "execution_count": null, "id": "e8288ab9", "metadata": {}, "outputs": [], "source": [ "def classical_energy(space, resolution):\n", " \"\"\"Assemble the same physical operator on a conforming P2 reference mesh.\"\"\"\n", " u, v = ufl.TrialFunction(space), ufl.TestFunction(space)\n", " x = ufl.SpatialCoordinate(space.ufl_domain())\n", " mu = ufl.exp(1.5 * ufl.sin(16 * np.pi * x[0]) * ufl.sin(16 * np.pi * x[1]))\n", " dx = ufl.Measure(\n", " \"dx\",\n", " domain=space.ufl_domain(),\n", " metadata={\"quadrature_degree\": 64 if resolution == 4 else 16},\n", " )\n", " return (\n", " 2 * mu * ufl.inner(ufl.sym(ufl.grad(u)), ufl.sym(ufl.grad(v)))\n", " + mu * ufl.div(u) * ufl.div(v)\n", " ) * dx\n", "\n", "\n", "references = {\n", " n: elasticity_reference(n, classical_energy, extension_boundary) for n in (4, 64, 128, 256)\n", "}\n", "coarse_evaluator = references[4][0]\n", "reference_evaluators = [references[n][0] for n in (64, 128, 256)]\n", "reference_rows = [references[n][1] for n in (64, 128, 256)]\n", "reference_rows\n" ] }, { "cell_type": "code", "execution_count": null, "id": "777cd395", "metadata": {}, "outputs": [], "source": [ "measure_error = partial(elasticity_errors, modulus=micro_modulus)\n", "reference_refinement = [\n", " measure_error(a, b, error_points, error_weights)\n", " for a, b in zip(reference_evaluators[:-1], reference_evaluators[1:])\n", "]\n", "reference = reference_evaluators[-1]\n", "mhm_errors = measure_error(mhm_evaluator, reference, error_points, error_weights)\n", "coarse_errors = measure_error(coarse_evaluator, reference, error_points, error_weights)\n", "print(\"Successive reference differences:\", reference_refinement)\n", "print(\"MHM versus fine reference:\", mhm_errors)\n", "print(\"Coarse Galerkin versus fine reference:\", coarse_errors)\n", "assert reference_refinement[-1][\"displacement_L2\"] < reference_refinement[0][\"displacement_L2\"]\n", "assert reference_refinement[-1][\"stress_L2\"] < reference_refinement[0][\"stress_L2\"]\n", "# Check error integration independently by increasing the common-grid Gauss rule.\n", "check_points, check_weights = triangle_grid_quadrature(256, order=7)\n", "quadrature_check = measure_error(mhm_evaluator, reference, check_points, check_weights)\n", "print(\"Higher-order norm quadrature:\", quadrature_check)\n" ] }, { "cell_type": "markdown", "id": "3c7a53c6", "metadata": {}, "source": [ "### Baseline refinement increments\n", "\n", "These measured displacement and stress differences compare successive conforming references. Their slopes describe refinement increments rather than errors against an analytical solution. Plotting both fields makes the reference uncertainty visible." ] }, { "cell_type": "code", "execution_count": null, "id": "411dca43", "metadata": {}, "outputs": [], "source": [ "figure = plot_convergence(\n", " [1 / 64, 1 / 128],\n", " {\n", " \"reference displacement L2 increment\": [r[\"displacement_L2\"] for r in reference_refinement],\n", " \"reference Cauchy stress L2 increment\": [r[\"stress_L2\"] for r in reference_refinement],\n", " },\n", ")\n", "for axis in figure.axes:\n", " axis.set_xlabel(\"Reference square width\")\n", "plt.show()\n", "print(\n", " \"Reference uncertainty/MHM difference ratios:\",\n", " {\n", " field: reference_refinement[-1][field] / mhm_errors[field]\n", " for field in (\"displacement_L2\", \"stress_L2\", \"energy\")\n", " },\n", ")\n" ] }, { "cell_type": "markdown", "id": "432a9112", "metadata": {}, "source": [ "## 7. See why a multiscale discretization helps\n", "\n", "The stiff regions resist strain; softer regions accommodate more displacement gradients. Compare the coarse classical stress with the fine reference and MHM stress, not just displacement magnitude. Each spatial panel overlays the actual MHM macro mesh. Fields on macroface intersections remain independent.\n", "\n", "The material has multiple oscillations inside a macro triangle. Its effect enters MHM through the local energy responses, without placing all reference-mesh coefficients in the global system." ] }, { "cell_type": "code", "execution_count": null, "id": "75447cbc", "metadata": {}, "outputs": [], "source": [ "plot_elasticity_fields(macro, mhm_evaluator, reference, coarse_evaluator, micro_modulus)\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "80f871c1", "metadata": {}, "source": [ "## 8. Scope of the verified result\n", "\n", "The baseline is finite, so its successive displacement and stress differences remain part of the evidence. Compare that uncertainty with MHM–reference differences. The higher-order norm quadrature checks the integration of fields that are discontinuous in gradient.\n", "\n", "This tutorial verifies a heterogeneous, compressible plane-strain application. It does not claim locking-free nearly incompressible behavior, fine-cell stress conservation, or reproduction of a particular published material experiment. Global dimensions alone do not establish a timing gain; local setup and response solves are part of the computational work.\n", "\n", "The bound coefficient convention, declared rigid basis, physical traction sign, moment rows, Dirichlet convention and approximation spaces make the numerical result auditable. Changing the material or boundary condition requires revisiting those choices before interpreting a small algebraic residual." ] }, { "cell_type": "markdown", "id": "references-multiscale_elasticity", "metadata": {}, "source": [ "## References\n", "\n", "- Christopher Harder, Alexandre L. Madureira, and Frédéric Valentin (2016). *A hybrid-mixed method for elasticity*, ESAIM: M2AN 50, 311–336. [DOI: 10.1051/m2an/2015046](https://doi.org/10.1051/m2an/2015046)." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python" }, "pymhm_editorial_revision": { "date": "2026-10-06", "scope": "Markdown citations and bibliography only", "preserved": "Code cells, execution counts, execution timestamps and numerical outputs", "preserved_code_cells_sha256": "f4359127d7c6a015351e3b37287006af600512a78a245e5698b6f3a389e47496" } }, "nbformat": 4, "nbformat_minor": 5 }