diff --git a/doc/manual/manual-structuralmechanicsmodel.tex b/doc/manual/manual-structuralmechanicsmodel.tex index 340012f5f..9796d047b 100644 --- a/doc/manual/manual-structuralmechanicsmodel.tex +++ b/doc/manual/manual-structuralmechanicsmodel.tex @@ -1,212 +1,224 @@ \chapter{Structural Mechanics model} Static structural mechanics problems can be handled using the \code{StructuralMechanicsModel}. So far, \akantu\ provides 2D and 3D Bernoulli beam elements \cite{frey2009}. Just as for the \code{SolidMechanicsModel}, the model is created for a given \code{Mesh}. The model will create its own \code{FEM} object to compute the interpolation, gradient, integration and assembly operations. The \code{StructuralMechanicsModel} constructor is used like \begin{cpp} StructuralMechanicsModel model(mesh, spatial_dimension); \end{cpp} where \code{mesh} is a \code{Mesh} object defining the structure for which the equations of statics are to be solved, and \code{spatial\_dimension} is the dimensionality of the problem. If \code{spatial\_dimension} is omitted, the problem is assumed to have the same dimensionality as the one specified by the mesh. \note[\ 1]{Dynamic computations are not supported to date.} \note[\ 2]{Structural meshes are created and loaded as described in Section~\ref{sect:common:mesh} with \code{MeshIOMSHStruct} instead of \code{MeshIOMSH}.} \vspace{1cm} This model contains at least the following \code{Arrays}: \begin{description} -\item[boundary] contains a boolean value for each degree of freedom +\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{boundary} value of a degree - of freedom to \code{true}. A Neumann boundary condition can be applied - by setting the corresponding \textbf{boundary} value to \code{false}. - The \textbf{displacement} is computed for all degrees of freedom for which the \textbf{boundary} value is set to \code{false}. For the remaining degrees of freedom, the imposed values (zero by default after initialization) are kept. + condition can be prescribed by setting the \textbf{blocked\_dofs} value of a +degree + of freedom to \code{true}. The \textbf{displacement} is 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\_rotation] contains the generalized displacements (displacements and rotations) 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[force\_moment] contains the generalized external forces (forces and moments) applied to the nodes ($\vec{f_{\st{ext}}}$ in the following). \item[residual] contains the difference between external and internal forces and moments. On the 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 the free degrees of freedom. \end{description} An example to help understand how to use this model will be presented in the next section. \section{Model setup} \label{sec:structMechMod:setup} \subsection{Initialization} The easiest way to initialize the structural mechanics model is: -%\begin{cpp} -% model.initModel(); -% -% StructuralMaterial mat1; -% mat1.E=2.05e11; -% mat1.I=0.00128; -% mat1.A=0.01; // for example -% -% model.addMaterial(mat1); ASK NICO ABOUT THIS CRAP -% -% model.initArrays(); -% model.initImplicitSolver(); -%\end{cpp} -%The method \code{initModel} computes the shape functions, \code{addMaterial} sets the material parameters to be used, \code{initArrays} initializes all the internal arrays mentioned before and \code{initImplicitSolver} creates the stiffness matrix. \begin{cpp} model.initFull(); \end{cpp} The method \code{initFull} computes the shape functions, initializes the internal vectors mentioned above and allocates the memory for the stiffness matrix. Material properties are defined using the \code{StructuralMaterial} structure described in Table~\ref{tab:structMechMod:strucMaterial}. Such a definition could for instance look like \begin{cpp} - StructuralMaterial mat; + StructuralMaterial mat1; mat.E=3e10; mat.I=0.0025; mat.A=0.01; \end{cpp} \begin{table}[htb] \centering \begin{tabular}{c|c} field & description \\\hline\hline \code{E} & Young's modulus \\\hline \code{A} & Cross section area \\\hline \code{I} & Second cross sectional moment of inertia (for 2D elements) \\\hline \code{Iy} & \code{I} around beam $y$--axis (for 3D elements) \\\hline \code{Iz} & \code{I} around beam $z$--axis (for 3D elements) - \\\hline \code{GJ} & Polar moment of inertia of beam cross section (for 3D elements) + \\\hline \code{GJ} & Polar moment of inertia of beam cross section +(for 3D elements) \end{tabular} \caption{Material properties for structural elements as defined by the structure \code{StructuralMaterial}.} \label{tab:structMechMod:strucMaterial} \end{table} Materials can be added to the model's \code{element\_material} vector using \begin{cpp} - model.addMaterial(material); + model.addMaterial(mat1); +\end{cpp} + +They are successively numbered and then assigned to specific elements. +\begin{cpp} +for (UInt i = 0; i < nb_element_mat_1; ++i) { + model.getElementMaterial(_bernoulli_beam_2)(i,0) = 1; + } \end{cpp} -They are used just like for the \code{SolidMechanicsModel} -see Section~\ref{sect:smm:CL}. \subsection{Setting boundary conditions}\label{sect:structMechMod:boundary} -Both Dirichlet and Neumann type boundary conditions are applied to nodes the -same exact way as for \code{SolidMechanicsModel}, see -Section~\ref{sect:smm:boundary}. The method \code{computeForcesFromFunction} -can still be used to apply Neumann boundary conditions. +As explained before, the Dirichlet boundary conditions are applied through the +array \textbf{blocked\_dofs}. Two options exist to define Neumann conditions. +If a +nodal force is applied, it has to be directly set in the array +\textbf{force\_momentum}. For loads distributed along the beam length, the +method +\code{computeForcesFromFunction} integrates them into nodal forces. +The method takes as input a function describing the distribution of loads +along the beam and a \code{BoundaryFunctionType} specifing if the function is +expressed in the local coordinates (\code{\_bft\_traction\_local}) +or in the global system of coordinates (\code{\_bft\_traction}). +\begin{cpp} + static void lin_load(double * position, double * load, + Real * normal, UInt surface_id){ + memset(load,0,sizeof(Real)*3); + load[1] = position[0]*position[0]-250; +} +int main(int argc, char *argv[]){ +... +model.computeForcesFromFunction<_bernoulli_beam_2>(lin_load, + _bft_traction_local); +...} +\end{cpp} + \section{Static analysis\label{sect:structMechMod:static}} The \code{StructuralMechanicsModel} class can perform static analyses of structures. In this case, the equation to solve is the same as for the \code{SolidMechanicsModel} used for static analyses \begin{equation}\label{eqn:structMechMod:static} \mat{K} \vec{u} = \vec{f_{\st{ext}}}~, \end{equation} where $\mat{K}$ is the global stiffness matrix, $\vec{u}$ the generalized displacement vector and $\vec{f_{\st{ext}}}$ the vector of generalized external forces applied to the system. - To solve such a problem, the static solver of the \code{StructuralMechanicsModel}\index{StructuralMechanicsModel} object is used. First a -model has to be created and initialized. To create the model, a mesh is required. -Once an instance of a \code{StructuralMechanicsModel} is obtained, the easiest way to -initialize it is to use the \code{initFull}\index{StructuralMechanicsModel!initFull} -function. +model has to be created and initialized. \begin{cpp} StructuralMechanicsModel model(mesh); model.initFull(); \end{cpp} \begin{itemize} \item \code{model.initFull} initializes all internal vectors to zero. \end{itemize} Once the model is created and initialized, the boundary conditions can be set as explained in Section \ref{sect:structMechMod:boundary}. Boundary conditions will prescribe the external forces or moments for the free degrees of freedom $\vec{f_{\st{ext}}}$ and displacements or rotations for the others. To completely define the system represented by equation (\ref{eqn:structMechMod:static}), the global stiffness matrix $\mat{K}$ must be assembled. \index{StructuralMechanicsModel!assembleStiffnessMatrix} \begin{cpp} model.assembleStiffnessMatrix(); \end{cpp} The computation of the static equilibrium is performed using the same Newton-Raphson algorithm as described in Section~\ref{sect:smm:static}. \note{To date, \code{StructuralMechanicsModel} handles only constitutively and geometrically linear problems, the algorithm is therefore guaranteed -to converge in one iteration.} +to converge in two iterations.} \begin{cpp} model.updateResidual(); model.solve(); \end{cpp} \index{StructuralMechanicsModel!updateResidual} \index{StructuralMechanicsModel!solve} \begin{itemize} \item \code{model.updateResidual} assembles the internal forces and removes them from the external forces. \item \code{model.solve} solves the equations - (\ref{eqn:smm:static-newton-raphson}). The \textbf{increment} vector of the + (\ref{eqn:structMechMod:static}). The \textbf{increment} vector of the model will contain the new increment of displacements, and the - \textbf{displacement} vector is also updated to the new displacements. + \textbf{displacement\_rotation} vector is also updated to the new +displacements. \end{itemize} At the end of the analysis, the final solution is stored in the -\textbf{displacement} vector. A full example of how to solve a structural +\textbf{displacement\_rotation} vector. A full example of how to solve a +structural mechanics problem is presented in the code -\shellcode{\examplesdir/WHERE-DO-EXAMPLES-GO?}\todo{set the right example file}. This example is composed of a 2D -beam, clamped at the left end and supported by two rollers as shown in Figure -\ref{fig:structMechMod:exem1_1}. The problem is defined by the applied load -$q=\unit{6}{\kilo\newton\per\metre}$, moment $\bar{M} = -\unit{3.6}{\kilo\newton\metre}$, moments of inertia $I_1 = -\unit{250\,000}{\power{\centi\metre}{4}}$ and $I_2 = -\unit{128\,000}{\power{\centi\metre}{4}}$ and lengths $L_1 = \unit{10}{\metre}$ -and $L_2 = \unit{8}{\metre}$. The resulting rotations at node two and three are -$ \varphi_2 = 0.001\,167\ \mbox{and}\ \varphi_3 = -0.000\,771.$ - +\shellcode{\examplesdir/structural\_mechanics/bernoulli\_beam\_2\_example.cc}. +This example is composed of a 2D beam, clamped at the left end and supported +by two rollers as shown in Figure \ref{fig:structMechMod:exem1_1}. The problem + is defined by the applied load $q=\unit{6}{\kilo\newton\per\metre}$, + moment $\bar{M} = \unit{3.6}{\kilo\newton\metre}$, moments + of inertia $I_1 = \unit{250\,000}{\power{\centi\metre}{4}}$ + and $I_2 = \unit{128\,000}{\power{\centi\metre}{4}}$ and + lengths $L_1 = \unit{10}{\metre}$ and $L_2 = \unit{8}{\metre}$. The resulting +rotations at node two and three are $ \varphi_2 = 0.001\,167\ \mbox{and}\ +\varphi_3 = -0.000\,771.$ \begin{figure}[htb] \centering \includegraphics[scale=1.1]{figures/beam_example} \caption{2D beam example} \label{fig:structMechMod:exem1_1} \end{figure} %%% Local Variables: %%% mode: latex %%% TeX-master: "manual" %%% End: diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc index 1c698d0ce..49d59e092 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.cc +++ b/src/model/structural_mechanics/structural_mechanics_model.cc @@ -1,673 +1,673 @@ /** * @file structural_mechanics_model.cc * @author Fabian Barras * @date Thu May 5 15:52:38 2011 * * @brief * * @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 "structural_mechanics_model.hh" #include "aka_math.hh" #include "integration_scheme_2nd_order.hh" #include "static_communicator.hh" #include "sparse_matrix.hh" #include "solver.hh" #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif #ifdef AKANTU_USE_IOHELPER # include "dumper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, dim, id, memory_id), Dumpable(), stress("stress", id, memory_id), element_material("element_material", id, memory_id), set_ID("beam sets", id, memory_id), stiffness_matrix(NULL), jacobian_matrix(NULL), solver(NULL), rotation_matrix("rotation_matices", id, memory_id) { AKANTU_DEBUG_IN(); registerFEEngineObject("StructuralMechanicsFEEngine", mesh, spatial_dimension); this->displacement_rotation = NULL; this->force_momentum = NULL; this->residual = NULL; this->blocked_dofs = NULL; this->increment = NULL; if(spatial_dimension == 2) nb_degree_of_freedom = 3; else if (spatial_dimension == 3) nb_degree_of_freedom = 6; else { AKANTU_DEBUG_TO_IMPLEMENT(); } this->registerDumper("paraview_all", id, true); this->addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::~StructuralMechanicsModel() { AKANTU_DEBUG_IN(); delete solver; delete stiffness_matrix; delete jacobian_matrix; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFull(std::string material) { initModel(); initArrays(); initSolver(); displacement_rotation->clear(); force_momentum ->clear(); residual ->clear(); blocked_dofs ->clear(); increment ->clear(); Mesh::type_iterator it = getFEEngine().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural); Mesh::type_iterator end = getFEEngine().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural); for (; it != end; ++it) { computeRotationMatrix(*it); } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = getFEEngine().getMesh().getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":blocked_dofs"; std::stringstream sstr_incr; sstr_incr << id << ":increment"; displacement_rotation = &(alloc(sstr_disp.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE)); force_momentum = &(alloc(sstr_forc.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE)); blocked_dofs = &(alloc(sstr_boun.str(), nb_nodes, nb_degree_of_freedom, false)); increment = &(alloc(sstr_incr.str(), nb_nodes, nb_degree_of_freedom, REAL_INIT_VALUE)); Mesh::type_iterator it = getFEEngine().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural); Mesh::type_iterator end = getFEEngine().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural); for (; it != end; ++it) { UInt nb_element = getFEEngine().getMesh().getNbElement(*it); UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(*it); element_material.alloc(nb_element, 1, *it, _not_ghost); set_ID.alloc(nb_element, 1, *it, _not_ghost); UInt size = getTangentStiffnessVoigtSize(*it); stress.alloc(nb_element * nb_quadrature_points, size , *it, _not_ghost); } dof_synchronizer = new DOFSynchronizer(getFEEngine().getMesh(), nb_degree_of_freedom); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initModel() { getFEEngine().initShapeFunctions(_not_ghost); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) { AKANTU_DEBUG_IN(); const Mesh & mesh = getFEEngine().getMesh(); #if !defined(AKANTU_USE_MUMPS) // or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake AKANTU_DEBUG_ERROR("You should at least activate one solver."); #else UInt nb_global_node = mesh.getNbGlobalNodes(); std::stringstream sstr; sstr << id << ":jacobian_matrix"; jacobian_matrix = new SparseMatrix(nb_global_node * nb_degree_of_freedom, _symmetric, nb_degree_of_freedom, sstr.str(), memory_id); jacobian_matrix->buildProfile(mesh, *dof_synchronizer); std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solv; sstr_solv << id << ":solver"; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS solver->initialize(options); #endif //AKANTU_HAS_SOLVER AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt StructuralMechanicsModel::getTangentStiffnessVoigtSize(const ElementType & type) { UInt size; #define GET_TANGENT_STIFFNESS_VOIGT_SIZE(type) \ size = getTangentStiffnessVoigtSize(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(GET_TANGENT_STIFFNESS_VOIGT_SIZE); #undef GET_TANGENT_STIFFNESS_VOIGT_SIZE return size; } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); stiffness_matrix->clear(); Mesh::type_iterator it = getFEEngine().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural); Mesh::type_iterator end = getFEEngine().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural); for (; it != end; ++it) { ElementType type = *it; #define ASSEMBLE_STIFFNESS_MATRIX(type) \ assembleStiffnessMatrix(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(ASSEMBLE_STIFFNESS_MATRIX); #undef ASSEMBLE_STIFFNESS_MATRIX } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::computeRotationMatrix<_bernoulli_beam_2>(Array & rotations){ ElementType type = _bernoulli_beam_2; Mesh & mesh = getFEEngine().getMesh(); UInt nb_element = mesh.getNbElement(type); Array::iterator< Vector > connec_it = mesh.getConnectivity(type).begin(2); Array::vector_iterator nodes_it = mesh.getNodes().begin(spatial_dimension); Array::matrix_iterator R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); for (UInt e = 0; e < nb_element; ++e, ++R_it, ++connec_it) { Matrix & R = *R_it; Vector & connec = *connec_it; Vector x2; x2 = nodes_it[connec(1)]; // X2 Vector x1; x1 = nodes_it[connec(0)]; // X1 Real le = x1.distance(x2); Real c = (x2(0) - x1(0)) / le; Real s = (x2(1) - x1(1)) / le; /// Definition of the rotation matrix R(0,0) = c; R(0,1) = s; R(0,2) = 0.; R(1,0) = -s; R(1,1) = c; R(1,2) = 0.; R(2,0) = 0.; R(2,1) = 0.; R(2,2) = 1.; } } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::computeRotationMatrix<_bernoulli_beam_3>(Array & rotations){ ElementType type = _bernoulli_beam_3; Mesh & mesh = getFEEngine().getMesh(); UInt nb_element = mesh.getNbElement(type); Array::vector_iterator n_it = mesh.getNormals(type).begin(spatial_dimension); Array::iterator< Vector > connec_it = mesh.getConnectivity(type).begin(2); Array::vector_iterator nodes_it = mesh.getNodes().begin(spatial_dimension); Matrix Pe (spatial_dimension, spatial_dimension); Matrix Pg (spatial_dimension, spatial_dimension); Matrix inv_Pg(spatial_dimension, spatial_dimension); Vector x_n(spatial_dimension); // x vect n Array::matrix_iterator R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); for (UInt e=0 ; e < nb_element; ++e, ++n_it, ++connec_it, ++R_it) { Vector & n = *n_it; Matrix & R = *R_it; Vector & connec = *connec_it; Vector x; x = nodes_it[connec(1)]; // X2 Vector y; y = nodes_it[connec(0)]; // X1 Real l = x.distance(y); x -= y; // X2 - X1 x_n.crossProduct(x, n); Pe.eye(); Pe(0, 0) *= l; Pe(1, 1) *= -l; Pg(0,0) = x(0); Pg(0,1) = x_n(0); Pg(0,2) = n(0); Pg(1,0) = x(1); Pg(1,1) = x_n(1); Pg(1,2) = n(1); Pg(2,0) = x(2); Pg(2,1) = x_n(2); Pg(2,2) = n(2); inv_Pg.inverse(Pg); Pe *= inv_Pg; for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { R(i, j) = Pe(i, j); R(i + spatial_dimension,j + spatial_dimension) = Pe(i, j); } } } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeRotationMatrix(const ElementType & type) { Mesh & mesh = getFEEngine().getMesh(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); if(!rotation_matrix.exists(type)) { rotation_matrix.alloc(nb_element, nb_degree_of_freedom*nb_nodes_per_element * nb_degree_of_freedom*nb_nodes_per_element, type); } else { rotation_matrix(type).resize(nb_element); } rotation_matrix(type).clear(); Arrayrotations(nb_element, nb_degree_of_freedom * nb_degree_of_freedom); rotations.clear(); #define COMPUTE_ROTATION_MATRIX(type) \ computeRotationMatrix(rotations); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_ROTATION_MATRIX); #undef COMPUTE_ROTATION_MATRIX Array::matrix_iterator R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); Array::matrix_iterator T_it = rotation_matrix(type).begin(nb_degree_of_freedom*nb_nodes_per_element, nb_degree_of_freedom*nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++R_it, ++T_it) { Matrix & T = *T_it; Matrix & R = *R_it; T.clear(); for (UInt k = 0; k < nb_nodes_per_element; ++k){ for (UInt i = 0; i < nb_degree_of_freedom; ++i) for (UInt j = 0; j < nb_degree_of_freedom; ++j) T(k*nb_degree_of_freedom + i, k*nb_degree_of_freedom + j) = R(i, j); } } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeStresses() { AKANTU_DEBUG_IN(); Mesh::type_iterator it = getFEEngine().getMesh().firstType(spatial_dimension, _not_ghost, _ek_structural); Mesh::type_iterator end = getFEEngine().getMesh().lastType(spatial_dimension, _not_ghost, _ek_structural); for (; it != end; ++it) { ElementType type = *it; #define COMPUTE_STRESS_ON_QUAD(type) \ computeStressOnQuad(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_STRESS_ON_QUAD); #undef COMPUTE_STRESS_ON_QUAD } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::updateResidual() { AKANTU_DEBUG_IN(); residual->copy(*force_momentum); Array ku(*displacement_rotation, true); ku *= *stiffness_matrix; *residual -= ku; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::solve() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving an implicit step."); UInt nb_nodes = displacement_rotation->getSize(); /// todo residual = force - Kxr * d_bloq jacobian_matrix->copyContent(*stiffness_matrix); jacobian_matrix->applyBoundary(*blocked_dofs); increment->clear(); solver->setRHS(*residual); solver->solve(*increment); Real * increment_val = increment->storage(); Real * displacement_val = displacement_rotation->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*blocked_dofs_val)) { *displacement_val += *increment_val; } else { *increment_val = 0.0; } displacement_val++; blocked_dofs_val++; increment_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool StructuralMechanicsModel::testConvergenceIncrement(Real tolerance) { Real error; bool tmp = testConvergenceIncrement(tolerance, error); AKANTU_DEBUG_INFO("Norm of increment : " << error); return tmp; } /* -------------------------------------------------------------------------- */ bool StructuralMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error) { AKANTU_DEBUG_IN(); Mesh & mesh= getFEEngine().getMesh(); UInt nb_nodes = displacement_rotation->getSize(); UInt nb_degree_of_freedom = displacement_rotation->getNbComponent(); Real norm = 0; Real * increment_val = increment->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt d = 0; d < nb_degree_of_freedom; ++d) { if(!(*blocked_dofs_val) && is_local_node) { norm += *increment_val * *increment_val; } blocked_dofs_val++; increment_val++; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); error = sqrt(norm); AKANTU_DEBUG_ASSERT(!isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (error < tolerance); } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::computeTangentModuli<_bernoulli_beam_2>(Array & tangent_moduli) { UInt nb_element = getFEEngine().getMesh().getNbElement(_bernoulli_beam_2); UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_bernoulli_beam_2); UInt tangent_size = 2; Array::matrix_iterator D_it = tangent_moduli.begin(tangent_size, tangent_size); Array & el_mat = element_material(_bernoulli_beam_2, _not_ghost); for (UInt e = 0; e < nb_element; ++e) { UInt mat = el_mat(e); Real E = materials[mat].E; Real A = materials[mat].A; Real I = materials[mat].I; for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it) { Matrix & D = *D_it; D(0,0) = E * A; D(1,1) = E * I; } } } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::computeTangentModuli<_bernoulli_beam_3>(Array & tangent_moduli) { UInt nb_element = getFEEngine().getMesh().getNbElement(_bernoulli_beam_3); UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_bernoulli_beam_3); UInt tangent_size = 4; Array::matrix_iterator D_it = tangent_moduli.begin(tangent_size, tangent_size); for (UInt e = 0; e < nb_element; ++e) { UInt mat = element_material(_bernoulli_beam_3, _not_ghost)(e); Real E = materials[mat].E; Real A = materials[mat].A; Real Iz = materials[mat].Iz; Real Iy = materials[mat].Iy; Real GJ = materials[mat].GJ; for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it) { Matrix & D = *D_it; D(0,0) = E * A; D(1,1) = E * Iz; D(2,2) = E * Iy; D(3,3) = GJ; } } } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::transferBMatrixToSymVoigtBMatrix<_bernoulli_beam_2>(Array & b, bool local) { UInt nb_element = getFEEngine().getMesh().getNbElement(_bernoulli_beam_2); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_bernoulli_beam_2); UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_bernoulli_beam_2); MyFEEngineType & fem = getFEEngineClass(); Array::const_vector_iterator shape_Np = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 0).begin(nb_nodes_per_element); Array::const_vector_iterator shape_Mpp = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 1).begin(nb_nodes_per_element); Array::const_vector_iterator shape_Lpp = fem.getShapesDerivatives(_bernoulli_beam_2, _not_ghost, 2).begin(nb_nodes_per_element); UInt tangent_size = getTangentStiffnessVoigtSize<_bernoulli_beam_2>(); UInt bt_d_b_size = nb_nodes_per_element * nb_degree_of_freedom; b.clear(); Array::matrix_iterator B_it = b.begin(tangent_size, bt_d_b_size); for (UInt e = 0; e < nb_element; ++e) { for (UInt q = 0; q < nb_quadrature_points; ++q) { Matrix & B = *B_it; const Vector & Np = *shape_Np; const Vector & Lpp = *shape_Lpp; const Vector & Mpp = *shape_Mpp; B(0,0) = Np(0); B(0,3) = Np(1); B(1,1) = Mpp(0); B(1,2) = Lpp(0); B(1,4) = Mpp(1); B(1,5) = Lpp(1); ++B_it; ++shape_Np; ++shape_Mpp; ++shape_Lpp; } // ++R_it; } } /* -------------------------------------------------------------------------- */ template<> void StructuralMechanicsModel::transferBMatrixToSymVoigtBMatrix<_bernoulli_beam_3>(Array & b, bool local) { MyFEEngineType & fem = getFEEngineClass(); UInt nb_element = getFEEngine().getMesh().getNbElement(_bernoulli_beam_3); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(_bernoulli_beam_3); UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(_bernoulli_beam_3); Array::const_vector_iterator shape_Np = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 0).begin(nb_nodes_per_element); Array::const_vector_iterator shape_Mpp = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 1).begin(nb_nodes_per_element); Array::const_vector_iterator shape_Lpp = fem.getShapesDerivatives(_bernoulli_beam_3, _not_ghost, 2).begin(nb_nodes_per_element); UInt tangent_size = getTangentStiffnessVoigtSize<_bernoulli_beam_3>(); UInt bt_d_b_size = nb_nodes_per_element * nb_degree_of_freedom; b.clear(); Array::matrix_iterator B_it = b.begin(tangent_size, bt_d_b_size); for (UInt e = 0; e < nb_element; ++e) { for (UInt q = 0; q < nb_quadrature_points; ++q) { Matrix & B = *B_it; const Vector & Np = *shape_Np; const Vector & Lpp = *shape_Lpp; const Vector & Mpp = *shape_Mpp; B(0,0) = Np(0); B(0,6) = Np(1); B(1,1) = Mpp(0); B(1,5) = Lpp(0); B(1,7) = Mpp(1); B(1,11) = Lpp(1); B(2,2) = Mpp(0); B(2,4) = -Lpp(0); B(2,8) = Mpp(1); B(2,10) = -Lpp(1); B(3,3) = Np(0); B(3,9) = Np(1); ++B_it; ++shape_Np; ++shape_Mpp; ++shape_Lpp; } } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(dumper_name, id, field, type, n, stride) \ internalAddDumpFieldToDumper(dumper_name, id, \ new DumperIOHelper::NodalField(*field, n, stride)) UInt n; if(spatial_dimension == 2) { n = 2; } else n = 3; if(field_id == "displacement" ) { ADD_FIELD(dumper_name, "displacement", displacement_rotation, Real, - n, nb_degree_of_freedom - n); } + n, 0); } else if(field_id == "rotation") { ADD_FIELD(dumper_name, "rotation", displacement_rotation, Real, nb_degree_of_freedom - n, n); } else if(field_id == "force" ) { ADD_FIELD(dumper_name, "force", force_momentum, Real, n, 0); } else if(field_id == "momentum") { ADD_FIELD(dumper_name, "momentum", force_momentum, Real, nb_degree_of_freedom - n, n); } else if(field_id == "residual") { ADD_FIELD(dumper_name, "residual", residual, Real, nb_degree_of_freedom, 0); } else if(field_id == "blocked_dofs") { ADD_FIELD(dumper_name, "blocked_dofs", blocked_dofs, bool, nb_degree_of_freedom, 0); } else if(field_id == "element_index_by_material") { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::ElementalField(element_material, spatial_dimension, _not_ghost, _ek_regular)); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(dumper_name, id, field, type, n, stride) \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field, n, stride); \ f->setPadding(3); \ internalAddDumpFieldToDumper(dumper_name, id, f) UInt n; if(spatial_dimension == 2) { n = 2; } else n = 3; if(field_id == "displacement" ) { ADD_FIELD(dumper_name, "displacement", displacement_rotation, Real, n, 0); } else if(field_id == "force" ) { ADD_FIELD(dumper_name, "force", force_momentum, Real, n, 0); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::addDumpFieldTensorToDumper(__attribute__((unused)) const std::string & dumper_name, __attribute__((unused)) const std::string & field_id) { } __END_AKANTU__ diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_complicated.cc b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_complicated.cc index 96ad57da8..5c49ad7dd 100644 --- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_complicated.cc +++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_complicated.cc @@ -1,165 +1,161 @@ /** * @file test_structural_mechanics_model_bernoulli_beam_2_complicated.cc * @author Fabian Barras * @date Wed Jun 1 16:06:45 2011 * * @brief A very complicated structure * * @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 #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh_struct.hh" #include "structural_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #define TYPE _bernoulli_beam_2 using namespace akantu; //Linear load function static void lin_load(double * position, double * load, __attribute__ ((unused)) Real * normal, __attribute__ ((unused)) UInt surface_id){ memset(load,0,sizeof(Real)*3); if(position[1]>=0.-Math::getTolerance()) { if ((position[0]<=10.)){ load[1]= -100; } else if (position[0]<=20.){ load[1]= -70; } } } int main(int argc, char *argv[]){ initialize(argc, argv); Mesh beams(2); debug::setDebugLevel(dblWarning); /* -------------------------------------------------------------------------- */ // Defining the mesh akantu::MeshIOMSHStruct mesh_io; mesh_io.read("complicated.msh", beams); - mesh_io.write("complicated_tata.msh", beams); /* -------------------------------------------------------------------------- */ // Defining the material const akantu::ElementType type = akantu::_bernoulli_beam_2; akantu::StructuralMechanicsModel model(beams); StructuralMaterial mat1; mat1.E=3e10; mat1.I=0.0025; mat1.A=0.01; model.addMaterial(mat1); StructuralMaterial mat2 ; mat2.E=3e10; mat2.I=0.003125; mat2.A=0.01; model.addMaterial(mat2); /* -------------------------------------------------------------------------- */ // Defining the forces model.initFull(); UInt nb_element = beams.getNbElement(type); for (unsigned int i = 0; i < nb_element; ++i) { model.getElementMaterial(type)(i,0) = beams.getData("tag_0", type)(i,0) - 1; } Array & forces = model.getForce(); Array & displacement = model.getDisplacement(); Array & boundary = model.getBlockedDOFs(); - // const Array & N_M = model.getStress(_bernoulli_beam_2); - - // Array & element_material = model.getElementMaterial(_bernoulli_beam_2); forces.clear(); displacement.clear(); model.computeForcesFromFunction<_bernoulli_beam_2>(lin_load, akantu::_bft_traction); /* -------------------------------------------------------------------------- */ // Defining the boundary conditions boundary(0,0) = true; boundary(0,1) = true; boundary(3,0) = true; boundary(3,1) = true; boundary(4,0) = true; boundary(4,1) = true; boundary(4,2) = true; boundary(5,0) = true; boundary(5,1) = true; boundary(5,2) = true; boundary(2,1) = true; boundary(2,0) = true; boundary(1,1) = true; boundary(1,0) = true; /* -------------------------------------------------------------------------- */ // Solve Real error; model.assembleStiffnessMatrix(); model.getStiffnessMatrix().saveMatrix("Kb.mtx"); UInt count = 0; model.addDumpFieldVector("displacement"); model.addDumpField("rotation"); model.addDumpField("force"); model.addDumpField("momentum"); do { if(count != 0) std::cerr << count << " - " << error << std::endl; model.updateResidual(); model.solve(); count++; } while (!model.testConvergenceIncrement(1e-10, error) && count < 10); std::cerr << count << " - " << error << std::endl; /* -------------------------------------------------------------------------- */ // Post-Processing model.computeStresses(); model.getStiffnessMatrix().saveMatrix("Ka.mtx"); std::cout<< " x1 = " << displacement(1,2) << std::endl; std::cout<< " x2 = " << displacement(2,2) << std::endl; model.dump(); } diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_exemple_1_1.cc b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_exemple_1_1.cc index 078b4eda5..34ab4d224 100644 --- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_exemple_1_1.cc +++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_exemple_1_1.cc @@ -1,183 +1,183 @@ /** * @file test_structural_mechanics_model_bernoulli_beam_2_exemple_1_1.cc * @author Fabian Barras * @date Tue May 31 19:10:26 2011 * * @brief Computation of the analytical exemple 1.1 in the TGC vol 6 * * @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 #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "structural_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #define TYPE _bernoulli_beam_2 using namespace akantu; //Linear load function static void lin_load(double * position, double * load, __attribute__ ((unused)) Real * normal, __attribute__ ((unused)) UInt surface_id){ memset(load,0,sizeof(Real)*3); if (position[0]<=10){ load[1]= -6000; } } int main(int argc, char *argv[]){ initialize(argc, argv); debug::setDebugLevel(dblWarning); /* -------------------------------------------------------------------------- */ // Defining the mesh Mesh beams(2); UInt nb_nodes=3; UInt nb_nodes_1=1; UInt nb_nodes_2=nb_nodes-nb_nodes_1 - 1; UInt nb_element=nb_nodes-1; Array & nodes = const_cast &>(beams.getNodes()); nodes.resize(nb_nodes); beams.addConnectivityType(_bernoulli_beam_2); Array & connectivity = const_cast &>(beams.getConnectivity(_bernoulli_beam_2)); connectivity.resize(nb_element); for(UInt i=0; i & forces = model.getForce(); Array & displacement = model.getDisplacement(); Array & boundary = model.getBlockedDOFs(); const Array & N_M = model.getStress(_bernoulli_beam_2); Array & element_material = model.getElementMaterial(_bernoulli_beam_2); forces.clear(); displacement.clear(); for (UInt i = 0; i < nb_nodes_2; ++i) { element_material(i+nb_nodes_1)=1; } forces(nb_nodes-1,2) += M; model.computeForcesFromFunction<_bernoulli_beam_2>(lin_load, akantu::_bft_traction); /* -------------------------------------------------------------------------- */ // Defining the boundary conditions boundary(0,0) = true; boundary(0,1) = true; boundary(0,2) = true; boundary(nb_nodes_1,1) = true; boundary(nb_nodes-1,1) = true; /* -------------------------------------------------------------------------- */ // Solve Real error; model.assembleStiffnessMatrix(); model.getStiffnessMatrix().saveMatrix("Kb.mtx"); UInt count = 0; - model.addDumpField("displacememt"); + model.addDumpField("displacement"); model.addDumpField("rotation"); model.addDumpField("force"); model.addDumpField("momentum"); do { if(count != 0) std::cerr << count << " - " << error << std::endl; model.updateResidual(); model.solve(); count++; } while (!model.testConvergenceIncrement(1e-10, error) && count < 10); std::cerr << count << " - " << error << std::endl; /* -------------------------------------------------------------------------- */ // Post-Processing model.computeStresses(); model.getStiffnessMatrix().saveMatrix("Ka.mtx"); std::cout<< " d1 = " << displacement(nb_nodes_1,2) << std::endl; std::cout<< " d2 = " << displacement(nb_nodes-1,2) << std::endl; std::cout<< " M1 = " << N_M(0,1) << std::endl; std::cout<< " M2 = " << N_M(2*(nb_nodes-2),1) << std::endl; model.dump(); finalize(); }