diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index 860ac272e..2918f5af2 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,693 +1,789 @@ /** * @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" #include "group_manager_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactMechanicsModel::ContactMechanicsModel( Mesh & mesh, 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) { 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, id + ":contact_detector"); registerFEEngineObject("ContactFacetsFEEngine", mesh, Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ 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() { AKANTU_DEBUG_IN(); getFEEngine("ContactMechanicsModel").initShapeFunctions(_not_ghost); getFEEngine("ContactMechanicsModel").initShapeFunctions(_ghost); getFEEngine("ContactFacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("ContactFacetsFEEngine").initShapeFunctions(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ FEEngine & ContactMechanicsModel::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::initSolver( TimeStepSolverType /*time_step_solver_type*/, NonLinearSolverType) { // for alloc type of solvers this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->displacement_increment, spatial_dimension, "displacement_increment"); this->allocNodalField(this->internal_force, spatial_dimension, "internal_force"); this->allocNodalField(this->external_force, spatial_dimension, "external_force"); this->allocNodalField(this->normal_force, spatial_dimension, "normal_force"); this->allocNodalField(this->tangential_force, spatial_dimension, "tangential_force"); this->allocNodalField(this->gaps, 1, "gaps"); this->allocNodalField(this->nodal_area, 1, "areas"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); this->allocNodalField(this->contact_state, 1, "contact_state"); this->allocNodalField(this->previous_master_elements, 1, "previous_master_elements"); this->allocNodalField(this->normals, spatial_dimension, "normals"); auto surface_dimension = spatial_dimension - 1; this->allocNodalField(this->tangents, surface_dimension*spatial_dimension, "tangents"); this->allocNodalField(this->projections, surface_dimension, "projections"); this->allocNodalField(this->previous_projections, surface_dimension, "previous_projections"); this->allocNodalField(this->previous_tangents, surface_dimension*spatial_dimension, "previous_tangents"); this->allocNodalField(this->tangential_tractions, surface_dimension, "tangential_tractions"); this->allocNodalField(this->previous_tangential_tractions, surface_dimension, "previous_tangential_tractions"); // todo register multipliers as dofs for lagrange multipliers } /* -------------------------------------------------------------------------- */ std::tuple ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _explicit_consistent_mass: { return std::make_tuple("explicit", TimeStepSolverType::_dynamic); } case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } case _implicit_dynamic: { return std::make_tuple("implicit", TimeStepSolverType::_dynamic); } default: return std::make_tuple("unknown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions ContactMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson_contact; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_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(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the contact forces"); UInt nb_nodes = mesh.getNbNodes(); this->internal_force->clear(); this->normal_force->clear(); this->tangential_force->clear(); internal_force->resize(nb_nodes, 0.); normal_force->resize(nb_nodes, 0.); tangential_force->resize(nb_nodes, 0.); // assemble the forces due to contact 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); // TODO : uncomment when developing code for parallelization, // currently it addes the force twice for not ghost elements // hence source of error AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { // save the previous state this->savePreviousState(); contact_elements.clear(); contact_elements.resize(0); // this needs to be resized if cohesive elements are added UInt nb_nodes = mesh.getNbNodes(); auto resize_arrays = [&](auto & internal_array) { internal_array->clear(); internal_array->resize(nb_nodes, 0.); }; resize_arrays(gaps); resize_arrays(normals); resize_arrays(tangents); resize_arrays(projections); resize_arrays(tangential_tractions); resize_arrays(contact_state); resize_arrays(nodal_area); resize_arrays(external_force); this->detector->search(contact_elements, *gaps, *normals, *tangents, *projections); // intepenetration value must be positive for contact mechanics // model to work by default the gap value from detector is negative std::for_each((*gaps).begin(), (*gaps).end(), [](Real & gap){ gap *= -1.; }); if (contact_elements.size() != 0) { this->computeNodalAreas(); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::savePreviousState() { AKANTU_DEBUG_IN(); // saving previous natural projections previous_projections->clear(); previous_projections->resize(projections->size(), 0.); (*previous_projections).copy(*projections); // saving previous tangents previous_tangents->clear(); previous_tangents->resize(tangents->size(), 0.); (*previous_tangents).copy(*tangents); // saving previous tangential traction previous_tangential_tractions->clear(); previous_tangential_tractions->resize(tangential_tractions->size(), 0.); (*previous_tangential_tractions).copy(*tangential_tractions); previous_master_elements->clear(); previous_master_elements->resize(projections->size()); for (auto & element : contact_elements) { (*previous_master_elements)[element.slave] = element.master; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void ContactMechanicsModel::computeNodalAreas() { +void ContactMechanicsModel::computeNodalAreas(GhostType ghost_type) { UInt nb_nodes = mesh.getNbNodes(); nodal_area->clear(); external_force->clear(); nodal_area->resize(nb_nodes, 0.); external_force->resize(nb_nodes, 0.); auto & fem_boundary = getFEEngineClassBoundary("ContactMechanicsModel"); fem_boundary.initShapeFunctions(getPositions(), _not_ghost); fem_boundary.initShapeFunctions(getPositions(), _ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); - switch (spatial_dimension) { + /*switch (spatial_dimension) { case 1: { std::fill((*nodal_area).begin(), (*nodal_area).end(), 1.); break; } case 2: case 3: { this->applyBC( BC::Neumann::FromHigherDim(Matrix::eye(spatial_dimension, 1)), mesh.getElementGroup("contact_surface")); for (auto && tuple : zip(*nodal_area, make_view(*external_force, spatial_dimension))) { auto & area = std::get<0>(tuple); Vector force(std::get<1>(tuple)); area = force.norm(); } break; } default: break; } + this->external_force->clear();*/ + + IntegrationPoint quad_point; + quad_point.ghost_type = ghost_type; + + auto & group = mesh.getElementGroup("contact_surface"); + UInt nb_degree_of_freedom = external_force->getNbComponent(); + + for (auto && type : group.elementTypes(spatial_dimension - 1, ghost_type)) { + const auto & element_ids = group.getElements(type, ghost_type); + + UInt nb_quad_points = + fem_boundary.getNbIntegrationPoints(type, ghost_type); + UInt nb_elements = element_ids.size(); + + UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); + + Array dual_before_integ(nb_elements * nb_quad_points, + nb_degree_of_freedom, 0.); + Array quad_coords(nb_elements * nb_quad_points, spatial_dimension); + + const auto & normals_on_quad = + fem_boundary.getNormalsOnIntegrationPoints(type, ghost_type); + + auto normals_begin = normals_on_quad.begin(spatial_dimension); + decltype(normals_begin) normals_iter; + auto quad_coords_iter = quad_coords.begin(spatial_dimension); + auto dual_iter = dual_before_integ.begin(nb_degree_of_freedom); + + quad_point.type = type; + + Element subelement; + subelement.type = type; + subelement.ghost_type = ghost_type; + for (auto el : element_ids) { + subelement.element = el; + const auto & element_to_subelement = + mesh.getElementToSubelement(type)(el); + + Vector outside(spatial_dimension); + mesh.getBarycenter(subelement, outside); + + Vector inside(spatial_dimension); + if (mesh.isMeshFacets()) { + mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside); + } + else { + mesh.getBarycenter(element_to_subelement[0], inside); + } + + Vector inside_to_outside(spatial_dimension); + inside_to_outside = outside - inside; + + normals_iter = normals_begin + el * nb_quad_points; + + quad_point.element = el; + for (auto q : arange(nb_quad_points)) { + quad_point.num_point = q; + auto ddot = inside_to_outside.dot(*normals_iter); + Vector normal(*normals_iter); + if (ddot < 0) + normal *= -1.0; + + (*dual_iter).mul(Matrix::eye(spatial_dimension, 1), + normal); + ++dual_iter; + ++quad_coords_iter; + ++normals_iter; + } + } + + Array dual_by_shapes(nb_elements * nb_quad_points, + nb_degree_of_freedom * nb_nodes_per_element); + + fem_boundary.computeNtb(dual_before_integ, dual_by_shapes, type, + ghost_type, element_ids); + + Array dual_by_shapes_integ(nb_elements, nb_degree_of_freedom * + nb_nodes_per_element); + fem_boundary.integrate(dual_by_shapes, dual_by_shapes_integ, + nb_degree_of_freedom * nb_nodes_per_element, type, + ghost_type, element_ids); + + this->getDOFManager().assembleElementalArrayLocalArray( + dual_by_shapes_integ, *external_force, type, ghost_type, 1., element_ids); + } + + for (auto && tuple : + zip(*nodal_area, + make_view(*external_force, spatial_dimension))) { + + auto & area = std::get<0>(tuple); + Vector force(std::get<1>(tuple)); + area = force.norm(); + } + this->external_force->clear(); } /* -------------------------------------------------------------------------- */ 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")); } for (auto & resolution : resolutions) { resolution->assembleStiffnessMatrix(_not_ghost); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleLumpedMatrix(const ID & /*matrix_id*/) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::beforeSolveStep() { for (auto & resolution : resolutions) resolution->beforeSolveStep(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::afterSolveStep(bool converged) { for (auto & resolution : resolutions) resolution->afterSolveStep(converged); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr 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->internal_force; real_nodal_fields["normal_force"] = this->normal_force; real_nodal_fields["tangential_force"] = this->tangential_force; real_nodal_fields["blocked_dofs"] = this->blocked_dofs; real_nodal_fields["normals"] = this->normals; real_nodal_fields["tangents"] = this->tangents; real_nodal_fields["gaps"] = this->gaps; real_nodal_fields["areas"] = this->nodal_area; real_nodal_fields["contact_state"] = this->contact_state; real_nodal_fields["tangential_traction"] = this->tangential_tractions; std::shared_ptr field; 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 /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr ContactMechanicsModel::createNodalFieldReal(const std::string & , const std::string & , bool) { 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; 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/contact_mechanics_model.hh b/src/model/contact_mechanics/contact_mechanics_model.hh index 26677e0d1..30561e567 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.hh +++ b/src/model/contact_mechanics/contact_mechanics_model.hh @@ -1,383 +1,383 @@ /** * @file contact_mechanics_model.hh * * @author Mohit Pundir * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Model of Contact Mechanics * * @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 "boundary_condition.hh" #include "contact_detector.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_MECHANICS_MODEL_HH__ #define __AKANTU_CONTACT_MECHANICS_MODEL_HH__ namespace akantu { class Resolution; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class ContactMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ using MyFEEngineType = FEEngineTemplate; public: ContactMechanicsModel( Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "contact_mechanics_model", const MemoryID & memory_id = 0, std::shared_ptr dof_manager = nullptr, const ModelType model_type = ModelType::_contact_mechanics_model); ~ContactMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl(const ModelOptions & options) override; /// allocate all vectors void initSolver(TimeStepSolverType, NonLinearSolverType) override; /// initialize all internal arrays for resolutions void initResolutions(); /// initialize the modelType void initModel() override; /// call back for the solver, computes the force residual void assembleResidual() override; /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep(bool converted = true) override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /* ------------------------------------------------------------------------ */ /* Contact Detection */ /* ------------------------------------------------------------------------ */ public: void search(); - void computeNodalAreas(); + void computeNodalAreas(GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Contact Resolution */ /* ------------------------------------------------------------------------ */ public: /// register an empty contact resolution of a given type Resolution & registerNewResolution(const ID & res_name, const ID & res_type, const ID & opt_param); protected: /// register a resolution in the dynamic database Resolution & registerNewResolution(const ParserSection & res_section); /// read the resolution files to instantiate all the resolutions void instantiateResolutions(); /// save the parameters from previous state void savePreviousState(); /* ------------------------------------------------------------------------ */ /* Solver Interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the contact stiffness matrix virtual void assembleStiffnessMatrix(); /// assembles the contant internal forces virtual void assembleInternalForces(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: /// contact detection class friend class ContactDetector; /// contact resolution class friend class Resolution; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// get the ContactMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the ContactMechanicsModel::increment vector \warn only consistent /// if ContactMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *displacement_increment, Array &); /// get the ContactMechanics::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the ContactMechanicsModel::external_force vector (external forces) AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the ContactMechanics::normal_force vector (normal forces) AKANTU_GET_MACRO(NormalForce, *normal_force, Array &); /// get the ContactMechanics::tangential_force vector (friction forces) AKANTU_GET_MACRO(TangentialForce, *tangential_force, Array &); /// get the ContactMechanics::traction vector (friction traction) AKANTU_GET_MACRO(TangentialTractions, *tangential_tractions, Array &); /// get the ContactMechanics::previous_tangential_tractions vector AKANTU_GET_MACRO(PreviousTangentialTractions, *previous_tangential_tractions, Array &); /// get the ContactMechanicsModel::force vector (external forces) Array & getForce() { AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, " "use getExternalForce instead"); return *external_force; } /// get the ContactMechanics::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the ContactMechanics::gaps (contact gaps) AKANTU_GET_MACRO(Gaps, *gaps, Array &); /// get the ContactMechanics::normals (normals on slave nodes) AKANTU_GET_MACRO(Normals, *normals, Array &); /// get the ContactMechanics::tangents (tangents on slave nodes) AKANTU_GET_MACRO(Tangents, *tangents, Array &); /// get the ContactMechanics::previous_tangents (tangents on slave nodes) AKANTU_GET_MACRO(PreviousTangents, *previous_tangents, Array &); /// get the ContactMechanics::areas (nodal areas) AKANTU_GET_MACRO(NodalArea, *nodal_area, Array &); /// get the ContactMechanics::previous_projections (previous_projections) AKANTU_GET_MACRO(PreviousProjections, *previous_projections, Array &); /// get the ContactMechanics::projections (projections) AKANTU_GET_MACRO(Projections, *projections, Array &); /// get the ContactMechanics::contact_state vector (no_contact/stick/slip /// state) AKANTU_GET_MACRO(ContactState, *contact_state, Array &); /// get the ContactMechanics::previous_master_elements AKANTU_GET_MACRO(PreviousMasterElements, *previous_master_elements, Array &); /// get contact detector AKANTU_GET_MACRO_NOT_CONST(ContactDetector, *detector, ContactDetector &); /// get the contact elements inline Array & getContactElements() { return contact_elements; } /// get the current positions of the nodes inline Array & getPositions() { return detector->getPositions(); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// tells if the resolutions are instantiated bool are_resolutions_instantiated; /// displacements array Array * displacement{nullptr}; /// increment of displacement Array * displacement_increment{nullptr}; /// contact forces array Array * internal_force{nullptr}; /// external forces array Array * external_force{nullptr}; /// normal force array Array * normal_force{nullptr}; /// friction force array Array * tangential_force{nullptr}; /// friction traction array Array * tangential_tractions{nullptr}; /// previous friction traction array Array * previous_tangential_tractions{nullptr}; /// boundary vector Array * blocked_dofs{nullptr}; /// array to store gap between slave and master Array * gaps{nullptr}; /// array to store normals from master to slave Array * normals{nullptr}; /// array to store tangents on the master element Array * tangents{nullptr}; /// array to store previous tangents on the master element Array * previous_tangents{nullptr}; /// array to store nodal areas Array * nodal_area{nullptr}; /// array to store stick/slip state : Array * contact_state{nullptr}; /// array to store previous projections in covariant basis Array * previous_projections{nullptr}; // array to store projections in covariant basis Array * projections{nullptr}; /// contact detection std::unique_ptr detector; /// list of contact resolutions std::vector> resolutions; /// mapping between resolution name and resolution internal id std::map resolutions_names_to_id; /// array to store contact elements Array contact_elements; /// array to store previous master elements Array * previous_master_elements{nullptr}; }; } // namespace akantu /* ------------------------------------------------------------------------ */ /* inline functions */ /* ------------------------------------------------------------------------ */ #include "parser.hh" #include "resolution.hh" /* ------------------------------------------------------------------------ */ #endif /* __AKANTU_CONTACT_MECHANICS_MODEL_HH__ */ diff --git a/src/model/contact_mechanics/geometry_utils.cc b/src/model/contact_mechanics/geometry_utils.cc index 935942276..4751dfeca 100644 --- a/src/model/contact_mechanics/geometry_utils.cc +++ b/src/model/contact_mechanics/geometry_utils.cc @@ -1,859 +1,847 @@ /** * @file geometry_utils.cc * * @author Mohit Pundir * * @date creation: Mon Mmay 20 2019 * @date last modification: Mon May 20 2019 * * @brief Implementation of various utilities needed for contact geometry * * @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 "geometry_utils.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void GeometryUtils::normal(const Mesh & mesh, const Array & positions, const Element & element, Vector & normal, bool outward) { UInt spatial_dimension = mesh.getSpatialDimension(); UInt surface_dimension = spatial_dimension - 1; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); Matrix coords(spatial_dimension, nb_nodes_per_element); UInt * elem_val = mesh.getConnectivity(element.type, _not_ghost).storage(); mesh.extractNodalValuesFromElement(positions, coords.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); Matrix vectors(spatial_dimension, surface_dimension); switch (spatial_dimension) { case 1: { normal[0] = 1; break; } case 2: { vectors(0) = Vector(coords(1)) - Vector(coords(0)); Math::normal2(vectors.storage(), normal.storage()); break; } case 3: { vectors(0) = Vector(coords(1)) - Vector(coords(0)); vectors(1) = Vector(coords(2)) - Vector(coords(0)); Math::normal3(vectors(0).storage(), vectors(1).storage(), normal.storage()); break; } default: { AKANTU_ERROR("Unknown dimension : " << spatial_dimension); } } // to ensure that normal is always outwards from master surface if (outward) { const auto & element_to_subelement = mesh.getElementToSubelement(element.type)(element.element); Vector outside(spatial_dimension); mesh.getBarycenter(element, outside); // check if mesh facets exists for cohesive elements contact Vector inside(spatial_dimension); if (mesh.isMeshFacets()) { mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside); } else { mesh.getBarycenter(element_to_subelement[0], inside); } Vector inside_to_outside = outside - inside; auto projection = inside_to_outside.dot(normal); if (projection < 0) { normal *= -1.0; } } } /* -------------------------------------------------------------------------- */ void GeometryUtils::normal(const Mesh & mesh, const Element & element, Matrix & tangents, Vector & normal, bool outward) { UInt spatial_dimension = mesh.getSpatialDimension(); // to ensure that normal is always outwards from master surface we // compute a direction vector form inside of element attached to the // suurface element Vector inside_to_outside(spatial_dimension); if (outward) { const auto & element_to_subelement = mesh.getElementToSubelement(element.type)(element.element); Vector outside(spatial_dimension); mesh.getBarycenter(element, outside); // check if mesh facets exists for cohesive elements contact Vector inside(spatial_dimension); if (mesh.isMeshFacets()) { mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside); } else { mesh.getBarycenter(element_to_subelement[0], inside); } inside_to_outside = outside - inside; //auto projection = inside_to_outside.dot(normal); //if (projection < 0) { // normal *= -1.0; //} } // to ensure that direction of tangents are correct, cross product // of tangents should give be in the same direction as of inside to outside switch (spatial_dimension) { case 2: { - /*Vector e_z(3); - e_z[0] = 0.; - e_z[1] = 0.; - e_z[2] = 1.; - - Vector tangent(3); - tangent[0] = tangents(0, 0); - tangent[1] = tangents(0, 1); - tangent[2] = 0.; - - normal = e_z.crossProduct(tangent);*/ - normal[0] = -tangents(0, 1); normal[1] = tangents(0, 0); auto ddot = inside_to_outside.dot(normal); if (ddot < 0) { tangents *= -1.0; normal *= -1.0; } break; } case 3: { auto tang_trans = tangents.transpose(); auto tang1 = Vector(tang_trans(0)); auto tang2 = Vector(tang_trans(1)); auto tang1_cross_tang2 = tang1.crossProduct(tang2); normal = tang1_cross_tang2 / tang1_cross_tang2.norm(); auto ddot = inside_to_outside.dot(normal); if (ddot < 0) { tang_trans(1) *= -1.0; normal *= -1.0; } tangents = tang_trans.transpose(); break; } default: break; } } /* -------------------------------------------------------------------------- */ void GeometryUtils::covariantBasis(const Mesh & mesh, const Array & positions, const Element & element, const Vector & normal, Vector & natural_coord, Matrix & tangents) { UInt spatial_dimension = mesh.getSpatialDimension(); const ElementType & type = element.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(positions, nodes_coord.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); UInt surface_dimension = spatial_dimension - 1; Matrix dnds(surface_dimension, nb_nodes_per_element); #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(natural_coord, dnds) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL tangents.mul(dnds, nodes_coord); auto temp_tangents = tangents.transpose(); for (UInt i = 0; i < spatial_dimension - 1; ++i) { auto temp = Vector(temp_tangents(i)); temp_tangents(i) = temp.normalize(); } tangents = temp_tangents.transpose(); // to ensure that direction of tangents are correct, cross product // of tangents should give the normal vector computed earlier switch (spatial_dimension) { case 2: { Vector e_z(3); e_z[0] = 0.; e_z[1] = 0.; e_z[2] = 1.; Vector tangent(3); tangent[0] = tangents(0, 0); tangent[1] = tangents(0, 1); tangent[2] = 0.; auto exp_normal = e_z.crossProduct(tangent); auto ddot = normal.dot(exp_normal); if (ddot < 0) { tangents *= -1.0; } break; } case 3: { auto tang_trans = tangents.transpose(); auto tang1 = Vector(tang_trans(0)); auto tang2 = Vector(tang_trans(1)); auto tang1_cross_tang2 = tang1.crossProduct(tang2); auto exp_normal = tang1_cross_tang2 / tang1_cross_tang2.norm(); auto ddot = normal.dot(exp_normal); if (ddot < 0) { tang_trans(1) *= -1.0; } tangents = tang_trans.transpose(); break; } default: break; } } /* -------------------------------------------------------------------------- */ void GeometryUtils::covariantBasis(const Mesh & mesh, const Array & positions, const Element & element, Vector & natural_coord, Matrix & tangents) { UInt spatial_dimension = mesh.getSpatialDimension(); const ElementType & type = element.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(positions, nodes_coord.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); UInt surface_dimension = spatial_dimension - 1; Matrix dnds(surface_dimension, nb_nodes_per_element); #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(natural_coord, dnds) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL tangents.mul(dnds, nodes_coord); auto temp_tangents = tangents.transpose(); for (UInt i = 0; i < spatial_dimension - 1; ++i) { auto temp = Vector(temp_tangents(i)); temp_tangents(i) = temp.normalize(); } tangents = temp_tangents.transpose(); } /* -------------------------------------------------------------------------- */ void GeometryUtils::curvature(const Mesh & mesh, const Array & positions, const Element & element, const Vector & natural_coord, Matrix & curvature) { UInt spatial_dimension = mesh.getSpatialDimension(); auto surface_dimension = spatial_dimension - 1; const ElementType & type = element.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); Matrix dn2ds2(surface_dimension * surface_dimension, Mesh::getNbNodesPerElement(type)); #define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \ ElementClass::computeDN2DS2(natural_coord, dn2ds2) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL); #undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL Matrix coords(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(positions, coords.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); curvature.mul(coords, dn2ds2); } /* -------------------------------------------------------------------------- */ UInt GeometryUtils::orthogonalProjection(const Mesh & mesh, const Array & positions, const Vector & slave, const Array & elements, Real & gap, Vector & natural_projection, Vector & normal, Real alpha, UInt max_iterations, Real projection_tolerance, Real extension_tolerance) { UInt index = UInt(-1); Real min_gap = std::numeric_limits::max(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt surface_dimension = spatial_dimension - 1; UInt nb_same_sides{0}; UInt nb_boundary_elements{0}; UInt counter{0}; auto & contact_group = mesh.getElementGroup("contact_surface"); for (auto & element : elements) { // filter out elements which are not there in the element group // contact surface created by the surface selector and is stored // in the mesh or mesh_facet, if a element is not there it // returnas UInt(-1) auto & elements_of_type = contact_group.getElements(element.type); if(elements_of_type.find(element.element) == UInt(-1)) continue; //if (!GeometryUtils::isBoundaryElement(mesh, element)) // continue; nb_boundary_elements++; // find the natural coordinate corresponding to the minimum gap // between slave node and master element /* Vector normal_ele(spatial_dimension); GeometryUtils::normal(mesh, positions, element, normal_ele); Vector master(spatial_dimension); GeometryUtils::realProjection(mesh, positions, slave, element, normal_ele, master); Vector xi(natural_projection.size()); GeometryUtils::naturalProjection(mesh, positions, element, master, xi, max_iterations, projection_tolerance);*/ Vector master(spatial_dimension); Vector xi(natural_projection.size()); GeometryUtils::naturalProjection(mesh, positions, element, slave, master, xi, max_iterations, projection_tolerance); Matrix tangent_ele(surface_dimension, spatial_dimension); GeometryUtils::covariantBasis(mesh, positions, element, xi, tangent_ele); Vector normal_ele(spatial_dimension); GeometryUtils::normal(mesh, element, tangent_ele, normal_ele); // if gap between master projection and slave point is zero, then // it means that slave point lies on the master element, hence the // normal from master to slave cannot be computed in that case auto master_to_slave = slave - master; Real temp_gap = master_to_slave.norm(); if (temp_gap != 0) master_to_slave /= temp_gap; // also the slave point should lie inside the master surface in // case of explicit or outside in case of implicit, one way to // check that is by checking the dot product of normal at each // master element, if the values of all dot product is same then // the slave point lies on the same side of each master element // A alpha parameter is introduced which is 1 in case of explicit // and -1 in case of implicit, therefor the variation (dot product // + alpha) should be close to zero (within tolerance) for both // cases Real direction_tolerance = 1e-8; auto product = master_to_slave.dot(normal_ele); auto variation = std::abs(product + alpha); if (variation <= direction_tolerance and temp_gap <= min_gap and GeometryUtils::isValidProjection(xi, extension_tolerance)) { gap = -temp_gap; min_gap = temp_gap; index = counter; natural_projection = xi; normal = normal_ele; } if(temp_gap == 0 or variation <= direction_tolerance) nb_same_sides++; counter++; } // if point is not on the same side of all the boundary elements // than it is not consider even if the closet master element is // found if(nb_same_sides != nb_boundary_elements) index = UInt(-1); return index; } /* -------------------------------------------------------------------------- */ UInt GeometryUtils::orthogonalProjection(const Mesh & mesh, const Array & positions, const Vector & slave, const Array & elements, Real & gap, Vector & natural_projection, Vector & normal, Matrix & tangent, Real alpha, UInt max_iterations, Real projection_tolerance, Real extension_tolerance) { UInt index = UInt(-1); Real min_gap = std::numeric_limits::max(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt surface_dimension = spatial_dimension - 1; UInt nb_same_sides{0}; UInt nb_boundary_elements{0}; auto & contact_group = mesh.getElementGroup("contact_surface"); UInt counter{0}; for (auto & element : elements) { // filter out elements which are not there in the element group // contact surface created by the surface selector and is stored // in the mesh or mesh_facet, if a element is not there it // returnas UInt(-1) auto & elements_of_type = contact_group.getElements(element.type); if(elements_of_type.find(element.element) == UInt(-1)) continue; //if (!GeometryUtils::isBoundaryElement(mesh, element)) // continue; nb_boundary_elements++; // find the natural coordinate corresponding to the minimum gap // between slave node and master element /* Vector normal_ele(spatial_dimension); GeometryUtils::normal(mesh, positions, element, normal_ele); Vector master(spatial_dimension); GeometryUtils::realProjection(mesh, positions, slave, element, normal_ele, master); Vector xi(natural_projection.size()); GeometryUtils::naturalProjection(mesh, positions, element, master, xi, max_iterations, projection_tolerance);*/ Vector master(spatial_dimension); Vector xi_ele(natural_projection.size()); GeometryUtils::naturalProjection(mesh, positions, element, slave, master, xi_ele, max_iterations, projection_tolerance); Matrix tangent_ele(surface_dimension, spatial_dimension); GeometryUtils::covariantBasis(mesh, positions, element, xi_ele, tangent_ele); Vector normal_ele(spatial_dimension); GeometryUtils::normal(mesh, element, tangent_ele, normal_ele); // if gap between master projection and slave point is zero, then // it means that slave point lies on the master element, hence the // normal from master to slave cannot be computed in that case auto master_to_slave = slave - master; Real temp_gap = master_to_slave.norm(); if (temp_gap != 0) master_to_slave /= temp_gap; // also the slave point should lie inside the master surface in // case of explicit or outside in case of implicit, one way to // check that is by checking the dot product of normal at each // master element, if the values of all dot product is same then // the slave point lies on the same side of each master element // A alpha parameter is introduced which is 1 in case of explicit // and -1 in case of implicit, therefor the variation (dot product // + alpha) should be close to zero (within tolerance) for both // cases Real direction_tolerance = 1e-8; auto product = master_to_slave.dot(normal_ele); auto variation = std::abs(product + alpha); if (variation <= direction_tolerance and temp_gap <= min_gap and GeometryUtils::isValidProjection(xi_ele, extension_tolerance)) { gap = -temp_gap; min_gap = temp_gap; index = counter; natural_projection = xi_ele; normal = normal_ele; tangent = tangent_ele; } if(temp_gap == 0 or variation <= direction_tolerance) nb_same_sides++; counter++; } // if point is not on the same side of all the boundary elements // than it is not consider even if the closet master element is // found if(nb_same_sides != nb_boundary_elements) index = UInt(-1); return index; } /* -------------------------------------------------------------------------- */ void GeometryUtils::realProjection(const Mesh & mesh, const Array & positions, const Vector & slave, const Element & element, const Vector & normal, Vector & projection) { UInt spatial_dimension = mesh.getSpatialDimension(); const ElementType & type = element.type; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(positions, nodes_coord.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); Vector point(nodes_coord(0)); Real alpha = (slave - point).dot(normal); projection = slave - alpha * normal; } /* -------------------------------------------------------------------------- */ void GeometryUtils::realProjection(const Mesh & mesh, const Array & positions, const Element & element, const Vector & natural_coord, Vector & projection) { UInt spatial_dimension = mesh.getSpatialDimension(); const ElementType & type = element.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element.type); Vector shapes(nb_nodes_per_element); #define GET_SHAPE_NATURAL(type) \ ElementClass::computeShapes(natural_coord, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); #undef GET_SHAPE_NATURAL Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); mesh.extractNodalValuesFromElement(positions, nodes_coord.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); projection.mul(nodes_coord, shapes); } /* -------------------------------------------------------------------------- */ void GeometryUtils::naturalProjection(const Mesh & mesh, const Array & positions, const Element & element, const Vector & slave_coords, Vector & master_coords, Vector & natural_projection, UInt max_iterations, Real projection_tolerance) { UInt spatial_dimension = mesh.getSpatialDimension(); UInt surface_dimension = spatial_dimension - 1; const ElementType & type = element.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(positions, nodes_coord.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); // initial guess natural_projection.clear(); // obhjective function computed on the natural_guess Matrix f(surface_dimension, 1); // jacobian matrix computed on the natural_guess Matrix J(surface_dimension, surface_dimension); // Jinv = J^{-1} Matrix Jinv(surface_dimension, surface_dimension); // dxi = \xi_{k+1} - \xi_{k} in the iterative process Matrix dxi(surface_dimension, 1); // gradient at natural projection Matrix gradient(surface_dimension, spatial_dimension); // second derivative at natural peojection Matrix double_gradient(surface_dimension, surface_dimension); // second derivative of shape function at natural projection Matrix dn2ds2(surface_dimension * surface_dimension, nb_nodes_per_element); auto compute_double_gradient = [&double_gradient, &dn2ds2, &nodes_coord, surface_dimension, spatial_dimension](UInt & alpha, UInt & beta) { auto index = alpha * surface_dimension + beta; Vector d_alpha_beta(spatial_dimension); auto dn2ds2_transpose = dn2ds2.transpose(); Vector dn2ds2_alpha_beta(dn2ds2_transpose(index)); - d_alpha_beta.mul(nodes_coord, dn2ds2_alpha_beta); + d_alpha_beta.mul(nodes_coord, dn2ds2_alpha_beta); return d_alpha_beta; }; /* --------------------------- */ /* init before iteration loop */ /* --------------------------- */ // do interpolation auto update_f = [&f, &master_coords, &natural_projection, &nodes_coord, &slave_coords, &gradient, surface_dimension, spatial_dimension, nb_nodes_per_element, type]() { // compute real coords on natural projection Vector shapes(nb_nodes_per_element); #define GET_SHAPE_NATURAL(type) \ ElementClass::computeShapes(natural_projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); #undef GET_SHAPE_NATURAL master_coords.mul(nodes_coord, shapes); auto distance = slave_coords - master_coords; // first derivative of shape function at natural projection Matrix dnds(surface_dimension, nb_nodes_per_element); // computing gradient at projection point #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(natural_projection, dnds) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL gradient.mul(dnds, nodes_coord); // gradient transpose at natural projection Matrix gradient_transpose(surface_dimension, spatial_dimension); gradient_transpose = gradient.transpose(); // loop over surface dimensions for(auto alpha : arange(surface_dimension)) { Vector gradient_alpha(gradient_transpose(alpha)); f(alpha, 0) = -2.*gradient_alpha.dot(distance); } // compute initial error auto error = f.norm(); return error; }; auto projection_error = update_f(); /* --------------------------- */ /* iteration loop */ /* --------------------------- */ UInt iterations{0}; while(projection_tolerance < projection_error and iterations < max_iterations) { // compute covariant components of metric tensor auto a = GeometryUtils::covariantMetricTensor(gradient); // computing second derivative at natural projection #define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \ ElementClass::computeDN2DS2(natural_projection, dn2ds2) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL); #undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL // real coord - physical guess auto distance = slave_coords - master_coords; // computing Jacobian J for(auto alpha : arange(surface_dimension)) { for(auto beta : arange(surface_dimension)) { auto dgrad_alpha_beta = compute_double_gradient(alpha, beta); J(alpha, beta) = 2.*(a(alpha, beta) - dgrad_alpha_beta.dot(distance)); } } Jinv.inverse(J); // compute increment dxi.mul(Jinv, f, -1.0); // update our guess natural_projection += Vector(dxi(0)); projection_error = update_f(); iterations++; } /* #define GET_NATURAL_COORDINATE(type) \ ElementClass::inverseMap(slave_coords, nodes_coord, \ natural_projection, max_iterations, projection_tolerance) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NATURAL_COORDINATE); #undef GET_NATURAL_COORDINATE*/ } /* -------------------------------------------------------------------------- */ void GeometryUtils::contravariantBasis(const Matrix & covariant, Matrix & contravariant) { auto inv_A = GeometryUtils::contravariantMetricTensor(covariant); contravariant.mul(inv_A, covariant); } /* -------------------------------------------------------------------------- */ Matrix GeometryUtils::covariantMetricTensor(const Matrix & covariant_bases) { Matrix A(covariant_bases.rows(), covariant_bases.rows()); A.mul(covariant_bases, covariant_bases); return A; } /* -------------------------------------------------------------------------- */ Matrix GeometryUtils::contravariantMetricTensor(const Matrix & covariant_bases) { auto A = GeometryUtils::covariantMetricTensor(covariant_bases); auto inv_A = A.inverse(); return inv_A; } /* -------------------------------------------------------------------------- */ Matrix GeometryUtils::covariantCurvatureTensor(const Mesh & mesh, const Array & positions, const Element & element, const Vector & natural_coord, const Vector & normal) { UInt spatial_dimension = mesh.getSpatialDimension(); auto surface_dimension = spatial_dimension - 1; const ElementType & type = element.type; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); Matrix dn2ds2(surface_dimension * surface_dimension, nb_nodes_per_element); #define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \ ElementClass::computeDN2DS2(natural_coord, dn2ds2) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL); #undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL Matrix coords(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(positions, coords.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); Matrix curvature(spatial_dimension, surface_dimension*surface_dimension); curvature.mul(coords, dn2ds2); Matrix H(surface_dimension, surface_dimension); UInt i = 0; for (auto alpha : arange(surface_dimension)) { for (auto beta : arange(surface_dimension)) { Vector temp(curvature(i)); H(alpha, beta) = temp.dot(normal); i++; } } return H; } }