{ "cells": [ { "cell_type": "markdown", "id": "7a4285c5", "metadata": {}, "source": [ "# Mixed stress MHM: write H(div) local and global equations\n", "\n", "This notebook uses row-wise BDM stress, discontinuous displacement and weak rotation. You write the compliance, force balance, symmetry and normal-traction equations before assembling. Basix-based package kernels handle finite-element integration and Piola maps. The method follows [Devloo et al. (2021)](https://doi.org/10.1051/m2an/2021013), with the stable stress spaces of [Arnold, Falk and Winther (2007)](https://doi.org/10.1090/S0025-5718-07-01998-9).\n", "\n", "Use the locked `introduction` Pixi environment for the independent classical UFL control. The multiscale equations themselves use the portable NumPy/SciPy backend." ] }, { "cell_type": "code", "execution_count": null, "id": "companion-bootstrap", "metadata": { "tags": [ "pymhm-companion-bootstrap" ] }, "outputs": [], "source": [ "import os\n", "import sys\n", "from pathlib import Path\n", "\n", "from pymhm.io.workspace import workspace_from_archive\n", "\n", "# Download checksum-verified support files; extraction does not execute them.\n", "COMPANION_URL = \"https://ipes-lncc.github.io/pymhm/downloads/17bd46e20ba59df4ea3f71dc01f26671558349962403336dc84074d614bcac52/mixed_stress_mhm_workflow-companion.zip\"\n", "COMPANION_SHA256 = \"17bd46e20ba59df4ea3f71dc01f26671558349962403336dc84074d614bcac52\"\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", "# Prepare only this notebook's declared inputs, without executing its equations.\n", "from scripts.notebook_reproduction import notebook_workspace\n", "\n", "ROOT = notebook_workspace(\"elasticity/mixed_stress_mhm_workflow.ipynb\", directory=ROOT)" ] }, { "cell_type": "code", "execution_count": null, "id": "ad83e55b", "metadata": {}, "outputs": [], "source": [ "import os\n", "import sys\n", "from pathlib import Path\n", "\n", "ROOT = Path(os.environ.get(\"PYMHM_WORKSPACE\", Path.cwd())).resolve()\n", "if not (ROOT / \"examples\").is_dir():\n", " ROOT = next(parent for parent in ROOT.parents if (parent / \"examples\").is_dir())\n", "os.environ[\"PYMHM_WORKSPACE\"] = str(ROOT)\n", "if str(ROOT) not in sys.path:\n", " sys.path.insert(0, str(ROOT))\n", "print(\"Workspace:\", ROOT)" ] }, { "cell_type": "code", "execution_count": null, "id": "d3776cf5", "metadata": {}, "outputs": [], "source": [ "from dataclasses import replace\n", "\n", "import basix.ufl\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import ufl\n", "from scipy import sparse\n", "from threadpoolctl import threadpool_limits\n", "\n", "from pymhm import (\n", " Equation,\n", " FaceSpace,\n", " MeshHierarchy,\n", " MultiscaleProblem,\n", " SkeletonSpace,\n", " TriangleMesh,\n", " assemble,\n", " bind_interface,\n", " bind_problem,\n", " solve,\n", ")\n", "from pymhm.core.equations import compile_form\n", "from pymhm.fem.hdiv.bdm_family import BDMFamily\n", "from pymhm.fem.scalar.triangle import multiindices, nodal_space\n", "from pymhm.fem.vector.stress import (\n", " displacement_rigid_moments,\n", " mixed_elasticity_operators,\n", " rigid_values,\n", " stress_trace_moments,\n", " traction_mapping,\n", ")\n", "from pymhm.postprocessing.nodal import nodal_field\n", "from pymhm.postprocessing.piola import hdiv_field\n", "from pymhm.postprocessing.stress import MixedElasticitySolution" ] }, { "cell_type": "markdown", "id": "f0014a91", "metadata": {}, "source": [ "## 1. State physical fields and boundary data\n", "\n", "On the unit square, choose a smooth solenoidal displacement:\n", "\n", "$$\n", "u=(\\sin(\\pi x)\\cos(\\pi y),-\\cos(\\pi x)\\sin(\\pi y)),\\qquad \\mu=1,\\quad\\lambda=1.\n", "$$\n", "\n", "Then $p=-\\lambda\\nabla\\cdot u=0$, $\\sigma=2\\mu\\varepsilon(u)$ and $f=2\\pi^2u$. The weak rotation convention is $r=(\\partial_yu_1-\\partial_xu_2)/2$. The exact boundary displacement is nonzero on parts of the exterior; it must enter the global functional." ] }, { "cell_type": "code", "execution_count": null, "id": "e0d7535d", "metadata": {}, "outputs": [], "source": [ "def smooth_fields(points):\n", " \"\"\"Define independent solenoidal displacement, physical stress, force and weak rotation.\"\"\"\n", " x, y = np.pi * points.T\n", " u = np.column_stack((np.sin(x) * np.cos(y), -np.cos(x) * np.sin(y)))\n", " sigma = np.zeros((len(points), 2, 2))\n", " sigma[:, 0, 0] = 2 * np.pi * np.cos(x) * np.cos(y)\n", " sigma[:, 1, 1] = -sigma[:, 0, 0]\n", " return u, sigma, 2 * np.pi**2 * u, -np.pi * np.sin(x) * np.sin(y)\n", "\n", "\n", "def exact_u(points):\n", " \"\"\"Return the nonhomogeneous prescribed displacement.\"\"\"\n", " return smooth_fields(points)[0]\n", "\n", "\n", "def exact_ufl(domain):\n", " \"\"\"State the same exact physical data for native classical assembly.\"\"\"\n", " x = ufl.SpatialCoordinate(domain)\n", " ue = ufl.as_vector(\n", " (\n", " ufl.sin(np.pi * x[0]) * ufl.cos(np.pi * x[1]),\n", " -ufl.cos(np.pi * x[0]) * ufl.sin(np.pi * x[1]),\n", " )\n", " )\n", " return ue, ufl.as_ufl(0), 2 * np.pi**2 * ue" ] }, { "cell_type": "markdown", "id": "1782cf69", "metadata": {}, "source": [ "## 2. Declare compatible spaces and independent local meshes\n", "\n", "Stress rows use BDM2; displacement and rotation use discontinuous $P_1$. They contain all three two-dimensional rigid modes. BDM1/P0/P0 omits the rigid rotation and is inadmissible for this local MHM construction.\n", "\n", "Interior macrofaces use vector $P_1$ negative traction. Exterior faces retain full $P_2$ moments on each fine boundary edge, as required to represent the prescribed normal-stress space in this comparison. The normal/interior degrees can be enriched independently through `BDMFamily`; its divergence and rotation degrees change with the selected interior family." ] }, { "cell_type": "code", "execution_count": null, "id": "8099ed19", "metadata": {}, "outputs": [], "source": [ "n, refinement = 8, 2\n", "lam, mu = 1.0, 1.0\n", "macro = TriangleMesh.unit_square(n)\n", "boundary = set(macro.boundary_faces)\n", "family = BDMFamily(2, 0)\n", "skeleton = SkeletonSpace(\n", " macro,\n", " tuple(\n", " FaceSpace.uniform(2, refinement) if f in boundary else FaceSpace.uniform(1)\n", " for f in range(len(macro.faces))\n", " ),\n", " 2,\n", ")\n", "hierarchy = MeshHierarchy(\n", " macro, tuple(macro.submesh(i, refinement) for i in range(len(macro.cells)))\n", ")\n", "interface = bind_interface(skeleton, convention=\"normal\")\n", "\n", "\n", "def source(points):\n", " \"\"\"Return the independently derived physical body force.\"\"\"\n", " return smooth_fields(points)[2]\n", "\n", "\n", "def displacement(points):\n", " \"\"\"Return the physical displacement including nonzero exterior values.\"\"\"\n", " return smooth_fields(points)[0]" ] }, { "cell_type": "markdown", "id": "8205d6a3", "metadata": {}, "source": [ "## 3. Identify each mathematical block of the local operator\n", "\n", "For $\\operatorname{asym}\\tau=\\tau_{12}-\\tau_{21}$, the mixed forms are\n", "\n", "$$\n", "\\begin{aligned}\n", "(\\mathcal A\\sigma,\\tau)_T+(u,\\nabla\\cdot\\tau)_T\n", " +(r,\\operatorname{asym}\\tau)_T-\\langle \\widehat u,\\tau n\\rangle_{\\partial T}&=0,\\\\\n", "(\\nabla\\cdot\\sigma,v)_T&=-(f,v)_T,\\\\\n", "(\\operatorname{asym}\\sigma,s)_T&=0,\\\\\n", "\\langle\\sigma n+\\lambda_H,w\\rangle_{\\partial T}&=0.\n", "\\end{aligned}\n", "$$\n", "\n", "Here $\\widehat u$ is an auxiliary local boundary displacement and $\\lambda_H=-\\sigma n$ is the global negative traction. The compliance in two dimensions is\n", "\n", "$$\n", "\\mathcal A\\sigma=\\frac{1}{2\\mu}\n", "\\left(\\sigma-\\frac{\\lambda}{2(\\mu+\\lambda)}\\operatorname{tr}(\\sigma)I\\right).\n", "$$\n", "\n", "`mixed_elasticity_operators` integrates precisely the compliance mass `mass`, row divergence `divergence`, weak symmetry `asymmetry` and volume force `force`. The following UFL expression is the same three volume equations; the actual portable solve below uses the shared kernels." ] }, { "cell_type": "code", "execution_count": null, "id": "df84820a", "metadata": {}, "outputs": [], "source": [ "def mixed_volume_form(space):\n", " \"\"\"State the independent row-wise stress/displacement/rotation UFL volume form.\"\"\"\n", " row0, row1, u, r = ufl.TrialFunctions(space)\n", " test0, test1, v, s = ufl.TestFunctions(space)\n", " sigma = ufl.as_matrix(((row0[0], row0[1]), (row1[0], row1[1])))\n", " tau = ufl.as_matrix(((test0[0], test0[1]), (test1[0], test1[1])))\n", " compliance = (sigma - lam / (2 * (mu + lam)) * ufl.tr(sigma) * ufl.Identity(2)) / (2 * mu)\n", " dx = ufl.Measure(\"dx\", domain=space.mesh, metadata={\"quadrature_degree\": 10})\n", " _, _, f = exact_ufl(space.mesh)\n", " a = (\n", " ufl.inner(compliance, tau)\n", " + ufl.inner(u, ufl.div(tau))\n", " + ufl.inner(ufl.div(sigma), v)\n", " + r * (tau[0, 1] - tau[1, 0])\n", " + s * (sigma[0, 1] - sigma[1, 0])\n", " ) * dx\n", " return Equation(a, -ufl.inner(f, v) * dx)" ] }, { "cell_type": "code", "execution_count": null, "id": "199ab253", "metadata": {}, "outputs": [], "source": [ "from pymhm.backends.spaces import bind_space\n", "\n", "native_volume = bind_space(\n", " hierarchy.local_mesh(0),\n", " basix.ufl.mixed_element(\n", " [\n", " basix.ufl.element(\"BDM\", \"triangle\", 2),\n", " basix.ufl.element(\"BDM\", \"triangle\", 2),\n", " basix.ufl.element(\n", " \"DG\", \"triangle\", 1, shape=(2,), lagrange_variant=basix.LagrangeVariant.equispaced\n", " ),\n", " basix.ufl.element(\n", " \"DG\", \"triangle\", 1, lagrange_variant=basix.LagrangeVariant.equispaced\n", " ),\n", " ]\n", " ),\n", ")\n", "try:\n", " native_equation = mixed_volume_form(native_volume.space)\n", " native_matrix = compile_form(native_equation.a)\n", " native_load = compile_form(native_equation.L)\n", " print(\"Independent UFL volume blocks:\", native_matrix.shape, native_load.shape)\n", "finally:\n", " native_volume.close()" ] }, { "cell_type": "markdown", "id": "744b4d0e", "metadata": {}, "source": [ "## 4. Include local boundary displacement and rigid moments\n", "\n", "With $S$ extracting the stress normal moments and $N$ mapping skeletal traction into the fine normal basis, the four local blocks are\n", "\n", "$$\n", "\\begin{pmatrix}\n", "M_\\mathcal A&D^T&W^T&-S\\\\\n", "D&0&0&0\\\\\n", "W&0&0&0\\\\\n", "-S^T&0&0&0\n", "\\end{pmatrix}\n", "\\begin{pmatrix}\\sigma\\\\u\\\\r\\\\\\widehat u\\end{pmatrix}\n", "+\\begin{pmatrix}0\\\\0\\\\0\\\\-N\\end{pmatrix}\\lambda_H\n", "=\\begin{pmatrix}0\\\\-F\\\\0\\\\0\\end{pmatrix}.\n", "$$\n", "\n", "The boundary selector records the declared BDM normal-moment ordering; `traction_mapping` owns geometric orientation and trace partitions. Its blocks already use global coordinates, so `coordinates=\"global\"` prevents applying incidence signs twice.\n", "\n", "The kernel comprises translations and the coupled rigid displacement/rotation/boundary trace. Declared displacement integrals fix their physical coordinates. The linear boundary rigid trace has a constant and first Legendre coefficient; higher coefficients vanish. These are definitions of the physical modes, not a numerical nullspace search." ] }, { "cell_type": "code", "execution_count": null, "id": "5541fbf9", "metadata": {}, "outputs": [], "source": [ "def local_equations(local):\n", " \"\"\"Declare mixed volume/normal-traction equations and moment-fixed physical rigid modes.\"\"\"\n", " fine = local.mesh\n", " mass, divergence, asymmetry, force = mixed_elasticity_operators(\n", " fine, lam, mu, source, 10, family\n", " )\n", " ns, nu, nr = mass.shape[0], divergence.shape[0], asymmetry.shape[0]\n", " normal_size = 2 * (family.degree + 1)\n", " nb = normal_size * len(fine.boundary_faces)\n", " normal_rows = (normal_size * fine.boundary_faces[:, None] + np.arange(normal_size)).ravel()\n", " S = sparse.coo_matrix((np.ones(nb), (normal_rows, np.arange(nb))), shape=(ns, nb)).tocsc()\n", " A = sparse.bmat(\n", " [\n", " [mass, divergence.T, asymmetry.T, -S],\n", " [divergence, None, None, None],\n", " [asymmetry, None, None, None],\n", " [-S.T, None, None, None],\n", " ],\n", " format=\"csc\",\n", " )\n", " L = np.r_[np.zeros(ns), -force, np.zeros(nr + nb)]\n", " N = traction_mapping(macro, local.cell, fine, skeleton, family)\n", " B = np.zeros((len(L), N.shape[1]))\n", " B[-nb:] = -N\n", " center = macro.points[macro.cells[local.cell]].mean(axis=0)\n", " nodes = np.einsum(\"qi,tij->tqj\", multiindices(1), fine.points[fine.cells])\n", " Z = np.zeros((len(L), 3))\n", " Z[ns : ns + nu] = rigid_values(nodes, center).reshape(-1, 3)\n", " Z[ns + nu : ns + nu + nr, 2] = -1\n", " rigid_edges = rigid_values(fine.points[fine.faces[fine.boundary_faces]], center)\n", " edge_coefficients = np.zeros((len(rigid_edges), family.degree + 1, 2, 3))\n", " edge_coefficients[:, 0] = rigid_edges.mean(axis=1)\n", " edge_coefficients[:, 1] = (rigid_edges[:, 1] - rigid_edges[:, 0]) / 2\n", " Z[ns + nu + nr :] = edge_coefficients.reshape(-1, 3)\n", " M = np.zeros_like(Z)\n", " M[ns : ns + nu] = displacement_rigid_moments(fine, center, family)\n", " plain, weighted, scale = stress_trace_moments(fine, ns, lam, mu, 10, family)\n", " pressure_weights = np.r_[weighted / scale, np.zeros(nu + nr + nb)]\n", " selector = sparse.eye(len(L), format=\"csr\")\n", " fields = (\n", " hdiv_field(\"stress\", fine, family, components=2, reconstruction=selector[:ns]),\n", " nodal_field(\n", " \"displacement\",\n", " fine,\n", " 1,\n", " components=2,\n", " discontinuous=True,\n", " reconstruction=selector[ns : ns + nu],\n", " ),\n", " nodal_field(\n", " \"rotation\", fine, 1, discontinuous=True, reconstruction=selector[ns + nu : ns + nu + nr]\n", " ),\n", " )\n", " equations = local.equations(\n", " a=A,\n", " L=L,\n", " b=B,\n", " c=-B.T,\n", " kernel=Z,\n", " moments=M,\n", " coordinates=\"global\",\n", " metadata=(pressure_weights, ns, nu, nr),\n", " )\n", " return replace(equations, field_data=fields)" ] }, { "cell_type": "markdown", "id": "c44e37a1", "metadata": {}, "source": [ "## 5. Declare the global boundary equation, assemble and solve\n", "\n", "The global equation enforces continuity of the weak boundary displacement and the prescribed exterior moments. The sign follows the mixed local equations above; it differs from a primal convention because the boundary auxiliary appears with $-S$.\n", "\n", "Finite $\\lambda$ also implies\n", "\n", "$$\n", "\\int_\\Omega\\frac{\\operatorname{tr}\\sigma}{2(\\mu+\\lambda)}\n", "=\\int_{\\partial\\Omega}g\\cdot n=0.\n", "$$\n", "\n", "The normalized weighted-stress row below states this identity. At incompressibility the trace integral becomes a pressure gauge; pure traction instead needs displacement rigid gauges. `MixedElasticitySolution` only interprets the executed blocks and supplies physical diagnostics." ] }, { "cell_type": "code", "execution_count": null, "id": "e86e113f", "metadata": {}, "outputs": [], "source": [ "def global_equation(context):\n", " \"\"\"Enforce assembled boundary displacement moments with the mixed-method sign.\"\"\"\n", " boundary, _ = context.boundary_data(displacement, order=10)\n", " return Equation(0, context.trace_load(boundary))\n", "\n", "\n", "with threadpool_limits(1):\n", " problem = bind_problem(\n", " hierarchy, interface, local_equations, global_equation=global_equation, retained=3\n", " )\n", " system = assemble(problem)\n", " # The chosen solenoidal data have zero integrated normal boundary displacement.\n", " constraints = [system.mean_constraint([m[0] for m in system.local_metadata], 0.0)]\n", " coefficients = solve(system, constraints=constraints)\n", " sf, uf, rf = (coefficients.field(name) for name in (\"stress\", \"displacement\", \"rotation\"))\n", " view = MixedElasticitySolution(\n", " skeleton,\n", " tuple(f.mesh for f in uf),\n", " tuple(f.coefficients.reshape(-1, 2) for f in sf),\n", " tuple(f.coefficients.reshape(len(f.mesh.cells), 3, 2) for f in uf),\n", " tuple(f.coefficients.reshape(len(f.mesh.cells), 3) for f in rf),\n", " coefficients,\n", " source,\n", " 10,\n", " family,\n", " )" ] }, { "cell_type": "markdown", "id": "2ba38683", "metadata": {}, "source": [ "## 6. Check physical errors, fine force balance and weak symmetry\n", "\n", "Do not symmetrize the computed stress before measuring its error. Symmetry is enforced in the rotation test space; its pointwise skew part may be nonzero. Fine-cell force moments, skeletal normal traction and full original equations are checked independently." ] }, { "cell_type": "code", "execution_count": null, "id": "65ab87c5", "metadata": {}, "outputs": [], "source": [ "from examples.core_elasticity_field_archive import ProductionObservation\n", "from examples.minimal_flow_originals import original_diagnostics\n", "\n", "errors = {\n", " \"displacement_l2\": view.l2_error(displacement, 12),\n", " \"stress_l2\": view.stress_l2_error(lambda x: smooth_fields(x)[1], 12),\n", " \"rotation_l2\": view.rotation_l2_error(lambda x: smooth_fields(x)[3], 12),\n", " \"stress_divergence_l2\": view.divergence_l2_error(lambda x: -source(x), 12),\n", "}\n", "invariants = {\n", " \"fine_force\": max(float(np.max(abs(r))) for r in view.fine_force_residuals()),\n", " \"weak_symmetry\": max(float(np.max(abs(r))) for r in view.weak_symmetry_residuals()),\n", " \"normal_traction\": max(float(np.max(abs(r))) for r in view.normal_traction_residuals()),\n", "}\n", "print(errors)\n", "print(invariants)\n", "assert max(invariants.values()) < 1e-10\n", "\n", "observation = ProductionObservation(\n", " system=system,\n", " applied_boundary=-system.global_load[: system.trace_size],\n", " constraints=constraints,\n", ")\n", "blocks = [\n", " {\n", " \"stress_constitutive\": ns,\n", " \"force_balance\": nu,\n", " \"weak_symmetry\": nr,\n", " \"normal_traction_compatibility\": len(field) - ns - nu - nr,\n", " }\n", " for field, (_, ns, nu, nr) in zip(coefficients.fields, system.local_metadata, strict=True)\n", "]\n", "physical = original_diagnostics(view, observation, blocks)\n", "print(\"Original physical residual:\", physical[\"full_uncondensed_relative_to_physical_rhs\"])\n", "assert physical[\"accepted\"]" ] }, { "cell_type": "markdown", "id": "659ec339", "metadata": {}, "source": [ "## 7. Build and refine the independent classical baseline\n", "\n", "The classical Taylor–Hood displacement–pressure system uses the same material, force and nonhomogeneous boundary values on global $32\\times32$ and $64\\times64$ meshes. These meshes are independent of the coarse $8\\times8$ macro partition. Analytical errors qualify the refined baseline before its field is plotted." ] }, { "cell_type": "code", "execution_count": null, "id": "b82c3e59", "metadata": {}, "outputs": [], "source": [ "from dolfinx import fem\n", "\n", "from pymhm.backends.spaces import bind_space, coefficient_map\n", "from pymhm.postprocessing.fields import DiscreteField, FieldDefinition\n", "\n", "\n", "def conforming_reference(resolution):\n", " \"\"\"Solve the notebook's unstabilized Taylor–Hood form on a fine global mesh.\"\"\"\n", " mesh = TriangleMesh.unit_square(resolution)\n", " element = basix.ufl.mixed_element(\n", " [\n", " basix.ufl.element(\n", " \"Lagrange\",\n", " \"triangle\",\n", " 2,\n", " shape=(2,),\n", " lagrange_variant=basix.LagrangeVariant.equispaced,\n", " ),\n", " basix.ufl.element(\n", " \"Lagrange\", \"triangle\", 1, lagrange_variant=basix.LagrangeVariant.equispaced\n", " ),\n", " ]\n", " )\n", " native = bind_space(mesh, element)\n", " try:\n", " u, p = ufl.TrialFunctions(native.space)\n", " v, q = ufl.TestFunctions(native.space)\n", " domain = native.space.mesh\n", " ue, pe, force = exact_ufl(domain)\n", " dx = ufl.Measure(\"dx\", domain=domain, metadata={\"quadrature_degree\": 12})\n", " a = (\n", " 2 * mu * ufl.inner(ufl.sym(ufl.grad(u)), ufl.sym(ufl.grad(v)))\n", " - p * ufl.div(v)\n", " - q * ufl.div(u)\n", " - p * q / lam\n", " ) * dx\n", " L = ufl.inner(force, v) * dx\n", " _, nodes = nodal_space(mesh, 2)\n", " boundary_nodes = np.flatnonzero(np.any(np.isclose(nodes, 0) | np.isclose(nodes, 1), axis=1))\n", " portable = (2 * boundary_nodes[:, None] + np.arange(2)).ravel()\n", " native_ids = coefficient_map(native, component=0)[portable]\n", " fixed = dict(zip(native_ids, exact_u(nodes[boundary_nodes]).ravel(), strict=True))\n", " reference_problem = MultiscaleProblem.from_global(\n", " Equation(a, L),\n", " native.size,\n", " fixed=fixed,\n", " )\n", " # Finite lambda determines this zero integral by the displacement data.\n", " pressure_integral = compile_form(q * dx)\n", " result = solve(assemble(reference_problem), constraints=[(pressure_integral, 0.0)])\n", " numerical = fem.Function(native.space)\n", " numerical.x.array[:] = result.trace\n", " uh, ph = ufl.split(numerical)\n", " sigma_h = 2 * mu * ufl.sym(ufl.grad(uh)) - ph * ufl.Identity(2)\n", " sigma_e = 2 * mu * ufl.sym(ufl.grad(ue)) - pe * ufl.Identity(2)\n", " differences = {\n", " \"displacement_l2\": uh - ue,\n", " \"pressure_l2\": ph - pe,\n", " \"stress_l2\": sigma_h - sigma_e,\n", " }\n", " errors = {\n", " name: float(np.sqrt(fem.assemble_scalar(fem.form(ufl.inner(delta, delta) * dx))))\n", " for name, delta in differences.items()\n", " }\n", " fields = {\n", " name: DiscreteField(\n", " FieldDefinition(name, descriptor=native.descriptor(component=i)), result.trace\n", " )\n", " for i, name in enumerate((\"displacement\", \"pressure\"))\n", " }\n", " return fields, {\n", " \"resolution\": resolution,\n", " \"triangles\": len(mesh.cells),\n", " \"unknowns\": native.size,\n", " \"errors\": errors,\n", " \"residual\": result.residual,\n", " }\n", " finally:\n", " native.close()\n", "\n", "\n", "with threadpool_limits(1):\n", " reference32, control32 = conforming_reference(32)\n", " reference64, control64 = conforming_reference(64)\n", "print(control32)\n", "print(control64)\n", "assert all(\n", " control64[\"errors\"][name] < control32[\"errors\"][name]\n", " for name in (\"displacement_l2\", \"stress_l2\")\n", ")" ] }, { "cell_type": "markdown", "id": "964e7ade", "metadata": {}, "source": [ "## 8. Plot stress and displacement with the actual macro mesh\n", "\n", "Independent fine-cell samples preserve the two traces at reconstructed interfaces; the display helper never averages them." ] }, { "cell_type": "code", "execution_count": null, "id": "7c2a72ee", "metadata": {}, "outputs": [], "source": [ "from IPython.display import SVG, display\n", "\n", "from examples.introduction.vector import plot_field_panels\n", "from examples.plot_mixed_elasticity import spatial_arrays\n", "\n", "triangulation, x, numerical = spatial_arrays(view)\n", "triangles = triangulation.triangles\n", "ue, se, _, re = smooth_fields(x)\n", "baseline = reference64[\"displacement\"].evaluate(x)\n", "panels = {\n", " \"Exact displacement\\ncomponent 1\": (x, triangles, ue[:, 0]),\n", " \"Mixed MHM displacement\\ncomponent 1\": (x, triangles, numerical[:, 0]),\n", " \"Fine Taylor–Hood displacement\\ncomponent 1\": (x, triangles, baseline[:, 0]),\n", " \"Exact Cauchy\\nstress xx\": (x, triangles, se[:, 0, 0]),\n", " \"Mixed MHM Cauchy\\nstress xx\": (x, triangles, numerical[:, 1]),\n", "}\n", "# Numerical x component and weak rotation are sampled independently in each fine cell.\n", "panels[\"Mixed displacement\\ncomponent error\"] = (x, triangles, numerical[:, 0] - ue[:, 0])\n", "figure = plot_field_panels(macro, panels, figsize=(13, 8))\n", "output = ROOT / \"build\" / \"tutorial-method-workflows\"\n", "output.mkdir(parents=True, exist_ok=True)\n", "figure.savefig(output / \"mixed_stress-elasticity-fields.png\", dpi=180, bbox_inches=\"tight\")\n", "figure.savefig(output / \"mixed_stress-elasticity-fields.svg\", bbox_inches=\"tight\")\n", "display(SVG(filename=str(output / \"mixed_stress-elasticity-fields.svg\")))\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "6178899f", "metadata": {}, "source": [ "## 9. Reach the smooth asymptotic regime\n", "\n", "The archived refinement series uses this same smooth field, BDM2/DG-P1/DG-P1 local spaces, refinement two, interior vector $P_1$ tractions and full exterior fine-edge $P_2$ tractions. It increases the macro resolution through $1,2,4,8,16,32$. The current explicit equations above run its $n=8$ level; reading the remaining attributed records is not another solve.\n", "\n", "The compatible smooth mixed stress/rotation estimates give order two. Displacement is shown separately. Rough coefficients, singular domains, different trace partitions or incomplete rigid spaces require their own regularity and stability assessment." ] }, { "cell_type": "code", "execution_count": null, "id": "d44ce6a8", "metadata": {}, "outputs": [], "source": [ "import json\n", "from examples.tutorial_convergence import refinement_series,plot_method_series,asymptotic_summary\n", "record_path=ROOT/'examples/results/tutorial-methods/mixed-elasticity-asymptotic.json'\n", "record=json.loads(record_path.read_text())\n", "rows=[dict(row,H=1/row['macro_resolution']) for row in record['rows']]\n", "current_row=next(row for row in rows if row['macro_resolution']==n)\n", "for name,value in errors.items():\n", " print(name, 'current/archive:', value, current_row[name])\n", " np.testing.assert_allclose(value,current_row[name],rtol=1e-8,atol=1e-12)\n", "series=refinement_series(\n", " record_path,rows,'H',\n", " {'displacement L2':'displacement_l2','stress L2':'stress_l2','rotation L2':'rotation_l2'},\n", " {'stress L2':2.0,'rotation L2':2.0},method='mixed-elasticity',\n", " spaces='BDM2 stress, DG P1 displacement and rotation, interior P1 and full exterior fine P2 traces',\n", " rate_provenance='Devloo et al. (2021), theorem 5.2, equation (5.5): conditional smooth two-dimensional compatible mixed estimate; stress and weak rotation target order two; displacement measured separately',\n", " refinement='macro grid spacing',root=ROOT,\n", ")\n", "print(asymptotic_summary(series))\n", "figure=plot_method_series(series)\n", "convergence_path=output/\"mixed-stress-convergence.svg\"\n", "figure.savefig(convergence_path,bbox_inches=\"tight\")\n", "display(SVG(filename=str(convergence_path)))\n", "plt.close(figure)\n" ] }, { "cell_type": "markdown", "id": "57a32935", "metadata": {}, "source": [ "This smooth verification uses independently stated analytical data and is not a claim of matching every original heterogeneous table.\n", "\n", "## References\n", "\n", "- [Devloo, Farias, Gomes, Pereira, dos Santos and Valentin (2021), *New H(div)-conforming multiscale hybrid-mixed methods for the elasticity problem on polygonal meshes*](https://doi.org/10.1051/m2an/2021013).\n", "- [Arnold, Falk and Winther (2007), *Mixed finite element methods for linear elasticity with weakly imposed symmetry*](https://doi.org/10.1090/S0025-5718-07-01998-9).\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python" }, "pymhm": { "problem": "elasticity", "title": "Mixed stress MHM through explicit local and global equations", "workflow": "mesh-first explicit local/global equations" } }, "nbformat": 4, "nbformat_minor": 5 }