-# -# -# .. _LinearElasticity2D: -# -# ========================= -# 2D linear elasticity -# ========================= -# -# -# Introduction -# ------------ -# -# In this first numerical tour, we will show how to compute a small strain solution for -# a 2D isotropic linear elastic medium, either in plane stress or in plane strain, -# in a tradtional displacement-based finite element formulation. The corresponding -# file can be obtained from :download:2D_elasticity.py. -# -# .. seealso:: -# -# Extension to 3D is straightforward and an example can be found in the :ref:ModalAnalysis example. -# -# We consider here the case of a cantilever beam modeled as a 2D medium of dimensions -# :math:L\times H. Geometrical parameters and mesh density are first defined -# and the rectangular domain is generated using the RectangleMesh function. -# We also choose a criss-crossed structured mesh:: - -from __future__ import print_function -from fenics import * - -L = 25. -H = 1. -Nx = 250 -Ny = 10 -mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, "crossed") - - -# Constitutive relation -# --------------------- -# -# We now define the material parameters which are here given in terms of a Young's -# modulus :math:E and a Poisson coefficient :math:\nu. In the following, we will -# need to define the constitutive relation between the stress tensor :math:\boldsymbol{\sigma} -# and the strain tensor :math:\boldsymbol{\varepsilon}. Let us recall -# that the general expression of the linear elastic isotropic constitutive relation -# for a 3D medium is given by: -# -# .. math:: -# \boldsymbol{\sigma} = \lambda \text{tr}(\boldsymbol{\varepsilon})\mathbf{1} + 2\mu\boldsymbol{\varepsilon} -# :label: constitutive_3D -# -# for a natural (no prestress) initial state where the Lamé coefficients are given by: -# -# .. math:: -# \lambda = \dfrac{E\nu}{(1+\nu)(1-2\nu)}, \quad \mu = \dfrac{E}{2(1+\nu)} -# :label: Lame_coeff -# -# In this demo, we consider a 2D model either in plane strain or in plane stress conditions. -# Irrespective of this choice, we will work only with a 2D displacement vector :math:\boldsymbol{u}=(u_x,u_y) -# and will subsequently define the strain operator eps as follows:: - -def eps(v): - return sym(grad(v)) - -# which computes the 2x2 plane components of the symmetrized gradient tensor of -# any 2D vectorial field. In the plane strain case, the full 3D strain tensor is defined as follows: -# -# .. math:: -# \boldsymbol{\varepsilon} = \begin{bmatrix} \varepsilon_{xx} & \varepsilon_{xy} & 0\\ -# \varepsilon_{xy} & \varepsilon_{yy} & 0 \\ 0 & 0 & 0\end{bmatrix} -# -# so that the 2x2 plane part of the stress tensor is defined in the same way as for the 3D case -# (the out-of-plane stress component being given by :math:\sigma_{zz}=\lambda(\varepsilon_{xx}+\varepsilon_{yy}). -# -# In the plane stress case, an out-of-plane strain component :math:\varepsilon_{zz} -# must be considered so that :math:\sigma_{zz}=0. Using this condition in the -# 3D constitutive relation, one has :math:\varepsilon_{zz}=-\dfrac{\lambda}{\lambda+2\mu}(\varepsilon_{xx}+\varepsilon_{yy}). -# Injecting into :eq:constitutive_3D, we have for the 2D plane stress relation: -# -# .. math:: -# \boldsymbol{\sigma} = \lambda^* \text{tr}(\boldsymbol{\varepsilon})\mathbf{1} + 2\mu\boldsymbol{\varepsilon} -# -# where :math:\boldsymbol{\sigma}, \boldsymbol{\varepsilon}, \mathbf{1} are 2D tensors and with -# :math:\lambda^* = \dfrac{2\lambda\mu}{\lambda+2\mu}. Hence, the 2D constitutive relation -# is identical to the plane strain case by changing only the value of the Lamé coefficient :math:\lambda. -# We can then have:: - -E = Constant(1e5) -nu = Constant(0.3) -model = "plane_stress" - -mu = E/2/(1+nu) -lmbda = E*nu/(1+nu)/(1-2*nu) -if model == "plane_stress": - lmbda = 2*mu*lmbda/(lmbda+2*mu) - -def sigma(v): - return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v) - -# .. note:: -# Note that we used the variable name lmbda to avoid any confusion with the -# lambda functions of Python -# -# We also used an intrinsic formulation of the constitutive relation. Example of -# constitutive relation implemented with a matrix/vector engineering notation -# will be provided in the :ref:OrthotropicElasticity example. -# -# -# Variational formulation -# ----------------------- -# -# For this example, we consider a continuous polynomial interpolation of degree 2 -# and a uniformly distributed loading :math:\boldsymbol{f}=(0,-f) corresponding -# to the beam self-weight. The continuum mechanics variational formulation (obtained -# from the virtual work principle) is given by: -# -# .. math:: -# \text{Find } \boldsymbol{u}\in V \text{ s.t. } \int_{\Omega} -# \boldsymbol{\sigma}(\boldsymbol{u}):\boldsymbol{\varepsilon}(\boldsymbol{v}) d\Omega -# = \int_{\Omega} \boldsymbol{f}\cdot\boldsymbol{v} d\Omega \quad \forall\boldsymbol{v} \in V -# -# which translates into the following FEniCS code:: - -rho_g = 1e-3 -f = Constant((0,-rho_g)) - -V = VectorFunctionSpace(mesh, 'Lagrange', degree=2) -du = TrialFunction(V) -u_ = TestFunction(V) -a = inner(sigma(du), eps(u_))*dx -l = inner(f, u_)*dx - - -# Resolution -# ---------- -# -# Fixed displacements are imposed on the left part of the beam, the solve -# function is then called and solution is plotted by deforming the mesh:: - -def left(x, on_boundary): - return near(x[0],0.) - -bc = DirichletBC(V, Constant((0.,0.)), left) - -u = Function(V, name="Displacement") -solve(a == l, u, bc) - -plot(1e3*u, mode="displacement") - -# The (amplified) solution should look like this: -# -# .. image:: cantilever_deformed.png -# :scale: 15% -# -# -# Validation and post-processing -# ------------------------------ -# -# The maximal deflection is compared against the analytical solution from -# Euler-Bernoulli beam theory which is here :math:w_{beam} = \dfrac{qL^4}{8EI}:: - -print("Maximal deflection:", -u(L,H/2.)[1]) -print("Beam theory deflection:", float(3*rho_g*L**4/2/E/H**3)) - -# One finds :math:w_{FE} = 5.8638\text{e-3} against :math:w_{beam} = 5.8594\text{e-3} -# that is a 0.07% difference. -# -# -# The stress tensor must be projected on an appropriate function space in order to -# evaluate pointwise values or export it for Paraview vizualisation. Here we choose -# to describe it as a (2D) tensor and project it onto a piecewise constant function -# space:: - -Vsig = TensorFunctionSpace(mesh, "DG", degree=0) -sig = Function(Vsig, name="Stress") -sig.assign(project(sigma(u), Vsig)) -print("Stress at (0,H):", sig(0, H)) - -# Fields can be exported in a suitable format for vizualisation using Paraview. -# VTK-based extensions (.pvd,.vtu) are not suited for multiple fields and parallel -# writing/reading. Prefered output format is now .xdmf:: - -file_results = XDMFFile("elasticity_results.xdmf") -file_results.parameters["flush_output"] = True -file_results.parameters["functions_share_mesh"] = True -file_results.write(u, 0.) -file_results.write(sig, 0.) \ No newline at end of file diff --git a/doc/_build/html/_downloads/axisymmetric_elasticity.ipynb b/doc/_build/html/_downloads/axisymmetric_elasticity.ipynb deleted file mode 100644 index 45745e4a2c775951548079e3980c2874a7ac97cb..0000000000000000000000000000000000000000 --- a/doc/_build/html/_downloads/axisymmetric_elasticity.ipynb +++ /dev/null @@ -1,2616 +0,0 @@ -{ - "cells": [ - { - "cell_type": "raw", - "metadata": { - "raw_mimetype": "text/restructuredtext" - }, - "source": [ - ".. AxisymmetricElasticity:" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Axisymmetric formulation for elastic structures of revolution\n", - "\n", - "In this numerical tour, we will deal with axisymmetric problems of elastic solids. We will consider a solid of revolution around a fixed axis $(Oz)$, the loading, boundary conditions and material properties being also invariant with respect to a rotation along the symmetry axis. The solid cross-section in a plane $\\theta=\\text{cst}$ will be represented by a two-dimensional domain $\\omega$ for which the first spatial variable (x[0] in FEniCS) will represent the radial coordinate $r$ whereas the second spatial variable will denote the axial variable $z$.\n", - "\n", - "## Problem position\n", - "\n", - "We will investigate here the case of a hollow hemisphere of inner (resp. outer) radius $R_i$ (resp. $R_e$). Due to the revolution symmetry, the 2D cross-section corresponds to a quarter of a hollow cylinder." - ] - }, - { - "cell_type": "code", - "execution_count": 51, - "metadata": {}, - "outputs": [ - { - "data": { - "application/javascript": [ - "/* Put everything inside the global mpl namespace */\n", - "window.mpl = {};\n", - "\n", - "\n", - "mpl.get_websocket_type = function() {\n", - " if (typeof(WebSocket) !== 'undefined') {\n", - " return WebSocket;\n", - " } else if (typeof(MozWebSocket) !== 'undefined') {\n", - " return MozWebSocket;\n", - " } else {\n", - " alert('Your browser does not have WebSocket support.' +\n", - " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", - " 'Firefox 4 and 5 are also supported but you ' +\n", - " 'have to enable WebSockets in about:config.');\n", - " };\n", - "}\n", - "\n", - "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", - " this.id = figure_id;\n", - "\n", - " this.ws = websocket;\n", - "\n", - " this.supports_binary = (this.ws.binaryType != undefined);\n", - "\n", - " if (!this.supports_binary) {\n", - " var warnings = document.getElementById(\"mpl-warnings\");\n", - " if (warnings) {\n", - " warnings.style.display = 'block';\n", - " warnings.textContent = (\n", - " \"This browser does not support binary websocket messages. \" +\n", - " \"Performance may be slow.\");\n", - " }\n", - " }\n", - "\n", - " this.imageObj = new Image();\n", - "\n", - " this.context = undefined;\n", - " this.message = undefined;\n", - " this.canvas = undefined;\n", - " this.rubberband_canvas = undefined;\n", - " this.rubberband_context = undefined;\n", - " this.format_dropdown = undefined;\n", - "\n", - " this.image_mode = 'full';\n", - "\n", - " this.root = $(' ');\n", - " this._root_extra_style(this.root)\n", - " this.root.attr('style', 'display: inline-block');\n", - "\n", - "$(parent_element).append(this.root);\n", - "\n", - " this._init_header(this);\n", - " this._init_canvas(this);\n", - " this._init_toolbar(this);\n", - "\n", - " var fig = this;\n", - "\n", - " this.waiting = false;\n", - "\n", - " this.ws.onopen = function () {\n", - " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", - " fig.send_message(\"send_image_mode\", {});\n", - " if (mpl.ratio != 1) {\n", - " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", - " }\n", - " fig.send_message(\"refresh\", {});\n", - " }\n", - "\n", - " this.imageObj.onload = function() {\n", - " if (fig.image_mode == 'full') {\n", - " // Full images could contain transparency (where diff images\n", - " // almost always do), so we need to clear the canvas so that\n", - " // there is no ghosting.\n", - " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", - " }\n", - " fig.context.drawImage(fig.imageObj, 0, 0);\n", - " };\n", - "\n", - " this.imageObj.onunload = function() {\n", - " this.ws.close();\n", - " }\n", - "\n", - " this.ws.onmessage = this._make_on_message_function(this);\n", - "\n", - " this.ondownload = ondownload;\n", - "}\n", - "\n", - "mpl.figure.prototype._init_header = function() {\n", - " var titlebar = $(\n", - " ' ');\n", - " var titletext =$(\n", - " '
');\n", - " titlebar.append(titletext)\n", - " this.root.append(titlebar);\n", - " this.header = titletext[0];\n", - "}\n", - "\n", - "\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._init_canvas = function() {\n", - " var fig = this;\n", - "\n", - " var canvas_div = $(' ');\n", - "\n", - " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", - "\n", - " function canvas_keyboard_event(event) {\n", - " return fig.key_event(event, event['data']);\n", - " }\n", - "\n", - " canvas_div.keydown('key_press', canvas_keyboard_event);\n", - " canvas_div.keyup('key_release', canvas_keyboard_event);\n", - " this.canvas_div = canvas_div\n", - " this._canvas_extra_style(canvas_div)\n", - " this.root.append(canvas_div);\n", - "\n", - " var canvas =$('');\n", - " canvas.addClass('mpl-canvas');\n", - " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", - "\n", - " this.canvas = canvas[0];\n", - " this.context = canvas[0].getContext(\"2d\");\n", - "\n", - " var backingStore = this.context.backingStorePixelRatio ||\n", - "\tthis.context.webkitBackingStorePixelRatio ||\n", - "\tthis.context.mozBackingStorePixelRatio ||\n", - "\tthis.context.msBackingStorePixelRatio ||\n", - "\tthis.context.oBackingStorePixelRatio ||\n", - "\tthis.context.backingStorePixelRatio || 1;\n", - "\n", - " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", - "\n", - " var rubberband = $('');\n", - " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", - "\n", - " var pass_mouse_events = true;\n", - "\n", - " canvas_div.resizable({\n", - " start: function(event, ui) {\n", - " pass_mouse_events = false;\n", - " },\n", - " resize: function(event, ui) {\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " stop: function(event, ui) {\n", - " pass_mouse_events = true;\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " });\n", - "\n", - " function mouse_event_fn(event) {\n", - " if (pass_mouse_events)\n", - " return fig.mouse_event(event, event['data']);\n", - " }\n", - "\n", - " rubberband.mousedown('button_press', mouse_event_fn);\n", - " rubberband.mouseup('button_release', mouse_event_fn);\n", - " // Throttle sequential mouse events to 1 every 20ms.\n", - " rubberband.mousemove('motion_notify', mouse_event_fn);\n", - "\n", - " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", - " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", - "\n", - " canvas_div.on(\"wheel\", function (event) {\n", - " event = event.originalEvent;\n", - " event['data'] = 'scroll'\n", - " if (event.deltaY < 0) {\n", - " event.step = 1;\n", - " } else {\n", - " event.step = -1;\n", - " }\n", - " mouse_event_fn(event);\n", - " });\n", - "\n", - " canvas_div.append(canvas);\n", - " canvas_div.append(rubberband);\n", - "\n", - " this.rubberband = rubberband;\n", - " this.rubberband_canvas = rubberband[0];\n", - " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", - " this.rubberband_context.strokeStyle = \"#000000\";\n", - "\n", - " this._resize_canvas = function(width, height) {\n", - " // Keep the size of the canvas, canvas container, and rubber band\n", - " // canvas in synch.\n", - " canvas_div.css('width', width)\n", - " canvas_div.css('height', height)\n", - "\n", - " canvas.attr('width', width * mpl.ratio);\n", - " canvas.attr('height', height * mpl.ratio);\n", - " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", - "\n", - " rubberband.attr('width', width);\n", - " rubberband.attr('height', height);\n", - " }\n", - "\n", - " // Set the figure to an initial 600x600px, this will subsequently be updated\n", - " // upon first draw.\n", - " this._resize_canvas(600, 600);\n", - "\n", - " // Disable right mouse context menu.\n", - "$(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", - " return false;\n", - " });\n", - "\n", - " function set_focus () {\n", - " canvas.focus();\n", - " canvas_div.focus();\n", - " }\n", - "\n", - " window.setTimeout(set_focus, 100);\n", - "}\n", - "\n", - "mpl.figure.prototype._init_toolbar = function() {\n", - " var fig = this;\n", - "\n", - " var nav_element = $(' ')\n", - " nav_element.attr('style', 'width: 100%');\n", - " this.root.append(nav_element);\n", - "\n", - " // Define a callback function for later on.\n", - " function toolbar_event(event) {\n", - " return fig.toolbar_button_onclick(event['data']);\n", - " }\n", - " function toolbar_mouse_event(event) {\n", - " return fig.toolbar_button_onmouseover(event['data']);\n", - " }\n", - "\n", - " for(var toolbar_ind in mpl.toolbar_items) {\n", - " var name = mpl.toolbar_items[toolbar_ind][0];\n", - " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", - " var image = mpl.toolbar_items[toolbar_ind][2];\n", - " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", - "\n", - " if (!name) {\n", - " // put a spacer in here.\n", - " continue;\n", - " }\n", - " var button =$('');\n", - " button.click(method_name, toolbar_event);\n", - " button.mouseover(tooltip, toolbar_mouse_event);\n", - " nav_element.append(button);\n", - " }\n", - "\n", - " // Add the status bar.\n", - " var status_bar = $('');\n", - " nav_element.append(status_bar);\n", - " this.message = status_bar[0];\n", - "\n", - " // Add the close button to the window.\n", - " var buttongrp =$('
');\n", - " var button = $('');\n", - " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", - " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", - " buttongrp.append(button);\n", - " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", - " titlebar.prepend(buttongrp);\n", - "}\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(el){\n", - " var fig = this\n", - " el.on(\"remove\", function(){\n", - "\tfig.close_ws(fig, {});\n", - " });\n", - "}\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(el){\n", - " // this is important to make the div 'focusable\n", - " el.attr('tabindex', 0)\n", - " // reach out to IPython and tell the keyboard manager to turn it's self\n", - " // off when our div gets focus\n", - "\n", - " // location in version 3\n", - " if (IPython.notebook.keyboard_manager) {\n", - " IPython.notebook.keyboard_manager.register_events(el);\n", - " }\n", - " else {\n", - " // location in version 2\n", - " IPython.keyboard_manager.register_events(el);\n", - " }\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._key_event_extra = function(event, name) {\n", - " var manager = IPython.notebook.keyboard_manager;\n", - " if (!manager)\n", - " manager = IPython.keyboard_manager;\n", - "\n", - " // Check for shift+enter\n", - " if (event.shiftKey && event.which == 13) {\n", - " this.canvas_div.blur();\n", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\n", - " }\n", - "}\n", - "\n", - "mpl.figure.prototype.handle_save = function(fig, msg) {\n", - " fig.ondownload(fig, null);\n", - "}\n", - "\n", - "\n", - "mpl.find_output_cell = function(html_output) {\n", - " // Return the cell and output element which can be found *uniquely* in the notebook.\n", - " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", - " // IPython event is triggered only after the cells have been serialised, which for\n", - " // our purposes (turning an active figure into a static one), is too late.\n", - " var cells = IPython.notebook.get_cells();\n", - " var ncells = cells.length;\n", - " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", - " data = data.data;\n", - " }\n", - " if (data['text/html'] == html_output) {\n", - " return [cell, data, j];\n", - " }\n", - " }\n", - " }\n", - " }\n", - "}\n", - "\n", - "// Register the function which deals with the matplotlib target/channel.\n", - "// The kernel may be null if the page has been refreshed.\n", - "if (IPython.notebook.kernel != null) {\n", - " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", - "}\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "plt.figure()\n", - "plot(100*u, mode=\"displacement\")\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The second loading case corresponds to a fully clamped condition on $z=0$, the vertical boundary remaining in smooth contact." - ] - }, - { - "cell_type": "code", - "execution_count": 56, - "metadata": {}, - "outputs": [ - { - "data": { - "application/javascript": [ - "/* Put everything inside the global mpl namespace */\n", - "window.mpl = {};\n", - "\n", - "\n", - "mpl.get_websocket_type = function() {\n", - " if (typeof(WebSocket) !== 'undefined') {\n", - " return WebSocket;\n", - " } else if (typeof(MozWebSocket) !== 'undefined') {\n", - " return MozWebSocket;\n", - " } else {\n", - " alert('Your browser does not have WebSocket support.' +\n", - " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", - " 'Firefox 4 and 5 are also supported but you ' +\n", - " 'have to enable WebSockets in about:config.');\n", - " };\n", - "}\n", - "\n", - "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", - " this.id = figure_id;\n", - "\n", - " this.ws = websocket;\n", - "\n", - " this.supports_binary = (this.ws.binaryType != undefined);\n", - "\n", - " if (!this.supports_binary) {\n", - " var warnings = document.getElementById(\"mpl-warnings\");\n", - " if (warnings) {\n", - " warnings.style.display = 'block';\n", - " warnings.textContent = (\n", - " \"This browser does not support binary websocket messages. \" +\n", - " \"Performance may be slow.\");\n", - " }\n", - " }\n", - "\n", - " this.imageObj = new Image();\n", - "\n", - " this.context = undefined;\n", - " this.message = undefined;\n", - " this.canvas = undefined;\n", - " this.rubberband_canvas = undefined;\n", - " this.rubberband_context = undefined;\n", - " this.format_dropdown = undefined;\n", - "\n", - " this.image_mode = 'full';\n", - "\n", - " this.root = $(' ');\n", - " this._root_extra_style(this.root)\n", - " this.root.attr('style', 'display: inline-block');\n", - "\n", - "$(parent_element).append(this.root);\n", - "\n", - " this._init_header(this);\n", - " this._init_canvas(this);\n", - " this._init_toolbar(this);\n", - "\n", - " var fig = this;\n", - "\n", - " this.waiting = false;\n", - "\n", - " this.ws.onopen = function () {\n", - " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", - " fig.send_message(\"send_image_mode\", {});\n", - " if (mpl.ratio != 1) {\n", - " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", - " }\n", - " fig.send_message(\"refresh\", {});\n", - " }\n", - "\n", - " this.imageObj.onload = function() {\n", - " if (fig.image_mode == 'full') {\n", - " // Full images could contain transparency (where diff images\n", - " // almost always do), so we need to clear the canvas so that\n", - " // there is no ghosting.\n", - " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", - " }\n", - " fig.context.drawImage(fig.imageObj, 0, 0);\n", - " };\n", - "\n", - " this.imageObj.onunload = function() {\n", - " this.ws.close();\n", - " }\n", - "\n", - " this.ws.onmessage = this._make_on_message_function(this);\n", - "\n", - " this.ondownload = ondownload;\n", - "}\n", - "\n", - "mpl.figure.prototype._init_header = function() {\n", - " var titlebar = $(\n", - " ' ');\n", - " var titletext =$(\n", - " '
');\n", - " titlebar.append(titletext)\n", - " this.root.append(titlebar);\n", - " this.header = titletext[0];\n", - "}\n", - "\n", - "\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._init_canvas = function() {\n", - " var fig = this;\n", - "\n", - " var canvas_div = $(' ');\n", - "\n", - " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", - "\n", - " function canvas_keyboard_event(event) {\n", - " return fig.key_event(event, event['data']);\n", - " }\n", - "\n", - " canvas_div.keydown('key_press', canvas_keyboard_event);\n", - " canvas_div.keyup('key_release', canvas_keyboard_event);\n", - " this.canvas_div = canvas_div\n", - " this._canvas_extra_style(canvas_div)\n", - " this.root.append(canvas_div);\n", - "\n", - " var canvas =$('');\n", - " canvas.addClass('mpl-canvas');\n", - " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", - "\n", - " this.canvas = canvas[0];\n", - " this.context = canvas[0].getContext(\"2d\");\n", - "\n", - " var backingStore = this.context.backingStorePixelRatio ||\n", - "\tthis.context.webkitBackingStorePixelRatio ||\n", - "\tthis.context.mozBackingStorePixelRatio ||\n", - "\tthis.context.msBackingStorePixelRatio ||\n", - "\tthis.context.oBackingStorePixelRatio ||\n", - "\tthis.context.backingStorePixelRatio || 1;\n", - "\n", - " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", - "\n", - " var rubberband = $('');\n", - " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", - "\n", - " var pass_mouse_events = true;\n", - "\n", - " canvas_div.resizable({\n", - " start: function(event, ui) {\n", - " pass_mouse_events = false;\n", - " },\n", - " resize: function(event, ui) {\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " stop: function(event, ui) {\n", - " pass_mouse_events = true;\n", - " fig.request_resize(ui.size.width, ui.size.height);\n", - " },\n", - " });\n", - "\n", - " function mouse_event_fn(event) {\n", - " if (pass_mouse_events)\n", - " return fig.mouse_event(event, event['data']);\n", - " }\n", - "\n", - " rubberband.mousedown('button_press', mouse_event_fn);\n", - " rubberband.mouseup('button_release', mouse_event_fn);\n", - " // Throttle sequential mouse events to 1 every 20ms.\n", - " rubberband.mousemove('motion_notify', mouse_event_fn);\n", - "\n", - " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", - " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", - "\n", - " canvas_div.on(\"wheel\", function (event) {\n", - " event = event.originalEvent;\n", - " event['data'] = 'scroll'\n", - " if (event.deltaY < 0) {\n", - " event.step = 1;\n", - " } else {\n", - " event.step = -1;\n", - " }\n", - " mouse_event_fn(event);\n", - " });\n", - "\n", - " canvas_div.append(canvas);\n", - " canvas_div.append(rubberband);\n", - "\n", - " this.rubberband = rubberband;\n", - " this.rubberband_canvas = rubberband[0];\n", - " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", - " this.rubberband_context.strokeStyle = \"#000000\";\n", - "\n", - " this._resize_canvas = function(width, height) {\n", - " // Keep the size of the canvas, canvas container, and rubber band\n", - " // canvas in synch.\n", - " canvas_div.css('width', width)\n", - " canvas_div.css('height', height)\n", - "\n", - " canvas.attr('width', width * mpl.ratio);\n", - " canvas.attr('height', height * mpl.ratio);\n", - " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", - "\n", - " rubberband.attr('width', width);\n", - " rubberband.attr('height', height);\n", - " }\n", - "\n", - " // Set the figure to an initial 600x600px, this will subsequently be updated\n", - " // upon first draw.\n", - " this._resize_canvas(600, 600);\n", - "\n", - " // Disable right mouse context menu.\n", - "$(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", - " return false;\n", - " });\n", - "\n", - " function set_focus () {\n", - " canvas.focus();\n", - " canvas_div.focus();\n", - " }\n", - "\n", - " window.setTimeout(set_focus, 100);\n", - "}\n", - "\n", - "mpl.figure.prototype._init_toolbar = function() {\n", - " var fig = this;\n", - "\n", - " var nav_element = $(' ')\n", - " nav_element.attr('style', 'width: 100%');\n", - " this.root.append(nav_element);\n", - "\n", - " // Define a callback function for later on.\n", - " function toolbar_event(event) {\n", - " return fig.toolbar_button_onclick(event['data']);\n", - " }\n", - " function toolbar_mouse_event(event) {\n", - " return fig.toolbar_button_onmouseover(event['data']);\n", - " }\n", - "\n", - " for(var toolbar_ind in mpl.toolbar_items) {\n", - " var name = mpl.toolbar_items[toolbar_ind][0];\n", - " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", - " var image = mpl.toolbar_items[toolbar_ind][2];\n", - " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", - "\n", - " if (!name) {\n", - " // put a spacer in here.\n", - " continue;\n", - " }\n", - " var button =$('');\n", - " button.click(method_name, toolbar_event);\n", - " button.mouseover(tooltip, toolbar_mouse_event);\n", - " nav_element.append(button);\n", - " }\n", - "\n", - " // Add the status bar.\n", - " var status_bar = $('');\n", - " nav_element.append(status_bar);\n", - " this.message = status_bar[0];\n", - "\n", - " // Add the close button to the window.\n", - " var buttongrp =$('
');\n", - " var button = $('');\n", - " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", - " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", - " buttongrp.append(button);\n", - " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", - " titlebar.prepend(buttongrp);\n", - "}\n", - "\n", - "mpl.figure.prototype._root_extra_style = function(el){\n", - " var fig = this\n", - " el.on(\"remove\", function(){\n", - "\tfig.close_ws(fig, {});\n", - " });\n", - "}\n", - "\n", - "mpl.figure.prototype._canvas_extra_style = function(el){\n", - " // this is important to make the div 'focusable\n", - " el.attr('tabindex', 0)\n", - " // reach out to IPython and tell the keyboard manager to turn it's self\n", - " // off when our div gets focus\n", - "\n", - " // location in version 3\n", - " if (IPython.notebook.keyboard_manager) {\n", - " IPython.notebook.keyboard_manager.register_events(el);\n", - " }\n", - " else {\n", - " // location in version 2\n", - " IPython.keyboard_manager.register_events(el);\n", - " }\n", - "\n", - "}\n", - "\n", - "mpl.figure.prototype._key_event_extra = function(event, name) {\n", - " var manager = IPython.notebook.keyboard_manager;\n", - " if (!manager)\n", - " manager = IPython.keyboard_manager;\n", - "\n", - " // Check for shift+enter\n", - " if (event.shiftKey && event.which == 13) {\n", - " this.canvas_div.blur();\n", - " // select the cell after this one\n", - " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", - " IPython.notebook.select(index + 1);\n", - " }\n", - "}\n", - "\n", - "mpl.figure.prototype.handle_save = function(fig, msg) {\n", - " fig.ondownload(fig, null);\n", - "}\n", - "\n", - "\n", - "mpl.find_output_cell = function(html_output) {\n", - " // Return the cell and output element which can be found *uniquely* in the notebook.\n", - " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", - " // IPython event is triggered only after the cells have been serialised, which for\n", - " // our purposes (turning an active figure into a static one), is too late.\n", - " var cells = IPython.notebook.get_cells();\n", - " var ncells = cells.length;\n", - " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", - " data = data.data;\n", - " }\n", - " if (data['text/html'] == html_output) {\n", - " return [cell, data, j];\n", - " }\n", - " }\n", - " }\n", - " }\n", - "}\n", - "\n", - "// Register the function which deals with the matplotlib target/channel.\n", - "// The kernel may be null if the page has been refreshed.\n", - "if (IPython.notebook.kernel != null) {\n", - " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", - "}\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "from __future__ import print_function\n", - "from dolfin import *\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "%matplotlib notebook\n", - "\n", - "a = 1. # unit cell width\n", - "b = sqrt(3)/2. # unit cell height\n", - "c = 0.5 # horizontal offset of top boundary\n", - "R = 0.2 # inclusion radius\n", - "vol = a*b # unit cell volume\n", - "# we define the unit cell vertices coordinates for later use\n", - "vertices = np.array([[0, 0.],\n", - " [a, 0.],\n", - " [a+c, b],\n", - " [c, b]])\n", - "mesh = Mesh(\"hexag_incl.xml\")\n", - "subdomains = MeshFunction(\"size_t\", mesh, \"hexag_incl_physical_region.xml\")\n", - "facets = MeshFunction(\"size_t\", mesh, \"hexag_incl_facet_region.xml\")\n", - "plt.figure()\n", - "plot(subdomains)\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**Remark**: mshr does not allow to generate a meshed domain with perfectly matching vertices on opposite boundaries as would be required when imposing periodic boundary conditions. For this reaseon, we used a Gmsh-generated mesh." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Periodic homogenization framework\n", - "\n", - "The goal of homogenization theory consists in computing the apparent elastic moduli of the homogenized medium associated with a given microstructure. In a linear elastic setting, this amounts to solving the following auxiliary problem defined on the unit cell $\\mathcal{A}$:\n", - "\n", - "\$$\\begin{cases}\\operatorname{div} \\boldsymbol{\\sigma} = \\boldsymbol{0} & \\text{in } \\mathcal{A} \\\\ \n", - "\\boldsymbol{\\sigma} = \\mathbb{C}(\\boldsymbol{y}):\\boldsymbol{\\varepsilon} & \\text{for }\\boldsymbol{y}\\in\\mathcal{A} \\\\\n", - "\\boldsymbol{\\varepsilon} = \\boldsymbol{E} + \\nabla^s \\boldsymbol{v} & \\text{in } \\mathcal{A} \\\\\n", - "\\boldsymbol{v} & \\text{is } \\mathcal{A}\\text{-periodic} \\\\\n", - "\\boldsymbol{T}=\\boldsymbol{\\sigma}\\cdot\\boldsymbol{n} & \\text{is } \\mathcal{A}\\text{-antiperiodic}\n", - "\\end{cases} \\label{auxiliary-problem}\n", - "\$$\n", - "\n", - "where $\\boldsymbol{E}$ is the **given** macroscopic strain, $\\boldsymbol{v}$ a periodic fluctuation and $\\mathbb{C}(\\boldsymbol{y})$ is the heterogeneous elasticity tensor depending on the microscopic space variable $\\boldsymbol{y}\\in\\mathcal{A}$. By construction, the local microscopic strain is equal on average to the macroscopic strain: $\\langle \\boldsymbol{\\varepsilon} \\rangle = \\boldsymbol{E}$. Upon defining the macroscopic stress $\\boldsymbol{\\Sigma}$ as the microscopic stress average: $\\langle \\boldsymbol{\\sigma} \\rangle = \\boldsymbol{\\Sigma}$, there will be a linear relationship between the auxiliary problem loading parameters $\\boldsymbol{E}$ and the resulting average stress:\n", - "\$$\\boldsymbol{\\Sigma} = \\mathbb{C}^{hom}:\\boldsymbol{E} \$$\n", - "where $\\mathbb{C}^{hom}$ represents the apparent elastic moduli of the homogenized medium. Hence, its components can be computed by solving elementary load cases corresponding to the different components of $\\boldsymbol{E}$ and performing a unit cell average of the resulting microscopic stress components.\n", - "\n", - "### Total displacement as the main unknown\n", - "\n", - "The previous problem can also be reformulated by using the total displacement $\\boldsymbol{u} = \\boldsymbol{E}\\cdot\\boldsymbol{y} + \\boldsymbol{v}$ as the main unknown with now $\\boldsymbol{\\varepsilon} = \\nabla^s \\boldsymbol{u}$. The periodicity condition is therefore equivalent to the following constraint: \n", - "\$$\n", - "\\boldsymbol{u}(\\boldsymbol{y}^+)-\\boldsymbol{u}(\\boldsymbol{y}^-) = \\boldsymbol{E}\\cdot(\\boldsymbol{y}^+-\\boldsymbol{y}^-)\$$\n", - "where $\\boldsymbol{y}^{\\pm}$ are opposite points on the unit cell boundary related by the periodicity condition. This formulation is widely used in solid mechanics FE software as it does not require specific change of the problem formulation but just adding tying constraints between some degrees of freedom.\n", - "\n", - "This formulation is however not easy to translate in FEniCS. It would indeed require introducing Lagrange multipliers defined on some part of the border only, a feature which does not seem to be available at the moment.\n", - "\n", - "### Periodic fluctuation as the main unknown\n", - "\n", - "Instead, we will keep the initial formulation and consider the periodic fluctuation $\\boldsymbol{v}$ as the main unknown. The periodicity constraint on $\\boldsymbol{v}$ will be imposed in the definition of the associated FunctionSpace using the constrained_domain optional keyword. To do so, one must define the periodic map linking the different unit cell boundaries. Here the unit cell is 2D and its boundary is represented by a parallelogram of vertices vertices and the corresponding base vectors a1 and a2 are computed. The right part is then mapped onto the left part, the top part onto the bottom part and the top-right corner onto the bottom-left one." - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "# class used to define the periodic boundary map\n", - "class PeriodicBoundary(SubDomain):\n", - " def __init__(self, vertices):\n", - " \"\"\" vertices stores the coordinates of the 4 unit cell corners\"\"\"\n", - " SubDomain.__init__(self)\n", - " self.vv = vertices\n", - " self.a1 = self.vv[1,:]-self.vv[0,:] # first vector generating periodicity\n", - " self.a2 = self.vv[3,:]-self.vv[0,:] # second vector generating periodicity\n", - " # check if UC vertices form indeed a parallelogram\n", - " assert np.linalg.norm(self.vv[2, :]-self.vv[3, :] - self.a1) <= 1e-8\n", - " assert np.linalg.norm(self.vv[2, :]-self.vv[1, :] - self.a2) <= 1e-8\n", - " \n", - " def inside(self, x, on_boundary):\n", - " # return True if on left or bottom boundary AND NOT on one of the \n", - " # bottom-right or top-left vertices\n", - " return bool((near(x[0], self.vv[0,0] + x[1]*self.a2[0]/self.vv[3,1]) or \n", - " near(x[1], self.vv[0,1] + x[0]*self.a1[1]/self.vv[1,0])) and \n", - " (not ((near(x[0], self.vv[1,0]) and near(x[1], self.vv[1,1])) or \n", - " (near(x[0], self.vv[3,0]) and near(x[1], self.vv[3,1])))) and on_boundary)\n", - "\n", - " def map(self, x, y):\n", - " if near(x[0], self.vv[2,0]) and near(x[1], self.vv[2,1]): # if on top-right corner\n", - " y[0] = x[0] - (self.a1[0]+self.a2[0])\n", - " y[1] = x[1] - (self.a1[1]+self.a2[1])\n", - " elif near(x[0], self.vv[1,0] + x[1]*self.a2[0]/self.vv[2,1]): # if on right boundary\n", - " y[0] = x[0] - self.a1[0]\n", - " y[1] = x[1] - self.a1[1]\n", - " else: # should be on top boundary\n", - " y[0] = x[0] - self.a2[0]\n", - " y[1] = x[1] - self.a2[1]\n", - "\n", - "V = VectorFunctionSpace(mesh, \"CG\", 2, constrained_domain=PeriodicBoundary(vertices))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We now define the constitutive law for both phases:" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "Em = 50e3\n", - "num = 0.2\n", - "Er = 210e3\n", - "nur = 0.3\n", - "material_parameters = [(Em, num), (Er, nur)]\n", - "nphases = len(material_parameters)\n", - "def eps(v):\n", - " return sym(grad(v))\n", - "def sigma(v, i, Eps):\n", - " E, nu = material_parameters[i]\n", - " lmbda = E*nu/(1+nu)/(1-2*nu)\n", - " mu = E/2/(1+nu)\n", - " return lmbda*tr(eps(v) + Eps)*Identity(2) + 2*mu*(eps(v)+Eps)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Variational formulation\n", - "\n", - "The previous problem is very similar to a standard linear elasticity problem, except for the periodicity constraint which has now been included in the FunctionSpace definition and for the presence of an eigenstrain term $\\boldsymbol{E}$. It can easily be shown that the variational formulation of the previous problem reads as: Find $\\boldsymbol{v}\\in V$ such that:\n", - "\\n", - "\\int_{\\mathcal{A}} (\