diff --git a/doc/dev-doc/manual/appendix/material-parameters.rst b/doc/dev-doc/manual/appendix/material-parameters.rst index 4ab43b393..622aef2cd 100644 --- a/doc/dev-doc/manual/appendix/material-parameters.rst +++ b/doc/dev-doc/manual/appendix/material-parameters.rst @@ -1,129 +1,150 @@ .. _app-material-parameters: Material Parameters =================== Linear elastic isotropic ------------------------ Keyword: :ref:`elastic ` Parameters: - ``rho``: (*Real*) Density - ``E``: (*Real*) Young's modulus - ``nu``: (*Real*) Poisson's ratio - ``Plane_stress``: (*bool*) Plane stress simplification (only 2D problems) +Energies: + +- ``potential``: elastic potential energy + Linear elastic anisotropic -------------------------- Keyword: :ref:`elastic_anisotropic ` Parameters: - ``rho``: (*Real*) Density - ``n1``: (*Vector*) Direction of main material axis - ``n2``: (*Vector*) Direction of second material axis - ``n3``: (*Vector*) Direction of third material axis - ``C..``: (*Real*) Coefficient ij of material tensor C (all the 36 values in Voigt notation can be entered) - ``alpha``: (*Real*) Viscous propertion (default is 0) -Linear elastic anisotropic ------------------------- +Linear elastic orthotropic +-------------------------- Keyword: :ref:`elastic_orthotropic ` Parameters: - ``rho``: (*Real*) Density - ``n1``: (*Vector*) Direction of main material axis - ``n2``: (*Vector*) Direction of second material axis (if applicable) - ``n3``: (*Vector*) Direction of third material axis (if applicable) - ``E1``: (*Real*) Young's modulus (n1) - ``E2``: (*Real*) Young's modulus (n2) - ``E3``: (*Real*) Young's modulus (n3) - ``nu1``: (*Real*) Poisson's ratio (n1) - ``nu2``: (*Real*) Poisson's ratio (n2) - ``nu3``: (*Real*) Poisson's ratio (n3) - ``G12``: (*Real*) Shear modulus (12) - ``G13``: (*Real*) Shear modulus (13) - ``G23``: (*Real*) Shear modulus (23) Neohookean (finite strains) --------------------------- Keyword: :ref:`neohookean ` Parameters: - ``rho``: (*Real*) Density - ``E``: (*Real*) Young's modulus - ``nu``: (*Real*) Poisson's ratio - ``Plane_stress``: (*bool*) Plane stress simplification (only 2D problems) Standard linear solid --------------------- Keyword: :ref:`sls_deviatoric ` Parameters: - ``rho``: (*Real*) Density - ``E``: (*Real*) Young's modulus - ``nu``: (*Real*) Poisson's ratio - ``Plane_stress``: (*bool*) Plane stress simplification (only 2D problems) - ``Eta``: (*Real*) Viscosity - ``Ev``: (*Real*) Stiffness of viscous element +Energies: + +- ``dissipated``: energy dissipated with viscosity + Elasto-plastic linear isotropic hardening ----------------------------------------- Keyword: :ref:`plastic_linear_isotropic_hardening ` Parameters: - ``rho``: (*Real*) Density - ``E``: (*Real*) Young's modulus - ``nu``: (*Real*) Poisson's ratio - ``h``: (*Real*) Hardening modulus - ``sigma_y``: (*Real*) Yield stress +Energies: + +- ``potential``: elastic part of the potential energy +- ``plastic``: dissipated plastic energy (integrated over time) + Marigo ------ Keyword: :ref:`marigo ` Parameters: - ``rho``: (*Real*) Density - ``E``: (*Real*) Young's modulus - ``nu``: (*Real*) Poisson's ratio - ``Plane_stress``: (*bool*) Plane stress simplification (only 2D problems) - ``Yd``: (*Random*) Hardening modulus - ``Sd``: (*Real*) Damage energy +Energies: + +- ``dissipated``: energy dissipated in damage + Mazars ------ Keyword: :ref:`mazars ` Parameters: - ``rho``: (*Real*) Density - ``E``: (*Real*) Young's modulus - ``nu``: (*Real*) Poisson's ratio - ``At``: (*Real*) Traction post-peak asymptotic value - ``Bt``: (*Real*) Traction decay shape - ``Ac``: (*Real*) Compression post-peak asymptotic value - ``Bc``: (*Real*) Compression decay shape - ``K0``: (*Real*) Damage threshold - ``beta``: (*Real*) Shear parameter + +Energies: + +- ``dissipated``: energy dissipated in damage diff --git a/doc/dev-doc/manual/solidmechanicsmodel.rst b/doc/dev-doc/manual/solidmechanicsmodel.rst index 374e930f3..63022b19c 100644 --- a/doc/dev-doc/manual/solidmechanicsmodel.rst +++ b/doc/dev-doc/manual/solidmechanicsmodel.rst @@ -1,880 +1,883 @@ .. _sect-smm: Solid Mechanics Model ===================== The solid mechanics model is a specific implementation of the :cpp:class:`Model ` interface dedicated to handle the equations of motion or equations of equilibrium. The model is created for a given mesh. It will create its own :cpp:class:`FEEngine ` object to compute the interpolation, gradient, integration and assembly operations. A :cpp:class:`SolidMechanicsModel ` object can simply be created like this:: SolidMechanicsModel model(mesh); where ``mesh`` is the mesh for which the equations are to be solved. A second parameter called ``spatial_dimension`` can be added after ``mesh`` if the spatial dimension of the problem is different than that of the mesh. This model contains at least the following six ``Arrays``: -:cpp:func:`blocked_dofs `: +:cpp:func:`blocked_dofs ` contains a Boolean value for each degree of freedom specifying whether that degree is blocked or not. A Dirichlet boundary condition can be prescribed by setting the **blocked_dofs** value of a degree of freedom to ``true``. A Neumann boundary condition can be applied by setting the **blocked_dofs** value of a degree of freedom to ``false``. The **displacement**, **velocity** and **acceleration** are computed for all degrees of freedom for which the **blocked_dofs** value is set to ``false``. For the remaining degrees of freedom, the imposed values (zero by default after initialization) are kept. -:cpp:func:`displacement `: +:cpp:func:`displacement ` contains the displacements of all degrees of freedom. It can be either a computed displacement for free degrees of freedom or an imposed displacement in case of blocked ones (:math:`\vec{u}` in the following). -:cpp:func:`velocity `: +:cpp:func:`velocity ` contains the velocities of all degrees of freedom. As **displacement**, it contains computed or imposed velocities depending on the nature of the degrees of freedom (:math:`\dot{\vec{u}}` in the following). -:cpp:func:`acceleration `: +:cpp:func:`acceleration ` contains the accelerations of all degrees of freedom. As **displacement**, it contains computed or imposed accelerations depending on the nature of the degrees of freedom (:math:`\ddot{\vec{u}}` in the following). -:cpp:func:`external_force `: +:cpp:func:`external_force ` contains the external forces applied on the nodes (:math:`\vec{f}_{\st{ext}}` in the following). -:cpp:func:`internal_force `: +:cpp:func:`internal_force ` contains the internal forces on the nodes (:math:`\vec{f}_{\mathrm{int}}` in the following). Some examples to help to understand how to use this model will be -presented in the next sections. +presented in the next sections. In addition to vector quantities, the solid +mechanics model can be queried for energies with the :cpp:func:`getEnergy +` which accepts an energy type as an +arguement (e.g. ``kinetic`` or ``potential``). Model Setup ----------- Setting Initial Conditions `````````````````````````` For a unique solution of the equations of motion, initial displacements and velocities for all degrees of freedom must be specified: .. math:: \vec{u}(t=0) & = \vec{u}_0\\ \dot{\vec u}(t=0) & = \vec{v}_0 The solid mechanics model can be initialized as follows:: model.initFull() This function initializes the internal arrays and sets them to zero. Initial displacements and velocities that are not equal to zero can be prescribed by running a loop over the total number of nodes. Here, the initial displacement in :math:`x`-direction and the initial velocity in :math:`y`-direction for all nodes is set to :math:`0.1` and :math:`1`, respectively:: auto & disp = model.getDisplacement(); auto & velo = model.getVelocity(); for (UInt node = 0; node < mesh.getNbNodes(); ++node) { disp(node, 0) = 0.1; velo(node, 1) = 1.; } .. _sect-smm-boundary: Setting Boundary Conditions ``````````````````````````` This section explains how to impose Dirichlet or Neumann boundary conditions. A Dirichlet boundary condition specifies the values that the displacement needs to take for every point :math:`x` at the boundary (:math:`\Gamma_u`) of the problem domain (:numref:`fig-smm-boundaries`): .. math:: \vec{u} = \bar{\vec u} \quad \forall \vec{x}\in \Gamma_{u} A Neumann boundary condition imposes the value of the gradient of the solution at the boundary :math:`\Gamma_t` of the problem domain (:numref:`fig-smm-boundaries`): .. math:: \vec{t} = \mat{\sigma} \vec{n} = \bar{\vec t} \quad \forall \vec{x}\in \Gamma_{t} .. _fig-smm-boundaries: .. figure:: figures/problem_domain.svg :align: center Problem domain :math:`\Omega` with boundary in three dimensions. The Dirchelet and the Neumann regions of the boundary are denoted with :math:`\Gamma_u` and :math:`\Gamma_t`, respecitvely. Different ways of imposing these boundary conditions exist. A basic way is to loop over nodes or elements at the boundary and apply local values. A more advanced method consists of using the notion of the boundary of the mesh. In the following both ways are presented. Starting with the basic approach, as mentioned, the Dirichlet boundary conditions can be applied by looping over the nodes and assigning the required values. :numref:`fig-smm-dirichlet_bc` shows a beam with a fixed support on the left side. On the right end of the beam, a load is applied. At the fixed support, the displacement has a given value. For this example, the displacements in both the :math:`x` and the :math:`y`-direction are set to zero. Implementing this displacement boundary condition is similar to the implementation of initial displacement conditions described above. However, in order to impose a displacement boundary condition for all time steps, the corresponding nodes need to be marked as boundary nodes using the function ``blocked``. While, in order to impose a load on the right side, the nodes are not marked. The detail codes are shown as follows .. code-block:: c++ auto & blocked = model.getBlockedDOFs(); const auto & pos = mesh.getNodes(); UInt nb_nodes = mesh.getNbNodes(); for (UInt node = 0; node < nb_nodes; ++node) { if(Math::are_float_equal(pos(node, _x), 0)) { blocked(node, _x) = true; // block dof in x-direction blocked(node, _y) = true; // block dof in y-direction disp(node, _x) = 0.; // fixed displacement in x-direction disp(node, _y) = 0.; // fixed displacement in y-direction } else if (Math::are_float_equal(pos(node, _y), 0)) { blocked(node, _x) = false; // unblock dof in x-direction forces(node, _x) = 10.; // force in x-direction } } .. _fig-smm-dirichlet_bc: .. figure:: figures/dirichlet.svg :align: center Beam with fixed support and load. For the more advanced approach, one needs the notion of a boundary in the mesh. Therefore, the boundary should be created before boundary condition functors can be applied. Generally the boundary can be specified from the mesh file or the geometry. For the first case, the function ``createGroupsFromMeshData`` is called. This function can read any types of mesh data which are provided in the mesh file. If the mesh file is created with Gmsh, the function takes one input strings which is either ``tag_0``, ``tag_1`` or ``physical_names``. The first two tags are assigned by Gmsh to each element which shows the physical group that they belong to. In Gmsh, it is also possible to consider strings for different groups of elements. These elements can be separated by giving a string ``physical_names`` to the function ``createGroupsFromMeshData`` .. code-block:: c++ mesh.createGroupsFromMeshData("physical_names"). Boundary conditions support can also be created from the geometry by calling ``createBoundaryGroupFromGeometry``. This function gathers all the elements on the boundary of the geometry. To apply the required boundary conditions, the function :cpp:func:`applyBC ` needs to be called on a :cpp:class:`SolidMechanicsModel `. This function gets a Dirichlet or Neumann functor and a string which specifies the desired boundary on which the boundary conditions is to be applied. The functors specify the type of conditions to apply. Three built-in functors for Dirichlet exist: :cpp:class:`FlagOnly `, :cpp:class:`FixedValue ` and :cpp:class:`IncrementValue `. The functor ``FlagOnly`` is used if a point is fixed in a given direction. Therefore, the input parameter to this functor is only the fixed direction. The ``FixedValue`` functor is used when a displacement value is applied in a fixed direction. The ``IncrementValue`` applies an increment to the displacement in a given direction. The following code shows the utilization of three functors for the top, bottom and side surface of the mesh which were already defined in the Gmsh .. code-block:: c++ model.applyBC(BC::Dirichlet::FixedValue(13.0, _y), "Top"); model.applyBC(BC::Dirichlet::FlagOnly(_x), "Bottom"); model.applyBC(BC::Dirichlet::IncrementValue(13.0, _x), "Side"); To apply a Neumann boundary condition, the applied traction or stress should be specified before. In case of specifying the traction on the surface, the functor :cpp:class:`FromTraction ` of Neumann boundary conditions is called. Otherwise, the functor :cpp:class:`FromStress ` should be called which gets the stress tensor as an input parameter .. code-block:: c++ Vector surface_traction{0., 0., 1.}; auto surface_stress(3, 3) = Matrix::eye(3); model.applyBC(BC::Neumann::FromTraction(surface_traction), "Bottom"); model.applyBC(BC::Neumann::FromStress(surface_stress), "Top"); If the boundary conditions need to be removed during the simulation, a functor is called from the Neumann boundary condition to free those boundary conditions from the desired boundary .. code-block:: c++ model.applyBC(BC::Neumann::FreeBoundary(), "Side"); User specified functors can also be implemented. A full example for setting both initial and boundary conditions can be found in ``examples/boundary_conditions.cc``. The problem solved in this example is shown in :numref:`fig-smm-bc_and_ic`. It consists of a plate that is fixed with movable supports on the left and bottom side. On the right side, a traction, which increases linearly with the number of time steps, is applied. The initial displacement and velocity in :math:`x`-direction at all free nodes is zero and two respectively. .. _fig-smm-bc_and_ic: .. figure:: figures/bc_and_ic_example.svg :align: center :width: 75% Plate on movable supports. .. \begin{figure}[!htb] \centering \includegraphics[scale=0.8]{figures/bc_and_ic_example} \caption{Plate on movable supports.\label{fig-smm-bc_and_ic}} \end{figure} As it is mentioned in Section \ref{sect:common:groups}, node and element groups can be used to assign the boundary conditions. A generic example is given below with a Dirichlet boundary condition:: // create a node group NodeGroup & node_group = mesh.createNodeGroup("nodes_fix"); /* fill the node group with the nodes you want */ // create an element group using the existing node group mesh.createElementGroupFromNodeGroup("el_fix", "nodes_fix", spatial_dimension-1); // boundary condition can be applied using the element group name model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "el_fix"); Material Selector ````````````````` If the user wants to assign different materials to different finite elements groups in ``Akantu``, a material selector has to be used. By default, ``Akantu`` assigns the first valid material in the material file to all elements present in the model (regular continuum materials are assigned to the regular elements and cohesive materials are assigned to cohesive elements or element facets). To assign different materials to specific elements, mesh data information such as tag information or specified physical names can be used. :cpp:class:`MeshDataMaterialSelector ` class uses this information to assign different materials. With the proper physical name or tag name and index, different materials can be assigned as demonstrated in the examples below:: auto mat_selector = std::make_shared>("physical_names", model); model.setMaterialSelector(mat_selector); In this example the physical names specified in a GMSH geometry file will by used to match the material names in the input file. Another example would be to use the first (``tag_0``) or the second (``tag_1``) tag associated to each elements in the mesh:: auto mat_selector = std::make_shared>( "tag_1", model, first_index); model.setMaterialSelector(*mat_selector); where ``first_index`` (default is 1) is the value of ``tag_1`` that will be associated to the first material in the material input file. The following values of the tag will be associated with the following materials. There are four different material selectors pre-defined in ``Akantu``. :cpp:class:`MaterialSelector ` and :cpp:class:`DefaultMaterialSelector ` is used to assign a material to regular elements by default. For the regular elements, as in the example above, :cpp:class:`MeshDataMaterialSelector ` can be used to assign different materials to different elements. Apart from the ``Akantu``'s default material selectors, users can always develop their own classes in the main code to tackle various multi-material assignment situations. For cohesive material, ``Akantu`` has a pre-defined material selector to assign the first cohesive material by default to the cohesive elements which is called :cpp:class:`DefaultMaterialCohesiveSelector ` and it inherits its properties from :cpp:class:`DefaultMaterialSelector `. Multiple cohesive materials can be assigned using mesh data information (for more details, see :ref:`sect-smm-intrinsic-insertion`). Insertion of Cohesive Elements `````````````````````````````` Cohesive elements are currently compatible only with static simulation and dynamic simulation with an explicit time integration scheme (see section :ref:`ssect-smm-expl-time-integration`). They do not have to be inserted when the mesh is generated (intrinsic) but can be added during the simulation (extrinsic). At any time during the simulation, it is possible to access the following energies with the relative function: .. code-block:: c++ Real Ed = model.getEnergy("dissipated"); Real Er = model.getEnergy("reversible"); Real Ec = model.getEnergy("contact"); A new model have to be call in a very similar way that the solid mechanics model: .. code-block:: c++ SolidMechanicsModelCohesive model(mesh); model.initFull(_analysis_method = _explicit_lumped_mass, _is_extrinsic = true); Cohesive element insertion can be either realized at the beginning of the simulation or it can be carried out dynamically during the simulation. The first approach is called ``intrinsic``, the second one ``extrinsic``. When an element is present from the beginning, a bi-linear or exponential cohesive law should be used instead of a linear one. A bi-linear law works exactly like a linear one except for an additional parameter :math:`\delta_0` separating an initial linear elastic part from the linear irreversible one. For additional details concerning cohesive laws see Section~\ref{sec:cohesive-laws}. .. _fig-smm-coh-insertion: .. figure:: figures/insertion.svg :align: center Insertion of a cohesive element. Extrinsic cohesive elements are dynamically inserted between two standard elements when .. math:: \sigma_\mathrm{eff} > \sigma_\mathrm{c} \quad\text {with} \quad \sigma_\mathrm{eff} = \sqrt{\sigma_\mathrm{n} ^ 2 + \frac{\tau ^ 2} {\beta ^ 2 }} in which :math:`\sigma_\mathrm { n } ` is the tensile normal traction and $\tau$ the resulting tangential one( :numref:`fig-smm-coh-insertion`). Extrinsic approach '''''''''''''''''' During the simulation, stress has to be checked along each facet in order to insert cohesive elements where the stress criterion is reached. This check is performed by calling the method :cpp:func:`checkCohesiveStress `, as example before each step resolution: .. code-block:: c++ model.checkCohesiveStress(); model.solveStep(); The area where stresses are checked and cohesive elements inserted can be limited using the method :cpp:func:`setLimit ` on the :cpp:class:`CohesiveInserter ` during initialization. As example, to limit insertion in the range :math:`[-1.5, 1.5]` in the $x$ direction: .. code-block:: c++ auto & inserter = model.getElementInserter(); inserter.setLimit(_x, -1.5, 1.5); Additional restrictions with respect to :math:`_y` and :math:`_z` directions can be added as well. .. _sect-smm-intrinsic-insertion: Intrinsic approach '''''''''''''''''' Intrinsic cohesive elements are inserted in the mesh with the method :cpp:func:`initFull `. Similarly, the range of insertion can be limited with :cpp:func:`setLimit ` before the :cpp:func:`initFull ` call. In both cases extrinsic and intrinsic the insertion can be restricted to surfaces or element groups. To do so the list of groups should be specified in the input file. .. code-block:: model solid_mechanics_model_cohesive [ cohesive_inserter [ cohesive_surfaces = [surface1, surface2, ...] cohesive_zones = [group1, group2, ...] ] material cohesive_linear [ name = insertion beta = 1 G_c = 10 sigma_c = 1e6 ] ] Static Analysis --------------- The :cpp:class:`SolidMechanicsModel ` class can handle different analysis methods, the first one being presented is the static case. In this case, the equation to solve is .. math:: \mat{K} \vec{u} = \vec{f}_{\mathrm{ext}} :label: eqn-smm-static where :math:`\mat{K}` is the global stiffness matrix, :math:`\vec{u}` the displacement vector and :math:`\vec{f}_{\st{ext}}` the vector of external forces applied to the system. To solve such a problem, the static solver of the :cpp:class:`SolidMechanicsModel ` object is used. First, a model has to be created and initialized. To create the model, a mesh (which can be read from a file) is needed, as explained in Section~\ref{sect:common:mesh}. Once an instance of a :cpp:class:`SolidMechanicsModel ` is obtained, the easiest way to initialize it is to use the :cpp:func:`initFull ` method by giving the :cpp:class:`SolidMechanicsModelOptions `. These options specify the type of analysis to be performed and whether the materials should be initialized with :cpp:func:`initMaterials ` or not .. code-block:: c++ SolidMechanicsModel model(mesh); model.initFull(_analysis_method = _static); Here, a static analysis is chosen by passing the argument :cpp:enumerator:`_static ` to the method. By default, the Boolean for no initialization of the materials is set to false, so that they are initialized during the ``initFull``. The method ``initFull`` also initializes all appropriate vectors to zero. Once the model is created and initialized, the boundary conditions can be set as explained in Section :ref:`sect-smm-boundary`. Boundary conditions will prescribe the external forces for some free degrees of freedom :math:`\vec{f}_{\st{ext}}` and displacements for some others. At this point of the analysis, the function :cpp:func:`solveStep ` can be called .. code-block:: c++ auto & solver = model.getNonLinearSolver(); solver.set("max_iterations", 1); solver.set("threshold", 1e-4); solver.set("convergence_type", SolveConvergenceCriteria::_residual); model.solveStep(); This function is templated by the solving method and the convergence criterion and takes two arguments: the tolerance and the maximum number of iterations (100 by default), which are :math:`10^{-4}` and :math:`1` for this example. The modified Newton-Raphson method is chosen to solve the system. In this method, the equilibrium equation (:eq:`eqn-smm-static`) is modified in order to apply a Newton-Raphson convergence algorithm: .. math:: \mat{K}^{i+1}\delta\vec{u}^{i+1} &= \vec{r} \\ &= \vec{f}_{\st{ext}} -\vec{f}_{\st{int}}\\ &= \vec{f}_{\st{ext}} - \mat{K}^{i} \vec{u}^{i}\\ \vec{u}^{i+1} &= \vec{u}^{i} + \delta\vec{u}^{i+1}~, where :math:`\delta\vec{u}` is the increment of displacement to be added from one iteration to the other, and :math:`i` is the Newton-Raphson iteration counter. By invoking the ``solveStep`` method in the first step, the global stiffness matrix :math:`\mat{K}` from (:eq:`eqn-smm-static`) is automatically assembled. A Newton-Raphson iteration is subsequently started, :math:`\mat{K}` is updated according to the displacement computed at the previous iteration and one loops until the forces are balanced (:cpp:enumerator:`SolveConvergenceCriteria::_residual `), i.e. :math:`||\vec{r}|| <` :cpp:member:`threshold `. One can also iterate until the increment of displacement is zero (:cpp:enumerator:`SolveConvergenceCriteria::_solution `) which also means that the equilibrium is found. For a linear elastic problem, the solution is obtained in one iteration and therefore the maximum number of iterations can be set to one. But for a non-linear case, one needs to iterate as long as the norm of the residual exceeds the tolerance threshold and therefore the maximum number of iterations has to be higher, e.g. :math:`100` .. code-block:: c++ solver.set("max_iterations", 100); model.solveStep(); At the end of the analysis, the final solution is stored in the **displacement** vector. A full example of how to solve a static problem is presented in the code ``examples/static/static.cc``. This example is composed of a 2D plate of steel, blocked with rollers on the left and bottom sides as shown in :numref:`fig-smm-static`. The nodes from the right side of the sample are displaced by :math:`0.01\%` of the length of the plate. .. _fig-smm-static: .. figure:: figures/static.svg :align: center :width: 75% Numerical setup. The results of this analysis is depicted in :numref:`fig-smm-implicit:static_solution`. .. _fig-smm-implicit:static_solution: .. figure:: figures/static_analysis.png :align: center :width: 75% Solution of the static analysis. Left: the initial condition, right: the solution (deformation magnified 50 times). Dynamic Methods --------------- Different ways to solve the equations of motion are implemented in the solid mechanics model. The complete equations that should be solved are: .. math:: \mat{M}\ddot{\vec{u}} + \mat{C}\dot{\vec{u}} + \mat{K}\vec{u} = \vec{f}_{\mathrm{ext}} :label: eqn-equation-motion where :math:`\mat{M}`, :math:`\mat{C}` and :math:`\mat{K}` are the mass, damping and stiffness matrices, respectively. In the previous section, it has already been discussed how to solve this equation in the static case, where :math:`\ddot{\vec{u}} = \dot{\vec{u}} = 0`. Here the method to solve this equation in the general case will be presented. For this purpose, a time discretization has to be specified. The most common discretization method in solid mechanics is the Newmark-:math:`\beta` method, which is also the default in ``Akantu``. For the Newmark-:math:`\beta` method, (:eq:`eqn-equation-motion`) becomes a system of three equations (see :cite:`curnier92a,hughes-83a` for more details): .. math:: \begin{eqnarray} \mat{M} \ddot{\vec{u}}_{n+1} + \mat{C}\dot{\vec{u}}_{n+1} + \mat{K} \vec{u}_{n+1} &={\vec{f}_{\st{ext}}}_{\, n+1}\\ \vec{u}_{n+1} &= \vec{u}_{n} + \left(1 - \alpha\right) \Delta t \dot{\vec{u}}_{n} + \alpha \Delta t \dot{\vec{u}}_{n+1} + \left(\frac{1}{2} - \alpha\right) \Delta t^2 \ddot{\vec{u}}_{n}\\ \dot{\vec{u}}_{n+1} &= \dot{\vec{u}}_{n} + \left(1 - \beta\right) \Delta t \ddot{\vec{u}}_{n} + \beta \Delta t \ddot{\vec{u}}_{n+1} \end{eqnarray} :label: eqn-equation-motion-discret In these new equations, :math:`\ddot{\vec{u}}_{n}`, :math:`\dot{\vec{u}}_{n}` and :math:`\vec{u}_{n}` are the approximations of :math:`\ddot{\vec{u}}(t_n)`, :math:`\dot{\vec{u}}(t_n)` and :math:`\vec{u}(t_n)`. Equation~(:eq:`eqn-equation-motion-discret`) is the equation of motion discretized in space (finite-element discretization), and the equations above are discretized in both space and time (Newmark discretization). The :math:`\alpha` and :math:`\beta` parameters determine the stability and the accuracy of the algorithm. Classical values for :math:`\alpha` and :math:`\beta` are usually :math:`\beta = 1/2` for no numerical damping and :math:`0 < \alpha < 1/2`. .. csv-table:: :header: ":math:`\\alpha`", "Method (:math:`\\beta = 1/2`)", "Type" ":math:`0`", "central difference", "explicit" ":math:`\frac{1}{6}`", "Fox-Goodwin(royal road)", "implicit" ":math:`\frac{1}{3}`", "Linear acceleration", "implicit" ":math:`\frac{1}{2}`", "Average acceleration (trapeziodal rule)", "implicit" The solution of this system of equations, (:eq:`eqn-equation-motion-discret`) is split into a predictor and a corrector system of equations. Moreover, in the case of a non-linear equations, an iterative algorithm such as the Newton-Raphson method is applied. The system of equations can be written as: - *Predictor*: .. math:: \vec{u}_{n+1}^{0} &= \vec{u}_{n} + \Delta t \dot{\vec{u}}_{n} + \frac{\Delta t^2}{2} \ddot{\vec{u}}_{n} \\ \dot{\vec{u}}_{n+1}^{0} &= \dot{\vec{u}}_{n} + \Delta t \ddot{\vec{u}}_{n} \\ \ddot{\vec{u}}_{n+1}^{0} &= \ddot{\vec{u}}_{n} - *Solve*: .. math:: \left(c \mat{M} + d \mat{C} + e \mat{K}_{n+1}^i\right) \vec{w} &= {\vec{f}_{\st{ext}}}_{\,n+1} - {\vec{f}_{\st{int}}}_{\,n+1}^i - \mat{C} \dot{\vec{u}}_{n+1}^i - \mat{M} \ddot{\vec{u}}_{n+1}^i\\ &= \vec{r}_{n+1}^i - *Corrector*: .. math:: \ddot{\vec{u}}_{n+1}^{i+1} &= \ddot{\vec{u}}_{n+1}^{i} +c \vec{w} \\ \dot{\vec{u}}_{n+1}^{i+1} &= \dot{\vec{u}}_{n+1}^{i} + d\vec{w} \\ \vec{u}_{n+1}^{i+1} &= \vec{u}_{n+1}^{i} + e \vec{w} where :math:`i` is the Newton-Raphson iteration counter and :math:`c`, :math:`d` and :math:`e` are parameters depending on the method used to solve the equations .. csv-table:: :header: "", ":math:`\\vec{w}`", ":math:`e`", ":math:`d`", ":math:`c`" "in acceleration", ":math:`\delta\ddot{\vec{u}}`", ":math:`\alpha\beta\Delta t^2`", ":math:`\beta\Delta t`", ":math:`1`" "in velocity", ":math:`\delta\dot{\vec{u}}`", ":math:`\alpha\Delta t`", ":math:`1`", ":math:`\frac{1}{\beta\Delta t}`" "in displacement", ":math:`\delta\vec{u}`", ":math:`1`", ":math:`\frac{1}{\alpha\Delta t}`", ":math:`\frac{1}{\alpha\beta \Delta t^2}`" .. note:: If you want to use the implicit solver ``Akantu`` should be compiled at least with one sparse matrix solver such as `Mumps `_ :cite:`mumps`. Implicit Time Integration ````````````````````````` To solve a problem with an implicit time integration scheme, first a :cpp:class:`SolidMechanicsModel ` object has to be created and initialized. Then the initial and boundary conditions have to be set. Everything is similar to the example in the static case (Section~\ref{sect:smm:static}), however, in this case the implicit dynamic scheme is selected at the initialization of the model:: SolidMechanicsModel model(mesh); model.initFull(_analysis_method = _implicit_dynamic); Because a dynamic simulation is conducted, an integration time step :math:`\Delta t` has to be specified. In the case of implicit simulations, ``Akantu`` implements a trapezoidal rule by default. That is to say :math:`\alpha = 1/2` and :math:`\beta = 1/2` which is unconditionally stable. Therefore the value of the time step can be chosen arbitrarily within reason:: model.setTimeStep(time_step); Since the system has to be solved for a given amount of time steps, the method ``solveStep()``, (which has already been used in the static example in Section~\ref{sect:smm:static}), is called inside a time loop:: /// time loop Real time = 0.; auto & solver = model.getNonLinearSolver(); solver.set("max_iterations", 100); solver.set("threshold", 1e-12); solver.set("convergence_type", SolveConvergenceCriteria::_solution); for (UInt s = 1; time ` class. In the initialization, the explicit scheme is selected using the ``_explicit_lumped_mass`` constant:: SolidMechanicsModel model(mesh); model.initFull(_analysis_method = _explicit_lumped_mass); .. note:: Writing ``model.initFull()`` or ``model.initFull();`` is equivalent to use the ``_explicit_lumped_mass`` keyword, as this is the default case. The explicit time integration scheme implemented in ``Akantu`` uses a lumped mass matrix :math:`\mat{M}` (reducing the computational cost). This matrix is assembled by distributing the mass of each element onto its nodes. The resulting :math:`\mat{M}` is therefore a diagonal matrix stored in the **mass** vector of the model. The explicit integration scheme is conditionally stable. The time step has to be smaller than the stable time step which is obtained in Akantu as follows:: critical_time_step = model.getStableTimeStep(); The stable time step corresponds to the time the fastest wave (the compressive wave) needs to travel the characteristic length of the mesh: .. math:: \Delta t_{\st{crit}} = \frac{\Delta x}{c} where :math:`\Delta x` is a characteristic length (\eg the inradius in the case of linear triangle element) and :math:`c` is the celerity of the fastest wave in the material. It is generally the compressive wave of celerity :math:`c = \sqrt{\frac{2 \mu + \lambda}{\rho}}`, :math:`\mu` and :math:`\lambda` are the first and second Lame's coefficients and :math:`\rho` is the density. However, it is recommended to impose a time step that is smaller than the stable time step, for instance, by multiplying the stable time step by a safety factor smaller than one:: const Real safety_time_factor = 0.8; Real applied_time_step = critical_time_step * safety_time_factor; model.setTimeStep(applied_time_step); The initial displacement and velocity fields are, by default, equal to zero if not given specifically by the user (see \ref{sect:smm:initial_condition}). Like in implicit dynamics, a time loop is used in which the displacement, velocity and acceleration fields are updated at each time step. The values of these fields are obtained from the Newmark:math:`-\beta` equations with :math:`\beta=1/2` and :math:`\alpha=0`. In ``Akantu`` these computations at each time step are invoked by calling the function ``solveStep``:: for (UInt s = 1; (s-1)*applied_time_step < total_time; ++s) { model.solveStep(); } The method ``solveStep`` wraps the four following functions: - ``model.explicitPred()`` allows to compute the displacement field at :math:`t+1` and a part of the velocity field at :math:`t+1`, denoted by :math:`\vec{\dot{u}^{\st{p}}}_{n+1}`, which will be used later in the method ``model.explicitCorr()``. The equations are: .. math:: \vec{u}_{n+1} &= \vec{u}_{n} + \Delta t \vec{\dot{u}}_{n} + \frac{\Delta t^2}{2} \vec{\ddot{u}}_{n}\\ \vec{\dot{u}^{\st{p}}}_{n+1} &= \vec{\dot{u}}_{n} + \Delta t \vec{\ddot{u}}_{n} - ``model.updateResidual()`` and ``model.updateAcceleration()`` compute the acceleration increment :math:`\delta \vec{\ddot{u}}`: .. math:: \left(\mat{M} + \frac{1}{2} \Delta t \mat{C}\right) \delta \vec{\ddot{u}} = \vec{f_{\st{ext}}} - \vec{f}_{\st{int}\, n+1} - \mat{C} \vec{\dot{u}^{\st{p}}}_{n+1} - \mat{M} \vec{\ddot{u}}_{n} The internal force :math:`\vec{f}_{\st{int}\, n+1}` is computed from the displacement :math:`\vec{u}_{n+1}` based on the constitutive law. - ``model.explicitCorr()`` computes the velocity and acceleration fields at :math:`t+1`: .. math:: \vec{\dot{u}}_{n+1} &= \vec{\dot{u}^{\st{p}}}_{n+1} + \frac{\Delta t}{2} \delta \vec{\ddot{u}} \\ \vec{\ddot{u}}_{n+1} &= \vec{\ddot{u}}_{n} + \delta \vec{\ddot{u}} The use of an explicit time integration scheme is illustrated by the example: ``examples/explicit/explicit_dynamic.cc``. This example models the propagation of a wave in a steel beam. The beam and the applied displacement in the :math:`x` direction are shown in :numref:`fig-smm-explicit`. .. _fig-smm-explicit: .. figure:: figures/explicit.svg :align: center :width: 90% Numerical setup. The length and height of the beam are :math:`L={10}\textrm{m}` and :math:`h = {1}\textrm{m}`, respectively. The material is linear elastic, homogeneous and isotropic (density: :math:`7800\mathrm{kg/m}^3`, Young's modulus: :math:`210\mathrm{GPa}` and Poisson's ratio: :math:`0.3`). The imposed displacement follow a Gaussian function with a maximum amplitude of :math:`A = {0.01}\textrm{m}`. The potential, kinetic and total energies are computed. The safety factor is equal to :math:`0.8`. .. include:: ./constitutive-laws.rst .. include:: ./new-constitutive-laws.rst diff --git a/src/model/solid_mechanics/materials/material_plastic/material_plastic.hh b/src/model/solid_mechanics/materials/material_plastic/material_plastic.hh index de3cad9df..c40cd1142 100644 --- a/src/model/solid_mechanics/materials/material_plastic/material_plastic.hh +++ b/src/model/solid_mechanics/materials/material_plastic/material_plastic.hh @@ -1,128 +1,132 @@ /** * @file material_plastic.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Dec 07 2017 * * @brief Common interface for plastic materials * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_elastic.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_PLASTIC_HH_ #define AKANTU_MATERIAL_PLASTIC_HH_ namespace akantu { /** * Parent class for the plastic constitutive laws * parameters in the material files : * - h : Hardening parameter (default: 0) * - sigmay : Yield stress */ template class MaterialPlastic : public MaterialElastic { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialPlastic(SolidMechanicsModel & model, const ID & id = ""); MaterialPlastic(SolidMechanicsModel & model, UInt a_dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id = ""); protected: void initialize(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - /// get the energy specifying the type for the time step + /** + * @brief Return potential or plastic energy + * + * Plastic dissipated energy is integrated over time. + */ Real getEnergy(const std::string & type) override; - /// Compute the plastic energy + /// Update the plastic energy for the current timestep void updateEnergies(ElementType el_type) override; - /// Compute the true potential energy + /// Compute the true potential energy (w/ elastic strain) void computePotentialEnergy(ElementType el_type) override; protected: /// compute the stress and inelastic strain for the quadrature point inline void computeStressAndInelasticStrainOnQuad( const Matrix & grad_u, const Matrix & previous_grad_u, Matrix & sigma, const Matrix & previous_sigma, Matrix & inelastic_strain, const Matrix & previous_inelastic_strain, const Matrix & delta_inelastic_strain) const; inline void computeStressAndInelasticStrainOnQuad( const Matrix & delta_grad_u, Matrix & sigma, const Matrix & previous_sigma, Matrix & inelastic_strain, const Matrix & previous_inelastic_strain, const Matrix & delta_inelastic_strain) const; - /// get the plastic energy for the time step + /// Get the integrated plastic energy for the time step Real getPlasticEnergy(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Yield stresss Real sigma_y; /// hardening modulus Real h; /// isotropic hardening, r InternalField iso_hardening; /// inelastic strain arrays ordered by element types (inelastic deformation) InternalField inelastic_strain; /// Plastic energy InternalField plastic_energy; /// @todo : add a coefficient beta that will multiply the plastic energy /// increment /// to compute the energy converted to heat /// Plastic energy increment InternalField d_plastic_energy; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ } // namespace akantu #include "material_plastic_inline_impl.hh" #endif /* AKANTU_MATERIAL_PLASTIC_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 48714facf..9ba719aea 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,580 +1,590 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Model of Solid Mechanics * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" #include "non_local_manager_callback.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLID_MECHANICS_MODEL_HH_ #define AKANTU_SOLID_MECHANICS_MODEL_HH_ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition, public NonLocalManagerCallback, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; using MyFEEngineType = FEEngineTemplate; protected: using EventManager = EventHandlerManager; public: SolidMechanicsModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "solid_mechanics_model", ModelType model_type = ModelType::_solid_mechanics_model); ~SolidMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl( const ModelOptions & options = SolidMechanicsModelOptions()) override; public: /// initialize all internal arrays for materials virtual void initMaterials(); protected: /// initialize the model void initModel() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this adds f_{ext} or f_{int} to the residual void assembleResidual(const ID & residual_part) override; bool canSplitResidual() override { return true; } /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep(bool converged = true) override; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } protected: /// update the current position vector void updateCurrentPosition(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register an empty material of a given type Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param); /// reassigns materials depending on the material selector virtual void reassignMaterial(); /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials virtual void assignMaterialToElements(const ElementTypeMapArray * filter = nullptr); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(ElementType type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* NonLocalManager inherited members */ /* ------------------------------------------------------------------------ */ protected: void initializeNonLocal() override; void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; void computeNonLocalStresses(GhostType ghost_type) override; void insertIntegrationPointsInNeighborhoods(GhostType ghost_type) override; /// update the values of the non local internal void updateLocalInternal(ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) override; /// copy the results of the averaging in the materials void updateNonLocalInternal(ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: void splitElementByMaterial(const Array & elements, std::vector> & elements_per_mat) const; template void splitByMaterial(const Array & elements, Operation && op) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array & element_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array & /*unused*/, const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const ChangedElementsEvent & /*unused*/) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); //! decide wether a field is a material internal or not bool isInternal(const std::string & field_name, ElementKind element_kind); //! give the amount of data per element virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, ElementKind kind); //! flatten a given material internal field ElementTypeMapArray & flattenInternal(const std::string & field_name, ElementKind kind, GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(ElementKind kind); std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; void dump(const std::string & dumper_name) override; void dump(const std::string & dumper_name, UInt step) override; void dump(const std::string & dumper_name, Real time, UInt step) override; void dump() override; void dump(UInt step) override; void dump(Real time, UInt step) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Displacement, displacement); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR(Displacement, displacement); /// get the SolidMechanicsModel::previous_displacement array AKANTU_GET_MACRO_DEREF_PTR(PreviousDisplacement, previous_displacement); /// get the SolidMechanicsModel::current_position array const Array & getCurrentPosition(); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR(Increment, displacement_increment); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Increment, displacement_increment); /// get the lumped SolidMechanicsModel::mass array AKANTU_GET_MACRO_DEREF_PTR(Mass, mass); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Velocity, velocity); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR(Velocity, velocity); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Acceleration, acceleration); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR(Acceleration, acceleration); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(ExternalForce, external_force); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR(ExternalForce, external_force); /// get the SolidMechanicsModel::force array (external forces) [[deprecated("Use getExternalForce instead of this function")]] Array & getForce() { return getExternalForce(); } /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(InternalForce, internal_force); /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR(InternalForce, internal_force); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(BlockedDOFs, blocked_dofs); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR(BlockedDOFs, blocked_dofs); /// get an iterable on the materials inline decltype(auto) getMaterials(); /// get an iterable on the materials inline decltype(auto) getMaterials() const; /// get a particular material (by numerical material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by numerical material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); } /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); - /// get the energies + /** + * @brief Returns the total energy for a given energy type + * + * Energy types of SolidMechanicsModel expected as argument are: + * - `kinetic` + * - `external work` + * + * Other energy types are passed on to the materials. All materials should + * define a `potential` energy type. For additional energy types, see material + * documentation. + */ Real getEnergy(const std::string & energy_id); - /// compute the energy for an element + /// Compute energy for an element type and material index Real getEnergy(const std::string & energy_id, ElementType type, UInt index); - /// compute the energy for an element + /// Compute energy for an individual element Real getEnergy(const std::string & energy_id, const Element & element) { return getEnergy(energy_id, element.type, element.element); } - /// compute the energy for an element group + /// Compute energy for an element group Real getEnergy(const ID & energy_id, const ID & group_id); AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); AKANTU_GET_MACRO(MaterialLocalNumbering, material_local_numbering, const ElementTypeMapArray &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, // material_local_numbering, UInt); AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, *material_selector, MaterialSelector &); void setMaterialSelector(std::shared_ptr material_selector) { this->material_selector = std::move(material_selector); } /// Access the non_local_manager interface AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); /// get the FEEngine object to integrate or interpolate on the boundary FEEngine & getFEEngineBoundary(const ID & name = "") override; protected: /// compute the stable time step Real getStableTimeStep(GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// release version of the displacement array UInt displacement_release{0}; /// release version of the current_position array UInt current_position_release{0}; /// Check if materials need to recompute the mass array bool need_to_reassemble_lumped_mass{true}; /// Check if materials need to recompute the mass matrix bool need_to_reassemble_mass{true}; /// mapping between material name and material internal id std::map materials_names_to_id; protected: /// conversion coefficient form force/mass to acceleration Real f_m2a{1.0}; /// displacements array std::unique_ptr> displacement; /// displacements array at the previous time step (used in finite deformation) std::unique_ptr> previous_displacement; /// increment of displacement std::unique_ptr> displacement_increment; /// lumped mass array std::unique_ptr> mass; /// velocities array std::unique_ptr> velocity; /// accelerations array std::unique_ptr> acceleration; /// external forces array std::unique_ptr> external_force; /// internal forces array std::unique_ptr> internal_force; /// array specifing if a degree of freedom is blocked or not std::unique_ptr> blocked_dofs; /// array of current position used during update residual std::unique_ptr> current_position; /// Arrays containing the material index for each element ElementTypeMapArray material_index; /// Arrays containing the position in the element filter of the material /// (material's local numbering) ElementTypeMapArray material_local_numbering; /// list of used materials std::vector> materials; /// class defining of to choose a material std::shared_ptr material_selector; using flatten_internal_map = std::map, std::unique_ptr>>; /// map a registered internals to be flattened for dump purposes flatten_internal_map registered_internals; /// non local manager std::unique_ptr non_local_manager; /// tells if the material are instantiated bool are_materials_instantiated{false}; friend class Material; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { using FromStress = FromHigherDim; using FromTraction = FromSameDim; } // namespace Neumann } // namespace BC } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" #include "solid_mechanics_model_inline_impl.hh" #include "solid_mechanics_model_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* AKANTU_SOLID_MECHANICS_MODEL_HH_ */