diff --git a/doc/manual/manual-contactmechanicsmodel.tex b/doc/manual/manual-contactmechanicsmodel.tex new file mode 100644 index 000000000..fec7da8b7 --- /dev/null +++ b/doc/manual/manual-contactmechanicsmodel.tex @@ -0,0 +1,9 @@ +\chapter{Contact Mechanics Model\index{ContactMechanicsModel}\label{sect:cmm}} + +The contact mechanics model is a specific implementation of the +\code{Model} interface dedicated to handle the contact between two deformable bodies. + +\section{The contact virtual work} +\begin{align} +G_c(\phi, \phi^{*}) = \int_{\Gamma_c^{(1)}}\Big[t_{N}\delta g + t_{T_{\alpha}}\delta\xi^{\alpha}\Big]d\Gamma +\end{align} diff --git a/src/model/contact_mechanics/resolution.cc b/src/model/contact_mechanics/resolution.cc index 36b6bbda7..b3f9e9ef3 100644 --- a/src/model/contact_mechanics/resolution.cc +++ b/src/model/contact_mechanics/resolution.cc @@ -1,424 +1,477 @@ /** * @file resolution.cc * * @author Mohit Pundir * * @date creation: Mon Jan 7 2019 * @date last modification: Mon Jan 7 2019 * * @brief Implementation of common part of the contact resolution class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "resolution.hh" #include "contact_mechanics_model.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Resolution::Resolution(ContactMechanicsModel & model, const ID & id) : Memory(id, model.getMemoryID()), Parsable(ParserType::_contact_resolution, id), fem(model.getFEEngine()), name(""), model(model), spatial_dimension(model.getMesh().getSpatialDimension()) { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Resolution::~Resolution() = default; /* -------------------------------------------------------------------------- */ void Resolution::initialize() { registerParam("name", name, std::string(), _pat_parsable | _pat_readable); registerParam("mu", mu, Real(0.), _pat_parsable | _pat_modifiable, "Friciton Coefficient"); } /* -------------------------------------------------------------------------- */ void Resolution::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; std::string type = getID().substr(getID().find_last_of(':') + 1); stream << space << "Contact Resolution " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void Resolution::assembleInternalForces(GhostType /*ghost_type*/) { AKANTU_DEBUG_IN(); auto & internal_force = const_cast &>(model.getInternalForce()); const auto local_nodes = model.getMesh().getElementGroup(name).getNodes(); auto & nodal_area = const_cast &>(model.getNodalArea()); auto & contact_map = model.getContactMap(); for (auto & slave : local_nodes) { if (contact_map.find(slave) == contact_map.end()) continue; auto & element = contact_map[slave]; const auto & conn = element.connectivity; const auto & type = element.master.type; auto nb_nodes_master = Mesh::getNbNodesPerElement(type); Vector shapes(nb_nodes_master); Matrix dnds(spatial_dimension - 1, nb_nodes_master); #define GET_SHAPES_NATURAL(type) \ - ElementClass::computeShapes(element.projection, shapes) + ElementClass::computeShapes(element.projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL); #undef GET_SHAPES_NATURAL #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDNDS(element.projection, dnds) + ElementClass::computeDNDS(element.projection, dnds) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL Vector fc(conn.size() * spatial_dimension); - Matrix m_alpha_beta(spatial_dimension - 1, spatial_dimension - 1); - computeMetricTensor(element.tangents, m_alpha_beta); - Vector n(conn.size() * spatial_dimension); - computeN(n, shapes, element.normal); - + computeNormalForce(fc, n, element.gap); - Array t_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); - Array n_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); - Array d_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); - - computeTalpha(t_alpha, shapes, element.tangents); - computeNalpha(n_alpha, dnds, element.normal); - computeDalpha(d_alpha, n_alpha, t_alpha, m_alpha_beta, element.gap); - computeFrictionForce(fc, d_alpha, m_alpha_beta, element.projection, - element.gap); + if(mu != 0) { + Matrix m_alpha_beta(spatial_dimension - 1, spatial_dimension - 1); + computeMetricTensor(element.tangents, m_alpha_beta); + + Array t_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); + Array n_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); + Array d_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); + computeTalpha(t_alpha, shapes, + element.tangents); + computeNalpha(n_alpha, dnds, + element.normal); + computeDalpha(d_alpha, n_alpha, t_alpha, + m_alpha_beta, element.gap); + + auto traction = computeFrictionalTraction(m_alpha_beta, + element.projection, + element.gap); + computeFrictionForce(fc, d_alpha, traction); + } + UInt nb_degree_of_freedom = internal_force.getNbComponent(); for (UInt i = 0; i < conn.size(); ++i) { UInt n = conn[i]; for (UInt j = 0; j < nb_degree_of_freedom; ++j) { UInt offset_node = n * nb_degree_of_freedom + j; internal_force[offset_node] += fc[i * nb_degree_of_freedom + j]; internal_force[offset_node] *= nodal_area[n]; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Resolution::assembleStiffnessMatrix(GhostType /*ghost_type*/) { AKANTU_DEBUG_IN(); auto & stiffness = const_cast(model.getDOFManager().getMatrix("K")); - const auto local_nodes = model.getMesh().getElementGroup(name).getNodes(); + const auto local_nodes = + model.getMesh().getElementGroup(name).getNodes(); - auto & nodal_area = const_cast &>(model.getNodalArea()); + auto & nodal_area = + const_cast &>(model.getNodalArea()); auto & contact_map = model.getContactMap(); for (auto & slave : local_nodes) { - if (contact_map.find(slave) == contact_map.end()) { + if (contact_map.find(slave) == contact_map.end()) continue; - } - - auto & master = contact_map[slave].master; - auto & gap = contact_map[slave].gap; - auto & projection = contact_map[slave].projection; - auto & normal = contact_map[slave].normal; - const auto & connectivity = contact_map[slave].connectivity; - const ElementType & type = master.type; - - UInt nb_nodes_master = Mesh::getNbNodesPerElement(master.type); + + auto & element = contact_map[slave]; + + const auto & conn = element.connectivity; + const auto & type = element.master.type; + + UInt nb_nodes_master = Mesh::getNbNodesPerElement(type); Vector shapes(nb_nodes_master); Matrix shapes_derivatives(spatial_dimension - 1, nb_nodes_master); Matrix shapes_second_derivatives((spatial_dimension-1)*(spatial_dimension-1) , nb_nodes_master); -#define GET_SHAPES_NATURAL(type) \ - ElementClass::computeShapes(projection, shapes) +#define GET_SHAPES_NATURAL(type) \ + ElementClass::computeShapes(element.projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL); #undef GET_SHAPES_NATURAL #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDNDS(projection, shapes_derivatives) + ElementClass::computeDNDS(element.projection, shapes_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL #define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDN2DS2(projection, shapes_second_derivatives) + ElementClass::computeDN2DS2(element.projection, shapes_second_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL); #undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL - - Matrix elementary_stiffness(connectivity.size() * spatial_dimension, - connectivity.size() * spatial_dimension); - Matrix tangents(spatial_dimension - 1, spatial_dimension); - Matrix global_coords(nb_nodes_master, spatial_dimension); + Matrix kc(conn.size() * spatial_dimension, + conn.size() * spatial_dimension); + + Matrix m_alpha_beta(spatial_dimension - 1, spatial_dimension - 1); + computeMetricTensor(element.tangents, m_alpha_beta); - computeCoordinates(master, global_coords); - computeTangents(shapes_derivatives, global_coords, tangents); + // normal tangent moduli + Vector n(conn.size() * spatial_dimension); + Array t_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); + Array n_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); + Array d_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); - Matrix m_alpha_beta(spatial_dimension - 1, spatial_dimension - 1); - computeMetricTensor(tangents, m_alpha_beta); - - Vector n(connectivity.size() * spatial_dimension); - Array t_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); - Array n_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); - Array d_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); - - computeN( n, shapes, normal); - computeTalpha( t_alpha, shapes, tangents); - computeNalpha( n_alpha, shapes_derivatives, normal); - computeDalpha( d_alpha, n_alpha, t_alpha, m_alpha_beta, gap); - - Array t_alpha_beta(connectivity.size() * spatial_dimension, (spatial_dimension - 1) * (spatial_dimension -1)); - computeTalphabeta(t_alpha_beta, shapes_derivatives, tangents); - - Array p_alpha(connectivity.size() * spatial_dimension, spatial_dimension - 1); - Vector p_t; - computePalpha(p_alpha, shapes_derivatives, p_t); - - Array n_alpha_beta(connectivity.size() * spatial_dimension, - (spatial_dimension - 1) * (spatial_dimension -1)); - computeNalphabeta(n_alpha_beta, shapes_second_derivatives, normal); - - Matrix kc(connectivity.size() * spatial_dimension, - connectivity.size() * spatial_dimension); - computeTangentModuli(kc, n, n_alpha, d_alpha, m_alpha_beta, gap); + computeN( n, shapes, element.normal); + computeTalpha( t_alpha, shapes, element.tangents); + computeNalpha( n_alpha, shapes_derivatives, element.normal); + computeDalpha( d_alpha, n_alpha, t_alpha, m_alpha_beta, element.gap); + + computeNormalModuli(kc, n, n_alpha, d_alpha, m_alpha_beta, element.gap); + + // frictional tangent moduli + if(mu != 0) { + Array t_alpha_beta(conn.size() * spatial_dimension, + (spatial_dimension - 1) * (spatial_dimension -1)); + Array p_alpha(conn.size() * spatial_dimension, + spatial_dimension - 1); + Array n_alpha_beta(conn.size() * spatial_dimension, + (spatial_dimension - 1) * (spatial_dimension -1)); + + auto traction = computeFrictionalTraction(m_alpha_beta, element.projection, + element.gap); + computeTalphabeta(t_alpha_beta, shapes_derivatives, element.tangents); + computeNalphabeta(n_alpha_beta, shapes_second_derivatives, element.normal); + computePalpha(p_alpha, shapes_derivatives, traction); + //computeGalpha(); + //computePbaralpha(); + //computeTbaralphabeta(); + + computeFrictionalModuli(kc, t_alpha_beta, n_alpha_beta, element.tangents, + shapes_second_derivatives, n, n_alpha, d_alpha, + element.gap); + } std::vector equations; UInt nb_degree_of_freedom = model.getSpatialDimension(); std::vector areas; for (UInt i : arange(connectivity.size())) { UInt n = connectivity[i]; for (UInt j : arange(nb_degree_of_freedom)) { equations.push_back(n * nb_degree_of_freedom + j); areas.push_back(nodal_area[n]); } } for (UInt i : arange(kc.rows())) { UInt row = equations[i]; for (UInt j : arange(kc.cols())) { UInt col = equations[j]; kc(i, j) *= areas[i]; stiffness.add(row, col, kc(i, j)); } } } AKANTU_DEBUG_OUT(); } -/* -------------------------------------------------------------------------- */ -void Resolution::computeTangents(Matrix & shapes_derivatives, - Matrix & global_coords, - Matrix & tangents) { - - tangents.mul(shapes_derivatives, global_coords); -} - /* -------------------------------------------------------------------------- */ void Resolution::computeMetricTensor(Matrix & tangents, Matrix & m_alpha_beta) { m_alpha_beta.mul(tangents, tangents); m_alpha_beta = m_alpha_beta.inverse(); } /* -------------------------------------------------------------------------- */ void Resolution::computeN(Vector & n, Vector & shapes, Vector & normal) { UInt dim = normal.size(); for (UInt i = 0; i < dim; ++i) { n[i] = normal[i]; for (UInt j = 0; j < shapes.size(); ++j) { n[(1 + j) * dim + i] = -normal[i] * shapes[j]; } } } /* -------------------------------------------------------------------------- */ void Resolution::computeTalpha(Array & t_alpha, Vector & shapes, Matrix & tangents) { t_alpha.clear(); for (auto && values : zip(tangents.transpose(), make_view(t_alpha, t_alpha.size()))) { auto & tangent = std::get<0>(values); auto & t_s = std::get<1>(values); for (UInt i : arange(spatial_dimension)) { t_s[i] = -tangent(i); for (UInt j : arange(shapes.size())) { t_s[(1 + j) * spatial_dimension + i] = -shapes[j] * tangent(i); } } } } /* -------------------------------------------------------------------------- */ void Resolution::computeNalpha(Array & n_alpha, Matrix & shapes_derivatives, Vector & normal) { n_alpha.clear(); for (auto && values : zip(shapes_derivatives.transpose(), make_view(n_alpha, n_alpha.size()))) { auto & dnds = std::get<0>(values); auto & n_s = std::get<1>(values); for (UInt i : arange(spatial_dimension)) { n_s[i] = 0; for (UInt j : arange(dnds.size())) { n_s[(1 + j) * spatial_dimension + i] = -dnds(j) * normal[i]; } } } } /* -------------------------------------------------------------------------- */ void Resolution::computeDalpha(Array & d_alpha, Array & n_alpha, Array & t_alpha, Matrix & m_alpha_beta, Real & gap) { d_alpha.clear(); for (auto && entry : zip(m_alpha_beta.transpose(), make_view(d_alpha, d_alpha.size()))) { auto & a_s = std::get<0>(entry); auto & d_s = std::get<1>(entry); for (auto && values : zip(arange(t_alpha.size()), make_view(t_alpha, t_alpha.size()), make_view(n_alpha, n_alpha.size()))) { auto & index = std::get<0>(values); auto & t_s = std::get<1>(values); auto & n_s = std::get<2>(values); d_s += (t_s + gap * n_s); d_s *= a_s(index); } } } /* -------------------------------------------------------------------------- */ void Resolution::computeTalphabeta(Array & t_alpha_beta, Matrix & shapes_derivatives, Matrix & tangents) { t_alpha_beta.clear(); auto t_alpha_size = t_alpha_beta.size() * (spatial_dimension - 1); for(auto && entry : zip(tangents.transpose(), make_view(t_alpha_beta, t_alpha_size))) { - auto & tangent = std::get<0>(entry); - auto & t_alpha = std::get<1>(entry); - //for(auto && values : zip(shapes_derivatives.transpose(), - // make_view(t_alpha, t_alpha_size))) { - // auto & dnds = std::get<0>(values); - //auto & t_alpha_beta = std::get<1>(values); - - //t_alpha_beta += dnds * tangent; - //} - + auto & tangent_s = std::get<0>(entry); + auto & t_alpha = std::get<1>(entry); + for(auto && values : zip(shapes_derivatives.transpose(), + make_view(t_alpha, t_alpha_beta.size()))) { + auto & dnds = std::get<0>(values); + auto & t_alpha_s = std::get<1>(values); + + for (UInt i : arange(spatial_dimension)) { + t_alpha_s[i] = 0; + for (UInt j : arange(dnds.size())) { + t_alpha_s[(1 + j) * spatial_dimension + i] = -dnds(j) * tangent_s[i]; + } + } + } } } /* -------------------------------------------------------------------------- */ void Resolution::computeNalphabeta(Array & n_alpha_beta, Matrix & shapes_second_derivatives, Vector & normal) { n_alpha_beta.clear(); for(auto && entry : zip(shapes_second_derivatives.transpose(), make_view(n_alpha_beta, n_alpha_beta.size()))) { - auto & dn2ds2 = std::get<0>(entry); - auto & n_alpha = std::get<1>(entry); - - + auto & dn2ds2 = std::get<0>(entry); + auto & n_alpha_s = std::get<1>(entry); + + for (UInt i : arange(spatial_dimension)) { + n_alpha_s[i] = 0; + for (UInt j : arange(dn2ds2.size())) { + n_alpha_s[(1 + j) * spatial_dimension + i] = -dn2ds2(j) * normal[i]; + } + } + } } /* -------------------------------------------------------------------------- */ void Resolution::computePalpha(Array & p_alpha, Matrix & shapes_derivatives, - Vector & p_t) { + Vector & traction) { + p_alpha.clear(); + auto normalized_traction = traction/traction.norm(); + + for(auto && entry : + zip(shapes_derivatives.transpose(), + make_view(p_alpha, p_alpha.size()))) { + auto & dnds = std::get<0>(entry); + auto & p_s = std::get<1>(entry); + + for(UInt i : arnage(spatial_dimension)) { + p_s[i] = 0; + for(UInt j : arange(dnds.size())){ + p_s[(1 + j) * spatial_dimension + i] = -dnds(j) * normalized_traction[i]; + } + } + } } /* -------------------------------------------------------------------------- */ void Resolution::computeGalpha(Array & /*t_alpha_beta*/, Array & /*d_alpha*/, Vector & /*tangential_gap*/) { Array g_alpha(d_alpha.size(), spatial_dimension - 1); + auto t_alpha_size = t_alpha_beta.size() * (spatial_dimension - 1); for(auto && value : - make_view(g_alpha, g_alpha.size())){ + zip(make_view(g_alpha, g_alpha.size()), + make_view(t_alpha_beta, t_alpha_size))){ auto & g_s = std::get<0>(value); - } - + } +} + +/* -------------------------------------------------------------------------- */ +void Resolution::computeTbaralphabeta() { + +} + +/* -------------------------------------------------------------------------- */ +void Resolution::computePbaralpha() { + } /* -------------------------------------------------------------------------- */ void Resolution::computeCoordinates(const Element & el, Matrix & coords) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); Vector connect = model.getMesh() .getConnectivity(el.type, _not_ghost) .begin(nb_nodes_per_element)[el.element]; - // change this to current position + // todo change this to current position auto & positions = model.getMesh().getNodes(); for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connect[n]; for (UInt s : arange(spatial_dimension)) { coords(n, s) = positions(node, s); } } } +/* -------------------------------------------------------------------------- */ +void Resolution::computeSecondDerivative(Matrix & dn2ds2, const Element & el) { + + UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); + Vector connect = model.getMesh().getConnectivity(el.type, _not_ghost).begin(nb_nodes_per_element)[el.element]; + + Matrix phi(dn2ds2.rows(), spatial_dimension); + Matrix nodal_values(nb_nodes_per_element, spatial_dimension); + + computeCoordinates(el, nodal_values); + phi.mul(dn2ds2, nodal_values); +} + } // namespace akantu diff --git a/src/model/contact_mechanics/resolution.hh b/src/model/contact_mechanics/resolution.hh index c999a4de6..c94ae6811 100644 --- a/src/model/contact_mechanics/resolution.hh +++ b/src/model/contact_mechanics/resolution.hh @@ -1,254 +1,263 @@ /** * @file resolution.hh * * @author Mohit Pundir * * @date creation: Mon Jan 7 2019 * @date last modification: Mon Jan 7 2019 * * @brief Mother class for all contact resolutions * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_factory.hh" #include "aka_memory.hh" #include "parsable.hh" #include "parser.hh" #include "fe_engine.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_RESOLUTION_HH__ #define __AKANTU_RESOLUTION_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class ContactMechanicsModel; } // namespace akantu namespace akantu { /** * Interface of all contact resolutions * Prerequisites for a new resolution * - inherit from this class * - implement the following methods: * \code * * virtual void computeNormalForce(); * virtual void computeFrictionForce(); * * \endcode * */ class Resolution : public Memory, public Parsable { /* ------------------------------------------------------------------------ */ /* Constructor/Destructor */ /* ------------------------------------------------------------------------ */ public: /// instantiate contact resolution with defaults Resolution(ContactMechanicsModel & model, const ID & id = ""); /// Destructor ~Resolution() override; protected: void initialize(); /// computes coordinates of a given element void computeCoordinates(const Element & , Matrix &); /// computes tangents void computeTangents(Matrix & , Matrix & , Matrix & ); /// computes metric tensor (@f$m_{\alpha\beta}@f$) where @f$\alpha, /// \beta@f$ are surface directions void computeMetricTensor(Matrix &, Matrix &); /// computes N array void computeN(Vector &, Vector &, Vector &); /// computes (@f$N_{\alpha}@f$) where \alpha is surface dimension void computeNalpha(Array &, Matrix &, Vector &); /// computes (@f$T_{\alpha}@f$) where @f$\alpha@f$ is surface dimension void computeTalpha(Array & , Vector & , Matrix & ); /// computes (@f$D_{\alpha}@f$) where @f$\alpha@f$ is surface dimension void computeDalpha(Array & , Array & , Array & , Matrix &, Real & ); /// computes @f$T_{\alpha\beta} @f$ void computeTalphabeta(Array &, Matrix &, Matrix &); /// computes @f$N_{\alpha\beta} @f$ void computeNalphabeta(Array &, Matrix &, Vector &); /// computes @f$P_{\alpha} @f$ void computePalpha(Array &, Matrix &, Vector &); /// computes @f$G_{\alpha}@f$ void computeGalpha(Array &, Array &, Array &, Vector &) /* ------------------------------------------------------------------------ */ /* Functions that resolutions can/should reimplement */ /* ------------------------------------------------------------------------ */ protected: - /// computes the normal force - virtual void computeNormalForce(__attribute__((unused)) Vector & /* force */, - __attribute__((unused)) Vector & /* n */, - __attribute__((unused)) Real & /* gap */) { + /// computes the force vector due to normal traction + virtual void computeNormalForce(__attribute__((unused)) Vector &, + __attribute__((unused)) Vector &, + __attribute__((unused)) Real &) { AKANTU_TO_IMPLEMENT(); } - /// computes the friction force - virtual void computeFrictionForce(__attribute__((unused)) Vector & /* force */, - __attribute__((unused)) Array & /* d_alpha */, - __attribute__((unused)) Matrix &, - __attribute__((unused)) Vector &, - __attribute__((unused)) Real & /* gap */) { + /// computes the force vector due to frictional traction + virtual void computeFrictionForce(__attribute__((unused)) Vector & , + __attribute__((unused)) Array & , + __attribute__((unused)) Vector &) { AKANTU_TO_IMPLEMENT(); } - /// compute the tangent moduli - virtual void computeTangentModuli(__attribute__((unused)) Matrix &, - __attribute__((unused)) Vector & /* n */, - __attribute__((unused)) Array & /* n_alpha */, - __attribute__((unused)) Array & /* d_alpha */, - __attribute__((unused)) Matrix & /* A */, - __attribute__((unused)) Real & /* gap */) { + /// compute the tangent moduli due to normal traction + virtual void computeNormalModuli(__attribute__((unused)) Matrix &, + __attribute__((unused)) Vector &, + __attribute__((unused)) Array &, + __attribute__((unused)) Array &, + __attribute__((unused)) Matrix &, + __attribute__((unused)) Real &) { + AKANTU_TO_IMPLEMENT(); + } + + + /// compute the tangent moduli due to frictional traction + virtual void computeFrictionalModuli(__attribute__((unused)) Matrix &, + __attribute__((unused)) Vector &, + __attribute__((unused)) Array &, + __attribute__((unused)) Array &, + __attribute__((unused)) Matrix &, + __attribute__((unused)) Real &) { AKANTU_TO_IMPLEMENT(); } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// assemble the residual for this resolution void assembleInternalForces(GhostType ghost_type); /// assemble the stiffness matrix for this resolution void assembleStiffnessMatrix(GhostType ghost_type); public: /// function to print the contain of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Link to the fem object in the model FEEngine & fem; /// resolution name std::string name; /// model to which the resolution belong ContactMechanicsModel & model; /// friciton coefficient : mu Real mu; /// spatial dimension UInt spatial_dimension; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Resolution & _this) { _this.printself(stream); return stream; } } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { using ResolutionFactory = Factory; /// macaulay bracket to convert positive gap to zero template T macaulay(T var) {return var < 0 ? 0 : var; } template T heaviside(T var) {return var < 0 ? 0 : 1; } } // namespace akantu #define INSTANTIATE_RESOLUTION_ONLY(res_name) \ class res_name #define RESOLUTION_DEFAULT_PER_DIM_ALLOCATOR(id, res_name) \ [](UInt dim, const ID &, ContactMechanicsModel & model, \ const ID & id) -> std::unique_ptr { \ switch (dim) { \ case 1: \ return std::make_unique(model, id); \ case 2: \ return std::make_unique(model, id); \ case 3: \ return std::make_unique(model, id); \ default: \ AKANTU_EXCEPTION("The dimension " \ << dim << "is not a valid dimension for the contact resolution " \ << #id); \ } \ } #define INSTANTIATE_RESOLUTION(id, res_name) \ INSTANTIATE_RESOLUTION_ONLY(res_name); \ static bool resolution_is_alocated_##id[[gnu::unused]] = \ ResolutionFactory::getInstance().registerAllocator( \ #id, RESOLUTION_DEFAULT_PER_DIM_ALLOCATOR(id, res_name)) #endif /* __AKANTU_RESOLUTION_HH__ */ diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.cc b/src/model/contact_mechanics/resolutions/resolution_penalty.cc index a1cb194be..501503d6f 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.cc +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.cc @@ -1,258 +1,319 @@ /** * @file resolution_penalty.cc * * @author Mohit Pundir * * @date creation: Mon Jan 7 2019 * @date last modification: Mon Jan 7 2019 * * @brief Specialization of the resolution class for the penalty method * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "resolution_penalty.hh" namespace akantu { /* -------------------------------------------------------------------------- */ ResolutionPenalty::ResolutionPenalty(ContactMechanicsModel & model, const ID & id) : Resolution(model, id) { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::initialize() { this->registerParam("epsilon", epsilon, Real(0.), _pat_parsable | _pat_modifiable, "Normal penalty parameter"); this->registerParam("epsilon_t", epsilon_t, Real(0.), _pat_parsable | _pat_modifiable, "Tangential penalty parameter"); } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeNormalForce(Vector & force, Vector & n, Real & gap) { force.clear(); Real tn = gap * epsilon; tn = macaulay(tn); for (UInt i : arange(force.size())) { force[i] += tn * n[i]; } } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeFrictionForce(Vector & force, Array & d_alpha, - Matrix & m_alpha_beta, - Vector & projection, - Real & gap) { + Vector & traction) { - auto traction = computeFrictionalTraction(m_alpha_beta, projection, gap); for (auto && values: zip(traction, make_view(d_alpha, d_alpha.size()))) { auto & t_s = std::get<0>(values); auto & d_s = std::get<1>(values); force += d_s * t_s; } } - -/* -------------------------------------------------------------------------- */ -Vector ResolutionPenalty::computeFrictionalTraction(Matrix& m_alpha_beta, - Vector& xi, - Real & gap) { - - Real tn = gap * epsilon; - tn = macaulay(tn); - - Vector previous_xi; - - auto delta_xi = xi - previous_xi; - - Vector previous_traction; - Vector trial_traction; - - trial_traction.mul(m_alpha_beta, delta_xi, epsilon); - //trial_traction += previous_traction; - - auto trial_slip_function = trial_traction.norm() - mu * tn; - - Vector traction; - if (trial_slip_function <= 0) { - traction = trial_traction; - } - else{ - traction = mu * tn * trial_traction / trial_traction.norm(); - } - - return traction; -} - - -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeTangentModuli(Matrix & kc, Vector & n, - Array & n_alpha, Array & d_alpha, - Matrix & m_alpha_beta, - Real & gap) { - computeNormalStiffness(kc, n, n_alpha, d_alpha, m_alpha_beta, gap); - computeFrictionalStiffness(n, n_alpha, d_alpha, gap); -} - - + /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeNormalStiffness(Matrix & ke, Vector & n, - Array & n_alpha, Array & d_alpha, - Matrix & /*m_alpha_beta*/, Real & gap) { +void ResolutionPenalty::computeNormalModuli(Matrix & ke, Vector & n, + Array & n_alpha, Array & d_alpha, + Matrix & /*m_alpha_beta*/, Real & gap) { Real tn = gap * epsilon; tn = macaulay(tn); Matrix n_mat(n.storage(), n.size(), 1); ke.mul(n_mat, n_mat); ke *= epsilon * heaviside(gap); for (auto && values: zip(make_view(n_alpha, n_alpha.size()), make_view(d_alpha, d_alpha.size()))) { auto & n_s = std::get<0>(values); auto & d_s = std::get<1>(values); Matrix ns_mat(n_s.storage(), n_s.size(), 1); Matrix ds_mat(d_s.storage(), d_s.size(), 1); Matrix tmp1(n_s.size(), n_s.size()); tmp1.mul(ns_mat, ds_mat); Matrix tmp2(n_s.size(), n_s.size()); tmp1.mul(ds_mat, ns_mat); ke -= (tmp1 + tmp2) * tn; } } /* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeFrictionalStiffness(Vector & /*n*/, - Array & /*n_alpha*/, - Array & /*d_alpha*/, - Real & gap) { +void ResolutionPenalty::computeFrictionalModuli(Matrix & ke, + Array & t_alpha_beta, + Array & n_alpha_beta, + Matrix & tangents, + Matrix & dn2ds2, + Vector & n, + Array & /*n_alpha*/, + Array & /*d_alpha*/, + Real & gap) { + + auto k_common = computeCommonModuli(t_alpha_beta, n_alpha_beta, d_alpha, + tangents, n, gap); + + bool stick; + Vector traction; + std::tie(traction, stick) = + computeFrictionalTraction(m_alpha_beta, projection, gap); + + Matrix k_t; + if(stick) + k_t = computeStickModuli(g_alpha, d_alpha, + m_alpha_beta); + else { + k_t = computeSlipModuli(g_alpha, d_alpha, + m_alpha_beta, traction); + } +} - computeCommonModuli(gap); - computeStickModuli(); - computeSlipModuli(); -} +/* -------------------------------------------------------------------------- */ +auto ResolutionPenalty::computeFrictionalTraction(Matrix& m_alpha_beta, + Vector& xi, + Real & gap) { + Real tn = gap * epsilon; + tn = macaulay(tn); + + Vector previous_xi; + + auto delta_xi = xi - previous_xi; + + Vector previous_traction; + Vector trial_traction; + + trial_traction.mul(m_alpha_beta, delta_xi, epsilon); + //trial_traction += previous_traction; + + auto trial_slip_function = trial_traction.norm() - mu * tn; + + bool stick = false; + Vector traction; + + if (trial_slip_function <= 0) { + traction = trial_traction; + stick = true; + } + else{ + traction = mu * tn * trial_traction / trial_traction.norm(); + } + + return std::make_pair(traction, stick); +} + /* -------------------------------------------------------------------------- */ -Matrix ResolutionPenalty::computeCommonModuli(Array & /*t_alpha_beta*/, - Array & /*n_alpha_beta*/, - Matrix & /*tangents*/ - Vector & /* n */ - Real & /*gap*/) { - - Array kt_alpha(spatial_dimension -1, d_alpha.size() * d_alpha.size(), - "k_T_alpha"); +Matrix ResolutionPenalty::computeCommonModuli(Array & t_alpha_beta, + Array & n_alpha_beta, + Array & d_alpha, + Matrix & tangents, + Vector & n + Real & gap) { + + Array kt_alpha(spatial_dimension -1, + d_alpha.size() * d_alpha.size(), "k_T_alpha"); + + auto t_alpha_beta_size = t_alpha_beta.size() + * (spatial_dimension - 1); + for(auto && values : zip(tangents.transpose(), make_view(kt_alpha, kt_alpha.size()), - make_view(t_alpha_beta, t_alpha_beta.size()), - make_view(n_alpha_beta, n_alpha_beta.size()))) { + make_view(t_alpha_beta, t_alpha_beta_size), + make_view(n_alpha_beta, n_alpha_beta_size))) { auto & tangent_s = std::get<0>(values); auto & kt_s = std::get<1>(values); auto & t_alpha_s = std::get<2>(values); auto & n_alpha_s = std::get<3>(values); - Matrix kt_s_mat(kt_s.storage(), d_alpha.size(), d_alpha.size()) + Matrix kt_s_mat(kt_s.storage(), d_alpha.size(), d_alpha.size()); + // loop over beta for(auto && tuple : - make_view(d_alpha, d_alpha.size())) { + zip(make_view(d_alpha, d_alpha.size()), + make_view(n_alpha_))) { auto & d_s = std::get<0>(tuple); + Matrix tmp(d_s.size(), d_s.size()); + // loop over gamma for(auto && entry : make_view(d_alpha, d_alpha.size())) { auto & d_t = std::get<0>(entry); + + // compute constant - Matrix tmp(d_t.size(), d_t.size()); - tmp.mul(d_s, d_t); - kt_s_mat += tmp; + Matrix tmp2(d_t.size(), d_t.size()); + tmp2.mul(d_s, d_t); + kt_s_mat += tmp2; } } } } /* -------------------------------------------------------------------------- */ Matrix ResolutionPenalty::computeStickModuli(Array & g_alpha, Array & d_alpha, Matrix & m_alpha_beta) { Matrix k_stick(d_alpha.size(), d_alpha.size()); for (auto && values : zip(arange(d_alpha.getNbComponent()) make_view(d_alpha, d_alpha.size()), make_view(g_alpha, g_alpha.size()))) { auto & s = std::get<0>(values); auto & d_s = std::get<1>(values); auto & g_s = std::get<2>(values); Matrix ds_mat(d_s.storage(), d_s.size(), 1); Matrix gs_mat(g_s.storage(), g_s.size(), 1); Matrix tmp1(d_s.size(), d_s.size()); tmp1.mul(ds_mat, gs_mat); k_stick += tmp1; for(auto && tuple : enumerate(make_view(d_alpha, d_alpha.size()))) { auto & t = std::get<0>(tuple); auto & d_t = std::get<1>(tuple); Matrix dt_mat(d_t.storage(), d_t.size(), 1); Matrix tmp2(d_t.size(), d_t.size()); tmp2.mul(ds_mat, dt_mat); k_stick += tmp2 * m_alpha_beta(s, t); } k_stick *= epsilon_t; return k_stick; } /* -------------------------------------------------------------------------- */ -Matrix ResolutionPenalty::computeSlipModuli() { +Matrix ResolutionPenalty::computeSlipModuli(Array & g_alpha, + Array & d_alpha, + Matrix & m_alpha_beta, + Vector & traction) { + + + Real tn = gap * epsilon; + tn = macaulay(tn); + + Real factor; + factor = epsilon_t * mu * tn; + factor /= traction.norm(); + + auto p_t = traction; + p_t /= p_t.norm(); + + Matrix k_slip(d_alpha.size(), d_alpha.size()); + // loop for alpha + for(auto && value : + make_view(d_alpha, d_alpha.size())) { + auto & d_s = std::get<0>(value); + + // loop for beta + for(auto && tuple : + zip(arange(spatial_dimension - 1), + make_view(d_alpha, d_alpha.size()), + make_view(g_alpha, g_alpha.size()))) { + auto & beta = std::get<0>(tuple); + auto & d_beta = std::get<1>(tuple); + auto & g_beta = std::get<2>(tuple); + + // loop for gamma + for(auto && entry : + zip(arange(spatial_dimension - 1), + make_view(d_alpha, d_alpha.size()))) { + auto & gamma = std::get<0>(entry); + auto & d_gamma = std::get<1>(entry); + + } + } + + + } + } INSTANTIATE_RESOLUTION(penalty, ResolutionPenalty); } // akantu diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.hh b/src/model/contact_mechanics/resolutions/resolution_penalty.hh index bb9dc0e2b..f6ac9ac73 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.hh +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.hh @@ -1,123 +1,131 @@ /** * @file contact_resolution_penalty.hh * * @author Mohit Pundir * * @date creation: Fri Jun 18 2010 * @date last modification: Mon Jan 29 2018 * * @brief Material isotropic thermo-elastic * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "resolution.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_RESOLUTION_PENALTY_HH__ #define __AKANTU_RESOLUTION_PENALTY_HH__ namespace akantu { class ResolutionPenalty : public Resolution { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ResolutionPenalty(ContactMechanicsModel & model, const ID & id = ""); ~ResolutionPenalty() override = default; protected: /// initialize the resolution void initialize(); /// local computation of stifnness matrix due to normal stress void computeNormalStiffness(Matrix & kc, Vector & n, Array & n_alpha, Array & d_alpha, Matrix & surface_matrix, Real &); /// local computation of stiffness matrix due to frictional stress void computeFrictionalStiffness(Vector &, Array &, Array &, Real & ); /// local computation of direct stiffness matrix due to friction, /// this matrix is common for both stick and slip part Matrix computeCommonModuli(Real &); /// local computaion of stiffness matrix due to stick state Matrix computeStickModuli(Array &, Array &, Matrix &); /// local computation of stiffness matrix due to slip state Matrix computeSlipModuli(); /// computes the tractions using return map algorithm Vector computeFrictionalTraction(Matrix &, Vector &, Real &); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - /// local computation of tangent moduli - void computeTangentModuli(Matrix &, - Vector & /* N */, - Array & /* N_alpha */, - Array & /* D_alpha */, - Matrix & /* A matrix */, - Real & /* gap */ - ) override; - + /// local computation of tangent moduli due to normal traction + void computeNormalModuli(Matrix & , + Vector & , + Array & , + Array & , + Matrix & , + Real & ) override; + + /// local computation of tangent moduli due to frictional traction + void computeFrictionalModuli(Matrix & , + Array & , + Array & , + Matrix & , + Matrix & , + Vector & , + Array & /*n_alpha*/, + Array & /*d_alpha*/, + Real & ) override; + /// local computation of normal force void computeNormalForce(Vector &, Vector &, Real & ) override; /// local computation of friction force void computeFrictionForce(Vector &, Array &, - Matrix &, - Vector &, - Real &) override; + Vector &) override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// penalty parameter for normal stress Real epsilon; /// penalty parameter for tangential stress Real epsilon_t; }; } // akantu #endif /* __AKANTU_RESOLUTION_PENALTY_HH__ */