diff --git a/doc/manual/manual-cohesive_elements_insertion.tex b/doc/manual/manual-cohesive_elements_insertion.tex index 9bc6e6f28..9ac6f3553 100644 --- a/doc/manual/manual-cohesive_elements_insertion.tex +++ b/doc/manual/manual-cohesive_elements_insertion.tex @@ -1,76 +1,82 @@ +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 +\code{DefaultMaterialCohesiveSelector} and it inherits its properties from +\code{DefaultMaterialSelector}. Multiple cohesive materials can be assigned +using mesh data information (for more details, see \ref{intrinsic_insertion}). + \subsection{Insertion of Cohesive Elements} Cohesive elements are currently compatible only with static simulation and dynamic simualtion with an explicit time integration scheme (see section~\ref{ssect:smm:expl-time-integr}). They do not have to be inserted when the mesh is generated but during the simulation. At any time during the simulation, it is possible to access the following energies with the relative function: \begin{cpp} Real Ed = model.getEnergy("dissipated"); Real Er = model.getEnergy("reversible"); Real Ec = model.getEnergy("contact"); \end{cpp} \subsubsection{Extrinsic approach} The dynamic insertion of extrinsic cohesive elements should be initialized the following way: \begin{cpp} SolidMechanicsModelCohesive model(mesh); model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); model.updateAutomaticInsertion(); \end{cpp} During the simulation, stress has to be checked along each facet in order to eventually insert cohesive elements. This check is performed by calling the method \code{checkCohesiveStress}, as example before each step resolution: \begin{cpp} model.checkCohesiveStress(); model.solveStep(); \end{cpp} The area where stress are checked and cohesive elements inserted can be limited using the method \code{limitInsertion} during initialization. As example, to limit insertion in the range $[-1.5, 1.5]$ in the $x$ direction: \begin{cpp} SolidMechanicsModelCohesive model(mesh); model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); model.limitInsertion(_x, -1.5, 1.5); model.updateAutomaticInsertion(); \end{cpp} Additional restrictions with respect to $y$ and $z$ directions can be added as well. \subsubsection{Intrinsic approach \label{intrinsic_insertion}} Intrinsic cohesive elements are inserted in the mesh with the method \code{insertIntrinsicElements}. Similarly, the range of insertion can me limited with \code{limitInsertion}. As example with a static simulation, \begin{cpp} SolidMechanicsModelCohesive model(mesh); model.initFull(SolidMechanicsModelCohesiveOptions(_static)); model.limitInsertion(_x, -1.5, 1.5); model.insertIntrinsicElements(); \end{cpp} Mesh data information becomes vital to the insertion of cohesive elements along surface with more sophisticated geometry or when multiple cohesive materials are wanted. At this purpose, cohesive elements can be easily inserted along a specific group of surface elements identified in a GMSH geometry file. This can be achieved, in the input file, by specifying the name of these physical groups in the \textit{mesh parameters} section as well as their corresponding cohesive materials. As example, with two physical surfaces named \textit{weak\_interface} and \textit{strong\_interface} defined in the GMSH geometry file: \begin{cpp} ... material %\emph{cohesive constitutive\_law}% [ name = weak_interface sigma_c = $value$ ... ] material %\emph{cohesive constitutive\_law}% [ name = strong_interface sigma_c = $value$ ... ] mesh parameters [ cohesive_surfaces = weak_interface,strong_interface ] \end{cpp} In this case, there is no need to call \code{insertIntrinsicElements} anymore since the insertion of cohesive elements along physical surfaces is performed automatically during \code{initFull} call. diff --git a/doc/manual/manual-io.tex b/doc/manual/manual-io.tex index 0142626bc..3ed3fa6eb 100644 --- a/doc/manual/manual-io.tex +++ b/doc/manual/manual-io.tex @@ -1,111 +1,154 @@ - \chapter{Input/Output}\index{I\/O} \section{Generic data} In this chapter, we address ways to get the internal data in human-readable formats. The models in \akantu handle data associated to the mesh, but this data can be split into several \code{Arrays}. For example, the data stored per element type in a \code{ElementTypeMapArray} is composed of as many \code{Array}s as types in the mesh. In order to get this data in a visualization software, the models contain a object to dump \code{VTK} files. These files can be visualized in software such as \code{ParaView}\cite{paraview}, \code{ViSit}\cite{visit} or \code{Mayavi}\cite{mayavi}. The internal dumper of the model can be configured to specify which data fields are to be written. This is done with the \code{addDumpField}\index{I\/O!addDumpField} method. By default all the files are generated in a folder called \code{paraview/} \begin{cpp} model.setBaseName("output"); // prefix for all generated files model.addDumpField("displacement"); model.addDumpField("stress"); ... model.dump() \end{cpp} The fields are dumped with the number of components of the memory. For example, in 2D, the memory has \code{Vector}s of 2 components, or the $2^{nd}$ order tensors with $2\times2$ components. This memory can be dealt with \code{addDumpFieldVector}\index{I\/O!addDumpFieldVector} which always dumps \code{Vector}s with 3 components or \code{addDumpFieldTensor}\index{I\/O!addDumpFieldTensor} which dumps $2^{nd}$ order tensors with $3\times3$ components respectively. The routines \code{addDumpFieldVector}\index{I\/O!addDumpFieldVector} and -\code{addDumpFieldTensor}\index{I\/O!addDumpFieldTensor} were introduced because of Paraview which mostly manipulate 3D data. +\code{addDumpFieldTensor}\index{I\/O!addDumpFieldTensor} were introduced because of \code{ParaView} which mostly manipulate 3D data. Those fields which are stored by quadrature point are modified to be seen in the \code{VTK} file as elemental data. To do this, the default is to average the values of all the quadrature points. The list of fields depends on the models (for \code{SolidMechanicsModel} see table~\ref{tab:io:smm_field_list}). \begin{table} \centering \begin{tabular}{llll} \toprule key & type & support \\ \midrule displacement & Vector & nodes \\ mass & Vector & nodes \\ velocity & Vector & nodes \\ acceleration & Vector & nodes \\ force & Vector & nodes \\ residual & Vector & nodes \\ increment & Vector & nodes \\ {blocked\_dofs} & Vector & nodes \\ partitions & Real & elements \\ material\_index & variable & elements\\ strain & Matrix & quadrature points \\ Green strain & Matrix & quadrature points \\ principal strain & Vector & quadrature points \\ principal Green strain & Vector & quadrature points \\ grad\_u & Matrix & quadrature points \\ stress & Matrix & quadrature points \\ Von Mises stress & Real & quadrature points \\ material\_index & variable & quadrature points \\ \bottomrule \end{tabular} \caption{List of dumpable fields for \code{SolidMechanicsModel}.} \label{tab:io:smm_field_list} \end{table} -The user can also register external fields which have the same mesh as the mesh from the model as support. To do this, an object of type \code{Field} has to be created.\index{I\/O!addDumpFieldExternal} - -\begin{itemize} -\item For nodal fields : -\begin{cpp} - Vector vect(nb_nodes, nb_component); - dumper::Field * field = new dumper::NodalField(vect); - model.addDumpFieldExternal("my_field", field); -\end{cpp} - -\item For elemental fields : -\begin{cpp} - ElementTypeMapArray arr; - dumper::Field * field = new dumper::ElementalField(arr, spatial_dimension); - model.addDumpFieldExternal("my_field", field); -\end{cpp} -\end{itemize} - \section{Cohesive elements' data} Cohesive elements and their relative data can be easily dumped thanks to a specific dumper contained in \code{SolidMechanicsModelCohesive}. In order to use it, one has just to add the string \code{"cohesive elements"} when calling each method already illustrated. Here is an example on how to dump displacement and damage: \begin{cpp} model.setBaseNameToDumper("cohesive elements", "cohesive_elements_output"); model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); model.addDumpFieldToDumper("cohesive elements", "damage"); ... model.dump("cohesive elements"); \end{cpp} -%%% Local Variables: + +\section{Advanced dumping} + +\subsection{Arbitrary fields} +In addition to the predetermined fields from the models and materials, the user +can add any data to a dumper as long as the support is the same. That is to say +data that have the size of the full mesh on if the dumper is dumping the mesh, +or of the size of an element group if it is a filtered dumper. + +For this the easiest is to use the ``external'' fields register functions\index{I\/O!addDumpFieldExternal} + +The simple case force nodal and elemental data are to pass directly the data +container itself if it as the good size. +\begin{itemize} +\item For nodal fields : +\begin{cpp} + Array nodal_data(nb_nodes, nb_component); + ... + model.addDumpFieldExternal("my_field", nodal_data); +\end{cpp} + +\item For elemental fields : +\begin{cpp} + ElementTypeMapArray elem_data; + ... + model.addDumpFieldExternal("my_field", elem_data); +\end{cpp} +\end{itemize} + +If some changes have to be applied on the data as for example a padding for +\code{ParaView} vectors, this can be done by using the +field interface. + +\begin{cpp} + model.addDumpFieldExternal(const std::string & field_name, + dumper::Field * field); +\end{cpp} + +An example of code presenting this interface is present in the +\shellcode{\examplesdir/io/dumper}. This interface is part of the +\code{Dumpable} class from which the \code{Mesh} inherits. + +\subsection{Creating a new dumper} + +You can also create you own dumpers, \akantu uses a third-party library in order +to write the output files, \code{IOHelper}. \akantu supports the \code{ParaView} +format and a Text format defined by \code{IOHelper}. + +This two files format are handled by the classes +\code{DumperParaview}\index{I\/O!DumperParaview} and +\code{DumperText}\index{I\/O!DumperText}. + +In order to use them you can instantiate on of this object in your code. This +dumper have a simple interface. You can register a mesh +\code{registerMesh}\index{I\/O!registerMesh}, +\code{registerFilteredMesh}\index{I\/O!registerFilteredMesh} or a field, +\code{registerField}\index{I\/O!registerField}. + +An example of code presenting this low level interface is present in the +\shellcode{\examplesdir/io/dumper}. The different types of \code{Field} that can +be created are present in the source folder \shellcode{src/io/dumper}. + +%%% Local Variables: %%% mode: latex %%% TeX-master: "manual" %%% End: diff --git a/doc/manual/manual-solidmechanicsmodel.tex b/doc/manual/manual-solidmechanicsmodel.tex index eb64db99b..883e80df6 100644 --- a/doc/manual/manual-solidmechanicsmodel.tex +++ b/doc/manual/manual-solidmechanicsmodel.tex @@ -1,1017 +1,1021 @@ \chapter{Solid Mechanics Model\index{SolidMechanicsModel}\label{sect:smm}} The solid mechanics model is a specific implementation of the \code{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 \code{FEEngine} object to compute the interpolation, gradient, integration and assembly operations. A \code{SolidMechanicsModel} object can simply be created like this: \begin{cpp} SolidMechanicsModel model(mesh); \end{cpp} where \code{mesh} is the mesh for which the equations are to be solved. A second parameter called \code{spatial\_dimension} can be added after \code{mesh} if the spatial dimension of the problem is different than that of the mesh. This model contains at least the following six \code{Arrays}: \begin{description} \item[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 \textbf{blocked\_dofs} value of a degree of freedom to \code{true}. A Neumann boundary condition can be applied by setting the \textbf{blocked\_dofs} value of a degree of freedom to \code{false}. The \textbf{displacement}, \textbf{velocity} and \textbf{acceleration} are computed for all degrees of freedom for which the \textbf{blocked\_dofs} value is set to \code{false}. For the remaining degrees of freedom, the imposed values (zero by default after initialization) are kept. \item[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 ($\vec{u}$ in the following). \item[velocity] contains the velocities of all degrees of freedom. As \textbf{displacement}, it contains computed or imposed velocities depending on the nature of the degrees of freedom ($\dot{\vec{u}}$ in the following). \item[acceleration] contains the accelerations of all degrees of freedom. As \textbf{displacement}, it contains computed or imposed accelerations depending on the nature of the degrees of freedom ($\ddot{\vec{u}}$ in the following). \item[force] contains the external forces applied on the nodes ($\vec{f}_{\st{ext}}$ in the following). \item[residual] contains the difference between external and internal forces. On blocked degrees of freedom, \textbf{residual} contains the support reactions. ($\vec{r}$ in the following). It should be mentioned that at equilibrium \textbf{residual} should be zero on free degrees of freedom. \end{description} Some examples to help to understand how to use this model will be presented in the next sections. \section{Model Setup} \subsection{Setting Initial Conditions \label{sect:smm:initial_condition}} For a unique solution of the equations of motion, initial displacements and velocities for all degrees of freedom must be specified: \begin{eqnarray} \vec{u}(t=0) = \vec{u}_0\\ \dot{\vec u}(t=0) =\vec{v}_0 \end{eqnarray} The solid mechanics model can be initialized as follows: \begin{cpp} model.initFull() \end{cpp} 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 $x$-direction and the initial velocity in $y$-direction for all nodes is set to $0.1$ and $1$, respectively. \begin{cpp} Array & disp = model.getDisplacement(); Array & velo = model.getVelocity(); for (UInt i = 0; i < mesh.getNbNodes(); ++i) { disp(i, 0) = 0.1; velo(i, 1) = 1.; } \end{cpp} \subsection{Setting Boundary Conditions\label{sect:smm:boundary}} 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 $x$ at the boundary ($\Gamma_u$) of the problem domain (Fig.~\ref{fig:smm:boundaries}): \begin{equation} \vec{u} = \bar{\vec u} \quad \forall \vec{x}\in \Gamma_{u} \end{equation} A Neumann boundary condition imposes the value of the gradient of the solution at the boundary $\Gamma_t$ of the problem domain (Fig.~\ref{fig:smm:boundaries}): \begin{equation} \vec{t} = \mat{\sigma} \vec{n} = \bar{\vec t} \quad \forall \vec{x}\in \Gamma_{t} \end{equation} \begin{figure} \centering \def\svgwidth{0.5\columnwidth} \input{figures/problemDomain.pdf_tex} \caption{Problem domain $\Omega$ with boundary in three dimensions. The Dirchelet and the Neumann regions of the boundary are denoted with $\Gamma_u$ and $\Gamma_t$, respecitvely.\label{fig:smm:boundaries}} \label{fig:problemDomain} \end{figure} 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. Figure~\ref{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 $x$ and the $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 as shown in the following code: \begin{cpp} Array & blocked = model.getBlockedDOFs(); const Array & pos = mesh.getNodes(); UInt nb_nodes = mesh.getNbNodes(); for (UInt i = 0; i < nb_nodes; ++i) { if(Math::are_float_equal(pos(i, 0), 0)) { blocked(i, 0) = true; //block displacement in x-direction blocked(i, 1) = true; //block displacement in y-direction disp(i, 0) = 0.; //fixed displacement in x-direction disp(i, 1)= 0.; //fixed displacement in y-direction } } \end{cpp} \begin{figure}[!htb] \centering \includegraphics[scale=0.4]{figures/dirichlet} \caption{Beam with fixed support.\label{fig:smm:dirichlet_bc}} \end{figure} 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 \code{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 \code{tag\_0}, \code{tag\_1} or \code{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 \code{physical\_names} to the function \code{createGroupsFromMeshData}: \begin{cpp} mesh.createGroupsFromMeshData("physical_names"). \end{cpp} Boundary conditions support can also be created from the geometry by calling \code{createBoundaryGroupFromGeometry}. This function gathers all the elements on the boundary of the geometry. To apply the required boundary conditions, the function \code{applyBC} needs to be called on a \code{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: \code{FlagOnly, FixedValue,} and \code{IncrementValue}. The functor \code{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 \code{FixedValue} functor is used when a displacement value is applied in a fixed direction. The \code{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 file: \begin{cpp} 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"); \end{cpp} 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 \code{FromTraction} of Neumann boundary conditions is called. Otherwise, the functor \code{FromStress} should be called which gets the stress tensor as an input parameter. \begin{cpp} Array surface_traction(3); surface_traction(0)=0.0; surface_traction(1)=0.0; surface_traction(2)=-1.0; Matrix surface_stress(3, 3, 0.0); surface_stress(0,0)=0.0; surface_stress(1,1)=0.0; surface_stress(2,2)=-1.0; model.applyBC(BC::Neumann::FromTraction(surface_traction), "Bottom"); model.applyBC(BC::Neumann::FromStress(surface_stress), "Top"); \end{cpp} 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. \begin{cpp} model.applyBC(BC::Neumann::FreeBoundary(), "Side"); \end{cpp} User specified functors can also be implemented. A full example for setting both initial and boundary conditions can be found in \shellcode{\examplesdir/boundary\_conditions.cc}. The problem solved in this example is shown in Fig.~\ref{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 $x$-direction at all free nodes is zero and two respectively. \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} \subsection{Material Selector\label{sect:smm:materialselector}} 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. \code{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. \begin{cpp} MeshDataMaterialSelector * mat_selector; mat_selector = new MeshDataMaterialSelector("physical_names", model); model.setMaterialSelector(*mat_selector); \end{cpp} 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 (\code{tag\_0}) or the second (\code{tag\_1}) tag associated to each elements in the mesh: +Another example would be to use the first (\code{tag\_0}) or the second +(\code{tag\_1}) tag associated to each elements in the mesh: \begin{cpp} MeshDataMaterialSelector * mat_selector; mat_selector = new MeshDataMaterialSelector("tag_1", model, first_index); model.setMaterialSelector(*mat_selector); \end{cpp} -where \code{first\_index} (default is 1) is the value of \code{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. +where \code{first\_index} (default is 1) is the value of \code{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. \code{MaterialSelector} and \code{DefaultMaterialSelector} is used to assign a material to regular elements by default. For the regular elements, as in the example above, \code{MeshDataMaterialSelector} can be used to assign different materials to different elements. -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 -\code{DefaultMaterialCohesiveSelector} and it inherits its properties from -\code{DefaultMaterialSelector}. Multiple cohesive materials can be assigned -using mesh data information (for more details, see \ref{intrinsic_insertion}). - 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. % An application of \code{DefaultMaterialCohesiveSelector} and usage in % a customly generated material selector class can be seen in % \shellcode{\examplesdir/cohesive\_element/cohesive\_extrinsic\_IG\_TG/cohesive\_extrinsic\_IG\_TG.cc}. \IfFileExists{manual-cohesive_elements_insertion.tex}{\input{manual-cohesive_elements_insertion}}{} \section{Static Analysis\label{sect:smm:static}} The \code{SolidMechanicsModel} class can handle different analysis methods, the first one being presented is the static case. In this case, the equation to solve is \begin{equation} \label{eqn:smm:static} \mat{K} \vec{u} = \vec{f}_{\st{ext}} \end{equation} where $\mat{K}$ is the global stiffness matrix, $\vec{u}$ the displacement vector and $\vec{f}_{\st{ext}}$ the vector of external forces applied to the system. To solve such a problem, the static solver of the \code{SolidMechanicsModel}\index{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 \code{SolidMechanicsModel} is obtained, the easiest way to initialize it is to use the \code{initFull}\index{SolidMechanicsModel!initFull} method by giving the \code{SolidMechanicsModelOptions}. These options specify the type of analysis to be performed and whether the materials should be initialized with \code{initMaterials} or not. \begin{cpp} SolidMechanicsModel model(mesh); model.initFull(SolidMechanicsModelOptions(_static, false)); \end{cpp} Here, a static analysis is chosen by passing the argument \code{\_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 \code{initFull}. The method \code{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 $\vec{f}_{\st{ext}}$ and displacements for some others. At this point of the analysis, the function \code{solveStep}\index{SolidMechanicsModel!solveStep} can be called: \begin{cpp} model.solveStep<_scm_newton_raphson_tangent_modified, _scc_residual>(1e-4, 1); \end{cpp} 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 $\num{1e-4}$ and $1$ for this example. The modified Newton-Raphson method is chosen to solve the system. In this method, the equilibrium equation (\ref{eqn:smm:static}) is modified in order to apply a Newton-Raphson convergence algorithm: \begin{align}\label{eqn:smm:static-newton-raphson} \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}~,\nonumber \end{align} where $\delta\vec{u}$ is the increment of displacement to be added from one iteration to the other, and $i$ is the Newton-Raphson iteration counter. By invoking the \code{solveStep} method in the first step, the global stiffness matrix $\mat{K}$ from Equation~(\ref{eqn:smm:static}) is automatically assembled. A Newton-Raphson iteration is subsequently started, $\mat{K}$ is updated according to the displacement computed at the previous iteration and one loops until the forces are balanced (\code{\_scc\_residual}), \ie $||\vec{r}|| < \mbox{\code{\_scc\_residual}}$. One can also iterate until the increment of displacement is zero (\code{\_scc\_increment}) 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. $100$: \begin{cpp} model.solveStep<_scm_newton_raphson_tangent_modified,_scc_residual>(1e-4, 100) \end{cpp} At the end of the analysis, the final solution is stored in the \textbf{displacement} vector. A full example of how to solve a static problem is presented in the code \code{\examplesdir/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 Figure \ref{fig:smm:static}. The nodes from the right side of the sample are displaced by $0.01\%$ of the length of the plate. \begin{figure}[!htb] \centering \includegraphics[scale=1.05]{figures/static} \caption{Numerical setup\label{fig:smm:static}} \end{figure} The results of this analysis is depicted in Figure~\ref{fig:smm:implicit:static_solution}. \begin{figure}[!htb] \centering \includegraphics[width=.7\linewidth]{figures/static_analysis} \caption{Solution of the static analysis. Left: the initial condition, right: the solution (deformation magnified 50 times)} \label{fig:smm:implicit:static_solution} \end{figure} \section{Dynamic Methods} \label{sect:smm: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: \begin{equation} \label{eqn:equation-motion} \mat{M}\ddot{\vec{u}} + \mat{C}\dot{\vec{u}} + \mat{K}\vec{u} = \vec{f}_{\st{ext}}~, \end{equation} where $\mat{M}$, $\mat{C}$ and $\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 $\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-$\beta$ method, which is also the default in \akantu. For the Newmark-$\beta$ method, (\ref{eqn:equation-motion}) becomes a system of three equations (see \cite{curnier92a} \cite{hughes-83a} for more details): \begin{align} \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} \label{eqn:equation-motion-discret} \\ \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} \label{eqn:finite-difference-1}\\ \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} \label{eqn:finite-difference-2} \end{align} In these new equations, $\ddot{\vec{u}}_{n}$, $\dot{\vec{u}}_{n}$ and $\vec{u}_{n}$ are the approximations of $\ddot{\vec{u}}(t_n)$, $\dot{\vec{u}}(t_n)$ and $\vec{u}(t_n)$. Equation~(\ref{eqn:equation-motion-discret}) is the equation of motion discretized in space (finite-element discretization), and equations (\ref{eqn:finite-difference-1}) and (\ref{eqn:finite-difference-2}) are discretized in both space and time (Newmark discretization). The $\alpha$ and $\beta$ parameters determine the stability and the accuracy of the algorithm. Classical values for $\alpha$ and $\beta$ are usually $\beta = 1/2$ for no numerical damping and $0 < \alpha < 1/2$. \begin{center} \begin{tabular}{cll} \toprule $\alpha$ & Method ($\beta = 1/2$) & Type\\ \midrule $0$ & central difference & explicit\\ $1/6$ & Fox-Goodwin (royal road) &implicit\\ $1/3$ & Linear acceleration &implicit\\ $1/2$ & Average acceleration (trapezoidal rule)& implicit\\ \bottomrule \end{tabular} \end{center} The solution of this system of equations, (\ref{eqn:equation-motion-discret})-(\ref{eqn:finite-difference-2}) 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: \begin{enumerate} \item \textit{Predictor:} \begin{align} \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} \end{align} \item \textit{Solve:} \begin{align} \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 \end{align} \item \textit{Corrector:} \begin{align} \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} \end{align} \end{enumerate} where $i$ is the Newton-Raphson iteration counter and $c$, $d$ and $e$ are parameters depending on the method used to solve the equations \begin{center} \begin{tabular}{lcccc} \toprule & $\vec{w}$ & $e$ & $d$ & $c$\\ \midrule in acceleration &$ \delta\ddot{\vec{u}}$ & $\alpha \beta\Delta t^2$ &$\beta \Delta t$ &$1$\\ in velocity & $ \delta\dot{\vec{u}}$& $\frac{1}{\beta} \Delta t$ & $1$ & $\alpha\Delta t$\\ in displacement &$\delta\vec{u}$ & $ 1$ & $\frac{1}{\alpha} \Delta t$ & $\frac{1}{\alpha \beta} \Delta t^2$\\ \bottomrule \end{tabular} \end{center} -%\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}.} +% \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}.} \subsection{Implicit Time Integration} To solve a problem with an implicit time integration scheme, first a \code{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. \begin{cpp} SolidMechanicsModel model(mesh); model.initFull(SolidMechanicsModelOptions(_implicit_dynamic)); /*Boundary conditions see Section~%\ref{sect:smm:boundary}% */ \end{cpp} Because a dynamic simulation is conducted, an integration time step $\Delta t$ has to be specified. In the case of implicit simulations, \akantu implements a trapezoidal rule by default. That is to say $\alpha = 1/2$ and $\beta = 1/2$ which is unconditionally stable. Therefore the value of the time step can be chosen arbitrarily within reason. \index{SolidMechanicsModel!setTimeStep} \begin{cpp} model.setTimeStep(time_step); \end{cpp} Since the system has to be solved for a given amount of time steps, the method \code{solveStep()}, (which has already been used in the static example in Section~\ref{sect:smm:static}), is called inside a time loop: \begin{cpp} /// time loop Real time = 0.; for (UInt s = 1; time (1e-12, 100); } \end{cpp} An example of solid mechanics with an implicit time integration scheme is presented in \shellcode{\examplesdir/implicit/implicit\_dynamic.cc}. This example consists of a 3D beam of $\SI{10}{\metre}\,\times\,\SI{1}{\metre}\,\times\,\SI{1}{\metre}$ blocked on one side and is on a roller on the other side. A constant force of \SI{5}{\kilo\newton} is applied in its middle. Figure~\ref{fig:smm:implicit:dynamic} presents the geometry of this case. The material used is a fictitious linear elastic material with a density of \SI{1000}{\kilo\gram\per\cubic\metre}, a Young's Modulus of \SI{120}{\mega\pascal} and Poisson's ratio of $0.3$. These values were chosen to simplify the analytical solution. An approximation of the dynamic response of the middle point of the beam is given by: \begin{equation} \label{eqn:smm:implicit} u\left(\frac{L}{2}, t\right) = \frac{1}{\pi^4} \left(1 - cos\left(\pi^2 t\right) + \frac{1}{81}\left(1 - cos\left(3^2 \pi^2 t\right)\right) + \frac{1}{625}\left(1 - cos\left(5^2 \pi^2 t\right)\right)\right) \end{equation} \begin{figure}[!htb] \centering \includegraphics[scale=.6]{figures/implicit_dynamic} \caption{Numerical setup} \label{fig:smm:implicit:dynamic} \end{figure} Figure \ref{fig:smm:implicit:dynamic_solution} presents the deformed beam at 3 different times during the simulation: time steps 0, 1000 and 2000. \begin{figure}[!htb] \centering \setlength{\unitlength}{0.1\textwidth} \begin{tikzpicture} \node[above right] (img) at (0,0) {\includegraphics[width=.6\linewidth]{figures/dynamic_analysis}}; \node[left] at (0pt,20pt) {$0$}; \node[left] at (0pt,60pt) {$1000$}; \node[left] at (0pt,100pt) {$2000$}; \end{tikzpicture} \caption{Deformed beam at 3 different times (displacement are magnified by a factor 10).} \label{fig:smm:implicit:dynamic_solution} \end{figure} \subsection{Explicit Time Integration} \label{ssect:smm:expl-time-integr} The explicit dynamic time integration scheme is based on the Newmark-$\beta$ scheme with $\alpha=0$ (see equations \ref{eqn:equation-motion-discret}-\ref{eqn:finite-difference-2}). In \akantu, $\beta$ is defaults to $\beta=1/2$, see section \ref{sect:smm:Dynamic_methods}. The initialization of the simulation is similar to the static and implicit dynamic version. The model is created from the \code{SolidMechanicsModel} class. In the initialization, the explicit scheme is selected using the \code{\_explicit\_lumped\_mass} constant. \begin{cpp} SolidMechanicsModel model(mesh); model.initFull(SolidMechanicsModelOptions(_explicit_lumped_mass)); \end{cpp} \index{SolidMechanicsModel!initFull} \note{Writing \code{model.initFull()} or \code{model.initFull(SolidMechanicsModelOptions());} is equivalent to use the \code{\_explicit\_lumped\_mass} keyword, as this is the default case.} The explicit time integration scheme implemented in \akantu uses a lumped mass matrix $\mat{M}$ (reducing the computational cost). This matrix is assembled by distributing the mass of each element onto its nodes. The resulting $\mat{M}$ is therefore a diagonal matrix stored in the \textbf{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: \begin{cpp} critical_time_step = model.getStableTimeStep(); \end{cpp} \index{SolidMechanicsModel!StableTimeStep} The stable time step corresponds to the time the fastest wave (the compressive wave) needs to travel the characteristic length of the mesh: \begin{equation} \label{eqn:smm:explicit:stabletime} \Delta t_{\st{crit}} = \frac{\Delta x}{c} \end{equation} -where $\Delta x$ is a characteristic length (\eg the inradius in the -case of linear triangle element) and $c$ is the celerity of the fastest wave in the material. It is generally the compressive -wave of celerity $c = \sqrt{\frac{2 \mu + \lambda}{\rho}}$, $\mu$ and $\lambda$ are the first and second Lame's coefficients and $\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. +where $\Delta x$ is a characteristic length (\eg the inradius in the case of +linear triangle element) and $c$ is the celerity of the fastest wave in the +material. It is generally the compressive wave of celerity +$c = \sqrt{\frac{2 \mu + \lambda}{\rho}}$, $\mu$ and $\lambda$ are the first and +second Lame's coefficients and $\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. \begin{cpp} const Real safety_time_factor = 0.8; Real applied_time_step = critical_time_step * safety_time_factor; model.setTimeStep(applied_time_step); \end{cpp} \index{SolidMechanicsModel!setTimeStep} 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$-\beta$ equations with $\beta=1/2$ and $\alpha=0$. In \akantu these computations at each time step are invoked by calling the function \code{solveStep}: \begin{cpp} for (UInt s = 1; (s-1)*applied_time_step < total_time; ++s) { model.solveStep(); } \end{cpp} \index{SolidMechanicsModel!solveStep} The method \code{solveStep} wraps the four following functions: \begin{itemize} \item \code{model.explicitPred()} allows to compute the displacement field at $t+1$ and a part of the velocity field at $t+1$, denoted by $\vec{\dot{u}^{\st{p}}}_{n+1}$, which will be used later in the method \code{model.explicitCorr()}. The equations are: \begin{align} \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} \label{eqn:smm:explicit:onehalfvelocity} \end{align} \item \code{model.updateResidual()} and \code{model.updateAcceleration()} compute the acceleration increment $\delta \vec{\ddot{u}}$: \begin{equation} \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}}_{n} - \mat{M} \vec{\ddot{u}}_{n} \end{equation} \note{The internal force $\vec{f}_{\st{int}\, n+1}$ is computed from the displacement $\vec{u}_{n+1}$ based on the constitutive law.} \item \code{model.explicitCorr()} computes the velocity and acceleration fields at $t+1$: \begin{align} \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}} \end{align} \end{itemize} The use of an explicit time integration scheme is illustrated by the example:\par \noindent \shellcode{\examplesdir/explicit/explicit\_dynamic.cc}\par \noindent This example models the propagation of a wave in a steel beam. The beam and the applied displacement in the $x$ direction are shown in Figure~\ref{fig:smm:explicit}. \begin{figure}[!htb] \centering \begin{tikzpicture} \coordinate (c) at (0,2); \draw[shift={(c)},thick, color=blue] plot [id=x, domain=-5:5, samples=50] ({\x, {(40 * sin(0.1*pi*3*\x) * exp(- (0.1*pi*3*\x)*(0.1*pi*3*\x) / 4))}}); \draw[shift={(c)},-latex] (-6,0) -- (6,0) node[right, below] {$x$}; \draw[shift={(c)},-latex] (0,-0.7) -- (0,1) node[right] {$u$}; \draw[shift={(c)}] (-0.1,0.6) node[left] {$A$}-- (1.5,0.6); \coordinate (l) at (0,0.6); \draw[shift={(0,-0.7)}] (-5, 0) -- (5,0) -- (5, 1) -- (-5, 1) -- cycle; \draw[shift={(l)}, latex-latex] (-5,0)-- (5,0) node [midway, above] {$L$}; \draw[shift={(l)}] (5,0.2)-- (5,-0.2); \draw[shift={(l)}] (-5,0.2)-- (-5,-0.2); \coordinate (h) at (5.3,-0.7); \draw[shift={(h)}, latex-latex] (0,0)-- (0,1) node [midway, right] {$h$}; \draw[shift={(h)}] (-0.2,1)-- (0.2,1); \draw[shift={(h)}] (-0.2,0)-- (0.2,0); \end{tikzpicture} \caption{Numerical setup \label{fig:smm:explicit}} \end{figure} The length and height of the beam are $L=\SI{10}{\metre}$ and $h = \SI{1}{\metre}$, respectively. The material is linear elastic, homogeneous and isotropic (density: \SI{7800}{\kilo\gram\per\cubic\metre}, Young's modulus: \SI{210}{\giga\pascal} and Poisson's ratio: $0.3$). The imposed displacement follow a Gaussian function with a maximum amplitude of $A = \SI{0.01}{\meter}$. The potential, kinetic and total energies are computed. The safety factor is equal to $0.8$. \input{manual-constitutive-laws} \section{Adding a New Constitutive Law}\index{Material!create a new material} There are several constitutive laws in \akantu as described in the previous Section~\ref{sect:smm:CL}. It is also possible to use a user-defined material for the simulation. These materials are referred to as local materials since they are local to the example of the user and not part of the \akantu library. To define a new local material, two files (\code {material\_XXX.hh} and \code{material\_XXX.cc}) have to be provided where \code{XXX} is the name of the new material. The header file \code {material\_XXX.hh} defines the interface of your custom material. Its implementation is provided in the \code{material\_XXX.cc}. The new law must inherit from the \code{Material} class or any other existing material class. It is therefore necessary to include the interface of the parent material in the header file of your local material and indicate the inheritance in the declaration of the class: \begin{cpp} /* ---------------------------------------------------------------------- */ #include "material.hh" /* ---------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_XXX_HH__ #define __AKANTU_MATERIAL_XXX_HH__ __BEGIN_AKANTU__ class MaterialXXX : public Material { /// declare here the interface of your material }; \end{cpp} In the header file the user also needs to declare all the members of the new material. These include the parameters that a read from the -material input file, as well as any other material parameters that will be computed during the simulation and internal variables.\\ -In the following the example of adding a new damage material will be presented. In this -case the parameters in the material will consist of the Young's modulus, the -Poisson coefficient, the resistance to damage and the damage threshold. The -material will then from these values compute its Lam\'{e} coefficients and its -bulk modulus. Furthermore, the user has to add a new internal variable -\code{damage} in order to store the amount of damage at each quadrature point in -each step of the simulation. For this specific material the member declaration -inside the class will look as follows: +material input file, as well as any other material parameters that will be +computed during the simulation and internal variables. + + +In the following the example of adding a new damage material will be +presented. In this case the parameters in the material will consist of the +Young's modulus, the Poisson coefficient, the resistance to damage and the +damage threshold. The material will then from these values compute its Lam\'{e} +coefficients and its bulk modulus. Furthermore, the user has to add a new +internal variable \code{damage} in order to store the amount of damage at each +quadrature point in each step of the simulation. For this specific material the +member declaration inside the class will look as follows: \begin{cpp} class LocalMaterialDamage : public Material { /// declare constructors/destructors here /// declare methods and accessors here /* -------------------------------------------------------------------- */ /* Class Members */ /* -------------------------------------------------------------------- */ AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); private: /// the young modulus Real E; /// Poisson coefficient Real nu; /// First Lame coefficient Real lambda; /// Second Lame coefficient (shear modulus) Real mu; /// resistance to damage Real Yd; /// damage threshold Real Sd; /// Bulk modulus Real kpa; /// damage internal variable InternalField damage; }; \end{cpp} In order to enable to print the material parameters at any point in the user's example file using the standard output stream by typing: \begin{cpp} for (UInt m = 0; m < model.getNbMaterials(); ++m) std::cout << model.getMaterial(m) << std::endl; \end{cpp} the standard output stream operator has to be redefined. This should be done at the end of the header file: \begin{cpp} class LocalMaterialDamage : public Material { /// declare here the interace of your material }: /* ---------------------------------------------------------------------- */ /* inline functions */ /* ---------------------------------------------------------------------- */ /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const LocalMaterialDamage & _this) { _this.printself(stream); return stream; } \end{cpp} However, the user still needs to register the material parameters that should be printed out. The registration is done during the call of the constructor. Like all definitions the implementation of the constructor has to be written in the \code{material\_XXX.cc} file. However, the declaration has to be provided in the \code{material\_XXX.hh} file: \begin{cpp} class LocalMaterialDamage : public Material { /* -------------------------------------------------------------------- */ /* Constructors/Destructors */ /* -------------------------------------------------------------------- */ public: LocalMaterialDamage(SolidMechanicsModel & model, const ID & id = ""); }; \end{cpp} The user can now define the implementation of the constructor in the \code{material\_XXX.cc} file: \begin{cpp} /* ---------------------------------------------------------------------- */ #include "local_material_damage.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* ---------------------------------------------------------------------- */ LocalMaterialDamage::LocalMaterialDamage(SolidMechanicsModel & model, const ID & id) : Material(model, id), damage("damage", *this) { AKANTU_DEBUG_IN(); this->registerParam("E", E, 0., _pat_parsable, "Young's modulus"); this->registerParam("nu", nu, 0.5, _pat_parsable, "Poisson's ratio"); this->registerParam("lambda", lambda, _pat_readable, "First Lame coefficient"); this->registerParam("mu", mu, _pat_readable, "Second Lame coefficient"); this->registerParam("kapa", kpa, _pat_readable, "Bulk coefficient"); this->registerParam("Yd", Yd, 50., _pat_parsmod); this->registerParam("Sd", Sd, 5000., _pat_parsmod); damage.initialize(1); AKANTU_DEBUG_OUT(); } \end{cpp} During the intializer list the reference to the model and the material id are assigned and the constructor of the internal field is called. Inside the scope of the constructor the internal values have to be initialized and the parameters, that should be printed out, are registered with the function: \code{registerParam}\index{Material!registerParam}: \begin{cpp} void registerParam(name of the parameter (key in the material file), member variable, default value (optional parameter), access permissions, description); \end{cpp} The available access permissions are as follows: \begin{itemize} \item \code{\_pat\_internal}: Parameter can only be output when the material is printed. \item \code{\_pat\_writable}: User can write into the parameter. The parameter is output when the material is printed. \item \code{\_pat\_readable}: User can read the parameter. The parameter is output when the material is printed. \item \code{\_pat\_modifiable}: Parameter is writable and readable. \item \code{\_pat\_parsable}: Parameter can be parsed, \textit{i.e.} read from the input file. \item \code{\_pat\_parsmod}: Parameter is modifiable and parsable. \end{itemize} In order to implement the new constitutive law the user needs to specify how the additional material parameters, that are not defined in the input material file, should be calculated. Furthermore, it has to be defined how stresses and the stable time step should be computed for the new local material. In the case of implicit simulations, in addition, the computation of the tangent stiffness needs to be defined. Therefore, the user needs to redefine the following functions of the parent material: \begin{cpp} void initMaterial(); // for explicit and implicit simulations void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); // for implicit simulations void computeTangentStiffness(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost); // for explicit and implicit simulations Real getStableTimeStep(Real h, const Element & element); \end{cpp} In the following a detailed description of these functions is provided: \begin{itemize} \item \code{initMaterial}:~ This method is called after the material file is fully read and the elements corresponding to each material are assigned. Some of the frequently used constant parameters are calculated in this method. For example, the Lam\'{e} constants of elastic materials can be considered as such parameters. \item \code{computeStress}:~ In this method, the stresses are computed based on the constitutive law as a function of the strains of the quadrature points. For example, the stresses for the elastic material are calculated based on the following formula: \begin{equation} \label{eqn:smm:constitutive_elastic} \mat{\sigma } =\lambda\mathrm{tr}(\mat{\varepsilon})\mat{I}+2 \mu \mat{\varepsilon} \end{equation} Therefore, this method contains a loop on all quadrature points assigned to the material using the two macros:\par \code{MATERIAL\_STRESS\_QUADRATURE\_POINT\_LOOP\_BEGIN}\par \code{MATERIAL\_STRESS\_QUADRATURE\_POINT\_LOOP\_END} \begin{cpp} MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(element_type); // sigma <- f(grad_u) MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; \end{cpp} \note{The strain vector in \akantu contains the values of $\nabla \vec{u}$, i.e. it is really the \emph{displacement gradient},} \item \code{computeTangentStiffness}:~ This method is called when the tangent to the stress-strain curve is desired (see Fig \ref {fig:smm:AL:K}). For example, it is called in the implicit solver when the stiffness matrix for the regular elements is assembled based on the following formula: \begin{equation} \label{eqn:smm:constitutive_elasc} \mat{K } =\int{\mat{B^T}\mat{D(\varepsilon)}\mat{B}} \end{equation} Therefore, in this method, the \code{tangent} matrix (\mat{D}) is computed for a given strain. \note{ The \code{tangent} matrix is a $4^{th}$ order tensor which is stored as a matrix in Voigt notation.} \begin{figure}[!htb] \begin{center} \includegraphics[width=0.4\textwidth,keepaspectratio=true]{figures/tangent.pdf} \caption{Tangent to the stress-strain curve.} \label{fig:smm:AL:K} \end{center} \end{figure} \item \code{getCelerity}:~The stability criterion of the explicit integration scheme depend on the fastest wave celerity~\eqref{eqn:smm:explicit:stabletime}. This celerity depend on the material, and therefore the value of this velocity should be defined in this method for each new material. By default, the fastest wave speed is the compressive wave whose celerity can be defined in~\code{getPushWaveSpeed}. \end{itemize} Once the declaration and implementation of the new material has been completed, this material can be used in the user's example by including the header file: \begin{cpp} #include "material_XXX.hh" \end{cpp} For existing materials, as mentioned in Section~\ref{sect:smm:CL}, by default, the materials are initialized inside the method \code{initFull}. If a local material should be used instead, the initialization of the material has to be postponed until the local material is registered in the model. Therefore, the model is initialized with the boolean for skipping the material initialization equal to true: \begin{cpp} /// model initialization model.initFull(SolidMechanicsModelOptions(_explicit_lumped_mass,true)); \end{cpp} Once the model has been initialized, the local material needs to be registered in the model: \begin{cpp} model.registerNewCustomMaterials("name_of_local_material"); \end{cpp} Only at this point the material can be initialized: \begin{cpp} model.initMaterials(); \end{cpp} A full example for adding a new damage law can be found in \shellcode{\examplesdir/new\_material}. %%% Local Variables: %%% mode: latex %%% TeX-master: "manual" %%% End: diff --git a/examples/io/dumper/dumpable_interface.cc b/examples/io/dumper/dumpable_interface.cc index a88f56a37..456bdd55d 100644 --- a/examples/io/dumper/dumpable_interface.cc +++ b/examples/io/dumper/dumpable_interface.cc @@ -1,179 +1,185 @@ /** * @file example_dumpable_interface.cc * @author Fabian Barras * @date Thu Jul 2 14:34:41 2015 * * @brief Example of dumper::Dumpable interface. * * @section LICENSE * * Copyright (©) 2010-2011 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 "mesh.hh" #include "element_group.hh" #include "dumper_paraview.hh" #include "dumpable_inline_impl.hh" #include "group_manager_inline_impl.cc" /* -------------------------------------------------------------------------- */ #include "locomotive_tools.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { /* In this example, we present dumper::Dumpable which is an interface for other classes who want to dump themselves. Several classes of Akantu inheritate from Dumpable (Model, Mesh, ...). In this example we reproduce the same tasks as example_dumper_low_level.cc using this time Dumpable interface inherted by Mesh, NodeGroup and ElementGroup. It is then advised to read first example_dumper_low_level.cc. */ initialize(argc, argv); // To start let us load the swiss train mesh and its mesh data information. UInt spatial_dimension = 2; Mesh mesh(spatial_dimension); mesh.read("swiss_train.msh"); mesh.createGroupsFromMeshData("physical_names"); /* swiss_train.msh has the following physical groups that can be viewed with GMSH: "$MeshFormat 2.2 0 8 $EndMeshFormat $PhysicalNames 6 2 1 "red" 2 2 "white" 2 3 "lwheel_1" 2 4 "lwheel_2" 2 5 "rwheel_2" 2 6 "rwheel_1" $EndPhysicalNames ..." */ // Grouping nodes and elements belonging to train wheels (=four mesh data). ElementGroup & wheels_elements = mesh.createElementGroup("wheels", spatial_dimension); wheels_elements.append(mesh.getElementGroup("lwheel_1")); wheels_elements.append(mesh.getElementGroup("lwheel_2")); wheels_elements.append(mesh.getElementGroup("rwheel_1")); wheels_elements.append(mesh.getElementGroup("rwheel_2")); const Array & lnode_1 = (mesh.getElementGroup("lwheel_1")).getNodes(); const Array & lnode_2 = (mesh.getElementGroup("lwheel_2")).getNodes(); const Array & rnode_1 = (mesh.getElementGroup("rwheel_1")).getNodes(); const Array & rnode_2 = (mesh.getElementGroup("rwheel_2")).getNodes(); Array & node = mesh.getNodes(); UInt nb_nodes = mesh.getNbNodes(); // This time a 2D Array is created and a padding size of 3 is passed to // NodalField in order to warp train deformation on Paraview. Array displacement(nb_nodes, spatial_dimension); // Create an ElementTypeMapArray for the colour ElementTypeMapArray colour("colour"); mesh.initElementTypeMapArray(colour, 1, spatial_dimension, false, _ek_regular, true); /* ------------------------------------------------------------------------ */ /* Creating dumpers */ /* ------------------------------------------------------------------------ */ // Create dumper for the complete mesh and register it as default dumper. DumperParaview dumper("train", "./paraview/dumpable", false); mesh.registerExternalDumper(dumper, "train", true); mesh.addDumpMesh(mesh); // The dumper for the filtered mesh can be directly taken from the // ElementGroup and then registered as "wheels_elements" dumper. DumperIOHelper & wheels = mesh.getGroupDumper("paraview_wheels", "wheels"); mesh.registerExternalDumper(wheels, "wheels"); mesh.setDirectoryToDumper("wheels", "./paraview/dumpable"); // Arrays and ElementTypeMapArrays can be added as external fields directly mesh.addDumpFieldExternal("displacement", displacement); - mesh.addDumpFieldExternal("color", colour); + + ElementTypeMapArrayFilter filtered_colour( + colour, wheels_elements.getElements()); + + dumper::Field * colour_field_wheel = + new dumper::ElementalField(filtered_colour); + mesh.addDumpFieldExternal("color", colour_field_wheel); mesh.addDumpFieldExternalToDumper("wheels", "displacement", displacement); mesh.addDumpFieldExternalToDumper("wheels", "colour", colour); // For some specific cases the Fields should be created, as when you want to // pad an array dumper::Field * displacement_vector_field = mesh.createNodalField(&displacement, "all", 3); mesh.addDumpFieldExternal("displacement_as_paraview_vector", displacement_vector_field); mesh.addDumpFieldExternalToDumper("wheels", "displacement_as_paraview_vector", displacement_vector_field); /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ // Fill the ElementTypeMapArray colour. fillColour(mesh, colour); /// Apply displacement and wheels rotation. Real tot_displacement = 50.; Real radius = 1.; UInt nb_steps = 100; Real theta = tot_displacement / radius; Vector l_center(spatial_dimension); Vector r_center(spatial_dimension); for (UInt i = 0; i < spatial_dimension; ++i) { l_center(i) = node(14, i); r_center(i) = node(2, i); } for (UInt i = 0; i < nb_steps; ++i) { displacement.clear(); Real step_ratio = Real(i) / Real(nb_steps); Real angle = step_ratio * theta; applyRotation(l_center, angle, node, displacement, lnode_1); applyRotation(l_center, angle, node, displacement, lnode_2); applyRotation(r_center, angle, node, displacement, rnode_1); applyRotation(r_center, angle, node, displacement, rnode_2); for (UInt j = 0; j < nb_nodes; ++j) { displacement(j, _x) += step_ratio * tot_displacement; } /// Dump call is finally made through Dumpable interface. mesh.dump(); mesh.dump("wheels"); } finalize(); return 0; }