diff --git a/examples/contact_mechanics/CMakeLists.txt b/examples/contact_mechanics/CMakeLists.txt index 816ca7ad4..6d6d6fa35 100644 --- a/examples/contact_mechanics/CMakeLists.txt +++ b/examples/contact_mechanics/CMakeLists.txt @@ -1,51 +1,39 @@ #=============================================================================== # @file CMakeLists.txt # # @author Mohit Pundir # # @date creation: Mon Jan 18 2016 # # @brief configuration for heat transfer example # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 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 . # # @section DESCRIPTION # #=============================================================================== -add_mesh(contact_hertz_2d contact_hertz_2d.geo 2 1) +add_mesh(hertz hertz.geo 2 1) register_example(contact_explicit_dynamic SOURCES contact_explicit_dynamic.cc - DEPENDS contact_hertz_2d - FILES_TO_COPY material.dat - ) - -register_example(solid_explicit_dynamic - SOURCES solid_explicit_dynamic.cc - DEPENDS contact_hertz_2d - FILES_TO_COPY material_solid.dat - ) - -register_example(node_selector - SOURCES test_node_selector.cc - DEPENDS contact_hertz_2d + DEPENDS hertz FILES_TO_COPY material.dat ) diff --git a/examples/contact_mechanics/contact_explicit_dynamic.cc b/examples/contact_mechanics/contact_explicit_dynamic.cc index 91023ec3b..74361dfa2 100644 --- a/examples/contact_mechanics/contact_explicit_dynamic.cc +++ b/examples/contact_mechanics/contact_explicit_dynamic.cc @@ -1,123 +1,123 @@ /** * @file test_contact_mechanics_model.cc * * @author Mohit Pundir * * @date creation: Tue Apr 30 2019 * @date last modification: Tue Apr 30 2019 * * @brief Test for contact mechanics model 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 "solid_mechanics_model.hh" #include "contact_mechanics_model.hh" #include "coupler_solid_contact.hh" #include "non_linear_solver.hh" +#include "surface_selector.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { const UInt spatial_dimension = 2; initialize("material.dat", argc, argv); Real time_step; Real time_factor = 0.8; UInt max_steps = 5000; Mesh mesh(spatial_dimension); - mesh.read("contact_hertz_2d.msh"); + mesh.read("hertz.msh"); CouplerSolidContact coupler(mesh); auto & solid = coupler.getSolidMechanicsModel(); auto & contact = coupler.getContactMechanicsModel(); auto && selector = std::make_shared>( "physical_names",solid); solid.setMaterialSelector(selector); solid.initFull( _analysis_method = _explicit_lumped_mass); contact.initFull(_analysis_method = _explicit_dynamic_contact); - + auto && surface_selector = std::make_shared(mesh); + contact.getContactDetector().setSurfaceSelector(surface_selector); + solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "bottom"); solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "bottom"); - solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "top"); - coupler.initFull(_analysis_method = _explicit_dynamic_contact); time_step = solid.getStableTimeStep(); std::cout << "Time Step = " << time_step * time_factor << "s (" << time_step << "s)" << std::endl; coupler.setTimeStep(time_step * time_factor); coupler.setBaseName("contact-explicit-dynamic"); coupler.addDumpFieldVector("displacement"); coupler.addDumpFieldVector("velocity"); coupler.addDumpFieldVector("normals"); - coupler.addDumpFieldVector("tangents"); coupler.addDumpFieldVector("contact_force"); coupler.addDumpFieldVector("external_force"); coupler.addDumpFieldVector("internal_force"); coupler.addDumpField("gaps"); coupler.addDumpField("areas"); coupler.addDumpField("blocked_dofs"); coupler.addDumpField("grad_u"); coupler.addDumpField("stress"); auto & velocity = solid.getVelocity(); Real damping_ratio = 0.7; for (auto i : arange(max_steps)) { Real factor = 5e-4; auto increment = time_step * factor; solid.applyBC(BC::Dirichlet::IncrementValue(-increment, _y), "top"); coupler.solveStep(); Real epot = solid.getEnergy("potential"); Real ekin = solid.getEnergy("kinetic"); std::cerr << i << "," << i * increment << "," << epot << "," << ekin << "," << epot + ekin << "," << std::endl; for (auto & v : make_view(velocity)) { v *= damping_ratio; } coupler.dump(); } finalize(); return EXIT_SUCCESS; } diff --git a/examples/contact_mechanics/contact_hertz_2d.geo b/examples/contact_mechanics/contact_hertz_2d.geo deleted file mode 100644 index 76e1e5abd..000000000 --- a/examples/contact_mechanics/contact_hertz_2d.geo +++ /dev/null @@ -1,33 +0,0 @@ -cl1 = 0.002; -cl2 = 0.0025; -cl3 = 0.0025; -Dy = 0.0; -radius = 0.1; -y = 0.1; -epsilon = -1e-4; -Point(1) = {0, y, 0, cl1}; -Point(2) = {radius, radius + y, 0, 5*cl2}; -Point(3) = {-radius, radius + y, 0, 5*cl2}; -Point(4) = {0.25, 0.1 + epsilon, 0, cl3}; -Point(5) = {-0.25, 0.1 + epsilon, 0, cl3}; -Point(6) = {-0.25, -0.1, 0, 10*cl2}; -Point(7) = {0.25, -0.1, 0, 10*cl2}; -Point(8) = {0, radius + y, 0, cl2}; -Circle(1) = {3, 8, 1}; -Circle(2) = {1, 8, 2}; -Line(3) = {2, 8}; -Line(13) = {8, 3}; -Line(4) = {6, 7}; -Line(5) = {7, 4}; -Line(6) = {4, 5}; -Line(7) = {5, 6}; -Line Loop(9) = {2, 3, 13, 1}; -Plane Surface(9) = {9}; -Line Loop(11) = {6, 7, 4, 5}; -Plane Surface(11) = {11}; -Physical Line("contact_bottom") = {6}; -Physical Line("contact_top") = {1, 2}; -Physical Line("top") = {3, 13}; -Physical Line("bottom") = {4}; -Physical Surface("top_body") = {9}; -Physical Surface("bot_body") = {11}; \ No newline at end of file diff --git a/examples/contact_mechanics/explicit.geo b/examples/contact_mechanics/explicit.geo deleted file mode 100644 index e2b8897e0..000000000 --- a/examples/contact_mechanics/explicit.geo +++ /dev/null @@ -1,51 +0,0 @@ -cl1 = 0.025; -cl2 = 0.025; -cl3 = 0.025; -Dy = 0.0; -radius = 0.1; -y = 0.1; -epsilon = 0.005; - -Point(1) = {0, y, 0, cl1}; -Point(2) = {radius, radius + y, 0, cl2}; -Point(3) = {-radius, radius + y, 0, cl2}; -Point(4) = {0.2, 0.1 + epsilon, 0, cl3}; -Point(5) = {-0.2, 0.1 + epsilon, 0, cl3}; -Point(6) = {-0.2, -0.15, 0, 10*cl3}; -Point(7) = {0.2, -0.15, 0, 10*cl3}; -Point(8) = {0, radius + y, 0, cl2}; -Point(9) = {0, y + epsilon, 0, cl3}; -Point(10) = {0, y-0.025, 0, cl3*0.5}; -Point(11) = {radius/3, y-0.0075, 0, cl3}; -Point(12) = {-radius/3, y-0.0075, 0, cl3}; - -Circle(1) = {3, 8, 1}; -Circle(2) = {1, 8, 2}; - -Line(3) = {2, 8}; -Line(13) = {8, 3}; -Line(4) = {6, 7}; -Line(5) = {7, 4}; -Line(6) = {4, 9}; -Line(7) = {5, 6}; -Line(8) = {9, 5}; - -Line Loop(9) = {2, 3, 13, 1}; -Plane Surface(9) = {9}; - -Line Loop(11) = {6, 8, 7, 4, 5}; -Plane Surface(11) = {11}; - -Point{10} In Surface{11}; -Point{11} In Surface{11}; -Point{12} In Surface{11}; - - -Physical Line("contact_bottom") = {6, 8}; -Physical Line("contact_top") = {1, 2}; - -Physical Line("top") = {3, 13}; -Physical Line("bottom") = {4}; - -Physical Surface("top_body") = {9}; -Physical Surface("bot_body") = {11}; \ No newline at end of file diff --git a/examples/contact_mechanics/explicit_detection.cc b/examples/contact_mechanics/explicit_detection.cc deleted file mode 100644 index 2dcede8ae..000000000 --- a/examples/contact_mechanics/explicit_detection.cc +++ /dev/null @@ -1,61 +0,0 @@ -/* -------------------------------------------------------------------------- */ -#include "contact_mechanics_model.hh" - -#include "dumper_elemental_field.hh" -#include "dumper_nodal_field.hh" - -#include "dumper_iohelper_paraview.hh" -/* -------------------------------------------------------------------------- */ - -using namespace akantu; - -/* -------------------------------------------------------------------------- */ -int main(int argc, char *argv[]) { - - std::string mesh_file = "explicit.msh"; - std::string material_file = "options_explicit.dat"; - - const UInt spatial_dimension = 2; - initialize(material_file, argc, argv); - - Mesh mesh(spatial_dimension); - mesh.read(mesh_file); - - ContactDetector detector(mesh); - - auto &&surface_selector = std::make_shared(mesh); - detector.setSurfaceSelector(surface_selector); - - std::map contact_map; - - detector.search(contact_map); - - UInt nb_nodes = mesh.getNbNodes(); - - Array gaps(nb_nodes, 1); - Array normals(nb_nodes, spatial_dimension); - - for (auto & entry: contact_map) { - const auto & node = entry.first; - const auto & element = entry.second; - - gaps[node] = element.gap; - for (auto i : arange(spatial_dimension)) { - normals(node, i) = element.normal[i]; - } - } - - DumperParaview dumper("explicit_detection", "./paraview", false); - dumper.registerMesh(mesh); - - auto gap_field = std::make_shared>(gaps); - auto normal_field = std::make_shared>(normals); - - dumper.registerField("gaps", gap_field); - dumper.registerField("normals", normal_field); - dumper.dump(); - - finalize(); - return EXIT_SUCCESS; -} - diff --git a/examples/contact_mechanics/explicit_dynamic_implicit_detection.cc b/examples/contact_mechanics/explicit_dynamic_implicit_detection.cc deleted file mode 100644 index 307024a8f..000000000 --- a/examples/contact_mechanics/explicit_dynamic_implicit_detection.cc +++ /dev/null @@ -1,102 +0,0 @@ -/* -------------------------------------------------------------------------- */ -#include "solid_mechanics_model.hh" -#include "contact_mechanics_model.hh" -#include "coupler_solid_contact.hh" -#include "non_linear_solver.hh" -#include "dumper_text.hh" -#include "dumper_variable.hh" -/* -------------------------------------------------------------------------- */ - -using namespace akantu; - -/* -------------------------------------------------------------------------- */ -int main(int argc, char *argv[]) { - - UInt max_steps = 50000; - UInt imposing_steps = 5000; - UInt damping_interval = 10; - Real max_weight = 0.004; - Real max_displacement = 0.01; - Real damping_ratio = 0.9; - std::string mesh_file = "hertz.msh"; - std::string material_file = "material_explicit.dat"; - - const UInt spatial_dimension = 2; - initialize(material_file, argc, argv); - - Real time_factor = 0.1; - - Mesh mesh(spatial_dimension); - mesh.read(mesh_file); - - CouplerSolidContact coupler(mesh); - - auto & solid = coupler.getSolidMechanicsModel(); - auto & contact = coupler.getContactMechanicsModel(); - - auto && material_selector = std::make_shared>( - "physical_names",solid); - solid.setMaterialSelector(material_selector); - - solid.initFull( _analysis_method = _explicit_lumped_mass); - contact.initFull(_analysis_method = _explicit_dynamic_contact); - - auto &&surface_selector = std::make_shared(mesh); - contact.getContactDetector().setSurfaceSelector(surface_selector); - - solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "bottom"); - solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "bottom"); - - //Vector weight = {0, -max_weight}; - //solid.applyBC(BC::Neumann::FromSameDim(weight), "top"); - - coupler.initFull(_analysis_method = _explicit_dynamic_contact); - - Real time_step = solid.getStableTimeStep(); - time_step *= time_factor; - - std::cout << "Time Step = " << time_step << "s (" << time_step - << "s)" << std::endl; - coupler.setTimeStep(time_step); - - coupler.setBaseName("explicit-dynamic-implicit"); - coupler.addDumpFieldVector("displacement"); - coupler.addDumpFieldVector("normals"); - coupler.addDumpFieldVector("tangents"); - coupler.addDumpFieldVector("contact_force"); - coupler.addDumpFieldVector("external_force"); - coupler.addDumpFieldVector("internal_force"); - coupler.addDumpField("gaps"); - coupler.addDumpField("areas"); - coupler.addDumpField("blocked_dofs"); - coupler.addDumpField("strain"); - coupler.addDumpField("stress"); - coupler.addDumpField("Von Mises stress"); - - auto & velocity = solid.getVelocity(); - - for (UInt s : arange(max_steps)) { - - std::cerr << "Step " << s << std::endl; - - Real increment = max_displacement / Real(imposing_steps / 2); - solid.applyBC(BC::Dirichlet::IncrementValue(-increment, _y), "top"); - - coupler.solveStep(); - - if (s % damping_interval == 0) { - for (auto & v : make_view(velocity)) - v *= damping_ratio; - } - - if (s % 100 == 0) { - coupler.dump(); - } - - } - - - finalize(); - return EXIT_SUCCESS; -} - diff --git a/examples/contact_mechanics/hertz.geo b/examples/contact_mechanics/hertz.geo index b2504e436..176cf0401 100644 --- a/examples/contact_mechanics/hertz.geo +++ b/examples/contact_mechanics/hertz.geo @@ -1,51 +1,50 @@ cl1 = 0.005; cl2 = 0.005; -cl3 = 0.0025; +cl3 = 0.001; Dy = 0.0; radius = 0.1; y = 0.1; -epsilon = 0; +epsilon = 1e-3; Point(1) = {0, y, 0, cl1}; Point(2) = {radius, radius + y, 0, cl2}; Point(3) = {-radius, radius + y, 0, cl2}; Point(4) = {0.2, 0.1 + epsilon, 0, 0.75*cl3}; Point(5) = {-0.2, 0.1 + epsilon, 0, 0.75*cl3}; Point(6) = {-0.2, -0.15, 0, 10*cl3}; Point(7) = {0.2, -0.15, 0, 10*cl3}; Point(8) = {0, radius + y, 0, cl2}; Point(9) = {0, y, 0, cl3*0.5}; Point(10) = {0, y-0.025, 0, cl3*0.5}; Point(11) = {radius/3, y-0.0075, 0, cl3*0.5}; Point(12) = {-radius/3, y-0.0075, 0, cl3*0.5}; Circle(1) = {3, 8, 1}; Circle(2) = {1, 8, 2}; Line(3) = {2, 8}; Line(13) = {8, 3}; Line(4) = {6, 7}; Line(5) = {7, 4}; Line(6) = {4, 9}; Line(7) = {5, 6}; Line(8) = {9, 5}; Line Loop(9) = {2, 3, 13, 1}; Plane Surface(9) = {9}; Line Loop(11) = {6, 8, 7, 4, 5}; Plane Surface(11) = {11}; Point{10} In Surface{11}; Point{11} In Surface{11}; Point{12} In Surface{11}; - Physical Line("contact_bottom") = {6, 8}; Physical Line("contact_top") = {1, 2}; Physical Line("top") = {3, 13}; Physical Line("bottom") = {4}; -Physical Surface("top_body") = {9}; -Physical Surface("bot_body") = {11}; \ No newline at end of file +Physical Surface("uppper") = {9}; +Physical Surface("lower") = {11}; \ No newline at end of file diff --git a/examples/contact_mechanics/implicit.geo b/examples/contact_mechanics/implicit.geo deleted file mode 100644 index be5487ac9..000000000 --- a/examples/contact_mechanics/implicit.geo +++ /dev/null @@ -1,40 +0,0 @@ -h = 1; -height = 1.0; -length = 1.0; -epsilon = 1e-3; - -Point(1) = {0, 0, 0, h*0.2}; -Point(2) = {length, 0, 0, h*0.2}; -Point(3) = {length, height/2-epsilon, 0, h*0.2}; -Point(4) = {length, height, 0, h*0.2}; -Point(5) = {0, height, 0, h*0.2}; -Point(6) = {0, height/2+epsilon, 0, h*0.2}; -Point(7) = {length, height/2+epsilon, 0, h*0.2}; -Point(8) = {0, height/2-epsilon, 0, h*0.2}; - -Line(1) = {1, 2}; -Line(2) = {2, 3}; -Line(3) = {7, 4}; -Line(4) = {4, 5}; -Line(5) = {5, 6}; -Line(6) = {6, 7}; -Line(7) = {8, 1}; -Line(8) = {3, 8}; - -Line Loop(1) = {1, 2, 8, 7}; -Line Loop(2) = {3, 4, 5, 6}; - -Plane Surface(1) = {1}; -Plane Surface(2) = {2}; - -Physical Line ("top") = {4}; -Physical Line ("bottom") = {1}; - -Physical Line ("contact_top") = {6}; -Physical Line ("contact_bottom") = {8}; - -Physical Surface("bot_body") = {1}; -Physical Surface("top_body") = {2}; - -Transfinite Surface{1, 2}; -Recombine Surface{1, 2}; \ No newline at end of file diff --git a/examples/contact_mechanics/implicit_detection.cc b/examples/contact_mechanics/implicit_detection.cc deleted file mode 100644 index 8e99b94f3..000000000 --- a/examples/contact_mechanics/implicit_detection.cc +++ /dev/null @@ -1,61 +0,0 @@ -/* -------------------------------------------------------------------------- */ -#include "contact_mechanics_model.hh" - -#include "dumper_elemental_field.hh" -#include "dumper_nodal_field.hh" - -#include "dumper_iohelper_paraview.hh" -/* -------------------------------------------------------------------------- */ - -using namespace akantu; - -/* -------------------------------------------------------------------------- */ -int main(int argc, char *argv[]) { - - std::string mesh_file = "implicit.msh"; - std::string material_file = "options_implicit.dat"; - - const UInt spatial_dimension = 2; - initialize(material_file, argc, argv); - - Mesh mesh(spatial_dimension); - mesh.read(mesh_file); - - ContactDetector detector(mesh); - - auto &&surface_selector = std::make_shared(mesh); - detector.setSurfaceSelector(surface_selector); - - std::map contact_map; - - detector.search(contact_map); - - UInt nb_nodes = mesh.getNbNodes(); - - Array gaps(nb_nodes, 1); - Array normals(nb_nodes, spatial_dimension); - - for (auto & entry: contact_map) { - const auto & node = entry.first; - const auto & element = entry.second; - - gaps[node] = element.gap; - for (auto i : arange(spatial_dimension)) { - normals(node, i) = element.normal[i]; - } - } - - DumperParaview dumper("implicit_detection", "./paraview", false); - dumper.registerMesh(mesh); - - auto gap_field = std::make_shared>(gaps); - auto normal_field = std::make_shared>(normals); - - dumper.registerField("gaps", gap_field); - dumper.registerField("normals", normal_field); - dumper.dump(); - - finalize(); - return EXIT_SUCCESS; -} - diff --git a/examples/contact_mechanics/material.dat b/examples/contact_mechanics/material.dat index af2dab78d..c8b95702d 100644 --- a/examples/contact_mechanics/material.dat +++ b/examples/contact_mechanics/material.dat @@ -1,30 +1,27 @@ material elastic [ - name = bot_body - rho = 7800 # density + name = upper + rho = 1.0 # density E = 1.0 # young's modulu nu = 0.0 # poisson's ratio ] material elastic [ - name = top_body - rho = 7800 # density + name = lower + rho = 1.0 # density E = 1.0 # young's modulu nu = 0.0 # poisson's ratio ] contact_detector [ - type = implicit - slave = contact_top - master = contact_bottom - two_pass_algorithm = false + type = explicit + master = contact_top + slave = contact_bottom ] contact_resolution penalty [ - name = contact_top - slave = contact_top - master = contact_bottom + name = contact_bottom mu = 0.0 - epsilon = 0.1 - epsilon_t = 100 - two_pass_algorithm = false -] \ No newline at end of file + epsilon_n = 50 + epsilon_t = 100.0 + is_master_deformable = false +] diff --git a/examples/contact_mechanics/material_explicit.dat b/examples/contact_mechanics/material_explicit.dat deleted file mode 100644 index 3512f0026..000000000 --- a/examples/contact_mechanics/material_explicit.dat +++ /dev/null @@ -1,27 +0,0 @@ -material elastic [ - name = bot_body - rho = 1.0 # density - E = 1.0 # young's modulu - nu = 0.0 # poisson's ratio -] - -material elastic [ - name = top_body - rho = 1.0 # density - E = 1e5 # young's modulu - nu = 0.0 # poisson's ratio -] - -contact_detector [ - type = implicit - slave = contact_bottom - master = contact_top - two_pass_algorithm = false -] - -contact_resolution penalty [ - name = contact_bottom - mu = 0.0 - epsilon_n = 100.0 - epsilon_t = 100 -] \ No newline at end of file diff --git a/examples/contact_mechanics/material_solid.dat b/examples/contact_mechanics/material_solid.dat deleted file mode 100644 index 5e3bf122e..000000000 --- a/examples/contact_mechanics/material_solid.dat +++ /dev/null @@ -1,13 +0,0 @@ -material elastic [ - name = bot_body - rho = 1 # density - E = 1.0 # young's modulu - nu = 0.3 # poisson's ratio -] - -material elastic [ - name = top_body - rho = 1 # density - E = 1.0 # young's modulu - nu = 0.3 # poisson's ratio -] \ No newline at end of file diff --git a/examples/contact_mechanics/options_explicit.dat b/examples/contact_mechanics/options_explicit.dat deleted file mode 100644 index 18e2068f7..000000000 --- a/examples/contact_mechanics/options_explicit.dat +++ /dev/null @@ -1,6 +0,0 @@ -contact_detector [ - type = explicit - slave = contact_bottom - master = contact_top - two_pass_algorithm = false -] diff --git a/examples/contact_mechanics/options_implicit.dat b/examples/contact_mechanics/options_implicit.dat deleted file mode 100644 index b3d60102b..000000000 --- a/examples/contact_mechanics/options_implicit.dat +++ /dev/null @@ -1,6 +0,0 @@ -contact_detector [ - type = implicit - slave = contact_bottom - master = contact_top - two_pass_algorithm = false -] diff --git a/examples/contact_mechanics/solid_explicit_dynamic.cc b/examples/contact_mechanics/solid_explicit_dynamic.cc deleted file mode 100644 index a4b639a83..000000000 --- a/examples/contact_mechanics/solid_explicit_dynamic.cc +++ /dev/null @@ -1,108 +0,0 @@ -/** - * @file test_contact_mechanics_model.cc - * - * @author Mohit Pundir - * - * @date creation: Tue Apr 30 2019 - * @date last modification: Tue Apr 30 2019 - * - * @brief Test for contact mechanics model 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 "solid_mechanics_model.hh" -#include "contact_mechanics_model.hh" -#include "coupler_solid_contact.hh" -#include "non_linear_solver.hh" -/* -------------------------------------------------------------------------- */ - -using namespace akantu; - -/* -------------------------------------------------------------------------- */ -int main(int argc, char *argv[]) { - - - const UInt spatial_dimension = 2; - initialize("material.dat", argc, argv); - - Real time_step; - Real time_factor = 0.8; - UInt max_steps = 5000; - - Mesh mesh(spatial_dimension); - mesh.read("contact_hertz_2d.msh"); - - SolidMechanicsModel solid(mesh); - - auto && selector = std::make_shared>( - "physical_names", solid); - solid.setMaterialSelector(selector); - - solid.initFull( _analysis_method = _explicit_lumped_mass); - - time_step = solid.getStableTimeStep(); - std::cout << "Time Step = " << time_step * time_factor << "s (" << time_step - << "s)" << std::endl; - solid.setTimeStep(time_step * time_factor); - - solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "bottom"); - solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "bottom"); - - solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "top"); - - solid.setBaseName("solid-explicit-dynamic"); - solid.addDumpFieldVector("displacement"); - solid.addDumpFieldVector("velocity"); - solid.addDumpFieldVector("external_force"); - solid.addDumpFieldVector("internal_force"); - solid.addDumpField("blocked_dofs"); - solid.addDumpField("grad_u"); - solid.addDumpField("stress"); - - auto & velocity = solid.getVelocity(); - - Real damping_ratio = 0.7; - - for (auto i : arange(max_steps)) { - auto increment = time_step * 5e-4; - - solid.applyBC(BC::Dirichlet::IncrementValue(-increment, _y), "top"); - - solid.solveStep(); - - Real epot = solid.getEnergy("potential"); - Real ekin = solid.getEnergy("kinetic"); - - std::cerr << i << "," << i * increment << "," << epot << "," << ekin << "," - << epot + ekin << "," << std::endl; - - for (auto & v : make_view(velocity)) { - v *= damping_ratio; - } - - solid.dump(); - } - - finalize(); - return EXIT_SUCCESS; -} - diff --git a/examples/contact_mechanics/test_node_selector.cc b/examples/contact_mechanics/test_node_selector.cc deleted file mode 100644 index 17e1805b7..000000000 --- a/examples/contact_mechanics/test_node_selector.cc +++ /dev/null @@ -1,32 +0,0 @@ -#include "aka_common.hh" -#include "surface_selector.hh" -#include "contact_mechanics_model.hh" -/* -------------------------------------------------------------------------- */ - -using namespace akantu; - -int main(int argc, char * argv[]) { - - const UInt spatial_dimension = 3; - - initialize("material.dat", argc, argv); - - Mesh mesh(spatial_dimension); - mesh.read("contact_hertz_2d.msh"); - - ContactMechanicsModel model(mesh); - - PhysicalSurfaceSelector selector(mesh); - auto & slave = selector.getSlaveList(); - auto & master = selector.getMasterList(); - - for (auto & s : slave) { - std::cerr << s << std::endl; - } - - for (auto m : master) { - std::cerr << m << std::endl; - } - - return 0; -} diff --git a/src/model/contact_mechanics/contact_mechanics_model.cc b/src/model/contact_mechanics/contact_mechanics_model.cc index fc197e906..6da43a0f1 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.cc +++ b/src/model/contact_mechanics/contact_mechanics_model.cc @@ -1,630 +1,602 @@ /** * @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.cc" #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"); 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() { 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) { // 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->stick_or_slip, 1, "stick_or_slip"); this->allocNodalField(this->normals, spatial_dimension, "normals"); this->allocNodalField(this->tangents, spatial_dimension, "tangents"); this->allocNodalField(this->projections, spatial_dimension - 1, "projections"); // todo register multipliers as dofs for lagrange multipliers } /* -------------------------------------------------------------------------- */ std::tuple ContactMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_contact: { return std::make_tuple("explicit_contact", TimeStepSolverType::_static); } case _implicit_contact: { return std::make_tuple("implicit_contact", TimeStepSolverType::_static); } case _explicit_dynamic_contact: { return std::make_tuple("explicit_dynamic_contact", TimeStepSolverType::_dynamic_lumped); break; } default: return std::make_tuple("unkown", 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); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::search() { UInt nb_nodes = mesh.getNbNodes(); contact_elements.clear(); contact_elements.resize(0); // this to resize if cohesive elements are added gaps->clear(); gaps->resize(nb_nodes, 0.); normals->clear(); normals->resize(nb_nodes, 0.); projections->clear(); projections->resize(nb_nodes, 0.); this->detector->search(contact_elements, *gaps, *normals, *projections); for (auto & gap : *gaps) { if (gap < 0) gap = std::abs(gap); else gap = -gap; } this->computeNodalAreas(); } -/* -------------------------------------------------------------------------- */ -void ContactMechanicsModel::assembleFieldsFromContactMap() { - - UInt nb_nodes = mesh.getNbNodes(); - - this->gaps->clear(); - - gaps->resize(nb_nodes, 0.); - normals->resize(nb_nodes, 0.); - tangents->resize(nb_nodes, 0.); - - if (this->contact_map.empty()) - return; - - 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]; - (*tangents)(node, i) = element.tangents(0, i); - } - } -} - /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::computeNodalAreas() { 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 = this->getFEEngineBoundary(); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); 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); auto & force = std::get<1>(tuple); for (auto & f : force) area += pow(f, 2); area = sqrt(area); } this->external_force->clear(); } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::beforeSolveStep() {} /* -------------------------------------------------------------------------- */ 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")); } for (auto & resolution : resolutions) { resolution->assembleStiffnessMatrix(_not_ghost); } } /* -------------------------------------------------------------------------- */ void ContactMechanicsModel::assembleLumpedMatrix(const ID & /*matrix_id*/) { AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ #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["stick_or_slip"] = this->stick_or_slip; 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 12eb8e8de..39e3a610a 100644 --- a/src/model/contact_mechanics/contact_mechanics_model.hh +++ b/src/model/contact_mechanics/contact_mechanics_model.hh @@ -1,363 +1,349 @@ /** * @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() 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 assembleFieldsFromContactMap(); - /* ------------------------------------------------------------------------ */ /* 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(); /* ------------------------------------------------------------------------ */ /* 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 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::areas (nodal areas) AKANTU_GET_MACRO(NodalArea, *nodal_area, Array &); /// get the ContactMechanics::areas (nodal areas) AKANTU_GET_MACRO(StickProjections, *stick_projections, Array &); /// get the ContactMechanics::areas (nodal areas) AKANTU_GET_MACRO(Projections, *projections, Array &); /// get the ContactMechanics::stick_or_slip vector (slip/stick /// state) AKANTU_GET_MACRO(StickSlip, *stick_or_slip, Array &); /// get contact detector AKANTU_GET_MACRO_NOT_CONST(ContactDetector, *detector, ContactDetector &); - /// get the contat map - inline std::map & getContactMap() { - return contact_map; - } - - /// get the contat map + /// get the contact elements inline Array & getContactElements() { return contact_elements; } - - /// - inline void setPositions(Array positions) { - detector->setPositions(positions); - } + /// 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}; /// 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 nodal areas Array * nodal_area{nullptr}; /// array to store stick/slip state : Array * stick_or_slip{nullptr}; /// array to store stick point projection in covariant basis Array * stick_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; - /// mapping between slave node its respective contact element - std::map contact_map; - /// Array contact_elements; }; } // namespace akantu /* ------------------------------------------------------------------------ */ /* inline functions */ /* ------------------------------------------------------------------------ */ #include "parser.hh" #include "resolution.hh" /* ------------------------------------------------------------------------ */ #endif /* __AKANTU_CONTACT_MECHANICS_MODEL_HH__ */ diff --git a/src/model/model_couplers/coupler_solid_cohesive_contact.cc b/src/model/model_couplers/coupler_solid_cohesive_contact.cc index ed632e883..25dfd0ec8 100644 --- a/src/model/model_couplers/coupler_solid_cohesive_contact.cc +++ b/src/model/model_couplers/coupler_solid_cohesive_contact.cc @@ -1,527 +1,550 @@ /** * @file coupler_solid_cohesive_contact.cc * * @author Mohit Pundir * * @date creation: Thu Jan 17 2019 * @date last modification: Thu May 22 2019 * * @brief class for coupling of solid mechanics and conatct 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 "coupler_solid_cohesive_contact.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { CouplerSolidCohesiveContact::CouplerSolidCohesiveContact( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, dof_manager, dim, id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject( "CohesiveFEEngine", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("coupler_solid_cohesive_contact", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_cohesive_contact", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); #endif this->registerDataAccessor(*this); solid = new SolidMechanicsModelCohesive(mesh, Model::spatial_dimension, "solid_mechanics_model_cohesive", 0, this->dof_manager); - contact = new ContactMechanicsModel(mesh, Model::spatial_dimension, + contact = new ContactMechanicsModel(mesh.getMeshFacets(), Model::spatial_dimension, "contact_mechanics_model", 0, this->dof_manager); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ CouplerSolidCohesiveContact::~CouplerSolidCohesiveContact() {} /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::initModel() { getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & CouplerSolidCohesiveContact::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::initSolver(TimeStepSolverType, NonLinearSolverType) {} /* -------------------------------------------------------------------------- */ std::tuple CouplerSolidCohesiveContact::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_contact: { return std::make_tuple("explicit_contact", TimeStepSolverType::_static); } case _implicit_contact: { return std::make_tuple("implicit_contact", TimeStepSolverType::_static); } case _explicit_dynamic_contact: { return std::make_tuple("explicit_dynamic_contact", TimeStepSolverType::_dynamic_lumped); break; } default: return std::make_tuple("unkown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ TimeStepSolverType CouplerSolidCohesiveContact::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions CouplerSolidCohesiveContact::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { 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::_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::_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 CouplerSolidCohesiveContact::assembleResidual() { // computes the internal forces this->assembleInternalForces(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); - auto & contact_map = contact->getContactMap(); + + auto get_connectivity = [&](auto & slave, auto & master) { + Vector master_conn(const_cast(mesh).getConnectivity(master)); + Vector elem_conn(master_conn.size() + 1); + + elem_conn[0] = slave; + for (UInt i = 1; i < elem_conn.size(); ++i) { + elem_conn[i] = master_conn[i - 1]; + } + return elem_conn; + }; + switch (method) { case _explicit_dynamic_contact: case _explicit_contact: { - for (auto & pair : contact_map) { - auto & connectivity = pair.second.connectivity; - for (auto node : connectivity) { - for (auto s : arange(spatial_dimension)) - external_force(node, s) = contact_force(node, s); + for (auto & element : contact->getContactElements()) { + for (auto & conn : get_connectivity(element.slave, element.master)) { + for (auto dim : arange(spatial_dimension)) { + external_force(conn, dim) = contact_force(conn, dim); + } } } - break; } default: break; } /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", external_force, 1); this->getDOFManager().assembleToResidual("displacement", internal_force, 1); switch (method) { case _implicit_contact: { this->getDOFManager().assembleToResidual("displacement", contact_force, 1); break; } default: break; } } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); auto & internal_force = solid->getInternalForce(); auto & external_force = solid->getExternalForce(); auto & contact_force = contact->getInternalForce(); - auto & contact_map = contact->getContactMap(); + auto get_connectivity = [&](auto & slave, auto & master) { + Vector master_conn(const_cast(mesh).getConnectivity(master)); + Vector elem_conn(master_conn.size() + 1); + + elem_conn[0] = slave; + for (UInt i = 1; i < elem_conn.size(); ++i) { + elem_conn[i] = master_conn[i - 1]; + } + return elem_conn; + }; + + switch (method) { case _explicit_dynamic_contact: case _explicit_contact: { - for (auto & pair : contact_map) { - auto & connectivity = pair.second.connectivity; - for (auto node : connectivity) { - for (auto s : arange(spatial_dimension)) - external_force(node, s) = contact_force(node, s); + for (auto & element : contact->getContactElements()) { + for (auto & conn : get_connectivity(element.slave, element.master)) { + for (auto dim : arange(spatial_dimension)) { + external_force(conn, dim) = contact_force(conn, dim); + } } } - break; } default: break; } + if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->getDOFManager().assembleToResidual("displacement", internal_force, 1); switch (method) { case _implicit_contact: { this->getDOFManager().assembleToResidual("displacement", contact_force, 1); break; } default: break; } AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::beforeSolveStep() {} /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::afterSolveStep() {} /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::predictor() { switch (method) { case _explicit_dynamic_contact: { Array displacement(0, Model::spatial_dimension); auto & current_positions = contact->getContactDetector().getPositions(); current_positions.copy(mesh.getNodes()); auto us = this->getDOFManager().getDOFs("displacement"); const auto blocked_dofs = this->getDOFManager().getBlockedDOFs("displacement"); for (auto && tuple : zip(make_view(us), make_view(blocked_dofs), make_view(current_positions))) { auto & u = std::get<0>(tuple); const auto & bld = std::get<1>(tuple); auto & cp = std::get<2>(tuple); if (not bld) cp += u; } contact->search(); break; } default: break; } } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::corrector() { switch (method) { case _implicit_contact: case _explicit_contact: { Array displacement(0, Model::spatial_dimension); auto & current_positions = contact->getContactDetector().getPositions(); current_positions.copy(mesh.getNodes()); auto us = this->getDOFManager().getDOFs("displacement"); const auto blocked_dofs = this->getDOFManager().getBlockedDOFs("displacement"); for (auto && tuple : zip(make_view(us), make_view(blocked_dofs), make_view(current_positions))) { auto & u = std::get<0>(tuple); const auto & bld = std::get<1>(tuple); auto & cp = std::get<2>(tuple); if (not bld) cp += u; } contact->search(); break; } default: break; } } /* -------------------------------------------------------------------------- */ MatrixType CouplerSolidCohesiveContact::getMatrixType(const ID & matrix_id) { if (matrix_id == "K") return _symmetric; if (matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { solid->assembleMass(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { solid->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); solid->assembleInternalForces(); contact->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); solid->assembleStiffnessMatrix(); switch (method) { case _implicit_contact: { contact->assembleStiffnessMatrix(); break; } default: break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMassLumped() { solid->assembleMassLumped(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMass() { solid->assembleMass(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMassLumped(GhostType ghost_type) { solid->assembleMassLumped(ghost_type); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::assembleMass(GhostType ghost_type) { solid->assembleMass(ghost_type); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) { std::shared_ptr field; field = solid->createElementalField(field_name, group_name, padding_flag, spatial_dimension, kind); return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldReal( const std::string & field_name, const std::string & group_name, bool padding_flag) { std::shared_ptr field; if (field_name == "contact_force" or field_name == "normals" or + field_name == "normal_force" or field_name == "tangential_force" or + field_name == "stick_or_slip" or field_name == "gaps" or field_name == "previous_gaps" or field_name == "areas" or field_name == "tangents") { - field = contact->createNodalFieldReal(field_name, group_name, padding_flag); + field = contact->createNodalFieldReal(field_name, group_name, padding_flag); } else { field = solid->createNodalFieldReal(field_name, group_name, padding_flag); } - + return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldBool( const std::string & field_name, const std::string & group_name, bool padding_flag) { std::shared_ptr field; field = solid->createNodalFieldBool(field_name, group_name, padding_flag); return field; } #else /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createElementalField(const std::string &, const std::string &, bool, const UInt &, const ElementKind &) { return nullptr; } /* ----------------------------------------------------------------------- */ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldReal(const std::string &, const std::string &, bool) { return nullptr; } /*-------------------------------------------------------------------*/ std::shared_ptr CouplerSolidCohesiveContact::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(const std::string & dumper_name) { solid->onDump(); mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(const std::string & dumper_name, UInt step) { solid->onDump(); mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(const std::string & dumper_name, Real time, UInt step) { solid->onDump(); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump() { solid->onDump(); mesh.dump(); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(UInt step) { solid->onDump(); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void CouplerSolidCohesiveContact::dump(Real time, UInt step) { solid->onDump(); mesh.dump(time, step); } } // namespace akantu