diff --git a/doc/conf.py b/doc/conf.py index 6ff4624c933a6ba00d1e84a0bc1122da739ebd7e..a1fea5e7981def738ffa1513814847a72cd862e1 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -24,8 +24,9 @@ sys.path.insert(0, os.path.abspath('..')) #sys.path.insert(0, os.path.abspath('../examples`_ -and the corresponding `issue `_ +and the corresponding `issue `_ for more information.:: from __future__ import print_function @@ -64,9 +64,9 @@ can also be related to a tangent elastic modulus :math:`E_t = \dfrac{EH}{E+H}`. The considered problem is that of a plane strain hollow cylinder of internal (resp. external) -radius :math:`R_i` (resp. :math:`R_e`) under internal uniform pressure :math:`q`. +radius :math:`R_i` (resp. :math:`R_e`) under internal uniform pressure :math:`q`. Only a quarter of cylinder is generated using ``Gmsh`` and converted to ``.xml`` format. :: - + # elastic parameters E = Constant(70e3) nu = Constant(0.3) @@ -80,7 +80,7 @@ Only a quarter of cylinder is generated using ``Gmsh`` and converted to ``.xml`` mesh = Mesh("thick_cylinder.xml") facets = MeshFunction("size_t", mesh, "thick_cylinder_facet_region.xml") ds = Measure('ds')[facets] - + Function spaces will involve a standard CG space for the displacement whereas internal state variables such as plastic strains will be represented using a **Quadrature** element. This choice will make it possible to express the complex non-linear material constitutive @@ -99,16 +99,16 @@ of Gauss points will be determined by the required degree of the Quadrature elem W = FunctionSpace(mesh, We) W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default') W0 = FunctionSpace(mesh, W0e) - + .. note:: In older versions, it was possible to define Quadrature function spaces directly using ``FunctionSpace(mesh, "Quadrature", 1)``. This is no longer the case since - FEniCS 2016.1 (see `this issue `_). Instead, Quadrature elements must first be defined + FEniCS 2016.1 (see `this issue `_). Instead, Quadrature elements must first be defined by specifying the associated degree and quadrature scheme before defining the associated FunctionSpace. - - - + + + Various functions are defined to keep track of the current internal state and currently computed increments:: @@ -124,21 +124,21 @@ currently computed increments:: u_ = TestFunction(V) Boundary conditions correspond to symmetry conditions on the bottom and left -parts (resp. numbered 1 and 3). Loading consists of a uniform pressure on the +parts (resp. numbered 1 and 3). Loading consists of a uniform pressure on the internal boundary (numbered 4). It will be progressively increased from 0 to :math:`q_{lim}=\dfrac{2}{\sqrt{3}}\sigma_0\log\left(\dfrac{R_e}{R_i}\right)` which is the analytical collapse load for a perfectly-plastic material (no hardening):: bc = [DirichletBC(V.sub(1), 0, facets, 1), DirichletBC(V.sub(0), 0, facets, 3)] - + n = FacetNormal(mesh) q_lim = float(2/sqrt(3)*ln(Re/Ri)*sig0) loading = Expression("-q*t", q=q_lim, t=0, degree=2) def F_ext(v): return loading*dot(n, v)*ds(4) - + ----------------------------- Constitutive relation update ----------------------------- @@ -147,7 +147,7 @@ Before writing the variational form, we now define some useful functions which will enable performing the constitutive relation update using a return mapping procedure. This step is quite classical in FEM plasticity for a von Mises criterion with isotropic hardening and follow notations from [BON2014]_. First, the strain -tensor will be represented in a 3D fashion by appending zeros on the out-of-plane +tensor will be represented in a 3D fashion by appending zeros on the out-of-plane components since, even if the problem is 2D, the plastic constitutive relation will involve out-of-plane plastic strains. The elastic consitutive relation is also defined and a function ``as_3D_tensor`` will enable to represent a 4 dimensional vector @@ -166,35 +166,35 @@ containing respectively :math:`xx, yy, zz` and :math:`xy` components as a 3D ten [0, 0, X[2]]]) The return mapping procedure consists in finding a new stress :math:`\boldsymbol{\sigma}_{n+1}` -and internal variable :math:`p_{n+1}` state verifying the current plasticity condition -from a previous stress :math:`\boldsymbol{\sigma}_{n}` and internal variable :math:`p_n` state and +and internal variable :math:`p_{n+1}` state verifying the current plasticity condition +from a previous stress :math:`\boldsymbol{\sigma}_{n}` and internal variable :math:`p_n` state and an increment of total deformation :math:`\Delta \boldsymbol{\varepsilon}`. An elastic -trial stress :math:`\boldsymbol{\sigma}_{\text{elas}} = \boldsymbol{\sigma}_{n} + \mathbf{C}\Delta \boldsymbol{\varepsilon}` -is first computed. The plasticity criterion is then evaluated with the previous plastic strain -:math:`f_{\text{elas}} = \sigma^{eq}_{\text{elas}} - \sigma_0 - H p_n` where +trial stress :math:`\boldsymbol{\sigma}_{\text{elas}} = \boldsymbol{\sigma}_{n} + \mathbf{C}\Delta \boldsymbol{\varepsilon}` +is first computed. The plasticity criterion is then evaluated with the previous plastic strain +:math:`f_{\text{elas}} = \sigma^{eq}_{\text{elas}} - \sigma_0 - H p_n` where :math:`\sigma^{eq}_{\text{elas}} = \sqrt{\frac{3}{2}\boldsymbol{s}:\boldsymbol{s}}` with the deviatoric elastic stress :math:`\boldsymbol{s} = \operatorname{dev}\boldsymbol{\sigma}_{\text{elas}}`. If :math:`f_{\text{elas}} < 0`, no plasticity occurs during this time increment and -:math:`\Delta p,\Delta \boldsymbol{\varepsilon}^p =0`. +:math:`\Delta p,\Delta \boldsymbol{\varepsilon}^p =0`. Otherwise, plasticity occurs and the increment of plastic strain is given by :math:`\Delta p = \dfrac{f_{\text{elas}}}{3\mu+H}`. Hence, both elastic and plastic evolution can be accounted for by defining the -plastic strain increment as follows: +plastic strain increment as follows: .. math:: \Delta p = \dfrac{\langle f_{\text{elas}}\rangle_+}{3\mu+H} where :math:`\langle \star \rangle_+` represents the positive part of :math:`\star` and is obtained by function ``ppos``. Plastic evolution also requires the computation -of the normal vector to the final yield surface given by +of the normal vector to the final yield surface given by :math:`\boldsymbol{n}_{\text{elas}} = \boldsymbol{s}/\sigma_{\text{elas}}^{eq}`. In the following, -this vector must be zero in case of elastic evolution. Hence, we multiply it by +this vector must be zero in case of elastic evolution. Hence, we multiply it by :math:`\dfrac{\langle f_{\text{elas}}\rangle_+}{ f_{\text{elas}}}` to tackle -both cases in a single expression. The final stress state is corrected by the -plastic strain as follows :math:`\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}_{\text{elas}} - -\beta \boldsymbol{s}` with :math:`\beta = \dfrac{3\mu}{\sigma_{\text{elas}}^{eq}}\Delta p`. -It can be observed that the last term vanishes in case of elastic evolution so +both cases in a single expression. The final stress state is corrected by the +plastic strain as follows :math:`\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}_{\text{elas}} - +\beta \boldsymbol{s}` with :math:`\beta = \dfrac{3\mu}{\sigma_{\text{elas}}^{eq}}\Delta p`. +It can be observed that the last term vanishes in case of elastic evolution so that the final stress is indeed the elastic predictor. :: ppos = lambda x: (x+abs(x))/2. @@ -202,7 +202,7 @@ that the final stress is indeed the elastic predictor. :: sig_n = as_3D_tensor(old_sig) sig_elas = sig_n + sigma(deps) s = dev(sig_elas) - sig_eq = sqrt(3/2.*inner(s, s)) + sig_eq = sqrt(3/2.*inner(s, s)) f_elas = sig_eq - sig0 - H*old_p dp = ppos(f_elas)/(3*mu+H) n_elas = s/sig_eq*ppos(f_elas)/f_elas @@ -212,7 +212,7 @@ that the final stress is indeed the elastic predictor. :: as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1]]), \ beta, dp -.. note:: +.. note:: We could have used conditionals to write more explicitly the difference between elastic and plastic evolution. @@ -224,10 +224,10 @@ need to derive the algorithmic consistent tangent matrix given by: \boldsymbol{n}_{\text{elas}} \otimes \boldsymbol{n}_{\text{elas}} - 2\mu\beta\mathbf{DEV} where :math:`\mathbf{DEV}` is the 4th-order tensor associated with the deviatoric -operator (note that :math:`\mathbf{C}_{\text{tang}}^{\text{alg}}=\mathbf{C}` for -elastic evolution). Contrary to what is done in the FEniCS book, we do not store it as the components -of a 4th-order tensor but it will suffice keeping track of the normal vector and -the :math:`\beta` parameter related to the plastic strains. We instead define a function +operator (note that :math:`\mathbf{C}_{\text{tang}}^{\text{alg}}=\mathbf{C}` for +elastic evolution). Contrary to what is done in the FEniCS book, we do not store it as the components +of a 4th-order tensor but it will suffice keeping track of the normal vector and +the :math:`\beta` parameter related to the plastic strains. We instead define a function computing the tangent stress :math:`\boldsymbol{\sigma}_\text{tang} = \mathbf{C}_{\text{tang}}^{\text{alg}} \boldsymbol{\varepsilon}` as follows:: @@ -242,26 +242,26 @@ Global problem and Newton-Raphson procedure We are now in position to derive the global problem with its associated Newton-Raphson procedure. Each iteration will require establishing equilibrium -by driving to zero the residual between the internal forces associated with the current +by driving to zero the residual between the internal forces associated with the current stress state ``sig`` and the external force vector. Because we use Quadrature elements a custom integration measure must be defined to match the quadrature degree and scheme used by the Quadrature elements:: metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"} dxm = dx(metadata=metadata) - + a_Newton = inner(eps(v), sigma_tang(eps(u_)))*dxm res = -inner(eps(u_), as_3D_tensor(sig))*dxm + F_ext(u_) -The consitutive update defined earlier will perform nonlinear operations on +The consitutive update defined earlier will perform nonlinear operations on the stress and strain tensors. These nonlinear expressions must then be projected back onto the associated Quadrature spaces. Since these fields are defined locally in each cell (in fact only at their associated Gauss point), this projection can be performed locally. For this reason, we define a ``local_project`` function that use the ``LocalSolver`` to gain in efficiency (see also :ref:`TipsTricksProjection`) for more details:: - + def local_project(v, V, u=None): dv = TrialFunction(V) v_ = TestFunction(V) @@ -279,17 +279,17 @@ for more details:: .. note:: We could have used the standard ``project`` if we are not interested in optimizing - the code. However, the use of Quadrature elements would have required telling + the code. However, the use of Quadrature elements would have required telling ``project`` to use an appropriate integration measure to solve the global :math:`L^2` - projection that occurs under the hood. This would have needed either redefining + projection that occurs under the hood. This would have needed either redefining explicitly the projection associated forms (as we just did) or specifiying the appropriate quadrature degree to the form compiler as follows :code:`project(sig_, W, form_compiler_parameters={"quadrature_degree":deg_stress})` Before defining the Newton-Raphson loop, we set up the output file and appropriate -FunctionSpace (here piecewise constant) and Function for output of the equivalent +FunctionSpace (here piecewise constant) and Function for output of the equivalent plastic strain since XDMF output does not handle Quadrature elements:: - + file_results = XDMFFile("plasticity_results.xdmf") file_results.parameters["flush_output"] = True file_results.parameters["functions_share_mesh"] = True @@ -303,7 +303,7 @@ the list of load steps). A nonlinear discretization is adopted to refine the load steps during the plastic evolution phase. At each time increment, the system is assembled and the residual norm is computed. The incremental displacement ``Du`` is initialized to zero and the inner iteration loop performing the constitutive -update is initiated. Inside this loop, corrections ``du`` to the displacement +update is initiated. Inside this loop, corrections ``du`` to the displacement increment are computed by solving the Newton system and the return mapping update is performed using the current total strain increment ``deps``. The resulting quantities are then projected onto their appropriate FunctionSpaces. The Newton @@ -338,21 +338,21 @@ the total displacement, stress and plastic strain states are updated :: u.assign(u+Du) sig_old.assign(sig) p.assign(p+local_project(dp_, W0)) - + ---------------- Post-processing ---------------- - + Inside the incremental loop, the displacement and plastic strains are exported at each time increment, the plastic strain must first be projected onto the -previously defined DG FunctionSpace. We also monitor the value of the cylinder +previously defined DG FunctionSpace. We also monitor the value of the cylinder displacement on the inner boundary. The load-displacement curve is then plotted:: file_results.write(u, t) p_avg.assign(project(p, P0)) file_results.write(p_avg, t) results[i+1, :] = (u(Ri, 0)[0], t) - + import matplotlib.pyplot as plt plt.plot(results[:, 0], results[:, 1], "-o") plt.xlabel("Displacement of inner boundary") @@ -364,13 +364,13 @@ The load-displacement curve looks as follows: .. image:: cylinder_expansion_load_displ.png :scale: 20% -It can also be checked that the analytical limit load is also well reproduced +It can also be checked that the analytical limit load is also well reproduced when considering a zero hardening modulus. ----------- References ----------- -.. [BON2014] Marc Bonnet, Attilio Frangi, Christian Rey. +.. [BON2014] Marc Bonnet, Attilio Frangi, Christian Rey. *The finite element method in solid mechanics.* McGraw Hill Education, pp.365, 2014 diff --git a/examples/elasticity/2D_elasticity.py.rst b/examples/elasticity/2D_elasticity.py.rst index cdf5a3f76f923f8e23f73431cee145faed1e7b12..c5d2dd17fc8fba6e6333c26de003a029d270b60e 100644 --- a/examples/elasticity/2D_elasticity.py.rst +++ b/examples/elasticity/2D_elasticity.py.rst @@ -12,10 +12,10 @@ Introduction ------------ -In this first numerical tour, we will show how to compute a small strain solution for +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`.Extension to 3D +in a tradtional displacement-based finite element formulation. The corresponding +file can be obtained from :download:`2D_elasticity.py`. 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 @@ -25,7 +25,7 @@ We also choose a criss-crossed structured mesh:: from __future__ import print_function from fenics import * - + L = 25. H = 1. Nx = 250 @@ -37,10 +37,10 @@ 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 +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:: @@ -52,61 +52,61 @@ for a natural (no prestress) initial state where the Lamé coefficients are give .. 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. + +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 +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 + +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}` +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})`. +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`. +: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 + 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 + 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 @@ -160,7 +160,7 @@ 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}` +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. @@ -173,11 +173,11 @@ space:: 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 diff --git a/examples/periodic_homog_elas/.ipynb_checkpoints/periodic_homog_elas-checkpoint.ipynb b/examples/periodic_homog_elas/.ipynb_checkpoints/periodic_homog_elas-checkpoint.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..d8169416cbebd5b4f95082cecdc7f03e6d8ea830 --- /dev/null +++ b/examples/periodic_homog_elas/.ipynb_checkpoints/periodic_homog_elas-checkpoint.ipynb @@ -0,0 +1,373 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Periodic homogenization of linear materials\n", + "\n", + "## Introduction\n", + "\n", + "In this tour will show how to perform periodic homogenization of linear elastic materials. The considered 2D plane strain problem deals with a skewed unit cell of dimensions $1\\times \\sqrt{3}/2$ consisting of circular inclusions (numbered $1$) of radius $R$ with elastic properties $(E_r, \\nu_r)$ and embedded in a matrix material (numbered $0$) of properties $(E_m, \\nu_m)$ following an hexagonal pattern. A classical result of homogenization theory ensures that the resulting overall behavior will be isotropic, a property that will be numerically verified later." + ] + }, + { + "cell_type": "raw", + "metadata": { + "raw_mimetype": "text/restructuredtext" + }, + "source": [ + "Complete sources including mesh files can be obtained from :download:`periodic_homog_elas.zip`.\n", + "\n", + "We suggest reading first the :ref:`LinearElasticity2D` tour if you are not familiar with implementation of elastic materials." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEclJREFUeJzt3UFuG1eex/HfP51dNmxPG+jOwgsG6MkugMy5gYJcQO4b\nSLmBjJygodzA9AkC+wKB1TeQBQSYVTfCwWyCBtJgtG0MkP8s+BiXyySriqx69d6r7wcQ7CJL9CtX\n1U8/PZZK5u4CAJTpo7EHAAAYDiEPAAUj5AGgYIQ8ABSMkAeAghHyAFCwViFvZm9arDM3s2szOw9/\nzk4fHgDgFHboOnkzO5c0l/TC3e3gC5m9dfen4e8zSS/d/VmfgwUAdHMw5H9bycwPhbyZnUm6cfcv\nK4/94u6/72eYAIBj9DUnP5f0UHtsHcIfADCSvkL+UU+vAwDo0cc9vc5aUv2N1r3Bb2ZXkq4k6ZNP\nPnn6+eef9zQMAMjf27dv/+Xuj/t4rb5CfqUdoe7u97tWdvelpKUkLRYLv7u762kYAJA/M/vfvl7r\n6OmacMnkTPowzM1sLun2xLEBAE50MOTN7MzMrsPfb8IllVs3kv5SWb7cXicv6ULSZe+jBQB00uoS\nyiExXQMA7ws/d7To47W4rQEAFIyQB4CCEfIAUDBCHgAKRsgDQMEIeQAoGCEPAAUj5AGgYIQ8ABSM\nkAeAghHyAFAwQh4ACkbIA0DBCHkAKBghDwAFI+QBoGCEPAAUjJAHgIIR8gBQMEIeAApGyANAwQh5\nACgYIQ8ABSPkAaBgH489AOTj13/+eewhdPbVp1+MPQQk4Puffhh7CKOhyQNAwQh5tJJjiwdy9NEf\n/97v6/X6agCApBDyaESLB+Lou8VLhDwAJGGIgJcIeTSgxQN5I+SxFwEPxDFUi5cIeQATMOWflyDk\nsVMJLX7KJzbyMWSLl1r8xKuZzSVdSLqXdCZp6e4Pe9Y9k7QIizNJr9191dNYAQAdtbmtwSt3fypJ\nZnYn6aWkZ3vWPXf3b7cLZvZC0tcnjxJRldDigRwM3eKlhuma0MzX2+XQ4M8PfMrXZjbraWwYAQGP\nUqU2fRcj4KXmJj+XVJ+aWZvZmbvf71j/RtL/mNnzsPx8xzrA4FI7oYGxNIX8oy4v5u7L0OS3UzS3\n+vCLhMzsStKVJD158qTLP4EB0eKBOGK1eKn56pq1Nm+gVu0NfjO7dvdvwxz+C0lvdq3n7kt3X7j7\n4vHjx50GDADHmuJ3eE0hv9KOUN81VWNm59pcgbNdZynpdZjXR+JKavFTPJGRj5gtXmoI+XqYh8sp\nb6vLlTda19pcYnnwNQBgqmIHvNTuEspLM7vWu+vkLyvP3WgzJbN09/sQ+lfhuZmk73odLQZBi8eU\nfPXpF5P6TVGNIR+a+LaN39aee1Zbft3f0BBDSQEPpGyMFi9xWwMUhBaPtqZ0rBDyE0aLB+IYq8VL\nhDwKMaVmhn5M5Zgh5CeKFg/EMWaLlwh5FGAqjQz9G/rYGTvgJUJ+kkpq8QQ8TlX6MUTIT0xJAQ+k\nLIUWLxHyyFjpDQzxlHwsEfITUlKLL/mkxDj6PKZSafESIY8MEfAYSonHFiE/EaW0+BJPQqTl1GMs\npRYvEfLICAGPWI491lILeImQn4QSWjwBj9hKOeYI+cIR8MDxuhx7KbZ4iZBH4gh4jC33Y5CQL1ju\nLT73kwvl+OrTLw4ej6m2eImQR4KaTihgLDkel4R8oXJt8TmeRJiW+jGacouX2v2OV2Tmy4+eScrr\n91gS7sjJ9nh98+urkUfSjJAv2PZATDnsCXdgWIR8YTYt/n0phj3hjtzl0OIlQn5Sxg57gh2Ij5Av\nyK4Wv0s1bIcOfIIdJcqlxUuE/OTtCuFTgp9QR+lyCniJkC9G2xbfBkENlIPr5AvQZ8AD2C+3Fi8R\n8gBQNEI+c7R4II4cW7xEyANA0Qj5jNHigThybfESIQ8AB+Uc8BIhny1aPIA2CPkMEfBAHLm3eImQ\nB4CiNf7Eq5nNJV1Iupd0Jmnp7g8H1r+oLrv761MHiXdo8UAcJbR4qd1tDV65+1NJMrM7SS8l7Uwa\nM7uWtHL312Y2k/Q3SYQ8AIzkYMib2Zmk9XbZ3R/M7HzPujNJ37j777frSnra41gnjxYPxFFKi5ea\n5+TnkupTM+sQ/nULSSszuzCzczO7DlM9AJCNkgJeap6uedThtebazNnfhsZ/J+mtpM/qK5rZlaQr\nSXry5EmHf2K6aPEAjtHU5NeSZrXH9gX/Spv5+Afpt+ma+a427+5Ld1+4++Lx48ddxzw5BDwQR2kt\nXmoO+ZV2hLq73+9Zt27vVTgAgOEdDPl6mIdWfltdDm+4yt1Xkh62y+HPVXgcR6LFA3GU2OKldpdQ\nXoZLI7fXyV9WnruR9EbSMiw/k/SNmf2ozVw8CQUgeaUGvNQi5EOb3zb629pzz2rLK0nPexvdxNHi\nAZyK2xokioAH4ii5xUuEPAAUjZBPEC0eiKP0Fi8R8gBQNEI+MbR4II4ptHiJkAeAohHyCaHFA3FM\npcVLhHwyCHggjikFvETIA0DRCPkE0OKBOKbW4iVCHgCKRsiPjBYPxDHFFi8R8gAmYKoBLxHyo6LF\nAxgaIT8SAh6IY8otXiLkAaBohPwIaPFAHFNv8RIhDwBFI+Qjo8UDcdDiNwj5iAh4IA4C/h1CHgAK\nRshHQosH4qDFv4+QB4CCEfIR0OKBOGjxHyLkARSBgN+NkB8YLR7AmAj5ARHwQBy0+P0IeQAoGCE/\nEFo8EAct/jBCHgAKRsgPgBYPxEGLb0bIA8gSAd9OY8ib2dzMrs3sPPw5a/PCZvbi9OHlhxYPICVt\nmvwrd//W3W8lLSW9bPoEMzuXdHXq4HJDwANx0OLbOxjyZnYmab1ddvcHSecNnzMLn/PQxwABAMdr\navJzfRjW6xD++5y7+/1pw8oPLR6IgxbfTVPIP+ryYmGa5vb44QAA+tQU8mtJ9Tdadwa/mc0lrcOU\nzkFmdmVmd2Z29/PPP7cbacJo8UActPjuPm54fqUdob5nOuZM0iMzW4TlmZldSbp191Xt85favImr\nxWLhnUcNYHII+OMcDHl3vzez35ZDW7+tLa/d/cHdX1c/18xehDAvGi0eQMraXEJ5ub1OXtKFpMvK\nczeS/lJd2cxmZnYd/n4dvhAUiYAH4qDFH69pumY7NbOdnrmtPfdByoU5+W/DBwBgRNzW4Ei0eCAO\nWvxpCHkAKBghfwRaPBAHLf50hHxHBDwQBwHfD0IeAApGyHdAiwfioMX3h5AHgIIR8i3R4oE4aPH9\nIuQBJIOA7x8h3wItHkCuCPkGBDwQBy1+GIQ8ABSMkD+AFg/EQYsfDiEPAAUj5PegxQNx0OKHRcgD\nGA0BPzxCfgdaPIBSEPI1BDwQBy0+DkIeAApGyFfQ4oE4aPHxEPIAUDBCPqDFA3HQ4uMi5AFEQ8DH\nR8iLFg+gXJMPeQIeiIMWP47JhzwAlGzSIU+LB+KgxY9n0iEPAKWbbMjT4oE4aPHjmmTIE/BAHAT8\n+CYZ8gAwFZMLeVo8EActPg2TC3kAmJJJhTwtHoiDFp+OSYU8gOER8Gn5uGkFM5tLupB0L+lM0tLd\nH/aseyZpIWkm6b8kPXf3VX/DPR4tHsAUNYa8pFfu/lSSzOxO0ktJHySmmc0kLdx9GZbPJb2R9Fl/\nwz0OAQ/EQYtPz8HpmtDM19vl0ODP96w+l/S8snwnaR7CHwAwgqYmP5dUn5pZm9mZu99XH3T3ezP7\nsvLQQtLDvqmdWGjxmKrvf/rhg8e++vSLwf49WnyamkL+UZcXq82/fy3pctd6ZnYl6UqSnjx50uWf\nABDsCvFTPmfILwAYT1PIr7V5E7WqMfhDiH/n7q93PR/m7ZeStFgsvMU4j0KLR2mOCfZjX7tL6NPi\n09UU8ivtCPX6VE1VeMN15e63J44NgIYN9jb/blPYE/BpOxjyYZ79t+VwOeVtbXm9nXffvlG7/SJg\nZhf72vzQaPHI2VjBvkt1LEzp5KfNJZSXZnatd9fJV+fZb7S5THIZAv+tJFW+MKwkRQ95Ah65Sinc\nd6m3e1p8+hpDPrTy7fTMbe25Z5W/rySZEvD9Tz/QOJCd1AO+inMsH22afFZ+/eefJbWfTwTGllO4\nV9Hi81D8vWtyPYFQvu9/+oHjE4MrKuS3Lb6Okwmpyf14/OiPfx97CGipqJBvkvuJhTLkfhwS8Hkp\nJuT3tfi63E8w5I3jD7EVEfJtA36LEw1jKOG4o8Xnp4iQP0YJJxzywfGGsWQf8l1bfBUnHmIo5Tij\nxecp+5A/VSknINJUyvFFwOcr65A/pcVXlXIiIi0cV0hBtiHfV8BvcUKiTyUdT7T4vGUb8kCqSgp4\n5C/LkO+7xW9xcgLvo8XnL8uQHxJBj1Nw/CA12YX8UC0eOFVpAU+LL0N2IR9DaScr0BUBX46sQp4W\nj1RRDJCqbEI+dsBz0mKqaPFlySbkgVRRCJCyLEJ+rGkaTl5MDS2+PFmEPJAqigBSl3zIj/1mKycx\npoIWX6bkQx7A8Aj4ciUd8mO3eOAQvstDDpIN+ZQCnpMZJaPFly3ZkAcAnC7JkE+pxQO7lPLdHS2+\nfEmGfIpKOakBTEtyIU+LB+KgxU9DUiFPwANxEPDTkVTIAwD6lUzI0+KRi9zfn6HFT0syIQ8A6N/H\nTSuY2VzShaR7SWeSlu7+cOq6VbR4IA5a/PQ0hrykV+7+VJLM7E7SS0nPelg3O9//9IO++vSLsYcB\nHIWAn6aD0zVmdiZpvV0Orfz81HXf83//3XKoAICumubk55Lq0y3rEOinrAsgIlr8dDVN1zzq8Fqt\n1zWzK0lXYfHfv/vTPzKq8/9ou+IfJP1rwIGMpdTtklpu2+/+FGEk/fqDZCXus5KPxf/s64WaQn4t\naVZ7bF+Yt17X3ZeSltJm7t7dFw3jyA7blZ9St43tyk94T7MXTdM1K+0Iane/P3FdAEAEB0O+HtDh\nEsnb6rKZzdqsCwCIr80llJdmdq13175fVp67kfRGYeqlYd19ls2rZIntyk+p28Z25ae3bTN37+u1\nAACJ4bYGAFAwQh7IlJm9abHO3Myuzew8/Dlr89yYWm7XmZldhXG/Cu8Bbp+7MTM3s1/M7E31ubG1\n3La94z9qn7l77x/a/GDUtTY/8XotaXbMul1eJ8ZHx+060+ZnAa4lvZI0rzx3I8kl/aLNexrzocfe\n43btHXtq++uIbfslbFv14zq1fRa25Wpz+jau+7by95k2tx5pfC7l7Qpjvap93o+V5au+xxZ5n+0d\n/zH7bKgNaj2QzA7CVuPJ7SDsuL96PQBT2bbw3Pm+bU1tn4UxecPzZ5Le1B77pem5sT9ablf1fJqF\nL8CzVPdV2207NP5j91nv0zV93e/m6HvhDKTjeOaSnleW7yTNU/l2uKqv/+fU9pfUfUzuXr08+EL5\nXwJ86FYj2d6GxDeXa39ZeWgh6cHf3fF2ZmYXYUrjJsXzrsG+8R+1z4aYk+/rfjepHYStx5PZQdj1\n/7nXA3BgXfbZb+uFbXrk7qvKKints7YO3Wqkyy1LklPbN1/r/cu1l+7+OnzR/k7S36IO7nT7xn/U\nPmtznXxXfd3vJrWDsNN4WhyED5JkZmttduLTk0d4nK7/z/vGntr+ko4f0zeS/lp7LKV91tahW410\nuWVJssJ9sL5z99fbx6pfsN39PrxJO/MWv9siBfvGryP32RBNvq/73aR2EB41njYHoaSzEZthp+06\nMPbU9pd0/JjO64GQ2D5r69CtRrK/DYmZnUtaVc+tEIhv6+vmEvAN4z9qnw0R8n3d7ya1g7DzeDI5\nCFtv1xAH4MCO3Wfr2mOp7bO92t5q5NBzKapuV1g+k7Tevo8S3kORNvv8RWW9c0mvlbDatu0d/7H7\nrPfpmvDtxd6BhOW1uz8cWrfpdWLrsl1heXsQ3oflixD2SR2EHbfr4AGY0v6Suu+z4EwfzuMntc/C\nsbW9QOFGmysuttvV5VYjx9yGZDBttyvst7dhve2nryS9dvcHM1uF76Al6TONvF1S+21rMf7O+2yQ\n2xpUNuiD3/VqZq+02cBli3X3PjeGttsVDsIfa5++cvfPwrrn2rwpKG124l9z2K6wvHfsqe2vpjHV\nty08di3pP9z9ee11ktpnQFvcuwYACsZtDQCgYIQ8ABSMkAeAghHyAFAwQh4ACkbIA0DBCHkAKBgh\nDwAFI+QBoGD/DzJW/F1/JEW9AAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from __future__ import print_function\n", + "from fenics import *\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\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", + "plot(subdomains)" + ] + }, + { + "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{equation}\\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", + "\\end{equation}\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", + "\\begin{equation} \\boldsymbol{\\Sigma} = \\mathbb{C}^{hom}:\\boldsymbol{E} \\end{equation}\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", + "\\begin{equation} \n", + "\\boldsymbol{u}(\\boldsymbol{y}^+)-\\boldsymbol{u}(\\boldsymbol{y}^-) = \\boldsymbol{E}\\cdot(\\boldsymbol{y}^+-\\boldsymbol{y}^-)\\end{equation}\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", + "\\begin{equation}\n", + "\\int_{\\mathcal{A}} (\\boldsymbol{E}+\\nabla^s\\boldsymbol{v}):\\mathbb{C}(\\boldsymbol{y}):\\nabla^s\\widehat{\\boldsymbol{v}}\\text{ d} \\Omega = 0 \\quad \\forall \\widehat{\\boldsymbol{v}}\\in V\n", + "\\end{equation}\n", + "\n", + "This readily translates into the following FEniCS code:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "v_ = TestFunction(V)\n", + "dv = TrialFunction(V)\n", + "v = Function(V)\n", + "dx = Measure('dx')(subdomain_data=subdomains)\n", + "\n", + "Eps = Constant(((0, 0), (0, 0)))\n", + "F = sum([inner(sigma(dv, i, Eps), eps(v_))*dx(i) for i in range(nphases)])\n", + "a, L = lhs(F), rhs(F)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have used a general implementation using a sum over the different phases for the functional `F`. We then used the `lhs` and `rhs` functions to respectively extract the corresponding bilinear and linear forms.\n", + "\n", + "## Resolution\n", + "\n", + "The resolution of the auxiliary problem is performed for elementary load cases consisting of uniaxial strain and pure shear sollicitations by assigning unit values of the correspnonding $E_{ij}$ components. For each load case, the average stress $\\boldsymbol{\\Sigma}$ is computed components by components and the macroscopic stiffness components $\\mathbb{C}^{hom}$ are then printed." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Solving Exx case...\n", + "Solving Eyy case...\n", + "Solving Exy case...\n", + "[[65267. 17338. 81.]\n", + " [17338. 65451. 79.]\n", + " [ 81. 79. 24009.]]\n" + ] + } + ], + "source": [ + "def macro_strain(i):\n", + " \"\"\"returns the macroscopic strain for the 3 elementary load cases\"\"\"\n", + " Eps_Voigt = np.zeros((3,))\n", + " Eps_Voigt[i] = 1\n", + " return np.array([[Eps_Voigt[0], Eps_Voigt[2]/2.], \n", + " [Eps_Voigt[2]/2., Eps_Voigt[1]]])\n", + "def stress2Voigt(s):\n", + " return as_vector([s[0,0], s[1,1], s[0,1]])\n", + "\n", + "Chom = np.zeros((3, 3))\n", + "for (j, case) in enumerate([\"Exx\", \"Eyy\", \"Exy\"]):\n", + " print(\"Solving {} case...\".format(case))\n", + " Eps.assign(Constant(macro_strain(j)))\n", + " solve(lhs(F) == rhs(F), v, [], solver_parameters={\"linear_solver\": \"cg\"})\n", + " Sigma = np.zeros((3,))\n", + " for k in range(3):\n", + " Sigma[k] = assemble(sum([stress2Voigt(sigma(v, i, Eps))[k]*dx(i) for i in range(nphases)]))/vol\n", + " Chom[j, :] = Sigma\n", + "\n", + "print(np.array_str(Chom, precision=0))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It can first be verified that the obtained macroscopic stiffness is indeed symmetric and that the corresponding behaviour is quasi-isotropic (up to the finite element discretization error). Indeed, if $\\lambda^{hom} = \\mathbb{C}_{xxyy}$ and $\\mu^{hom} = \\mathbb{C}_{xyxy}$ we have that $\\mathbb{C}_{xxxx}\\approx\\mathbb{C}_{yyyy}\\approx \\mathbb{C}_{xxyy}+2\\mathbb{C}_{xyxy} = \\lambda^{hom}+2\\mu^{hom}$." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "65266.53664945247 65355.28660576194\n" + ] + } + ], + "source": [ + "lmbda_hom = Chom[0, 1]\n", + "mu_hom = Chom[2, 2]\n", + "print(Chom[0, 0], lmbda_hom + 2*mu_hom)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We thus deduce that $E^{hom} = \\mu^{hom}\\dfrac{3\\lambda^{hom}+2\\mu^{hom}}{\\lambda^{hom}+\\mu^{hom}}$ and $\\nu^{hom} = \\dfrac{\\lambda^{hom}}{2(\\lambda^{hom}+\\mu^{hom})}$ that is:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Apparent Young modulus: 58085.12241295444\n", + "Apparent Poisson ratio: 0.20966347372127778\n" + ] + } + ], + "source": [ + "E_hom = mu_hom*(3*lmbda_hom + 2*mu_hom)/(lmbda_hom + mu_hom)\n", + "nu_hom = lmbda_hom/(lmbda_hom + mu_hom)/2\n", + "print(\"Apparent Young modulus:\", E_hom)\n", + "print(\"Apparent Poisson ratio:\", nu_hom)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAATkAAAEJCAYAAADvm1BcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvX2wZdlVH/Zb+9zX06DgNJOoYoQ9Gt7IQmCLjzcPFMsE\nCPRoNDgxVaQHEbDigvL0xBTEUYJGYAHClkBqfczEIXEyPRiC44+SR644dirC9AjHjk2c0NOFCXFK\ngmkJG5NU4Wp1Ukmqp+89e+WPvdc+a3+cc+9973a/d+9bv6pX755z9tln39vv/vq3PjcxMwwGg2FX\n4U56AQaDwXAvYSRnMBh2GkZyBoNhp2EkZzAYdhpGcgaDYadhJGcwGHYaRnIGg2GnMTvpBRgMLRDR\nPoArAPYBfDCe3gfwI8z8xSe2MMPWwUjOcCrBzDeJ6BqAR5n5E3KeiG4Q0T4z3zzB5Rm2CGauGrYG\nRHSBmV8CcOGk12LYHhjJGbYC0XzdBwBmvkFEF4noFSK6REQX1OuLRPR5IjqI9z1PRM+c6OINJwoz\nVw2nHftEdAnA0wDeIyeZ+SUiegzBb3cbway9DQBEdBXAg3Hoi1H9Gc4oTMkZTjtuRp/ckwhkJqoO\n0S93DcCTQnARHwTwHiK6AMB8d2ccRnKGrQAz31bBBjFFhcRuE9FFPTaev2wBCoORnOFUIqq1xwAc\nRl/bJSK6TEQvA7gZTdhPAbgO4HkALxLRZTXF8zAVZwBAm+4nR0TXmPmxJWMOABwiRMm+DsB75H9c\nIroC4BkE0+Q6gKftf2PDqpD0EiK6aL44A7DBwEM0F/YBXFwy7gKAQ2a+qu67BuCROOQVZqZNrctw\n5nCJiF6CqThDxL1QcjxFUlHFvcjMj8TjCwA+D+CLmfk2EV0WAjQYDIbj4r775Jj5BoKvRXAI4LaK\njl1Q+U5XIgkaDAbDkXAieXKFj+1pAE+p46sq3+kWgnP50fu4PIPBsEO47+ZqMfYygFu6NrE1H6Ip\n27j3MgC85jWvefRNb3rTMVZtMBi2DS+//PK/YObXLht3YiQneU06Ahb9dS8w86Pq3NL5Dg8P+fr1\n68dYtcFgOEks/s83VOdc4U1zv/cz2TERvczMh8vmvi/masx5uqXM0IN4fCMeX4pq7iZCfpPcdxHA\nqMozGAynG//37zyUXncUtIomLwdS18dDBCXBrYNNppAcIKaPxFy3a0qlXUFIE7kaCe/lOE5uvwng\nEzG6elMldT6C3F9nMBhOGL/7O1+aXmckhdzgcmgbYB4+EZ0Hp7E9+4zo9LjjYGMkF1XZDQAfblx7\nUr2+CYy8e+Tmq8FgOD34Z//8S9CBsEerpbF68CjRrYPjqDjAupAYDIYl+I3ffh26qLg6AnrmZHpq\n9OBKzY1hTKWVau64BAcYyRkMhhH8yj99GOdpgXMAPAh7xOjB2MtMysHcbEGruTFyLMdtGkZyBoMh\nwy//1j72yGfn9mjIwpgiq14U3whh6XuXqbnZ7/3NI62/hJGcwWDAL33uywEAHRh7BMzZ4YvcXQA5\nwY2R16aV2KYIDjCSMxjONP7bm1+DPVrgfIOf5uywRx7/L3c4Tx57YNxlxnki9Cw+uuXENkWAZaT1\nXpisRnIGwxnDi688inPUp+MOPv4eFNt5WgAAPBOcVnIFB+Xm52okNWXuAsADX7LZBjJGcgbDGcDP\nfuYbsBeJ67xyg+3RAj0cHDN6oCC/vBpqDqrOabLSBNdjuYnbwqYJDjCSMxh2Fs9/+pvQxQCCA+Dh\nsKdITCBKLowLr+fcoYukeE5d7zmoOZeOgyrrmVP267omp4es8d40RTKSMxh2CM/9728DAHyhezU7\n36lo6ZxnSbFpgktjowq7wzOcpwXucAcw8EVukUZ71ES3Dlr3fMGXfHatOVaFkZzBsMX44D/5tkyd\ndRhM0HNRibVwx+/hvJvjDp+L44HzNMccHc6hVnvn4zP2aEj61USn0cqda507CjkeBUZyBsOW4c/9\n+h9Lr8+7dhchz4S7mOEcLZomqvjnxDzdU8R2lzt05DP/2x3u0j3ZcxCIrkVYWSJwIsaB7PT413zJ\nb0295WPBSM5g2AK899e+I71+IMqnjjzm3GEPgKPa7ASAuxy+4proztPd9DqQWX6vmLYlOd5hFyoe\nwFWUFTh6ude9JDjASM5gOJV4169+V3r9gFugU9fmUWmBXUZ0mpQk7aM0WXsmzDFLqmzOM3R0Fx6U\nntHHee/wDHvkcZ4W6MCZshsjrimiEyRFx6vl2R0XRnIGwynB99/44+jAUZWF/DQhrh6uqbicMivn\n3GGP+pASosxPB07E5RHUmINDh2FuYDBZdZBCntmr9JFWKsmqKInxXqs4wEjOYDgxfO+vfC8ApDpR\nR4SOcvLwBWEJXCLD4jypFBD4QHhK4e1Rjw4enh1AfeZnE5KcgpflFQGIsKb18Hte90/XvONoMJIz\nGO4jvud/fioSFANwFVH1IPQF2Sy4CwECCnluQZFRICoE539QcBTJTUqufLxOcCPKa87BdO3gIwFy\nOhfW49CzSxUQQKhllXw5jVYQYlmXkvsBIzmD4R7iyV/+U3AUVJf4yRzVhKOvCzznys4zYY4ukZdW\nXj0okVQ5X89BzfXskjnbwac0kjnPAFqgZ4fXxKDEnGf4wligX5rJWrGt00NO48Lr/tna9xwVRnIG\nw4bx7f/gB+CIMVNE4yPRiKLyTJi5Qa3JdYEQ2Nw7PODK4IHL/GZ6DlF1aWwktFLNtXx8e4VZ7OP+\nUaLY7sKhA+NVDmouBSoi0WXkVwQg7mW/uGUwkjMYNoA/+vf/AwDAzIl/jaNJSZhFQlr4LhBb5Bo5\nFgJK1xsIvrlgUgaFl1/XCq4HYU8RWEc+EVyv1KErqiCy6Kwu5SpM1mBSU1wVcC6ZpnJvG2K6/iuv\n++2REfcGRnIGwxHw2P/wrvT6XKbIarUi54TAfFQ0MwoBgI4C0c2oVFLicxt+y1xz3+EBt0iqTsxV\ngY6wagXnVVRVo2eCJ8Jd7nCO+hiF9ZNR1HlSeu1i/FY6yf0mOMBIzmBYGd/4qXfDgUExaCCqbcEu\nM00B4G7f4VxXk5Y2SaWN0dx3wYcmZfQNnx2AmgQjoXSNsRJdFUUpaI0FAlEltRfH742MDfOslkJy\nv0q3pmAkZzBM4A//4g+jc4EshASYCXtdrYbK3mvlcXlujMw0HPloGgJABxASec19hz3XZ2kje6nG\ntG/OL8EHnUYimKPDeZqn4zu8hy9Shf6eKXUaGcxj4A6HvDmH6J8jYK/hf3vt6/750vd7L2AkZzAU\n+Ppf+DNwJIotoCSshQ8mpKg18au1iE3Ds4PnYLrqsYFwfHoWgCpnrmfCjMRHFsiuZ4e5j6VbRYBC\n1FzJN73y6Yl5qu/Zi5HYc9Tj//PnUpS1j++xA0fzeFhf6Ycri/dPiuAAIzmDAQBw8N//KIA8twxR\nsTETiDiZpSX5CVqqzTNlPrvsuop47om/jofcub6IWDbnUOoKGEhK/HkdhVQRnV8H5OQmStEh+vd4\neOoeLTIynkcfYni2vG/gVQ6kt0eh51za0IaAOR+tOmJTMJIznGl89X/3YyAKX1QhOO1CqhQcD40n\nHQ3KzTNhAZeZtUKAQixORUAl6iqpJD0TelCVLiJ+vKoNuXQPidHZB9w8u0/KsEK+WzA5y9STVsXE\nXZ5VJutrVEF/OBcSjucILZi84rA516ru93/p/1E9537CSM5w5vDmv/W+9LqLnvhASA5EDBeVB1Hw\nv3kAYEIXAw2i7HrvMt9cO8lXMv45kZ6LJKpTRwDggW4Rqxs4ESEIcMp8XXCHWUzoPe/mSXXN/Qx7\nblEFBJwU8CtfnXQInvMs60ji4dKxRGyBGIhQFQ59DJCI2ToVhDhpggOM5AxnCF/5N38CAEBElc+t\nJKipgKAHQdq4DcSVq72hdGtQg5rUznV1bzYxG8tcOW3C6vnmkRz3XJ8Umi4Lc6qNUj9i+ErTzODf\nI/TksEe5KhwqImbVfqzyzKAag8kq/riHf9/JExxwD0iOiK4x82NLxuwDuATgBoADAFeZ+fayawbD\nUfDGv/F+EDE6V5OX9wTntANdIqhA54Zoqqi4dB8oM/1a+XFy3jPBdTmJLiRtJBJhWf2Q1hNrUkuK\n0tUKfSLVfNetUAHhms0ux+YCgomqa1mBocZ1XtSxzuNae+IUYQWA/d/3O5PPvJ/YGMkR0UUA+wAu\nrjD8RWZ+NN53HcALAJ5c4ZrBsDLe8Nc/kJELM4EZQcFFUzReCVUEish8wwwryS6YijG9hBgLP+TO\nlWVUpeIr/WFCbGPNL318loaoM0n6lShr5+bV2j2HZOGhXjWqwCL3bs4zOPgUgOgRTFftp5N7g09u\nMdTIxkqIN//++5/wO4WNkRwzvwQARPT81DgiOgBwS913OxLk5DWDYRU88vGfTK+JEEqoxAcWzTxm\ngufaRO29AyJJdcToffSFueB/m3V5AMAhzIVGtFWCCkARoEipJi6rgNCEGCKYjFkRcJh7hz1Xm59Z\n/zeVM+cKciox7Mw1y9ReqHYYSPVcQZhinjpi3EWHczwEYk4bwQEn45PbB1Can7ciwY1eY+Yb92V1\nhlOH/b/6U1Wul0Zmgpbj1HEiJ8cZwYnJqs/prAdfBB08kOpRgUCOnjipuDEkkzaSr07/KIMSQoTS\nFDMl/HqHB5wvyK/LorJz7vBAVHN3eC/luVXrgQPYo6MhJ09UnMZd7mI9bPD9ab+c3vRmzGQ/aZwE\nyT14xGsZiOgygMsA8NBDDx13TYb7jC/7yx8ML0a+F1SorCkEE7Q1B8CeQC5XXKzUncD72m8mc4R7\nKCg9BELQaSFUBhhiOslsxPRMzyzKvIZzHnN2TSd/WS0h/rRX/SwRnYs1sdKR5FW/V6WYlLjLUvDV\n4TX0anW9A8faWYeePc7HZ92NhWIdPA4euj9NMNfFSZDcLQAXinMPrnAtAzNfBXAVAA4PD08229Aw\niof/6w+l103iEsJpEI+kalTj0Bor5mm4QE5y3hp+MQC9D2ZodV3UiPMoG01m62nAM2V+OUkeLsfo\nDiXwedAhi6LGgIMEHsQv5ylUTTh4eIRIquTLdQi+uc7NK39b2UcuPMPhLgPgDufdHOdiaklL0WnM\n2WV+v697/ekkOGD9jsWbwE00iCuao1PXDDsAbpk0Le5bcVz7GeF3RpAR2iTt/aDOSvKS52uzVcbI\nPetCAhTpWOXOeVDW7bcrFJuOuM69jr4OJqvs5QAEYrzDg4Z51e+pcjGfggeCc7FK4o7fCx1JFDXM\n0eEOz1LeXIm3vP6z63wM9x33heSIaJ+ILgA1YcWUkZeWXTNsHx7+Sx8KBCNck4iDxklMjc8ureDv\nYc5JKUGZnHoebYr6JcRV+pt8NF89aPz9IKg5UXiLFcmx7AY8+Ot8RpSaMEXJhft9RooAmmkkfUGo\n1fUUgID6rVopgU49wQEbJDkiOiCiZ+LrK0VU9AqA71THTxHRM3HMJQBPrXjNsK1Yl7iEGD1l5Jhd\nX0KK3ocfgSQAh59hrGeCawQNRF211ju19oV3Q36c9t8p351ek1Q8AIEU574LiblMsaDfJXLRpCfo\nERSYrFljHoMGd3gvtT7v2SWC080xS/QFUQ6pIx3u8B7e+vqbo/eeJhCfcPHsJnB4eMjXr18/6WUY\nFB7+S1fiq8bfVzNIMPF3uGpwQoKXjrPrzoW8OCIGxfNyTZPdLBbjz7rcnxeCC2GMEGVy/Ls8fUOr\nK+lQon17M5ersVncVlDfJ348HTEVv532sw2BhuHe88oXN9TY5ptF67SSPbVjV0ccC/JdmiusX9Jq\n8pSUb3n40zhJENHLzHy4bNxJ+OQMO46B4IDAPKulFoyqo1X/H27xaRmgQG626v/j5Zzkx8lxMEkH\nMmmts1RXEmVtjRPz1TNhweNfweSnS80xOTNDU8dg8bVJCVdjzvLeFnqmpAiX4aQJbh0YyRlODiP+\numVjl885RElbwYdV0lMyX5h0zPXhXBmRlY4iXJiS2g/nlbmpIflwC98NQYg4j5isOkAxtV7dAaWP\nJm6vghkC8bX1cJFEh+vnlPq7y3Uf4Tu8VwUtTjusQN+wUTz881fyE1S+4KVklQcHisEq1aSV0kET\n/23LeM+UVqPvl+vye6xoXxKCs4TitCgCIC2UfH5NQXbzqs4XH47k5DmEtI1AZDXEFJ1zl21cs0qb\n8l6SjolVqZhOZQmb6PTs8G/v/9rS+U4bTMkZNoaK4JpYHmxIIyPhVOquYWJWkVNVoZDd2iDQlsnK\nPKSYzH2XEZVWVcxthSYQs1VM05aik+ul4pKAhKg9H3PTSpPVqyiohhCcBC3EN6fN254d5tzhVd7D\nHR9+pDV6SPylRHh+S+liO1dt2C1oU1T9XhrN5HxcqGyQVBJJ7chvScSm5g194wZf3EpLVhUQYygj\noTq6KqSnfXIlEQJD+yXBPJqhc+6yfLkFd0nFCXGJWSl93/Q5WVePeg/XUMYViPMOn0tE+u37v7ry\n53OaYCRn2Age/q+utM3QkrwADMGIkaBE855CrXFj/ASYCVyQmOTGDcnDtSLTJNlH0szOpTliyggG\n9alVXOnf04EJHVwYtiEckoNf9bNgsjby3wSz1Mp8aHtekl3LN6fRUbvK4S53Tf/ctsBIznBsPPxz\nHx4OWvlr3Lh2DyAkVqaQtFRgSyVKBLVMDhZlKMm/JRnquXREdN7nqimNidHVqgIivtZ5c0C+/0M5\nn6g5zxTNTAlahLVKTpwOFohKawUQUsVEZk47PPnIy9XYbYGRnOHkkL7fq6eZAMurH9irGtTGtZLg\nNKlJBLWE7vohik4wZraWnYdb5qhWeJlpC045dPVaXPqtN6B+lffCptJpnaVyq/1+PqaN3OVZGp/8\nfDG6+443/Erz/W0LjOQMx0JScVwQS8us1CjVHYAm2a2h/IiGyK32y2XXG6/zOcLvVt7b1HO9kKU4\n9aUdeVRsY/659rN88teJ818rK0kMLmtJS/+dRlb7yl12POzE1VWR4Hf+gX80Oue2wEjOcCxUBTMl\n2QFHME8b/rqC7Kqo68QzOBEfJd9cmdPWwqLvqppWjn43MVmF3HQ6iz4nxzoHLs2fBQ6cCky46j5J\nJUndSyKkG2/feFYJibZKfzpJHi7vmXOHnt1OEBxgJGc4Bl7/s0HFjRbGj5HQWqQ3reyaEVhV97oM\niaiKubIVKEITf1d43jCHPD8FHFAkIxfPlKiqrnMtgxQ+jelSPWuJYT8HSoX2QakNe1DIWvri95zz\nNFmpa3VgfN8b/8Hk57ZNMJIzHAmv/4sfbpJVU9ml11hOcOtUNqxynQnsXUbE3AiAyCq9H/xtybfV\nu6zQH6hNw6YfLzr+tT9PzzulJMeiqFK0D+QEp+FiUu8dvzcQnCr0DxHbunxL5tm2ioZlMJIzbAZT\npKNNWK22ygjs6BzjqSaTJVoMcPTbB6ILXUmqLsKlOct595LyGWKOlo59UWHzvpsksIV31QbUUyhN\nVA/CXT9TnUlk79X66yzR1aTqJtJIpGri6S//e0vXtE0wkjOsjdf/xUbKCHKlNGm+VkGK1lg0iK+d\nCpL8XzxET9uVEsP2g9o/NwYdaQ3VF+F81nhTRVq1uVqa0ZJOIvdLYEHvC9Eq1i83uUlzQG9W7dO8\nWoVJSgmA0cJ7zw5znqWUkl0jOMBIzrAmXv8zHxlIqqVCMn9Z/lOPVfO08usac4Znqp8Vcu+yTWmi\n6anN0nq8BBKi6ecdFn3XrKDQyb3pGSPO/9b5u7EwX86PlX4BwTcnfjdJFJ5HM1SK/PVzdHWDdAOW\nwn1HPovOenb4gTf9UvsD2XIYyRk2h00n+ZbBimr+/ERSZXFs3dJcXlDVOLOcXpTbsC2hH0xZRUop\nebbRoaRXHUjKziU6EbhUdDIGyBODHXGVTqJR5tTp3nNitraaBQDAu77iF5vndwFGcoaV8fqf+Uh9\ncixlZJWgxBQYuXXavFfU3OCzYz+YoKmSoaHatN+NmeB7Nzy24efLUkRG/HWBFGWT6MHnVrZQAoao\nKjDs0Zr9qHQReV1tSK2ipTrRVwrpxwIIOvAwtm/DLsFIzrASHn7hIyBG+qkwRXYjJmzz/ikTmAF4\nAnpC6gRU+exUOoeubmDdKqleDxBMWSnfKhtrrhIgkPmnCvfLetZUBZH8ayoKq4ITDqFN+l0/SwnC\nQCjMH+auqyPmsdgeyDeeDvMGM3eXVRxgJGdYFSundiz31aVTq6aTtAiQI9EtW1cWFFDBkZHSL/aU\n1bz6RHj11LrAXqNzg8Nft4IqqyaWpZIE9RZaK931Iact9YlTc0nDzVf9DIuYT6dV3BRBP/OVnxy9\ntiswkjMsxcNXazNVq7pRZVedQ5OUKmWnSI17Hb3M7xkWU85L+bUJcNGJRC8VGN+qUJu6ZRRXlKDG\ngt2kwluwyxKDpyBdg4FAelrBOfLN+1/1syzlpWe3cqvzbYeRnGE9jBDVpAk7oexqX1n87ZFy3ORY\nrldqbCINJIxRvrp4b5Urh2GDm7KIX6DbKpWBiT6akLI+7YfTc/U+J7uyBVO5beFYkjGAqpnnXpEL\np9ssAUP7JelK8uN/6G/VH8IOwkjOMImWihvDkZSdxkICBQBp4hKT1U8QWkm8ilxT0ICh8ulqBclM\n43nHK/jkhKjKIv+qRfuEmSo+u7t9V+XHpQAEhu7APRMW3MWmmaH3HFC3a/IIHY4F73/zf7P0/ewK\njOQMo3j4v/zo9IAJf1iT7EaCE9QTqCcwxdfSF87T8KMJ1A+EV5m64qeTdYh/rSTHJpHl43QBvlZu\n0pa9hE9Bj3bXEq+CItJTrhojEVfnh/rWuAeDBCBkXGamxvN7KTE436NB4ywRHGAkZxhBIrgp9SRY\nQnbtewLhUT8QFvUigaZVWzanMkfTvGphyUQlpeIKgmIVTdV1rVluWyQkosHk1FsYJnWG4Vz5nJay\nG/O/CQEu2ClTtsuacnp26NQOXXoTag3x35Xm61mBkZyhwsP/+cdCmkavvjAsqRsrkF4DWaBi7kAL\nAi2KuTzCM1jdpH2ALcLUWw6qdJGg6MT2bEdR8xNxhy7XZmVRb2WhvlRGtIIKKWG4aKkkyHLpij1W\ny44kupoBGMq65r5LJqruM6e3ORyqJBw+9FWfaL6/XYZtSWjI8GU//THAKbW0oNqd1vFAThLZhPot\n4x0DCwLNOCmuBEIgIoeGPy3+pG3ra+IhjhwmiyMOwQkqggosa6S0gxegGmwW7y0r4Fdqrtz+UOpT\nS39bq84UGJpnalO3j/fPYspJq1OJbtt0zi2qeeW+OYdytRn1RYmZgyQVfvSrP95c267DlJwh4ct+\n+mMAgCynVNSXL8774qf13Y6BBMxdoQoRFVuhFmXekjS1qhSfXFPVTb8/9pRIOe0HIdPyoMpWQb7H\nw3C/jnbqQn2nyFIrNo0pH1243jY1U4G+DkjEKgjDhkmOiPaJ6Bkiuhh/X5gY+3ki4uLnmXjtSjz+\nPBFdI6L9Ta7T0IZ2+kOILQYG0hgGaB79ZxxNzjheyIfmlMboe5PvrTD5aJGP1WsYi5imtbC6lu6h\nIRqrH0UA95qcCL4njFVApNQTNX7oYqKipIqcZLObVJLVyKPTcMRZiVdZASHnKrMVNfnNG11MHHnM\nvTuzKg7YvLn6IjM/CgBEdB3ACwCeLAdF8nuSmV9S5y4z89V4+AqvErM3bAz7/8mz4YWQGkHZhHlK\nB8tuWFHZucykrf/ZMgUYlRTx8DrdoghVfsuzklmrJ4rH5GJ6iJ5I4JtLGoZlpi2pO+sCfmBIAnaO\nkwkKDN2DgUCGra39wrhohqontcxUICg7KdwXs9WBcbefBf8bIQs6OPKYs8MMQ2E+APz5r/1rIx/A\n2cDGlBwRHQC4JcfMfBvAxbHxBcFdAvDS2FjDvcUjzz6bzNHKLG1AEx7FhFxqKLQ0pvS3la+norMy\np5jFGowhOKEJqQyYrBIoKVNR9CWv61jlN7Kk3/Le0hzNOpco5adz18pIq1aDpSk7tbGOBBwAIzhg\ns+bqPoDbxblbkfwyRAIEkFTdg8x8Uw25QESXotl7pWX2EtFlIrpORNd/93d/d1Pv4czhkY89C+pR\nRzA92r4vUU+eRs3KlNtWmqtyXUValyELLqjnV2ZseR9jqJgoS7OqyGocVgQ4kjGR/8qgyUbv+cCR\nmLRZCwSCI7UpTUceC+/SPqy6HrbcplDOexDu9rNhXHMfCYe/cPCXGys+e9ikufrgEe/7EQAfLM5d\nFSIkolsAPgXgUT0gmrZXAeDw8HCFr4uhxBs+/BxQ+rIjYekIY7rURb5QyilLu/IAxVy0lNnhEPxw\nOvKqXhMAOK55SoiH4xza5yaTyfmo5tLYYt2DFRvHUgwyOM5bJfUEcgDF3nFV3aosX4IJihS9d3BF\nTzghs0Xv0Kln9d4BcWxXxl2ictNmaolQr5qfX7DDDO1OJGcdm1RytwCUimsV4ruolR2QKz1mvgHg\nYCqIYTgiODdNy6gqkIsgtwCoj+PiT+rLKOP7/J5VzF8JRowpu9H7y/FUiDaKik7y+1hODoRZFtgH\nBRj7somZKoStgx7EyT/XCjAsK7KXZGJdNF8W6Jc5c2luUCzo78KPmkNy4q4e/vzk888SNklyN9Eg\ntUhSTRDRRSg/Xjx3QEQvN+YpTWHDMfAHPvRceCGVB4sh8jn8DONb/rYQeR1+GmlcgzL0w/hwc2ts\no4QrPT8fN2yvhSEdRZ7XQukYFILTJjlpBZvvB1EtNSq9Vh84PUbSUqZIT4rmBdpslUirRGEFs4L5\n5bhnMoIrsDGSK8kspn3o4MJ+Q40doPbj3QTwvLrvIoCzl6Z9DyEEp0lnSMUYflKaSHL+iyKK51KN\nqZix0Q8n1QwjJEWiAHnsZyCtlrpL9aslypSRNZwYpZ+ulThMxLkrMJaB6Ry4spZWR1k1WY6RIpB3\nMNGmqz4X5vDhB0Mjzp/7up9b/U2fEWw6heSpmOt2A4HAnlLXrgC4huhHU9ABBzDzbSK6SUSX46lH\ninkMx8AbfzL44Sgm3bJTCswNgki+gzJOW3npPAD0BO646fynBWXVCtQrn5l6Bhd+uoToQ0t8MJUO\nku4R39wbOThfAAAgAElEQVSg7KhTD+Lc4SgVCENlQ7jM0W9G4r/Tvju1jNBJGOji+xTfXDsYQFmz\nSw0X1+KZcK4bUkBEwekKCCCWh0ldq6uVnWEA8VqN908nDg8P+fr16ye9jFOPL39/UHBN57wQnJCc\nL3irQS7kw1zNcS3SIq7HMjJylbWk+xVJpjGkjnUQooyOOqRAQ+I1x9WaKI0RZo/j43HX+ey5Mq6L\nCkvn0zniFIAQ4iPijNwccbpXrmuS0tfOdX0VaMhIEIxzXfhf6q+85QWcJRDRy8x8uGyc1a6eEQjB\nAbnpKAFH9Ag1q9Kc0tXRSopKKt2D4ZwWKLr6KM0XU0GoIBgAoX+cnHNhLSMbyGfI1qdTTSRBV13P\nBFw1EXJ5VozXSi79xmBejm1w3cc62FK9eaYQhcZQ+rWAG/xq3iWi89H/6JnSdb0Bjhz/tX+9NJAM\nAqtdPSMok32FfMaiq+m6CiqILy2ZsCqoUD7LlVFWHUjgeKIo0xIk83XM96bNXDnWSf46MFFiYr4p\n6DZNVbaLjkBnCm24V3LmBGm7QkWgZd1ruU1h6CtX586N1bQaAozkzgDe9BPPNb/IiZy4+K2ur5q0\ni4LQsvOM+hk6HWNZRYLUvBbzjgYg4nWmketFpJbVT3MqIdUyMZfr/Vt137kxSLCitbmNoFWsr6Hv\ne/Gt/8Xk2LMOI7kdx5ve91xOMgKlwlbKZSvubSpC/SPPEBQBjfyaOlBz15Kp5RjM700qkVYkZygz\ntPqMcqWZWqir5UuPuWxrw6jadORVEoNbLvAW2Yn6k58Si9gIwAhuOcwnt8P4ih97LvtvTNJFWBz7\nhSM/KKP8nKCKiqrAQpMgxc/WCBaUsYnkqysDE6WfTAcWWkpJ+tN5FdWl4ri6V4d2i4XpaPISxtR5\nc5qwSr+cmKQUCVD8b8l8Jcae6+sHyFtU8/+Nt/6FyTUZAkzJ7Si+8r3PZb4xraC0f00gvreElukq\nY1rmqB7L42ZuljaSTrbfQ/L9pTU22jTp/VN1EKKct9wPIq2bRwmMRRUiN1V1Tl0r+NBqr6TbpSfC\nW1IVUVU6FPlzhtVgSm6XUUimJum0VNiIj053Epe4QXZ5SuwoVUYxkitzpkBDVIUpIuqK+4GQk6eP\nx/6bjoECKupiV22KqaEDCACyQv5CbAJA6gA89jwh2q7R8bjlFywrIGbk8bf/jZ9e702cYZiS20F8\n5XulZAupvrRSYBGillLlg/jYxghLVSvowMSYTy8Lbox1EJZzqoJCfmdzu2Ft1fq01Sm/CxJJhNPY\np0G2I9Svw8VxP1p6LFMKOLT2kZVzrc1udGmYjGkFOPQYI7j1YEpuByEkJF1DuFB04Fx9iO8skQvC\nPTqDgVWCrlZeet5U0VAqPEksFkW0Sg6cR+UHZA4mq/a3pffHmP4vm3jIfZMgQ3Ia5iZrVt1QflbZ\ntbCAUnuFDadddg9F/xvQTiLWc865Q+d80zf3yW/88xNv0tCCkdyO4Q/+8HOJVKhHHWVsp6bVfrXi\nsiYdKslPm4OFj45d47wipOT20uuSbGNFpHlAIr8eXlIuo1x+f5qHAUIxf7w/9QUugy5So+ram9eU\nJmvLDG0h7OPqsgoHYDzIYX64o8HM1R3CH3r3c7mp1+qk2zArVx1XmqiAIq9GACLN1cJYUMIvyf2o\n0koUAao9VqvnxCinLtcqKx3K7r8AqkaaKdWEBzJK3BxN3ZKMZGNqYMiP69kps3VoseTVHK/2uQa5\n9s3PwbA+jOR2BG/+oecG/1UWkYwv5ByF6Kr+AQZSHN2RqzV34Ysj8e0pwivHTImRjKziWpu7dxVI\nnEjKtmxVUrTIXIiuVKCNJGLv831XW2NkXO8pJQtPbVdS+vB0+smMQmrJwjsjuGPASG7X0CACt1AB\ngoKE0j2i3gpnf5a8W74WIvOYFF8aZdAgHTeCCRWBTilGAKPlXJNlG62ognqpUlRaZqTezSs9rqiC\nGMPYHq2GzcJIbgfwVe/K/5fXEc+sBrVRYwpunJfo6UhEtkSLnMaeVRFla45lSb+rIKuEGEvEKw51\nW6hWBcQIWmauRFV771KOnCg8gezX0DNlm9uU1//ut3xstYUYmjCS23J81buUmapVWKnUInSibpXn\nNqLIUvffBiHK/NUzgbZfT+XZteYJ97WjI02u0vMRq5IyypVdIiJ1rZinlQeX1KoyO6fqXFtqT/xw\nZVNMgSQHi29Oz/H3v/UjzecYVoeR3Bbjq//0c1nLIwnSMXIfW7LWxvx1UOSlFRahqWaqTiQofi/z\nnQmWBSxiQi8wBFQB5EpPro+ov7KyojYjY26LrINyEqPUqklMVjW3Vm7xt1eVDWGMehLlRfnMYUcu\nnUsnqo+ZjOA2BCO5XQLnZKeVWZmUm9JLhAiLlKymaiuCACsV9euC/ZE1Z88cI8gsIEJpfKbUhFwm\nnpcIpWyeCaRNqocTdWCgJK3yPbSWrzedBoBFr6KpGH6A6f1UDUeDkdyW4qv/dPTDlQpMdwUh9RqK\n/Hj4nX5QK7uWGszIRhKIVeAgCxBolGSrMeF3y1xrrb9W9V6SDd5SoMrHRssiqowqIRgj6qyxlKwC\novmeiDOfXTlv7x3+4WNXRu83rAcjuS3E1/zAc1kgoSyurwrt1c9U8XxrLkCNLeaqyrS4vm+U9BrI\nXGU6n02O9WumupIDaJqyo3M0Bma+Nj1VJLrBL1esvZEbp++VEq4WAZY+vn/0eLkNseE4MJLbMnzt\n9z+XkZRTEdAsqqlUnUbGG0qFVRwkpFj63lCTZL5dYP3M5E8rf4ox2VyM0DlEqa4s2tuMkBTHjrOd\nxIJCm7hnjUiubHCTyIlrsmodL5vTCG7zMJLbUpAH3ByJiJppHIy2z02RjDZlowt+Wu2NkEKm9rwi\nR2C5yekb70EW03pm9WYLBVjM0bxWbl9Yzlcps/BcvavXmNlatlpipqrUq6UId2BPqVMJI7ktwsHT\nzw1qrfCrpdfSdUT9RkFACaWPri/GKgyKKPxohZUdpxvG50mvNZEycjXHyCOmek4AYBo1uVdCFmQY\nCI/UDlt6/czBZ5dFXidUYEl0Oj9urDb1V574qaO8E8MSGMltCQ6eVgm/yWQbTtGYGkKuzDJ1x8jN\nwxKlP02Nr5J9y+dyPqYiQX2+BBXzMQZFRvVNTKjaKgEYalll2lFrsbbXS1M03c/5Hq3ZGkdIV+aS\nUq8WbnzbB8YWZzgmjOS2AId/8tnhQH+ZVF1pE8WXsBVUYEKm7sZcXc3zBamWJm62nSHLRCNrXfKX\nONS1Dg/S7ZvK6HLT1K3ewPQzqzVkpqWouvZ8RByWpIv948uyusEI7t7CSO6UQwiOGHA9p5/Kx6ZT\nOYqUj1J9aTN2zJ9XIvsul+SpiDcLFIypQ6/GAs2/Qi7VnEbZAl0CCyWJkuoKXLJx9abiqUL9tVD6\n66oKiBiEoGrM9LyGewMjuVOOIVct/4YQR7IrTMgUZCginrrgvlXSVUYwqzQSGSPzFf607LfcwCPP\nKyOpOmFYSFJ8Y4x87wYho4mggX6d8uegggVTeXmrOPpaZWHAoCZXUIiSSvKP/633Lx9sOBY2SnJE\ntE9EzxDRxfj7wsTYK0TERPR5IrpGRPtHmWeX8XXfF1VcWfCuRUkkvzLVI9sisLiPqSC9lj+vMWf5\nbLm/PC6VXjOfbQzJFA9EkvHJiE+vip5WBF4HFprjgMoPN4oiD67yzyFvW+6KxpieCf/rH/uzSx5i\n2AQ23Rn4RWZ+FACI6DqAFwA8OTL2FR5PHFpnnp3E1/+JZwG1MbomOnZRMMSqBkIwy9jJNxm5WVhu\nCBOHuR7NlA7SZBg3nUm3tay5SJrcDfdVQzzq1ugj//rlmqgf5gYDPGswkN6CUB3roARpshz5773q\n+sukIq4DmZVBDL3ZjY7Mkgu+ubID8P/27T/RXoBh49iYkiOiAwC35JiZbwO4eFLz7AKa3X2BPHIp\nXzoOJix5rk1bGTuxEcxUKVbTtNV+tdLpD3Uuzq2DEHq+tA4hiZLgtPlHADqurwODn04HB8qghD5Z\nfq5Tyi3en6u1qDRVvlvWcqmlOuN1I7j7i02aq/sAbhfnbkXSauECEV2KJukVZZKuO8/O4S3vfHbI\ncyt20UrBgpaJmfnWOf3I/a3+cFkQgnNS1QSTvUYxJp4sa2LlfGlO6utVlxDOn5XeV4sY9XteNWku\nFeqi/iyOlXjXmIfyIAUzGcGdADZprj645virUaWBiG4B+BSAR1edh4guA7gMAA899NCajz69eMs7\nn63OpY4hkfAqtVOmTUCNK80qH7/bK6RYiErL9klNEyGRTxZYEPJyao6JuTPyimapfn+1iizCrq33\nkfndOLVLqu7X0ZQ100laSLuBITdb5fVmKNSwLjap5G4BKAMEo4QlBBdf3wBwENXcSvMw81VmPmTm\nw9e+9rVHX/UpRVVaVZBHGufbrwUcv3XEHNUa10GM4hlNjF1jRVhyasTPV3YsaRGUTvlg7WdLkUt1\nrhWUILWmYm4Ate9Ov4+SUVcKlFCVQjLmav70d/z4ChMaNo1NktxNtMnoRnmOiA6I6OXG2NvrzLNr\n+MPf/TG4BcMtlF9N/GyRGFxfE2AiD6qvSV6dEJygFY3N5kPtL5vK6s/m1ceKQFPyLqG5VSGrQEsz\nYlqYsBmXOMQecVwr2sk11zlvo9BE3phzal8HI7iTw8ZIriShmBLykj5WfrebAJ5X1y4C+MQq8+wq\n3vqOjw01p8n/NvjTNEhIS/vpeIJgUJ7PiTPDRABCormMfIzuXScR2cqvJuMKIs3WyQPpkRTQF+Qm\nKSOZWhRiE8e/UoNQJqQct99ceV6NSEW1aBNdo9he+/c+8+/82PjkhnuOTaeQPEVEzwC4AeAAwFPq\n2hUA1xB9cUR0M/rVAOCRYuzUPDuHt77jY8FXpv1TC+QpJP3wpWGHoPSIwnkqeqvJd1LSNlh9PxMZ\naeZA5mPLvsyaZEjNpW93+XPTuZaPsHid0tKq44GQ0jUVWW0mNJdoRUiyz2lNL1kVqBjm1vWs0oaJ\nwPjMJSO4k8ZGSS6qMFFiLxXXniyOR9XZ1Dy7CCEq+eKm757OD9PjI3lRz9mu9ule5auSoEWzUW3h\nx0pBBAyvoY6bqSStHLyC9FIKSYOUiAEf/wonidHxQGw0nAu/G++tCD6k95nONd4LYa0Iawo0FMEQ\nOW8EdzpgZV0njD9y6aPZMXkGmDOSauXLrbS/QgFiPV/jy8zDONcPAYUyNaSqXfX5mDI/TRNxy6dX\nlpAlOAxlXWOlVOVb0E0BCEOFQ3XreiouM1kjyk1tmi2YDCeOTZurhjXwDd/x0fBF9oD+0iVzFIi2\nj1wIyk6U3EAe+ReZHWU+MV3GJcfLoMVJlvBbfNczvutkPW0FOvr8QsFVfeRGJ0Jm0nPXIm4lT7Pn\nFlGMRurHsrQSnTJS4jee/NHpmw33DabkTgjf8B0fjdUJ9bVmpQO3i/HdYkQJFd/hcv7wux3YaJml\nmXorr5MiwnicghFleog6pwMIVSpIqzOJlGo5Xk6EOtiQHTfWXz6n8b/ApBkbr8kYI7jTBSO5E4J0\nEHF9XopVRzs5C9tpknMLPW4wG0MrJkVgNFyrTU9O/rdmaZeGpH+MkERWalbM0yyyV+eydBWlHFM9\nbCK4YpyUl3WM0p5uFc2Pvq+py/ofRbv6qlQVwm9+pxHcaYOR3AngG789bhocv8ikCEoTiNM7cukv\n1xL1l9xHI11/624j+Ze4TECulFo5pyYoOeVQN/SkSFpjhOO4eKMYmmSWkFMdp6AEGn6z0fuWnVPn\nl21AI3jlHe9daZzh/sJI7j7jm/7ohytnfXL4S0NMUXYsai9e9+0C/JVR+tjUuWzYSBS0OZ2s3bfv\nSwNL1dOIdFKfy8TmfHohHdeEJnPpoMMo2tJ11QirNcHcDhjJ3Ud807dFgusZ8AxEkxUoTMkRpTaa\nqKvP5xZbCjqMdTMJYwamyCy+0p8GjO5Oz91g7mqTUkc7Zf48V22YO+tUos3S8j4gBBxKFlwlYLFO\n5LNJsupleh9kKu4Uw0juPqLZPUP5z+R8Vr5VpJQAQlxR6bVaK+m9WKcCBvq6Jkrl96vqUCUSSvV0\nElgou6OUhDW2nhRQ0AGKkgAdqr/alQIRafDIApbcztx4wxGvfNefWe3ZhhOBkdx9wje//QrID3Wk\nghQxFWVXEN5AepyrrAKB9ALBNfc2LccrZaRrWlmSeYGYyoJkXmfrak6aH1brKO8rkoWp0RNukKSN\n57iiiF8uHzNPbZXoquXCbQ+M5O4D/s3HrxSBg7putNV+HFAkqCHkV+zbUN1bqLyUdiIpJ8U9iZQm\noqRQ9avpUpGHV1U66GdQPgYYKh4AhC7Es3y8VEJkeXC6DnXUF6iZT58fGb8Mjc/DVNzph5HcfYSQ\njhO1xoPJmcAA+oaqU+OzOZvkhqW+p1aLojGyBAqyEkWpHpNapascuGx9I+vJAxCc+Rir5N4y767w\n9wGFwmp+OHqO9putVVt7nBHcdsBI7h7jW771Q5WiYiK4hc8CAuJjGwMtVJDCqy4iaR+HfEd5nSA8\nRGuRqqPSfa3KBIm+Fn6xdE391s0vs/QQ9ZdF+rpKb9FEy47TdSnlygIoySxl1RNuDUm2bGiKytbu\ngqykS973UdWg4b7DSO4e4lu/+aeGJF9RZh5Dom5MxBVIoCEjRc/hB4HohsFVJuryBTVUEI2YrojK\nLNWv6u7DSwIIaT79g/x13ipp8MdxIkOux+nSNKlnLR+7zF9YBR7yOVbNibv575qK2xYYyd0jfOs3\n/1R6TXORTfW4qZy3VtrH4E/LFUciSEUkrXKvpvgR4tEqb8Q/qNfFBUfIfdxScSXBTpnHLSJ0AOIu\nXSvyUPHQ0UHZUWWaJuU2zGUEt10wkrtHcKK6pOa0H/xv5HlQdsjVmyv8cMu6jdR1pjxuxXFjfKkk\nNXc6RWiuvrd8TorMFg76zCQuiDELOmhfoh8GpS0IPdXNMsu3Q+paWT1xXLQCMYZTDyO5e4DHvuEn\nASBvY47CFG2YscO4QZmlYw/AUWqqoa9rZGpLfyllrwdNJEphZT5Dh8k0FDFhGQA8cr9c4b8Lfe+W\nu89Sfpy6Nz1L0I10Mx5dZD5XY9AKE+W4+d2m4rYNRnIbxmN/5ANtx1Dr3IipKpvODL4sRZQ8kB5P\nmHxV/SoUAbZIZ0rpjKRfEDD0fNNpJCrXTtaSVTUUfrqq+B6An6kE6KJN1BCwWMMndxRUbk9LjttG\nGMltEI+99QPDgU73UNuuV2khUdUNgYY4tLyXuSLFoNbqb3VtkhbO9UbqRZhPq04V2Wx9t2mYJ/ut\nI6ZCSGPpI7L5jEzJ4Zg7ZXITB19cx+lwJSXXiJIeF5/9nh85/iSG+w4juXsECTaIOUo6qloSndyz\niJtBF0SXfakjKVb+MXH6l+qjG8zkMlhRzZ0mw2AuugaxjAkaMVtF2bUCCGpsMHMpRFZ1E1BNcJPP\nW5G5WikwjYmroIP533YCRnIbwtu+/s+BvM97w/U+/IgPrvDPOZ0SUib6tlosjZBjWdu67Ltf7q4V\nXuj58rnK9JG0QfWSbsO+G8a35p4ykbmxf0PaqYswGnRollvx9LNWgam47YWR3Abw+KPvGxRYrFYA\nEI7L6KiYpwvO25wDzW8oNVQboFRh2eV3UqGppFb9uwUeUWCcb4qT+ciOVN8azNMyjQVA+Ossduga\nGgrUgYWV/HE8csANn1s8NILbbhjJHROPf+37MvOSvJipeZJaUHntfJAs+TdGTYU0mfJKhjHyKCse\ngDx1g91AcFUDzFa+mzxLk93Is6tW5zRCkNl6lVobI9tGe6llvv8m0Y2aq8vPG8FtP4zkjoHHv/Z9\nw8HCh29Yr4hOfGqLyBpekR1FH5wyQamPcxBSj7ehOiKOUWSoyWTwf5GaD4PvTfnkdHulstdclj5S\nkpUbfHQaVfJvvLd8LZHVKslXntUVUdZZToKTTTllmok8wNGbWrBA6s7ASO6IePvX/HjWogjAYKZK\n3pv3QdExgxY+H6/8RHUhviKzqncaDcQl9azxm52RAHNWK1r2o2spLU2QqS+cL0zIEo3zVfJveQsj\nj6C20Gp3tMQEHvXHHQVsKm5XYCR3TFDfp59JMIO7NlPo82WCL7e+uTr6ScgUmpi+7ChGdDH4u7Ln\n6GcgP6/nLgMH+vmS8wbkyrJh5o6Zr9L/ThMjeWS+uPaNdaR4k/jsHzeC2xUYyR0Bb3/ze4G+Dz8K\n5Au1JoGH0mwVRDJobT8I5uFfJ6aGJMe7Iq1ASkUe3IRZp/u5tfq9NRFJrYrgFlHLVMyvP5aG768s\nrm+Wrq3a6VfmReGPW6lu1XAWsFGSI6J9InqGiC7G3xcmxh4Q0eU47kUi2lfXrhARE9HnieiavnbS\nePub3zvkvsV6VPHFAQgmKnP6SZDgxMIH39yiHYjIfWoMtxD/XjznedzMlfuSfw55gq66NvjqItnF\nTayrPDsaxmXnoVRb+R66nDyXma/SdSSte6/9vu5ZwUExr6m43cJs+ZC18CIzPwoARHQdwAsAniwH\nRfI7ZOar8fgigGsAHolDXuFTWkNDCz/i/EEIPgAA90DXBV8cAPRCIgw4B2IeAgvMiKNAPasoaCQd\nV6g0fewo3hOvdZQ1q0xUoRVXob7CnPLm4o82PQuCnPTPpYHqWSUSAU/Iq/JSY81j19Ou9kuL88sQ\ns9x/Kv/sDMfAxpQcER0AuCXHzHwbwMWR4fsA3qOOrwPYn1J+pwFPvPE9QO+BRZ8UHPVcm6FEgC/I\nsJn1P3wzsyjrpENeq0MGz0ipvHyo3ny6JEtgCVkViq/ymSmEQII8dFCHle8Qg8pL/eOyNJc4vtUR\neGrxVUR1GcGVeSnDy8+984dbNxi2GJs0V/cB3C7O3Yrkl4GZbwB4TJ06BHA7EiMAXCCiS9HsvXIa\nyO+JN74nP9H7gfA8p9SQygT1Pg8e+JAmUtWmRkjQoIRWUCnqOStVnkwSzvtCp1c+OB1YmFJ35Vqq\niG9OfL5rEGhD3aVqCuJ8zhH2XacR8DQU+6q1GcHtJjZprj64zmBmvqkOnwbwlDq+KoRHRLcAfArA\no/p+IroM4DIAPPTQQ0dZ78p44g3vBpwLBOXcQEzNyCclZScbsoQ0EjdEN8WmSkSH/DvXBxXIHSVz\ntox6MiFUTcQ5k5XZ8KsJJHghr/Xv9GzZhlDGcizPkutiBqsARrX+xhrIK1+drEH/9ZUmabZYVIGP\nyV5y+t5sQH1KYAS3u9gkyd0CUCqupcQXyerjzPwJOacUHZj5RgxSXCjOXwVwFQAODw/vWfzsif0f\nArxsYU9D/amjmuR8ZAV9rWegi+Zrp20zDo0htd8rlV1RU80BuaJLpFkRCsPPaCAshC0NvfjsGn41\n8uE6oIiIhvvhkXQ/xVgLSsUm1yOhgYdzzb0korIsc/jSQ44DJnzue5853hyGncAmSe4mGqQWTdMm\nYsDhJjO/pM4dAHhBAhhqntIUvud44qH/MBCWgAo7zfuB9Pbyj5KJhogrEAkyvnYAu2EurdzSOa9U\nWlR8vkjWTc9qBSc4JxbZkFqXd1XrzeYYxhDnpm96V1T89khVEWOOEDFr/Tn5z2LgMy++uCUE91t/\n8t2T1w0GjY2RXFRc6TimfbxUHN9SZuhBPL4Rjy9FNXcTwPPqvosAkso7MfjoQBKiWywA15InGBRZ\nynMLUdXsXDlWDpNZN5xnotrUEh6gxrnyvD7NscC+5Z8rxqWecAVZjkK3SyqTfJUSbNXC6nk/d9lI\nzLA5bDqF5CkiegbADQAHyP1sVxDSRK5GwnsZABQx3gTwCWa+TUQ3oxkLhLQSPc99wRNf+oPBvOyV\nR33RA+iBrmAqRyExWIgsBha0WgOQp4tEdSbHVWVD9NuNVUmEOYSwOJElx+NAMoooHTXN1HQNQ5R0\nzKem78/MT6Uayec+vGGCYc3sOJE59YDfC9c+9/0/NPpeDYajYqMkF1WZmKcvFdeeVK9vYsINrM3X\nk8ATX/qD4YUQV6pAcIVfDgB7APFbLYGERQ/MOtDCg2cF0fXRDJUgBnKCKzdxDjc1FinBB9/2zeUE\nN5i/WqFpv77u4Jsl8o40vJQbCTUhuqji/FRgAQhdVvaAz/7gf9x4iMGwGWxayW09nvjX/tQQ+aSc\njDKyA4D5HJh1gQzJhQCDoPeD4ovzsAQePIaxWvUxJ3OYGCEQ4AFCnlYSfGRKuaF24Jf+srFARphv\nUJhSliUEFVJVlnxoRUZGlWKirmmV98oP/UdLJjYYjg8jOYUnXvvvx1cxQqqJDqijo8BgznZA8KLr\n4ABlY6nvA9EJMWVjw3hRekFp5SYtxQqIqqqBRZrlz01RzrRmVD68ZtRT7nU5YWm1l+XctfyFjdQS\nwW++28jNcP9gJBfx9gefAmkS85xHVoHcfG2h90HZAXV6CTCUdIEqfx0ApE4lyseWXdYE14f1kTI9\ngwIbCM7PUBGNVnxDf7dhnqZ5WpIlMGxLKKkiQEaakiAsPr5P/9i7GhMbDPceRnKCvs8rGR0BfTif\nkR8wKLtZ/PjEjC2DBL0HKN5bkIeuX83OS/mrqkcN5+v0j7E60tZuXcTtgEBmBhfPDusciKoZUFD3\nVkTojNwMJw8jOQBv/5e/L7zQlQwtJaehKyCAIRihy7VmBcH1fTJ562gq0vO4Yc5mEdDWshr96kri\nKe/VycLVdJro+oHoymtpvIqu/pMPGLEZTg/OPMklgtNI5mgknfkCJMm+ouISualKBukvJwqPKHfq\nq5QQ6nvwnjJtCYmoOFZMaPWmm2O2HF2Vr07nohUEmMiuSggeCRq0qhKKIIcRm+G04kyTXEVwWpkB\n2RaCrM1WokyVVZAIbEwlSfe0oDJQErKC/vh7wlQE2mkfZcH+VCKvLv6vkoSVX09vYfjrHzFiM5x+\nnGmSywIIQixaqWUBBgpENytbe8T7StNWF+GX5WACqQLgyCJEbd/WTAUGOhqetYI/LlUglPNWfrv6\ndfQ/3igAAAlxSURBVLmOX3vWSM2wfTizJPf23/O9+Ymys0hllhYpI0I0zFHVxbR9Umkko7sdx2sp\nuFDaiDycd/k9YWPnoUmAkJ/r68DAYLLmeXZpk+bCBxieoVRdD/zqf2bEZthunFmSY9k2sJXKsQpa\n+XCN66kjMACezfJE4Ab8bCBGYg5dhacCIBGSRqLTQzLTszQ3ZX4g+AxVOyMjNsMu4UyS3OP/0p9I\nryuym1J0RGCO4+EGv5xzg8Kb5fMQM9gT0NEQbMhqTd1o0iwAeFUWJi2Uwnrk/tx0lS0EgSH9pJlq\noo5vPG+kZthdnDmS0wSnwd7nqk4n5GYEOKGqxOzUJV2C2FeurGfN8uWYQax7vsn5IlqqVZkU8hc+\ntaHjbktlAtd/xqoODGcDZ47kprr6VkTXghBYKuuKaSPSsUQrOaKidKuudKj7uA3ENhZRHfrMyRzR\nJzejOvAQj/+XnzdSM5xNnCmSe/w1/97SMU2ia6SWMFSktRVhjekltOjBDzQistrvhmC2hn1Ulfry\nALmhY2/KpwMGNaeCB1kumwP+p79q3T0MhjNDco9/wTvDCyfVDOriWARUm6ya6HRkVZd2oUF2QDBV\nHWJ01IPhALUJTTO6qqKvWetyN7RL8sVGNr/8143UDIYSZ4bkElSC70bm6ntgb284NxEJTabqzAUz\ntbXjVqHQwrmGaR39cf/j37QuugbDFM4EyT1+/nvCi7Gk3KLrx6hvTtRcK7AAhPOyabITNTYosmSa\nir9O/HexTjavK2X48hnRbP17f9uIzWBYFTtPcm87992gZF7GNA9yuaKT3nFAbbqWbY9SbznXVm3K\nlOVU7wpgr44ipEhrvMcthhSRUOUQIq9/9xffU91rMBhWw86T3EbQ6u9WBh+A6TrVmFOXmZ5jScFE\n+KVrtg+owbAJ7DTJve3cdwMA2POg5lqQtkqrmKxACkAkgpvPgXPngrnae+CBc+G694MfTtJJZEpd\niRBV5LV/+KPrvkWDwbAEO01ya8HHppe6GkGqIXS3EZUYzHfnoC84P1zTPjTn6iiuqiUVGLEZDPcW\nO0tyb9v7rvAiBhsyNad9cwLtl1OYypmjrgPuzoHzD+T39kVeHZCpxL/z8p9d/w0ZDIYjYSdJLhHc\nMrAfj7imIdFsLYMPQmKdC+aqpJF0HdLuXOeG1JJf+MfvX+ctGAyGDWEnSS6DUm0cI6pN/9xIu/Os\ncL8kOjFPZ7Nm3zha9Pjkpz907LdgMBiOjp0juce6dwAYIbISpbkKjFc/lHAOvFjk0dU4xyd/48Or\nzWEwGO45dorkhOBWRpk3V+TLLYuulvjkb35kvecbDIZ7jp0iOY2ptJHsmlZzhcm6tPIhmr+f/O3/\ndGPrNhgMm8VGSY6I9gFcAnADwAGAq8x8e92x68wjeMw9GSceSGlpflyJxtgpovuF/+tnV5/bYDCc\nCDat5F5k5kcBgIiuA3gBwJNHGLvOPJuHTvf4f37+vj3WYDBsHhsjOSI6AHBLjpn5NhFdXHfsOvMI\nPvPyTbyFvizesDwthBzh79z5K8veksFg2AFsUsntAyhNyltEdMDMN1Ydu+o8RHQZwGUAOI8vBABc\n8y8e7x0YDIadwxG3qmriwQ2NXWkeZr7KzIfMfPgHH/0KIziDwdDEJknuFoALxbkxwpoau848BoPB\nMIlNmqs30SCjhqk6OZaC03/VeQwGg2ESG1NyJQnFNJCX9DERXVg2dtk8BoPBsA42nULyFBE9gyG/\n7Sl17QqAawCurjB26prBYDCsDOJGe6Ftw+HhIV+/fv2kl2EwGO4jiOhlZj5cNm6TgQeDwWA4dTCS\nMxgMOw0jOYPBsNMwkjMYDDsNIzmDwbDTMJIzGAw7DSM5g8Gw0zCSMxgMOw0jOYPBsNMwkjMYDDsN\nIzmDwbDTMJIzGAw7DSM5g8Gw0zCSMxgMOw0jOYPBsNMwkjMYDDsNIzmDwbDTMJIzGAw7DSM5g8Gw\n0zCSMxgMOw0jOYPBsNMwkjMYDDsNIzmDwbDTMJIzGAw7DSM5g8Gw0zCSMxgMOw0jOYPBsNPYGMkR\n0T4RPUNEF+PvCxNjD4jochz3IhHtq2tXiIiJ6PNEdE1fMxgMhnUx2+BcLzLzowBARNcBvADgyXJQ\nJL9DZr4ajy8CuAbgkTjkFWamDa7LYDCcYWyE5IjoAMAtOWbm25G8WtgH8B4AV+PxdQD7RHSBmW+v\n8czLAC7Hw1eJ6NfXX/mZwr8K4F+c9CJOMezzWY7T9hm9fpVBm1Jy+wBKgrpFRAfMfEOfZOYbRPSY\nOnUI4LYiuAtEdCnO9xiAD7bILypBUYPXmflwQ+9lJ2Gf0TTs81mObf2MNkVyD64zmJlvqsOnATyl\njq8KqRHRLQCfAvDosVdoMBjOJCZJLpqEj0wMucbMLyGYqmWgYSnxxfk/zsyfkHNatUXVd7CuKWsw\nGAyCSZKT4MAKuIkGqZWmqkb02d2MJCnnDgC8IAEMNc8yglt1nWcZ9hlNwz6f5djKz2gjKSQlmcW0\nD01e+zqlRAIVQnDRBwcEsnxejbsIIKm8iedv5Yd/P2Gf0TTs81mObf2MiJk3M1EgrosAbgA4QO5b\nexHBtL0aCfCV4vabzPxIHHsRIZABBFO5GXgwGAyGVbAxkjMYDIbTCCvrMpwpENG1FcasXL2za1jx\n89mqqqRNVjzcU8QP8hIa5vBxxu4K1vx8rgB4BiEX8TqAp4u0np2DcoOMJalrrFS9s0tY8/PZqqqk\nrTFXiehl9Yd3ASEK2/zDW2fsrmDNz+fytjqRjwsi4qkvaPQtX2Hmx9S5zzPzF9+XBZ4wln0+ccxW\n/f1shbnaKhvDyP8464zdFZzF93wPMVq9cxKLOaW4QESXojl/5bSb81tBcljvD+8s/pGu+5636o/0\nPmOt6p0ziqvM/ImYAvZxhKqkU4tt8cmt84d3Fv9I133PVjo3jiNV75wlbFtV0rYouXX+8M7iH+la\n77n8IwVwYGouYe3qnbOESGgvl+dPK8EB20Ny6/zhncU/0pXf8zb+kd5r6IqcZdU7ZxFFxdKRqpJO\nEltBcuuUjZ3FP9I1y+q27o90E4jk/kx8faXod3gFwHeq46ckTw4hLUd3ydlJrPr5xP8Mb8bO3pcR\n2qGd6s9nm1JIViobWzZ2V7Hm52Olc4Yzg60hOYPBYDgKtsJcNRgMhqPCSM5gMOw0jOQMBsNOw0jO\nYDDsNIzkDAbDTsNIzmAw7DSM5AwGw07DSM5gMOw0/n+byNOviH3a6QAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# plotting deformed unit cell with total displacement u = Eps*y + v\n", + "plot(0.5*(dot(Eps, Expression((\"x[0]\",\"x[1]\"), degree=1))+v), mode=\"displacement\", title=case)" + ] + } + ], + "metadata": { + "celltoolbar": "Raw Cell Format", + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/periodic_homog_elas/hexag_incl.geo b/examples/periodic_homog_elas/hexag_incl.geo new file mode 100644 index 0000000000000000000000000000000000000000..10a46a29e55360c2b9013a250a20de18ca2b376b --- /dev/null +++ b/examples/periodic_homog_elas/hexag_incl.geo @@ -0,0 +1,75 @@ +// Gmsh project created on Mon Sep 14 09:33:57 2015 +a = 1.; +R = 0.2; +t = Pi/3; +b = Sin(t)*a; +bR = Sin(t)*R; +c = a*Cos(t); +cR = R*Cos(t); + +d = 1.; + +nx = 80; +nR = R*nx; +ns = nx-2*nR; + +Point(1) = {0, 0, 0, d}; +Point(2) = {a, 0, 0, d}; +Point(3) = {a+c, b, 0, d}; +Point(4) = {c, b, 0, d}; +Point(5) = {R,0,0,d}; +Point(6) = {a-R,0,0,d}; +Point(7) = {a+cR,bR,0,d}; +Point(8) = {a+c-cR,b-bR,0,d}; +Point(9) = {a+c-R,b,0,d}; +Point(10) = {c+R,b,0,d}; +Point(11) = {c-cR,b-bR,0,d}; +Point(12) = {cR,bR,0,d}; + +Line(1) = {1, 5}; +Line(2) = {5, 6}; +Line(3) = {6, 2}; +Line(4) = {2, 7}; +Line(5) = {7, 8}; +Line(6) = {8, 3}; +Line(7) = {3, 9}; +Line(8) = {9, 10}; +Line(9) = {10, 4}; +Line(10) = {4, 11}; +Line(11) = {11, 12}; +Line(12) = {12, 1}; +Circle(13) = {5, 1, 12}; +Circle(14) = {7, 2, 6}; +Circle(15) = {9, 3, 8}; +Circle(16) = {11, 4, 10}; +Line Loop(17) = {13, -11, 16, -8, 15, -5, 14, -2}; +Plane Surface(18) = {17}; +Line Loop(19) = {1, 13, 12}; +Plane Surface(20) = {19}; +Line Loop(21) = {3, 4, 14}; +Plane Surface(22) = {21}; +Line Loop(23) = {7, 15, 6}; +Plane Surface(24) = {23}; +Line Loop(25) = {10, 16, 9}; +Plane Surface(26) = {25}; + +Transfinite Line{1,3,4,6,7,9,10,12} = nR; +Transfinite Line{2,5,8,11} = ns; +Transfinite Line{14,16} = 2*nR; +Transfinite Line{13,15} = nR; + +Periodic Line{1} = {9}; +Periodic Line{2} = {8}; +Periodic Line{3} = {7}; +Periodic Line{4} = {12}; +Periodic Line{5} = {11}; +Periodic Line{6} = {10}; + +Physical Line(1) = {1,2,3}; +Physical Line(2) = {4,5,6}; +Physical Line(3) = {7,8,9}; +Physical Line(4) = {10,11,12}; + +Physical Surface(0) = {18}; +Physical Surface(1) = {20,22,24,26}; + diff --git a/examples/periodic_homog_elas/hexag_incl.xml b/examples/periodic_homog_elas/hexag_incl.xml new file mode 100644 index 0000000000000000000000000000000000000000..e71a76b3e04407f197da224ca1e58960096883b1 --- /dev/null +++ b/examples/periodic_homog_elas/hexag_incl.xml @@ -0,0 +1,20013 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +