test_LinearElasticity_fenics_2023.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "40873d58-713a-467b-b4ec-8c1807e41431",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyvista\n",
    "from dolfinx import mesh, fem, plot, io, default_scalar_type\n",
    "from dolfinx.fem.petsc import LinearProblem\n",
    "from mpi4py import MPI\n",
    "import ufl\n",
    "import numpy as np\n",
    "import gmsh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9b5b56ee-9638-43d5-afa2-583d05cdbc5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#from fenics import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "593b9f69-5396-4f8b-bee1-3fded4b8528d",
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 1\n",
    "W = 0.2\n",
    "mu = 1\n",
    "rho = 1\n",
    "delta = W / L\n",
    "gamma = 0.4 * delta**2\n",
    "beta = 1.25\n",
    "lambda_ = beta\n",
    "g = gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c03fea11-7dd9-45c1-a082-1359f840a66c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],\n",
    "#                         [30, 9, 9], cell_type=mesh.CellType.hexahedron)\n",
    "#V = fem.VectorFunctionSpace(domain, (\"Lagrange\", 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ac243274-5ea4-454a-900f-38712281c1c7",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:py.warnings:/tmp/ipykernel_4120907/210129278.py:14: DeprecationWarning: This method is deprecated. Use FunctionSpace with an element shape argument instead\n",
      "  V = fem.VectorFunctionSpace(domain, (\"Lagrange\", 1))\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info    : Meshing 1D...\n",
      "Info    : [  0%] Meshing curve 1 (Line)\n",
      "Info    : [ 10%] Meshing curve 2 (Line)\n",
      "Info    : [ 20%] Meshing curve 3 (Line)\n",
      "Info    : [ 30%] Meshing curve 4 (Line)\n",
      "Info    : [ 40%] Meshing curve 5 (Line)\n",
      "Info    : [ 50%] Meshing curve 6 (Line)\n",
      "Info    : [ 50%] Meshing curve 7 (Line)\n",
      "Info    : [ 60%] Meshing curve 8 (Line)\n",
      "Info    : [ 70%] Meshing curve 9 (Line)\n",
      "Info    : [ 80%] Meshing curve 10 (Line)\n",
      "Info    : [ 90%] Meshing curve 11 (Line)\n",
      "Info    : [100%] Meshing curve 12 (Line)\n",
      "Info    : Done meshing 1D (Wall 0.000572084s, CPU 0s)\n",
      "Info    : Meshing 2D...\n",
      "Info    : [  0%] Meshing surface 1 (Plane, Frontal-Delaunay)\n",
      "Info    : [ 20%] Meshing surface 2 (Plane, Frontal-Delaunay)\n",
      "Info    : [ 40%] Meshing surface 3 (Plane, Frontal-Delaunay)\n",
      "Info    : [ 50%] Meshing surface 4 (Plane, Frontal-Delaunay)\n",
      "Info    : [ 70%] Meshing surface 5 (Plane, Frontal-Delaunay)\n",
      "Info    : [ 90%] Meshing surface 6 (Plane, Frontal-Delaunay)\n",
      "Info    : Done meshing 2D (Wall 0.0229622s, CPU 0.020074s)\n",
      "Info    : Meshing 3D...\n",
      "Info    : 3D Meshing 1 volume with 1 connected component\n",
      "Info    : Tetrahedrizing 1220 nodes...\n",
      "Info    : Done tetrahedrizing 1228 nodes (Wall 0.0104322s, CPU 0.007015s)\n",
      "Info    : Reconstructing mesh...\n",
      "Info    :  - Creating surface mesh\n",
      "Info    :  - Identifying boundary edges\n",
      "Info    :  - Recovering boundary\n",
      "Info    : Done reconstructing mesh (Wall 0.0226027s, CPU 0.019179s)\n",
      "Info    : Found volume 1\n",
      "Info    : It. 0 - 0 nodes created - worst tet radius 3.43958 (nodes removed 0 0)\n",
      "Info    : It. 500 - 500 nodes created - worst tet radius 1.105 (nodes removed 0 0)\n",
      "Info    : 3D refinement terminated (1875 nodes total):\n",
      "Info    :  - 0 Delaunay cavities modified for star shapeness\n",
      "Info    :  - 0 nodes could not be inserted\n",
      "Info    :  - 7624 tetrahedra created in 0.0219555 sec. (347247 tets/s)\n",
      "Info    : 0 node relocations\n",
      "Info    : Done meshing 3D (Wall 0.0694416s, CPU 0.065978s)\n",
      "Info    : Optimizing mesh...\n",
      "Info    : Optimizing volume 1\n",
      "Info    : Optimization starts (volume = 0.04) with worst = 0.0123988 / average = 0.782072:\n",
      "Info    : 0.00 < quality < 0.10 :        17 elements\n",
      "Info    : 0.10 < quality < 0.20 :        40 elements\n",
      "Info    : 0.20 < quality < 0.30 :        75 elements\n",
      "Info    : 0.30 < quality < 0.40 :        95 elements\n",
      "Info    : 0.40 < quality < 0.50 :       165 elements\n",
      "Info    : 0.50 < quality < 0.60 :       296 elements\n",
      "Info    : 0.60 < quality < 0.70 :       935 elements\n",
      "Info    : 0.70 < quality < 0.80 :      1876 elements\n",
      "Info    : 0.80 < quality < 0.90 :      2803 elements\n",
      "Info    : 0.90 < quality < 1.00 :      1322 elements\n",
      "Info    : 132 edge swaps, 5 node relocations (volume = 0.04): worst = 0.216892 / average = 0.793341 (Wall 0.00182112s, CPU 0.001618s)\n",
      "Info    : 134 edge swaps, 5 node relocations (volume = 0.04): worst = 0.258376 / average = 0.793498 (Wall 0.00230595s, CPU 0.001948s)\n",
      "Info    : No ill-shaped tets in the mesh :-)\n",
      "Info    : 0.00 < quality < 0.10 :         0 elements\n",
      "Info    : 0.10 < quality < 0.20 :         0 elements\n",
      "Info    : 0.20 < quality < 0.30 :         2 elements\n",
      "Info    : 0.30 < quality < 0.40 :        96 elements\n",
      "Info    : 0.40 < quality < 0.50 :       139 elements\n",
      "Info    : 0.50 < quality < 0.60 :       289 elements\n",
      "Info    : 0.60 < quality < 0.70 :       950 elements\n",
      "Info    : 0.70 < quality < 0.80 :      1880 elements\n",
      "Info    : 0.80 < quality < 0.90 :      2816 elements\n",
      "Info    : 0.90 < quality < 1.00 :      1331 elements\n",
      "Info    : Done optimizing mesh (Wall 0.00736043s, CPU 0.006866s)\n",
      "Info    : 1875 nodes 10139 elements\n"
     ]
    }
   ],
   "source": [
    "from dolfinx.io import gmshio\n",
    "gmsh.initialize()\n",
    "box=gmsh.model.occ.addBox(0,0,0,L,W,W)\n",
    "gmsh.model.occ.synchronize()\n",
    "gdim = 3\n",
    "gmsh.model.addPhysicalGroup(gdim, [box], 1)\n",
    "gmsh.option.setNumber(\"Mesh.CharacteristicLengthMin\", 0.03)\n",
    "gmsh.option.setNumber(\"Mesh.CharacteristicLengthMax\", 0.03)\n",
    "#gmsh.option.setNumber(\"Mesh.SubdivisionAlgorithm\", 1)\n",
    "gmsh.model.mesh.generate(gdim)\n",
    "gmsh_model_rank = 0\n",
    "mesh_comm = MPI.COMM_WORLD\n",
    "domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)\n",
    "V = fem.VectorFunctionSpace(domain, (\"Lagrange\", 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1b7ef132-16b2-4418-9106-3f65b147e0fd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6077dc802aa84b9da243efe09be1eb46",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Widget(value='<iframe src=\"http://localhost:44587/index.html?ui=P_0x79cfd2d10ce0_0&reconnect=auto\" class=\"pyvi…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = pyvista.Plotter()\n",
    "topology, cell_types, geometry = plot.vtk_mesh(V)\n",
    "grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)\n",
    "\n",
    "p.add_mesh(grid, show_edges=True)\n",
    "p.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "55a22075-2f0b-48d1-be19-d0df063c0dc9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "FunctionSpace(Mesh(blocked element (Basix element (P, tetrahedron, 1, equispaced, unset, False), (3,)), 0), blocked element (Basix element (P, tetrahedron, 1, gll_warped, unset, False), (3,)))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9ab40e4b-b83e-4e4f-b263-732f5775995b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#def clamped_boundary(x):\n",
    "#    return np.isclose(x[0], 0)\n",
    "#\n",
    "#\n",
    "#fdim = domain.topology.dim - 1\n",
    "#boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)\n",
    "#\n",
    "#u_D = np.array([0, 0, 0], dtype=default_scalar_type)\n",
    "#bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5a2bee50-a47a-4a0b-96fd-a50c39d3b84b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def clamped_boundary(x):\n",
    "    return np.isclose(x[0], 0)\n",
    "\n",
    "def clamped_boundary1(x):\n",
    "    return np.isclose(x[0], 1)\n",
    "\n",
    "boundary_dofs = fem.locate_dofs_geometrical(V, clamped_boundary)\n",
    "boundary_dofs1 = fem.locate_dofs_geometrical(V, clamped_boundary1)\n",
    "u_D = np.array([0, 0, 0], dtype=default_scalar_type)\n",
    "#bc = fem.dirichletbc(u_D, boundary_dofs, V)\n",
    "bc_0 = fem.dirichletbc(u_D, boundary_dofs, V)\n",
    "bc_1 = fem.dirichletbc(u_D, boundary_dofs1, V)\n",
    "#bc = [bc_0, bc_1]\n",
    "bc = [bc_0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4b0977b1-678b-441e-ac77-d3cb9f0ec0fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.02166747, 0.17854794, 0.        ],\n",
       "       [0.        , 0.2       , 0.        ],\n",
       "       [0.        , 0.17142857, 0.        ],\n",
       "       ...,\n",
       "       [0.97899825, 0.02076676, 0.2       ],\n",
       "       [1.        , 0.02857143, 0.2       ],\n",
       "       [1.        , 0.        , 0.2       ]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "domain.geometry.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "3596941a-ee07-465d-9c6a-3a2a25b33fb8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1875"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(domain.geometry.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "289ec5f5-cca6-4248-a381-afb9fbf8906a",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = ufl.Measure(\"ds\", domain=domain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "3fde53bf-4abf-46a3-95cd-6e121cbe31c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def epsilon(u):\n",
    "    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)\n",
    "\n",
    "\n",
    "def sigma(u):\n",
    "    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)\n",
    "\n",
    "T = fem.Constant(domain, default_scalar_type((0, 0, 0)))\n",
    "u = ufl.TrialFunction(V)\n",
    "v = ufl.TestFunction(V)\n",
    "f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))\n",
    "a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx\n",
    "L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "41277096-d572-440c-ac35-cac2ce37ada7",
   "metadata": {},
   "outputs": [],
   "source": [
    "problem = LinearProblem(a, L, bcs=bc, petsc_options={\"ksp_type\": \"preonly\", \"pc_type\": \"lu\"})\n",
    "uh = problem.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "e1bde932-e786-464e-94d7-7f020bbff05d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.75752411e-03,  5.29270075e-04, -1.03437145e-03, ...,\n",
       "        2.89496660e-02,  2.19170604e-04, -2.24796679e-01])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "uh.x.array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "3793a384-09d6-495b-a863-a33064b22dc8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8cc506e424f44f9a9628f683277f9241",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Widget(value='<iframe src=\"http://localhost:44587/index.html?ui=P_0x79cf13f21e50_1&reconnect=auto\" class=\"pyvi…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#pyvista.start_xvfb()\n",
    "\n",
    "# Create plotter and pyvista grid\n",
    "p = pyvista.Plotter()\n",
    "topology, cell_types, geometry = plot.vtk_mesh(V)\n",
    "grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)\n",
    "\n",
    "# Attach vector values to grid and warp grid by vector\n",
    "grid[\"u\"] = uh.x.array.reshape((geometry.shape[0], 3))\n",
    "#actor_0 = p.add_mesh(grid, style=\"wireframe\", color=\"k\")\n",
    "warped = grid.warp_by_vector(\"u\", factor=1.0)\n",
    "actor_1 = p.add_mesh(warped, show_edges=True)\n",
    "p.show_axes()\n",
    "if not pyvista.OFF_SCREEN:\n",
    "    p.show()\n",
    "else:\n",
    "    figure_as_array = p.screenshot(\"deflection.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b37d047d-bb68-439c-96f8-0614eabed69c",
   "metadata": {},
   "outputs": [],
   "source": [
    "s = sigma(uh) - 1. / 3 * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))\n",
    "von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))\n",
    "V_von_mises = fem.FunctionSpace(domain, (\"DG\", 0))\n",
    "stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())\n",
    "stresses = fem.Function(V_von_mises)\n",
    "stresses.interpolate(stress_expr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "9f40aaf0-d282-400a-b3cf-7ec948e47d6e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f38187498fda418e82a21b8c4fd3eec3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Widget(value='<iframe src=\"http://localhost:44587/index.html?ui=P_0x79cf230ab5c0_2&reconnect=auto\" class=\"pyvi…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "warped.cell_data[\"VonMises\"] = stresses.vector.array\n",
    "warped.set_active_scalars(\"VonMises\")\n",
    "p = pyvista.Plotter()\n",
    "p.add_mesh(warped, show_edges=True)\n",
    "p.show_axes()\n",
    "if not pyvista.OFF_SCREEN:\n",
    "    p.show()\n",
    "else:\n",
    "    stress_figure = p.screenshot(f\"stresses.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f63354f-91b1-4998-b17d-078ec27541ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "uh.vector.array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3f72c4d-5d2c-45bc-92d4-c927cd57e8a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(uh.vector.array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51979b4b-d8cb-49cd-b153-229af7387c53",
   "metadata": {},
   "outputs": [],
   "source": [
    "V_von_mises.mesh.geometry.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a205fc6f-b6a9-4dd8-81f6-c9f3c0a3bbd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(V_von_mises.mesh.geometry.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "472aa7b5-1eb6-4e11-8545-3a23a0cc4a8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "stresses.vector.array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecff6004-59e3-4fe7-b552-3a6035618b89",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(stresses.vector.array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a90d287-8917-4ab0-ab83-b311a74a9355",
   "metadata": {},
   "outputs": [],
   "source": [
    "uh.x.array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3f4e5fe-63f5-4da1-bb3a-1c780924aecb",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(uh.x.array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93dcc2ad-8a1d-4794-969d-80dd336bdf30",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}