diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index 6f2038b74..fd4533968 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,588 +1,591 @@ /** * @file coontact_mechanics_model.cc * * @author Mohit Pundir * * @date creation: Tue May 08 2012 * @date last modification: Wed Feb 21 2018 * * @brief Contact mechanics model * * @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 "contact_mechanics_model.hh" #include "boundary_condition_functor.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, Array & positions, UInt dim, const ID & id, const MemoryID & memory_id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, dof_manager, dim, id, memory_id), current_positions(positions) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("ContactMechanicsModel", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("contact_mechanics", id, true); this->mesh.addDumpMeshToDumper("contact_mechanics", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->registerDataAccessor(*this); this->detector = std::make_unique(this->mesh, current_positions, id + ":contact_detector"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, std::shared_ptr dof_manager, const ModelType model_type) : ContactMechanicsModel(mesh, mesh.getNodes(), dim, id, memory_id, dof_manager, model_type) {} /* -------------------------------------------------------------------------- */ ContactMechanicsModel::~ContactMechanicsModel() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); // initalize the resolutions if (this->parser.getLastParsedFile() != "") { this->instantiateResolutions(); } this->initResolutions(); this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::instantiateResolutions() { ParserSection model_section; bool is_empty; std::tie(model_section, is_empty) = this->getParserSection(); if (not is_empty) { auto model_resolutions = model_section.getSubSections(ParserType::_contact_resolution); for (const auto & section : model_resolutions) { this->registerNewResolution(section); } } auto sub_sections = this->parser.getSubSections(ParserType::_contact_resolution); for (const auto & section : sub_sections) { this->registerNewResolution(section); } if (resolutions.empty()) AKANTU_EXCEPTION("No contact resolutions where instantiated for the model" << getID()); are_resolutions_instantiated = true; } /* -------------------------------------------------------------------------- */ Resolution & ContactMechanicsModel::registerNewResolution(const ParserSection & section) { std::string res_name; std::string res_type = section.getName(); std::string opt_param = section.getOption(); try { std::string tmp = section.getParameter("name"); res_name = tmp; /** this can seem weird, but there is an ambiguous operator * overload that i couldn't solve. @todo remove the * weirdness of this code */ } catch (debug::Exception &) { AKANTU_ERROR("A contact resolution of type \'" << res_type << "\' in the input file has been defined without a name!"); } Resolution & res = this->registerNewResolution(res_name, res_type, opt_param); res.parseSection(section); return res; } /* -------------------------------------------------------------------------- */ Resolution & ContactMechanicsModel::registerNewResolution( const ID & res_name, const ID & res_type, const ID & opt_param) { AKANTU_DEBUG_ASSERT(resolutions_names_to_id.find(res_name) == resolutions_names_to_id.end(), "A resolution with this name '" << res_name << "' has already been registered. " << "Please use unique names for resolutions"); UInt res_count = resolutions.size(); resolutions_names_to_id[res_name] = res_count; std::stringstream sstr_res; sstr_res << this->id << ":" << res_count << ":" << res_type; ID res_id = sstr_res.str(); std::unique_ptr resolution = ResolutionFactory::getInstance().allocate(res_type, spatial_dimension, opt_param, *this, res_id); resolutions.push_back(std::move(resolution)); return *(resolutions.back()); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initResolutions() { AKANTU_DEBUG_ASSERT(resolutions.size() != 0, "No resolutions to initialize !"); if (!are_resolutions_instantiated) instantiateResolutions(); // \TODO check if each resolution needs a initResolution() method } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initModel() { getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & ContactMechanicsModel::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initSolver( TimeStepSolverType /*time_step_solver_type*/, NonLinearSolverType) { this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->displacement_increment, spatial_dimension, "displacement_increment"); this->allocNodalField(this->contact_force, spatial_dimension, "contact_force"); this->allocNodalField(this->external_force, spatial_dimension, "external_force"); this->allocNodalField(this->normals, spatial_dimension, "normals"); this->allocNodalField(this->gaps, 1, "gaps"); this->allocNodalField(this->nodal_area, 1, "areas"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); } /* -------------------------------------------------------------------------- */ std::tuple ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_contact: { return std::make_tuple("explicit_contact", _tsst_dynamic); } case _implicit_contact: { return std::make_tuple("implicit_contact", _tsst_static); } default: return std::make_tuple("unkown", _tsst_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions ContactMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case _tsst_dynamic: { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } case _tsst_static: { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); break; } return options; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the contact forces"); // assemble the forces due to local stresses auto assemble = [&](auto && ghost_type) { for (auto & resolution : resolutions) { resolution->assembleInternalForces(ghost_type); } }; AKANTU_DEBUG_INFO("Assemble residual for local elements"); assemble(_not_ghost); // assemble the stresses due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); assemble(_ghost); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->contact_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { this->detector->search(this->contact_map); this->assembleFieldsFromContactMap(); this->computeNodalAreas(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search(Array & increment) { this->detector->search(this->contact_map); for (auto & entry : contact_map) { auto & element = entry.second; const auto & connectivity = element.connectivity; auto master_node = connectivity[1]; Vector u(spatial_dimension); for (UInt s : arange(spatial_dimension)) { u(s) = increment(master_node, s); } u *= -1.0; const auto & normal = element.normal; Real uv = Math::vectorDot(u.storage(), normal.storage(), spatial_dimension); if (uv + element.gap <= 0) { element.gap = abs(uv + element.gap); } else { element.gap = 0.0; } } this->assembleFieldsFromContactMap(); this->computeNodalAreas(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleFieldsFromContactMap() { if (this->contact_map.empty()) AKANTU_ERROR( "Contact map is empty, Please run search before assembling the fields"); for (auto & entry : contact_map) { const auto & element = entry.second; auto connectivity = element.connectivity; auto node = connectivity(0); (*gaps)[node] = element.gap; for (UInt i = 0; i < spatial_dimension; ++i) (*normals)(node, i) = element.normal[i]; } } /* -------------------------------------------------------------------------- */ template void ContactMechanicsModel::computeNodalAreas() { this->nodal_area->clear(); this->external_force->clear(); this->applyBC( BC::Neumann::FromHigherDim(Matrix::eye(spatial_dimension, 1)), this->detector->getSurfaceId()); for (auto && tuple : zip(*nodal_area, make_view(*external_force, spatial_dimension))) { auto & area = std::get<0>(tuple); auto & force = std::get<1>(tuple); for (auto & f : force) area += pow(f, 2); area = sqrt(area); } this->external_force->clear(); } /* -------------------------------------------------------------------------- */ -void ContactMechanicsModel::beforeSolveStep() { this->search(); } +void ContactMechanicsModel::beforeSolveStep() { + this->search(); + +} /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::afterSolveStep() {} /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Contact Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << Model::spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + resolutions [" << std::endl; for (auto & resolution : resolutions) { resolution->printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ MatrixType ContactMechanicsModel::getMatrixType(const ID & matrix_id) { if (matrix_id == "K") return _symmetric; return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); if (!this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", getMatrixType("K")); } this->getDOFManager().clearMatrix("K"); for (auto & resolution : resolutions) { resolution->assembleStiffnessMatrix(_not_ghost); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleLumpedMatrix(const ID & /*matrix_id*/) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER dumper::Field * ContactMechanicsModel::createNodalFieldBool(const std::string & /*field_name*/, const std::string & /*group_name*/, bool /*padding_flag*/) { dumper::Field * field = nullptr; return field; } /* -------------------------------------------------------------------------- */ dumper::Field * ContactMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map *> real_nodal_fields; real_nodal_fields["contact_force"] = this->contact_force; real_nodal_fields["external_force"] = this->external_force; real_nodal_fields["blocked_dofs"] = this->blocked_dofs; real_nodal_fields["normals"] = this->normals; real_nodal_fields["gaps"] = this->gaps; real_nodal_fields["areas"] = this->nodal_area; dumper::Field * field = nullptr; if (padding_flag) field = this->mesh.createNodalField(real_nodal_fields[field_name], group_name, 3); else field = this->mesh.createNodalField(real_nodal_fields[field_name], group_name); return field; } #else /* -------------------------------------------------------------------------- */ dumper::Field * ContactMechanicsModel::createNodalFieldReal( __attribute__((unused)) const std::string & field_name, __attribute__((unused)) const std::string & group_name, __attribute__((unused)) bool padding_flag) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name) { mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name, UInt step) { mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(const std::string & dumper_name, Real time, UInt step) { mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump() { mesh.dump(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(UInt step) { mesh.dump(step); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::dump(Real time, UInt step) { mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ UInt ContactMechanicsModel::getNbData( const Array & elements, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; for (const Element & el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::packData(CommunicationBuffer & /*buffer*/, const Array & /*elements*/, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::unpackData(CommunicationBuffer & /*buffer*/, const Array & /*elements*/, const SynchronizationTag & /*tag*/) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt ContactMechanicsModel::getNbData( const Array & dofs, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); UInt size = 0; // UInt nb_nodes = mesh.getNbNodes(); AKANTU_DEBUG_OUT(); return size * dofs.size(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::packData(CommunicationBuffer & /*buffer*/, const Array & /*dofs*/, const SynchronizationTag & /*tag*/) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::unpackData(CommunicationBuffer & /*buffer*/, const Array & /*dofs*/, const SynchronizationTag & /*tag*/) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/contact_mechanics/resolution.cc b/src/model/contact_mechanics/resolution.cc index 69569c5d8..85e718cc3 100644 --- a/src/model/contact_mechanics/resolution.cc +++ b/src/model/contact_mechanics/resolution.cc @@ -1,350 +1,356 @@ /** * @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) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL); #undef GET_SHAPES_NATURAL #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ 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 surface_matrix(spatial_dimension - 1, spatial_dimension - 1); computeSurfaceMatrix(element.tangents, surface_matrix); 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, surface_matrix, element.gap); - //computeFrictionForce(fc, d_alpha, gap); 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(); 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 & 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); Vector shapes(nb_nodes_master); Matrix shapes_derivatives(spatial_dimension - 1, nb_nodes_master); #define GET_SHAPES_NATURAL(type) \ ElementClass::computeShapes(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) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_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); computeCoordinates(master, global_coords); computeTangents(shapes_derivatives, global_coords, tangents); Matrix surface_matrix(spatial_dimension - 1, spatial_dimension - 1); computeSurfaceMatrix(tangents, surface_matrix); 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, surface_matrix, gap); + Array t_alpha_beta(conn.size() * spatial_dimension, (spatial_dimension - 1) * (spatial_dimension -1)); + Array n_alpha_beta(conn.size() * spatial_dimension, (spatial_dimension - 1) * (spatial_dimension -1)); + Array p_alpha(conn.size() * spatial_dimension, spatial_dimension - 1); + + + + Matrix kc(connectivity.size() * spatial_dimension, connectivity.size() * spatial_dimension); computeTangentModuli(kc, n, n_alpha, d_alpha, surface_matrix, 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::computeSurfaceMatrix(Matrix & tangents, Matrix & surface_matrix) { surface_matrix.mul(tangents, tangents); surface_matrix = surface_matrix.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 & surface_matrix, Real & gap) { d_alpha.clear(); for (auto && entry : zip(surface_matrix.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::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 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); } } } } // akantu diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.cc b/src/model/contact_mechanics/resolutions/resolution_penalty.cc index ec2ef7117..18e3d19a2 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.cc +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.cc @@ -1,139 +1,154 @@ /** * @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, - "Penalty parameter"); + "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, Real & gap) { Vector tractions(d_alpha.getNbComponent()); computeFrictionalTraction(tractions); for (auto && values: zip(tractions, 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; } } -/* -------------------------------------------------------------------------- */ -void ResolutionPenalty::computeFrictionalTraction(Vector & tractions) { - -} /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeTangentModuli(Matrix & kc, Vector & n, Array & n_alpha, Array & d_alpha, Matrix & surface_matrix, Real & gap) { computeNormalStiffness(kc, n, n_alpha, d_alpha, surface_matrix, gap); computeFrictionalStiffness(n, n_alpha, d_alpha, gap); } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeNormalStiffness(Matrix & ke, Vector & n, Array & n_alpha, Array & d_alpha, Matrix & surface_matrix, 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::computeCommonModuli(Real & gap) { + + +} + +/* -------------------------------------------------------------------------- */ +void ResolutionPenalty::computeStickModuli() { + +} + +/* -------------------------------------------------------------------------- */ +void ResolutionPenalty::computeSlipModuli() { + } - 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 352b1bb86..21031c0f5 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.hh +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.hh @@ -1,105 +1,113 @@ /** * @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 & gap); /// local computation of stiffness matrix due to frictional stress void computeFrictionalStiffness(Vector & n, Array & n_alpha, Array & d_alpha, Real & gap); - /// computation of frictional tractions - void computeFrictionalTraction(Vector & tractions); + /// local computation of direct stiffness matrix due to friction, + /// this matrix is common for both stick and slip part + void computeCommonModuli(Real & gap); + + /// local computaion of stiffness matrix due to stick state + void computeStickModuli(); + + /// local computation of stiffness matrix due to slip state + void computeSlipModuli(); /* ------------------------------------------------------------------------ */ /* 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 normal force void computeNormalForce(Vector & /* force vector */, Vector & /* n */, Real & /* gap */) override; /// local computation of friction force void computeFrictionForce(Vector & /* force vector */, Array & /* D_alpha */, Real & /* gap */) override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: - /// penalty parameter + /// penalty parameter for normal stress Real epsilon; - + /// penalty parameter for tangential stress + Real epsilon_t; }; } // akantu #endif /* __AKANTU_RESOLUTION_PENALTY_HH__ */