diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_2d/test_contact_2d_expli/test_contact_2d_expli.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_2d/test_contact_2d_expli/test_contact_2d_expli.cc index 7ca3b9427..8af374dc9 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_2d/test_contact_2d_expli/test_contact_2d_expli.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_2d/test_contact_2d_expli/test_contact_2d_expli.cc @@ -1,265 +1,265 @@ /** * @file test_contact_2d_expli.cc * @author Leonardo Snozzi * @date Fri Nov 26 07:43:47 2010 * * @brief test explicit DCR contact algorithm for 2d * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "io_helper.h" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "io_helper.h" #endif //AKANTU_USE_IOHELPER #define NORMAL_PRESSURE -1.e6 using namespace akantu; static void reduceGap(const SolidMechanicsModel & model, const Real threshold, const Real gap); static void setBoundaryConditions(SolidMechanicsModel & model); void my_force(double * coord, double *T); static void reduceVelocities(const SolidMechanicsModel & model, const Real ratio); static void initParaview(SolidMechanicsModel & model); Real y_min, y_max; DumperParaview dumper; int main(int argc, char *argv[]) { UInt spatial_dimension = 2; UInt max_steps = 30000; Real time_factor = 0.2; Mesh mesh(spatial_dimension); MeshIOMSH mesh_io; mesh_io.read("squares.msh", mesh); SolidMechanicsModel * model = new SolidMechanicsModel(mesh); /// get two squares closer // reduceGap(*model, 0.05, 1.e-6); UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); UInt nb_elements = model->getFEM().getMesh().getNbElement(_triangle_3); /// model initialization model->initVectors(); model->readMaterials("materials.dat"); model->initMaterials(); model->initModel(); std::cout << model->getMaterial(0) << std::endl; model->assembleMassLumped(); /// set vectors to zero memset(model->getForce().values, 0, spatial_dimension*nb_nodes*sizeof(Real)); memset(model->getVelocity().values, 0, spatial_dimension*nb_nodes*sizeof(Real)); memset(model->getAcceleration().values, 0, spatial_dimension*nb_nodes*sizeof(Real)); memset(model->getDisplacement().values, 0, spatial_dimension*nb_nodes*sizeof(Real)); memset(model->getResidual().values, 0, spatial_dimension*nb_nodes*sizeof(Real)); memset(model->getMaterial(0).getStrain(_triangle_3).values, 0, spatial_dimension*spatial_dimension*nb_elements*sizeof(Real)); memset(model->getMaterial(0).getStress(_triangle_3).values, 0, spatial_dimension*spatial_dimension*nb_elements*sizeof(Real)); /// Paraview Helper #ifdef AKANTU_USE_IOHELPER initParaview(*model); #endif //AKANTU_USE_IOHELPER Real time_step = model->getStableTimeStep() * time_factor; std::cout << "Time Step = " << time_step << "s" << std::endl; model->setTimeStep(time_step); /// set boundary conditions setBoundaryConditions(*model); /// define and initialize contact Contact * my_contact = Contact::newContact(*model, _ct_2d_expli, _cst_2d_expli, _cnst_2d_grid); my_contact->initContact(true); my_contact->setFrictionCoefficient(0.); my_contact->initNeighborStructure(); my_contact->initSearch(); for (UInt s = 0; s < max_steps; ++s) { model->explicitPred(); model->updateCurrentPosition(); my_contact->solveContact(); model->updateResidual(); model->updateAcceleration(); model->explicitCorr(); if(s % 200 == 0) dumper.Dump(); if(s%100 == 0 && s>499) reduceVelocities(*model, 0.95); if(s % 500 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; } delete my_contact; delete model; finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ static void reduceGap(const SolidMechanicsModel & model, const Real threshold, const Real gap) { UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); Real * coord = model.getFEM().getMesh().getNodes().values; Real y_top = HUGE_VAL, y_bot = -HUGE_VAL; for (UInt n = 0; n < nb_nodes; ++n) { if (coord[2*n+1] > threshold) { if(coord[2*n+1] < y_top) y_top = coord[2*n+1]; } else { if (coord[2*n+1] > y_bot) y_bot = coord[2*n+1]; } } Real delta = y_top - y_bot - gap; /// move all nodes belonging to the top cube for (UInt n = 0; n < nb_nodes; ++n) { if (coord[2*n+1] > threshold) coord[2*n+1] -= delta; } } /* -------------------------------------------------------------------------- */ static void setBoundaryConditions(SolidMechanicsModel & model) { UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); Real * coord = model.getFEM().getMesh().getNodes().values; for (UInt n = 0; n < nb_nodes; ++n) { if (coord[2*n+1] > y_max) y_max = coord[2*n+1]; if (coord[2*n+1] < y_min) y_min = coord[2*n+1]; } FEM & b_fem = model.getFEMBoundary(); b_fem.initShapeFunctions(); - b_fem.computeNormalsOnQuadPoints(); + b_fem.computeNormalsOnControlPoints(); bool * id = model.getBoundary().values; memset(id, 0, 2*nb_nodes*sizeof(bool)); std::cout << "Nodes "; for (UInt i = 0; i < nb_nodes; ++i) { if (coord[2*i+1] < y_min + 1.e-5) { id[2*i+1] = true; std::cout << " " << i << " "; } } std::cout << "are blocked" << std::endl; model.computeForcesFromFunction(my_force, _bft_stress); } /* -------------------------------------------------------------------------- */ void my_force(double * coord, double *T) { memset(T, 0, 4*sizeof(double)); if(*(coord+1) > y_max-1.e-5) T[3] = NORMAL_PRESSURE; } /* -------------------------------------------------------------------------- */ /// artificial damping of velocities in order to reach a global static equilibrium static void reduceVelocities(const SolidMechanicsModel & model, const Real ratio) { UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); Real * velocities = model.getVelocity().values; if(ratio>1.) { fprintf(stderr,"**error** in Reduce_Velocities ratio bigger than 1!\n"); exit(-1); } for(UInt i =0; i # @date Wed Jan 19 12:37:24 2011 # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== #=============================================================================== add_mesh(test_squares_mesh squares.geo 2 1) register_test(test_contact_rigid_no_friction_triangle_3 test_contact_rigid_no_friction_triangle_3.cc) add_dependencies(test_contact_rigid_no_friction_triangle_3 test_squares_mesh) #=============================================================================== add_mesh(test_cubes_mesh cubes.geo 3 1) register_test(test_contact_rigid_no_friction_tetrahedron_4 test_contact_rigid_no_friction_tetrahedron_4.cc) add_dependencies(test_contact_rigid_no_friction_tetrahedron_4 test_cubes_mesh) #=============================================================================== add_mesh(test_hertz_2d_mesh hertz_2d.geo 2 1) register_test(test_contact_rigid_no_friction_hertz_2d test_contact_rigid_no_friction_hertz_2d.cc) add_dependencies(test_contact_rigid_no_friction_hertz_2d test_hertz_2d_mesh) #=============================================================================== add_mesh(test_hertz_3d_mesh hertz_3d.geo 3 1) register_test(test_contact_rigid_no_friction_hertz_3d test_contact_rigid_no_friction_hertz_3d.cc) add_dependencies(test_contact_rigid_no_friction_hertz_3d test_hertz_3d_mesh) #=============================================================================== add_mesh(test_hertz_3d_full_mesh hertz_3d_full.geo 3 1) register_test(test_contact_rigid_no_friction_hertz_3d_full test_contact_rigid_no_friction_hertz_3d_full.cc) add_dependencies(test_contact_rigid_no_friction_hertz_3d_full test_hertz_3d_full_mesh) #=============================================================================== add_mesh(test_force_2d_mesh force_2d.geo 2 1) register_test(test_contact_rigid_no_friction_force_2d test_contact_rigid_no_friction_force_2d.cc) add_dependencies(test_contact_rigid_no_friction_force_2d test_force_2d_mesh) #=============================================================================== add_mesh(test_force_3d_mesh force_3d.geo 3 1) register_test(test_contact_rigid_no_friction_force_3d test_contact_rigid_no_friction_force_3d.cc) add_dependencies(test_contact_rigid_no_friction_force_3d test_force_3d_mesh) #=============================================================================== add_mesh(test_sliding_cube_2d_mesh sliding_cube_2d.geo 2 1) register_test(test_contact_rigid_sliding_cube_2d test_contact_rigid_sliding_cube_2d.cc) add_dependencies(test_contact_rigid_sliding_cube_2d test_sliding_cube_2d_mesh) +#=============================================================================== +add_mesh(test_sliding_cube_3d_mesh sliding_cube_3d.geo 3 1) + +register_test(test_contact_rigid_sliding_cube_3d test_contact_rigid_sliding_cube_3d.cc) +add_dependencies(test_contact_rigid_sliding_cube_3d test_sliding_cube_3d_mesh) + #=============================================================================== configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/material.dat ${CMAKE_CURRENT_BINARY_DIR}/material.dat COPYONLY ) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview/hertz_2d) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview/hertz_3d) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview/hertz_3d_full) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview/force_2d) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview/force_3d) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview/sliding_cube_2d) +file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview/sliding_cube_3d) diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/sliding_cube_3d.geo b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/sliding_cube_3d.geo new file mode 100644 index 000000000..f56fe047f --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/sliding_cube_3d.geo @@ -0,0 +1,81 @@ +h_contact = 0.05; +h_top = 0.4; +h_plate = 0.2; + +Point(1) = {0, 0, 0, h_top}; +Point(2) = {0, -1, 0, h_contact}; +Point(3) = {1, -1, 0, h_contact}; +Point(4) = {1, 0, 0, h_top}; + +Point(5) = {0, 0, 1, h_top}; +Point(6) = {0, -1, 1, h_contact}; +Point(7) = {1, -1, 1, h_contact}; +Point(8) = {1, 0, 1, h_top}; + +Point(9) = {-0.25, -1, -0.25, h_plate}; +Point(10) = {-0.25, -1.5, -0.25, h_plate}; +Point(11) = { 2.25, -1.5, -0.25, h_plate}; +Point(12) = { 2.25, -1, -0.25, h_plate}; + +Point(13) = {-0.25, -1, 1.25, h_plate}; +Point(14) = {-0.25, -1.5, 1.25, h_plate}; +Point(15) = { 2.25, -1.5, 1.25, h_plate}; +Point(16) = { 2.25, -1, 1.25, h_plate}; + +Line(1) = {1, 5}; +Line(2) = {5, 8}; +Line(3) = {8, 4}; +Line(4) = {4, 1}; +Line(5) = {6, 2}; +Line(6) = {2, 3}; +Line(7) = {3, 7}; +Line(8) = {7, 6}; +Line(9) = {5, 6}; +Line(10) = {1, 2}; +Line(11) = {4, 3}; +Line(12) = {8, 7}; + +Line(13) = {9, 13}; +Line(14) = {13, 16}; +Line(15) = {16, 12}; +Line(16) = {12, 9}; +Line(17) = {14, 10}; +Line(18) = {10, 11}; +Line(19) = {11, 15}; +Line(20) = {15, 14}; +Line(21) = {13, 14}; +Line(22) = {9, 10}; +Line(23) = {12, 11}; +Line(24) = {16, 15}; + +Line Loop(26) = {4, 1, 2, 3}; +Plane Surface(26) = {26}; +Line Loop(28) = {9, -8, -12, -2}; +Plane Surface(28) = {28}; +Line Loop(30) = {-7, -11, -3, 12}; +Plane Surface(30) = {30}; +Line Loop(32) = {-10, -4, 11, -6}; +Plane Surface(32) = {32}; +Line Loop(34) = {10, -5, -9, -1}; +Plane Surface(34) = {34}; +Line Loop(36) = {8, 5, 6, 7}; +Plane Surface(36) = {36}; + +Line Loop(38) = {16, 13, 14, 15}; +Plane Surface(38) = {38}; +Line Loop(40) = {21, -20, -24, -14}; +Plane Surface(40) = {40}; +Line Loop(42) = {-19, -23, -15, 24}; +Plane Surface(42) = {42}; +Line Loop(44) = {-22, -16, 23, -18}; +Plane Surface(44) = {44}; +Line Loop(46) = {22, -17, -21, -13}; +Plane Surface(46) = {46}; +Line Loop(48) = {20, 17, 18, 19}; +Plane Surface(48) = {48}; + +Surface Loop(50) = {26, 28, 30, 32, 34, 36}; +Volume(50) = {50}; + +Surface Loop(52) = {38, 40, 42, 44, 46, 48}; +Volume(52) = {52}; diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_force_2d.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_force_2d.cc index 3d23bb3d2..44ad29641 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_force_2d.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_force_2d.cc @@ -1,239 +1,230 @@ /** * @file test_contact_rigid_no_friction_force_2d.cc * @author David Kammer * @date Mon Jan 24 10:04:42 2011 * * @brief test for force in 2d rigid contact in explicit * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" +#include "contact_rigid.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" #include "contact_search.hh" #include "contact_search_explicit.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 2; const ElementType element_type = _triangle_3; const UInt paraview_type = TRIANGLE1; UInt max_steps = 200000; UInt imposing_steps = 100000; Real max_displacement = -0.1; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; mesh_io.read("force_2d.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); UInt nb_nodes = my_mesh.getNbNodes(); /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors memset(my_model.getForce().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getBoundary().values, false, dim*nb_nodes*sizeof(bool)); + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); UInt nb_element = my_model.getFEM().getMesh().getNbElement(element_type); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); my_model.assembleMassLumped(); /// contact declaration - Contact * my_contact = Contact::newContact(my_model, - _ct_rigid, - _cst_expli, - _cnst_regular_grid); + Contact * contact = Contact::newContact(my_model, + _ct_rigid, + _cst_expli, + _cnst_regular_grid); + + ContactRigid * my_contact = dynamic_cast(contact); my_contact->initContact(false); Surface master = 1; my_contact->addMasterSurface(master); - /*const RegularGridNeighborStructure<2> & my_rgns = dynamic_cast &>(my_contact->getContactSearch().getContactNeighborStructure(master)); - const_cast&>(my_rgns).setGridSpacing(0.075, 0); - const_cast&>(my_rgns).setGridSpacing(0.075, 1);*/ - my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); my_contact->initSearch(); // does nothing so far // boundary conditions Surface impactor = 0; Vector * top_nodes = new Vector(0, 1); - UInt middle_point; Real * coordinates = my_mesh.getNodes().values; Real * displacement = my_model.getDisplacement().values; bool * boundary = my_model.getBoundary().values; UInt * surface_to_nodes_offset = my_contact->getSurfaceToNodesOffset().values; UInt * surface_to_nodes = my_contact->getSurfaceToNodes().values; // symetry boundary conditions for(UInt n = surface_to_nodes_offset[impactor]; n < surface_to_nodes_offset[impactor+1]; ++n) { UInt node = surface_to_nodes[n]; - Real x_coord = coordinates[node*dim]; Real y_coord = coordinates[node*dim + 1]; - if (x_coord < 0.00001) - boundary[node*dim] = true; if (y_coord > -0.00001) { boundary[node*dim + 1] = true; top_nodes->push_back(node); } - if (x_coord < 0.00001 && y_coord > -0.00001) - middle_point = node; } // ground boundary conditions for(UInt n = surface_to_nodes_offset[master]; n < surface_to_nodes_offset[master+1]; ++n) { UInt node = surface_to_nodes[n]; Real y_coord = coordinates[node*dim + 1]; if (y_coord < -1.2) boundary[node*dim] = true; boundary[node*dim + 1] = true; } UInt * top_nodes_val = top_nodes->values; #ifdef AKANTU_USE_IOHELPER /// initialize the paraview output DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(my_model.getFEM().getMesh().getNodes().values, dim, nb_nodes, "coordinates_force_2d"); dumper.SetConnectivity((int *)my_model.getFEM().getMesh().getConnectivity(element_type).values, paraview_type, nb_element, C_MODE); dumper.AddNodeDataField(my_model.getDisplacement().values, dim, "displacements"); dumper.AddNodeDataField(my_model.getVelocity().values, dim, "velocity"); dumper.AddNodeDataField(my_model.getResidual().values, dim, "force"); dumper.AddNodeDataField(my_model.getForce().values, dim, "applied_force"); dumper.AddElemDataField(my_model.getMaterial(0).getStrain(element_type).values, dim*dim, "strain"); dumper.AddElemDataField(my_model.getMaterial(0).getStress(element_type).values, dim*dim, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/force_2d/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER std::ofstream force_out; force_out.open("force_2d.csv"); force_out << "%id,ftop,fcont,zone" << std::endl; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for(UInt s = 1; s <= max_steps; ++s) { if(s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; if(s == imposing_steps){ my_model.updateCurrentPosition(); my_contact->updateContact(); } if(s <= imposing_steps) { Real current_displacement = max_displacement/(static_cast(imposing_steps))*s; for(UInt n=0; ngetSize(); ++n) { UInt node = top_nodes_val[n]; displacement[node*dim + 1] = current_displacement; } } my_model.explicitPred(); - my_model.updateCurrentPosition(); - - /// compute the penetration list - PenetrationList * my_penetration_list = new PenetrationList(); - const_cast(my_contact->getContactSearch()).findPenetration(master, *my_penetration_list); - UInt nb_nodes_pen = my_penetration_list->penetrating_nodes.getSize(); - Vector pen_nodes = my_penetration_list->penetrating_nodes; - UInt * pen_nodes_val = pen_nodes.values; - + my_model.initializeUpdateResidualData(); + my_contact->solveContact(); my_model.updateResidual(false); + my_contact->avoidAdhesion(); + + // find the total force applied at the imposed displacement surface (top) Real * residual = my_model.getResidual().values; Real top_force = 0.; for(UInt n=0; ngetSize(); ++n) { UInt node = top_nodes_val[n]; top_force += residual[node*dim + 1]; } - my_model.updateCurrentPosition(); + + // find the total contact force and contact area + ContactRigid::ImpactorNodesInfoPerMaster * imp_info = my_contact->getImpactorsInformation().at(master); + UInt * active_imp_nodes_val = imp_info->active_impactor_nodes->values; Real * current_position = my_model.getCurrentPosition().values; Real contact_force = 0.; Real contact_zone = 0.; - for (UInt i = 0; i < nb_nodes_pen; ++i) { - UInt node = pen_nodes_val[i]; + for (UInt i = 0; i < imp_info->active_impactor_nodes->getSize(); ++i) { + UInt node = active_imp_nodes_val[i]; contact_force += residual[node*dim + 1]; contact_zone = std::max(contact_zone, current_position[node*dim]); } - delete my_penetration_list; force_out << s << "," << top_force << "," << contact_force << "," << contact_zone << std::endl; my_model.updateAcceleration(); my_model.explicitCorr(); #ifdef AKANTU_USE_IOHELPER if(s % 1000 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER } force_out.close(); delete my_contact; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_force_3d.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_force_3d.cc index 098fb1b4b..1982b7c85 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_force_3d.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_force_3d.cc @@ -1,241 +1,241 @@ /** * @file test_contact_rigid_no_friction_force_3d.cc * @author David Kammer * @date Mon Jan 24 11:56:42 2011 * * @brief test force for rigid contact 3d in explicit * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" #include "contact_search.hh" #include "contact_search_explicit.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 3; const ElementType element_type = _tetrahedron_4; const UInt paraview_type = TETRA1; UInt max_steps = 200000; UInt imposing_steps = 100000; Real max_displacement = -0.1; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; mesh_io.read("force_3d.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); UInt nb_nodes = my_mesh.getNbNodes(); /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors memset(my_model.getForce().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getBoundary().values, false, dim*nb_nodes*sizeof(bool)); + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); UInt nb_element = my_model.getFEM().getMesh().getNbElement(element_type); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); my_model.assembleMassLumped(); /// contact declaration Contact * my_contact = Contact::newContact(my_model, _ct_rigid, _cst_expli, _cnst_regular_grid); my_contact->initContact(false); Surface master = 1; my_contact->addMasterSurface(master); /*const RegularGridNeighborStructure<3> & my_rgns = dynamic_cast &>(my_contact->getContactSearch().getContactNeighborStructure(master)); const_cast&>(my_rgns).setGridSpacing(0.075, 0); const_cast&>(my_rgns).setGridSpacing(0.075, 1); const_cast&>(my_rgns).setGridSpacing(0.075, 2);*/ my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); my_contact->initSearch(); // does nothing so far // boundary conditions Surface impactor = 0; Vector * top_nodes = new Vector(0, 1); Real * coordinates = my_mesh.getNodes().values; Real * displacement = my_model.getDisplacement().values; bool * boundary = my_model.getBoundary().values; UInt * surface_to_nodes_offset = my_contact->getSurfaceToNodesOffset().values; UInt * surface_to_nodes = my_contact->getSurfaceToNodes().values; // symetry boundary conditions for(UInt n = surface_to_nodes_offset[impactor]; n < surface_to_nodes_offset[impactor+1]; ++n) { UInt node = surface_to_nodes[n]; Real x_coord = coordinates[node*dim]; Real y_coord = coordinates[node*dim + 1]; Real z_coord = coordinates[node*dim + 2]; /*if (x_coord < 0.00001) boundary[node*dim] = true; if (z_coord < 0.00001) boundary[node*dim+2] = true;*/ if (y_coord > -0.00001) { boundary[node*dim + 1] = true; top_nodes->push_back(node); } } // ground boundary conditions for(UInt n = surface_to_nodes_offset[master]; n < surface_to_nodes_offset[master+1]; ++n) { UInt node = surface_to_nodes[n]; Real y_coord = coordinates[node*dim + 1]; if (y_coord < -1.2) boundary[node*dim] = true; boundary[node*dim + 1] = true; boundary[node*dim + 2] = true; } UInt * top_nodes_val = top_nodes->values; #ifdef AKANTU_USE_IOHELPER /// initialize the paraview output DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(my_model.getFEM().getMesh().getNodes().values, dim, nb_nodes, "coordinates_force_3d"); dumper.SetConnectivity((int *)my_model.getFEM().getMesh().getConnectivity(element_type).values, paraview_type, nb_element, C_MODE); dumper.AddNodeDataField(my_model.getDisplacement().values, dim, "displacements"); dumper.AddNodeDataField(my_model.getVelocity().values, dim, "velocity"); dumper.AddNodeDataField(my_model.getResidual().values, dim, "force"); dumper.AddNodeDataField(my_model.getForce().values, dim, "applied_force"); dumper.AddElemDataField(my_model.getMaterial(0).getStrain(element_type).values, dim*dim, "strain"); dumper.AddElemDataField(my_model.getMaterial(0).getStress(element_type).values, dim*dim, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/force_3d/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER std::ofstream force_out; force_out.open("force_3d.csv"); force_out << "%id,ftop,fcont,zone" << std::endl; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for(UInt s = 1; s <= max_steps; ++s) { if(s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; if(s == imposing_steps){ my_model.updateCurrentPosition(); my_contact->updateContact(); } if(s <= imposing_steps) { Real current_displacement = max_displacement/(static_cast(imposing_steps))*s; for(UInt n=0; ngetSize(); ++n) { UInt node = top_nodes_val[n]; displacement[node*dim + 1] = current_displacement; } } my_model.explicitPred(); - my_model.updateCurrentPosition(); + my_model.initializeUpdateResidualData(); /// compute the penetration list PenetrationList * my_penetration_list = new PenetrationList(); const_cast(my_contact->getContactSearch()).findPenetration(master, *my_penetration_list); UInt nb_nodes_pen = my_penetration_list->penetrating_nodes.getSize(); Vector pen_nodes = my_penetration_list->penetrating_nodes; UInt * pen_nodes_val = pen_nodes.values; my_contact->solveContact(); my_model.updateResidual(false); Real * residual = my_model.getResidual().values; Real top_force = 0.; for(UInt n=0; ngetSize(); ++n) { UInt node = top_nodes_val[n]; top_force += residual[node*dim + 1]; } my_model.updateCurrentPosition(); Real * current_position = my_model.getCurrentPosition().values; Real contact_force = 0.; Real contact_zone = 0.; for (UInt i = 0; i < nb_nodes_pen; ++i) { UInt node = pen_nodes_val[i]; contact_force += residual[node*dim + 1]; contact_zone = std::max(contact_zone, current_position[node*dim]); } delete my_penetration_list; force_out << s << "," << top_force << "," << contact_force << "," << contact_zone << std::endl; my_model.updateAcceleration(); my_model.explicitCorr(); #ifdef AKANTU_USE_IOHELPER if(s % 1000 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER } force_out.close(); delete my_contact; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_hertz_2d.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_hertz_2d.cc index d2c1fa7ec..986221368 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_hertz_2d.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_hertz_2d.cc @@ -1,258 +1,267 @@ /** * @file test_contact_rigid_no_friction_hertz_2d.cc * @author David Kammer * @date Wed Jan 19 15:04:42 2011 * * @brief test contact search for 2d hertz in explicit * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_rigid.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" #include "contact_search.hh" #include "contact_search_explicit.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 2; const ElementType element_type = _triangle_3; const UInt paraview_type = TRIANGLE1; UInt max_steps = 200000; UInt imposing_steps = 100000; Real max_displacement = -0.1; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; mesh_io.read("hertz_2d.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); UInt nb_nodes = my_mesh.getNbNodes(); /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors memset(my_model.getForce().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getBoundary().values, false, dim*nb_nodes*sizeof(bool)); + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); UInt nb_element = my_model.getFEM().getMesh().getNbElement(element_type); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); my_model.assembleMassLumped(); /// contact declaration Contact * contact = Contact::newContact(my_model, _ct_rigid, _cst_expli, _cnst_regular_grid); ContactRigid * my_contact = dynamic_cast(contact); my_contact->initContact(false); Surface master = 1; my_contact->addMasterSurface(master); /*const RegularGridNeighborStructure<2> & my_rgns = dynamic_cast &>(my_contact->getContactSearch().getContactNeighborStructure(master)); const_cast&>(my_rgns).setGridSpacing(0.075, 0); const_cast&>(my_rgns).setGridSpacing(0.075, 1); */ my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); my_contact->initSearch(); // does nothing so far // boundary conditions Surface impactor = 0; Vector * top_nodes = new Vector(0, 1); UInt middle_point; Real * coordinates = my_mesh.getNodes().values; Real * displacement = my_model.getDisplacement().values; bool * boundary = my_model.getBoundary().values; UInt * surface_to_nodes_offset = my_contact->getSurfaceToNodesOffset().values; UInt * surface_to_nodes = my_contact->getSurfaceToNodes().values; // symetry boundary conditions for(UInt n = surface_to_nodes_offset[impactor]; n < surface_to_nodes_offset[impactor+1]; ++n) { UInt node = surface_to_nodes[n]; Real x_coord = coordinates[node*dim]; Real y_coord = coordinates[node*dim + 1]; if (x_coord < 0.00001) boundary[node*dim] = true; if (y_coord > -0.00001) { boundary[node*dim + 1] = true; top_nodes->push_back(node); } if (x_coord < 0.00001 && y_coord > -0.00001) middle_point = node; } // ground boundary conditions for(UInt n = surface_to_nodes_offset[master]; n < surface_to_nodes_offset[master+1]; ++n) { UInt node = surface_to_nodes[n]; Real y_coord = coordinates[node*dim + 1]; if (y_coord < -1.2) boundary[node*dim] = true; boundary[node*dim + 1] = true; } UInt * top_nodes_val = top_nodes->values; #ifdef AKANTU_USE_IOHELPER /// initialize the paraview output DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(my_model.getFEM().getMesh().getNodes().values, dim, nb_nodes, "coordinates_2d"); dumper.SetConnectivity((int *)my_model.getFEM().getMesh().getConnectivity(element_type).values, paraview_type, nb_element, C_MODE); dumper.AddNodeDataField(my_model.getDisplacement().values, dim, "displacements"); dumper.AddNodeDataField(my_model.getVelocity().values, dim, "velocity"); dumper.AddNodeDataField(my_model.getResidual().values, dim, "force"); dumper.AddNodeDataField(my_model.getForce().values, dim, "applied_force"); dumper.AddElemDataField(my_model.getMaterial(0).getStrain(element_type).values, dim*dim, "strain"); dumper.AddElemDataField(my_model.getMaterial(0).getStress(element_type).values, dim*dim, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/hertz_2d/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER std::ofstream hertz; hertz.open("hertz_2d.csv"); hertz << "%id,ftop,fcont,zone" << std::endl; std::ofstream energy; energy.open("hertz_2d_energy.csv"); - energy << "%id,kin,pot" << std::endl; + energy << "%id,kin,pot,tot" << std::endl; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for(UInt s = 1; s <= max_steps; ++s) { if(s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; if(s % 200 == 0){ my_model.updateCurrentPosition(); my_contact->updateContact(); } if(s <= imposing_steps) { Real current_displacement = max_displacement/(static_cast(imposing_steps))*s; for(UInt n=0; ngetSize(); ++n) { UInt node = top_nodes_val[n]; displacement[node*dim + 1] = current_displacement; } } my_model.explicitPred(); my_model.initializeUpdateResidualData(); /// compute the penetration list /*PenetrationList * my_penetration_list = new PenetrationList(); const_cast(my_contact->getContactSearch()).findPenetration(master, *my_penetration_list); UInt nb_nodes_pen = my_penetration_list->penetrating_nodes.getSize(); Vector pen_nodes = my_penetration_list->penetrating_nodes; UInt * pen_nodes_val = pen_nodes.values;*/ my_contact->solveContact(); my_model.updateResidual(false); my_contact->avoidAdhesion(); // find the total force applied at the imposed displacement surface (top) Real * residual = my_model.getResidual().values; Real top_force = 0.; for(UInt n=0; ngetSize(); ++n) { UInt node = top_nodes_val[n]; top_force += residual[node*dim + 1]; } + // find index of master surface in impactors_information + Int master_index = -1; + for (UInt i=0; i < my_contact->getImpactorsInformation().size(); ++i) { + if (my_contact->getImpactorsInformation().at(i)->master_id == master) { + master_index = i; + break; + } + } + // find the total contact force and contact area - const Vector * active_imp_nodes = my_contact->getActiveImpactorNodes(); - UInt * active_imp_nodes_val = active_imp_nodes->values; + ContactRigid::ImpactorNodesInfoPerMaster * imp_info = my_contact->getImpactorsInformation().at(master_index); + UInt * active_imp_nodes_val = imp_info->active_impactor_nodes->values; Real * current_position = my_model.getCurrentPosition().values; Real contact_force = 0.; Real contact_zone = 0.; - for (UInt i = 0; i < active_imp_nodes->getSize(); ++i) { + for (UInt i = 0; i < imp_info->active_impactor_nodes->getSize(); ++i) { UInt node = active_imp_nodes_val[i]; contact_force += residual[node*dim + 1]; contact_zone = std::max(contact_zone, current_position[node*dim]); } //delete my_penetration_list; hertz << s << "," << top_force << "," << contact_force << "," << contact_zone << std::endl; my_model.updateAcceleration(); my_model.explicitCorr(); Real epot = my_model.getPotentialEnergy(); Real ekin = my_model.getKineticEnergy(); - energy << s << "," << ekin << "," << epot << std::endl; + energy << s << "," << ekin << "," << epot << "," << ekin+epot << std::endl; #ifdef AKANTU_USE_IOHELPER if(s % 100 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER } hertz.close(); energy.close(); delete my_contact; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_hertz_3d.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_hertz_3d.cc index 8b9bda6ae..85d5f1af2 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_hertz_3d.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_hertz_3d.cc @@ -1,244 +1,244 @@ /** * @file test_contact_rigid_no_friction_hertz_3d.cc * @author David Kammer * @date Thu Jan 20 15:50:42 2011 * * @brief test rigid contact for 3d hertz in explicit * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" #include "contact_search.hh" #include "contact_search_explicit.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 3; const ElementType element_type = _tetrahedron_4; const UInt paraview_type = TETRA1; UInt max_steps = 200000; UInt imposing_steps = 100000; Real max_displacement = -0.1; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; mesh_io.read("hertz_3d.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); UInt nb_nodes = my_mesh.getNbNodes(); /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors memset(my_model.getForce().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getBoundary().values, false, dim*nb_nodes*sizeof(bool)); + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); UInt nb_element = my_model.getFEM().getMesh().getNbElement(element_type); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); my_model.assembleMassLumped(); /// contact declaration Contact * my_contact = Contact::newContact(my_model, _ct_rigid, _cst_expli, _cnst_regular_grid); my_contact->initContact(false); Surface master = 1; my_contact->addMasterSurface(master); /* const RegularGridNeighborStructure<3> & my_rgns = dynamic_cast &>(my_contact->getContactSearch().getContactNeighborStructure(master)); const_cast&>(my_rgns).setGridSpacing(0.075, 0); const_cast&>(my_rgns).setGridSpacing(0.075, 1); const_cast&>(my_rgns).setGridSpacing(0.075, 2);*/ my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); my_contact->initSearch(); // does nothing so far // boundary conditions Surface impactor = 0; Vector * top_nodes = new Vector(0, 1); UInt middle_point; Real * coordinates = my_mesh.getNodes().values; Real * displacement = my_model.getDisplacement().values; bool * boundary = my_model.getBoundary().values; UInt * surface_to_nodes_offset = my_contact->getSurfaceToNodesOffset().values; UInt * surface_to_nodes = my_contact->getSurfaceToNodes().values; // symetry boundary conditions for(UInt n = surface_to_nodes_offset[impactor]; n < surface_to_nodes_offset[impactor+1]; ++n) { UInt node = surface_to_nodes[n]; Real x_coord = coordinates[node*dim]; Real y_coord = coordinates[node*dim + 1]; Real z_coord = coordinates[node*dim + 2]; if (x_coord < 0.00001) boundary[node*dim] = true; if (z_coord < 0.00001) boundary[node*dim+2] = true; if (y_coord > -0.00001) { boundary[node*dim + 1] = true; top_nodes->push_back(node); } if (x_coord < 0.00001 && y_coord > -0.00001 && z_coord < 0.00001) middle_point = node; } // ground boundary conditions for(UInt n = surface_to_nodes_offset[master]; n < surface_to_nodes_offset[master+1]; ++n) { UInt node = surface_to_nodes[n]; Real y_coord = coordinates[node*dim + 1]; if (y_coord < -1.2) boundary[node*dim] = true; boundary[node*dim + 1] = true; boundary[node*dim + 2] = true; } UInt * top_nodes_val = top_nodes->values; #ifdef AKANTU_USE_IOHELPER /// initialize the paraview output DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(my_model.getFEM().getMesh().getNodes().values, dim, nb_nodes, "coordinates_3d"); dumper.SetConnectivity((int *)my_model.getFEM().getMesh().getConnectivity(element_type).values, paraview_type, nb_element, C_MODE); dumper.AddNodeDataField(my_model.getDisplacement().values, dim, "displacements"); dumper.AddNodeDataField(my_model.getVelocity().values, dim, "velocity"); dumper.AddNodeDataField(my_model.getResidual().values, dim, "force"); dumper.AddNodeDataField(my_model.getForce().values, dim, "applied_force"); dumper.AddElemDataField(my_model.getMaterial(0).getStrain(element_type).values, dim*dim, "strain"); dumper.AddElemDataField(my_model.getMaterial(0).getStress(element_type).values, dim*dim, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/hertz_3d/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER std::ofstream hertz; hertz.open("hertz_3d.csv"); hertz << "%id,ftop,fcont,zone" << std::endl; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for(UInt s = 1; s <= max_steps; ++s) { if(s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; if(s % 250 == 0){ my_model.updateCurrentPosition(); my_contact->updateContact(); } if(s <= imposing_steps) { Real current_displacement = max_displacement/(static_cast(imposing_steps))*s; for(UInt n=0; ngetSize(); ++n) { UInt node = top_nodes_val[n]; displacement[node*dim + 1] = current_displacement; } } my_model.explicitPred(); - my_model.updateCurrentPosition(); + my_model.initializeUpdateResidualData(); /// compute the penetration list PenetrationList * my_penetration_list = new PenetrationList(); const_cast(my_contact->getContactSearch()).findPenetration(master, *my_penetration_list); UInt nb_nodes_pen = my_penetration_list->penetrating_nodes.getSize(); Vector pen_nodes = my_penetration_list->penetrating_nodes; UInt * pen_nodes_val = pen_nodes.values; my_contact->solveContact(); - my_model.updateResidual(); + my_model.updateResidual(false); Real * residual = my_model.getResidual().values; Real top_force = 0.; for(UInt n=0; ngetSize(); ++n) { UInt node = top_nodes_val[n]; top_force += residual[node*dim + 1]; } my_model.updateCurrentPosition(); Real * current_position = my_model.getCurrentPosition().values; Real contact_force = 0.; Real contact_zone = 0.; for (UInt i = 0; i < nb_nodes_pen; ++i) { UInt node = pen_nodes_val[i]; contact_force += residual[node*dim + 1]; contact_zone = std::max(contact_zone, current_position[node*dim]); } delete my_penetration_list; hertz << s << "," << top_force << "," << contact_force << "," << contact_zone << std::endl; my_model.updateAcceleration(); my_model.explicitCorr(); #ifdef AKANTU_USE_IOHELPER if(s % 100 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER } hertz.close(); delete my_contact; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_hertz_3d_full.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_hertz_3d_full.cc index 5416caf26..a06d96ef3 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_hertz_3d_full.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_hertz_3d_full.cc @@ -1,242 +1,242 @@ /** * @file test_contact_rigid_no_friction_hertz_3d_full.cc * @author David Kammer * @date Mon Jan 24 15:53:42 2011 * * @brief test rigid contact without friction for full 3d hertz in explicit * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" #include "contact_search.hh" #include "contact_search_explicit.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 3; const ElementType element_type = _tetrahedron_4; const UInt paraview_type = TETRA1; UInt max_steps = 200000; UInt imposing_steps = 100000; Real max_displacement = -0.1; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; mesh_io.read("hertz_3d_full.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); UInt nb_nodes = my_mesh.getNbNodes(); /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors memset(my_model.getForce().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getBoundary().values, false, dim*nb_nodes*sizeof(bool)); + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); UInt nb_element = my_model.getFEM().getMesh().getNbElement(element_type); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); my_model.assembleMassLumped(); /// contact declaration Contact * my_contact = Contact::newContact(my_model, _ct_rigid, _cst_expli, _cnst_regular_grid); my_contact->initContact(false); Surface master = 1; my_contact->addMasterSurface(master); /*const RegularGridNeighborStructure<3> & my_rgns = dynamic_cast &>(my_contact->getContactSearch().getContactNeighborStructure(master)); const_cast&>(my_rgns).setGridSpacing(0.075, 0); const_cast&>(my_rgns).setGridSpacing(0.075, 1); const_cast&>(my_rgns).setGridSpacing(0.075, 2);*/ my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); my_contact->initSearch(); // does nothing so far // boundary conditions Surface impactor = 0; Vector * top_nodes = new Vector(0, 1); Real * coordinates = my_mesh.getNodes().values; Real * displacement = my_model.getDisplacement().values; bool * boundary = my_model.getBoundary().values; UInt * surface_to_nodes_offset = my_contact->getSurfaceToNodesOffset().values; UInt * surface_to_nodes = my_contact->getSurfaceToNodes().values; // symetry boundary conditions for(UInt n = surface_to_nodes_offset[impactor]; n < surface_to_nodes_offset[impactor+1]; ++n) { UInt node = surface_to_nodes[n]; Real x_coord = coordinates[node*dim]; Real y_coord = coordinates[node*dim + 1]; Real z_coord = coordinates[node*dim + 2]; /*if (x_coord < 0.00001) boundary[node*dim] = true; if (z_coord < 0.00001) boundary[node*dim+2] = true;*/ if (y_coord > -0.00001) { boundary[node*dim + 1] = true; top_nodes->push_back(node); } } // ground boundary conditions for(UInt n = surface_to_nodes_offset[master]; n < surface_to_nodes_offset[master+1]; ++n) { UInt node = surface_to_nodes[n]; Real y_coord = coordinates[node*dim + 1]; if (y_coord < -1.2) boundary[node*dim] = true; boundary[node*dim + 1] = true; boundary[node*dim + 2] = true; } UInt * top_nodes_val = top_nodes->values; #ifdef AKANTU_USE_IOHELPER /// initialize the paraview output DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(my_model.getFEM().getMesh().getNodes().values, dim, nb_nodes, "coordinates_3d_full"); dumper.SetConnectivity((int *)my_model.getFEM().getMesh().getConnectivity(element_type).values, paraview_type, nb_element, C_MODE); dumper.AddNodeDataField(my_model.getDisplacement().values, dim, "displacements"); dumper.AddNodeDataField(my_model.getVelocity().values, dim, "velocity"); dumper.AddNodeDataField(my_model.getResidual().values, dim, "force"); dumper.AddNodeDataField(my_model.getForce().values, dim, "applied_force"); dumper.AddElemDataField(my_model.getMaterial(0).getStrain(element_type).values, dim*dim, "strain"); dumper.AddElemDataField(my_model.getMaterial(0).getStress(element_type).values, dim*dim, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/hertz_3d_full/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER std::ofstream hertz; hertz.open("hertz_3d_full.csv"); hertz << "%id,ftop,fcont,zone" << std::endl; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for(UInt s = 1; s <= max_steps; ++s) { if(s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; if(s % 200 == 0){ my_model.updateCurrentPosition(); my_contact->updateContact(); } if(s <= imposing_steps) { Real current_displacement = max_displacement/(static_cast(imposing_steps))*s; for(UInt n=0; ngetSize(); ++n) { UInt node = top_nodes_val[n]; displacement[node*dim + 1] = current_displacement; } } my_model.explicitPred(); - my_model.updateCurrentPosition(); + my_model.initializeUpdateResidualData(); /// compute the penetration list PenetrationList * my_penetration_list = new PenetrationList(); const_cast(my_contact->getContactSearch()).findPenetration(master, *my_penetration_list); UInt nb_nodes_pen = my_penetration_list->penetrating_nodes.getSize(); Vector pen_nodes = my_penetration_list->penetrating_nodes; UInt * pen_nodes_val = pen_nodes.values; my_contact->solveContact(); my_model.updateResidual(false); Real * residual = my_model.getResidual().values; Real top_force = 0.; for(UInt n=0; ngetSize(); ++n) { UInt node = top_nodes_val[n]; top_force += residual[node*dim + 1]; } my_model.updateCurrentPosition(); Real * current_position = my_model.getCurrentPosition().values; Real contact_force = 0.; Real contact_zone = 0.; for (UInt i = 0; i < nb_nodes_pen; ++i) { UInt node = pen_nodes_val[i]; contact_force += residual[node*dim + 1]; contact_zone = std::max(contact_zone, current_position[node*dim]); } delete my_penetration_list; hertz << s << "," << top_force << "," << contact_force << "," << contact_zone << std::endl; my_model.updateAcceleration(); my_model.explicitCorr(); #ifdef AKANTU_USE_IOHELPER if(s % 500 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER } hertz.close(); delete my_contact; delete top_nodes; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_tetrahedron_4.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_tetrahedron_4.cc index d7a42468b..e2b8b68f3 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_tetrahedron_4.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_tetrahedron_4.cc @@ -1,227 +1,227 @@ /** * @file test_contact_rigid_no_friction_tetrahedron_4.cc * @author David Kammer * @date Wed Jan 19 12:40:42 2011 * * @brief test rigid contact without friction for 3d case in explicit * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" #include "contact_search.hh" #include "contact_search_explicit.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 3; const ElementType element_type = _tetrahedron_4; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; mesh_io.read("cubes.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); UInt max_steps = 3; unsigned int nb_nodes = my_mesh.getNbNodes(); /// dump facet and surface information to paraview #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(my_mesh.getNodes().values, dim, nb_nodes, "tetrahedron_4_nodes_test-surface-extraction"); dumper.SetConnectivity((int*)my_mesh.getConnectivity(_tetrahedron_4).values, TETRA1, my_mesh.getNbElement(_tetrahedron_4), C_MODE); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors memset(my_model.getForce().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, dim*nb_nodes*sizeof(Real)); Real * displacement = my_model.getDisplacement().values; - + + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); my_model.assembleMassLumped(); /// contact declaration Contact * my_contact = Contact::newContact(my_model, _ct_rigid, _cst_expli, _cnst_regular_grid); my_contact->initContact(false); Surface master = 0; Surface impactor = 1; my_contact->addMasterSurface(master); my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); const NodesNeighborList & my_neighbor_list = dynamic_cast(my_contact->getContactSearch().getContactNeighborStructure(master).getNeighborList()); UInt nb_nodes_neigh = my_neighbor_list.impactor_nodes.getSize(); Vector impact_nodes = my_neighbor_list.impactor_nodes; UInt * impact_nodes_val = impact_nodes.values; /// print impactor nodes std::cout << "we have " << nb_nodes_neigh << " impactor nodes:"; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " " << impact_nodes_val[i]; } std::cout << std::endl; UInt * master_nodes_offset_val = my_neighbor_list.master_nodes_offset.values; UInt * master_nodes_val = my_neighbor_list.master_nodes.values; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " Impactor node: " << impact_nodes_val[i] << " has master nodes:"; for(UInt mn = master_nodes_offset_val[i]; mn < master_nodes_offset_val[i+1]; ++mn) { std::cout << " " << master_nodes_val[mn]; } std::cout << std::endl; } my_contact->initSearch(); // does nothing so far std::cout << std::endl << "epsilon = " << std::numeric_limits::epsilon() << std::endl; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for(UInt s = 1; s <= max_steps; ++s) { std::cout << std::endl << "passing step " << s << "/" << max_steps << std::endl; /// apply a displacement to the slave body if(s == 2) { Real * coord = my_mesh.getNodes().values; for(UInt n = 0; n < nb_nodes; ++n) { if(coord[n*dim + 2] > 1.0) { displacement[n*dim+2] = -0.01; } } /* UInt nb_elements = my_mesh.getNbElement(element_type); UInt nb_nodes_element = my_mesh.getNbNodesPerElement(element_type); Vector element_mat = my_model.getElementMaterial(element_type); UInt * element_mat_val = element_mat.values; UInt * connectivity = my_mesh.getConnectivity(element_type).values; for(UInt el = 0; el < nb_elements; ++el) { std::cout << "element: " << el << " with mat: " << element_mat_val[el] << std::endl; if(element_mat_val[el] == impactor) { for(UInt n = 0; n < nb_nodes_element; ++n) { displacement[connectivity[el * nb_nodes_element + n]+2] = -0.2; } } }*/ } /// central difference predictor my_model.explicitPred(); /// update current positions - my_model.updateCurrentPosition(); + my_model.initializeUpdateResidualData(); /// compute the penetration list std::cout << "Before solveContact" << std::endl; PenetrationList * my_penetration_list = new PenetrationList(); const_cast(my_contact->getContactSearch()).findPenetration(master, *my_penetration_list); UInt nb_nodes_pen = my_penetration_list->penetrating_nodes.getSize(); Vector pen_nodes = my_penetration_list->penetrating_nodes; UInt * pen_nodes_val = pen_nodes.values; std::cout << "we have " << nb_nodes_pen << " penetrating nodes:"; for (UInt i = 0; i < nb_nodes_pen; ++i) std::cout << " " << pen_nodes_val[i]; std::cout << std::endl; delete my_penetration_list; /// solve contact my_contact->solveContact(); /// compute the penetration list std::cout << "After solveContact" << std::endl; PenetrationList * my_penetration_list_2 = new PenetrationList(); const_cast(my_contact->getContactSearch()).findPenetration(master, *my_penetration_list_2); UInt nb_nodes_pen_2 = my_penetration_list_2->penetrating_nodes.getSize(); Vector pen_nodes_2 = my_penetration_list_2->penetrating_nodes; UInt * pen_nodes_2_val = pen_nodes_2.values; std::cout << "we have " << nb_nodes_pen_2 << " penetrating nodes:"; for (UInt i = 0; i < nb_nodes_pen_2; ++i) std::cout << " " << pen_nodes_2_val[i]; std::cout << std::endl; delete my_penetration_list_2; /// compute the residual my_model.updateResidual(false); /// compute the acceleration my_model.updateAcceleration(); /// central difference corrector my_model.explicitCorr(); } delete my_contact; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_triangle_3.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_triangle_3.cc index cccd7ab8a..a287803cf 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_triangle_3.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_no_friction_triangle_3.cc @@ -1,226 +1,226 @@ /** * @file test_contact_rigid_no_friction_triangle_3.cc * @author David Kammer * @date Mon Jan 17 11:13:42 2011 * * @brief test contact search for 2d case in explicit * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" #include "contact_search.hh" #include "contact_search_explicit.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 2; const ElementType element_type = _triangle_3; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; mesh_io.read("squares.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); UInt max_steps = 3; unsigned int nb_nodes = my_mesh.getNbNodes(); /// dump facet and surface information to paraview #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(my_mesh.getNodes().values, dim, nb_nodes, "triangle_3_nodes_test-surface-extraction"); dumper.SetConnectivity((int*)my_mesh.getConnectivity(_triangle_3).values, TETRA1, my_mesh.getNbElement(_triangle_3), C_MODE); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors memset(my_model.getForce().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, dim*nb_nodes*sizeof(Real)); Real * displacement = my_model.getDisplacement().values; - + + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); my_model.assembleMassLumped(); /// contact declaration Contact * my_contact = Contact::newContact(my_model, _ct_rigid, _cst_expli, _cnst_regular_grid); my_contact->initContact(false); Surface master = 0; Surface impactor = 1; my_contact->addMasterSurface(master); my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); const NodesNeighborList & my_neighbor_list = dynamic_cast(my_contact->getContactSearch().getContactNeighborStructure(master).getNeighborList()); UInt nb_nodes_neigh = my_neighbor_list.impactor_nodes.getSize(); Vector impact_nodes = my_neighbor_list.impactor_nodes; UInt * impact_nodes_val = impact_nodes.values; /// print impactor nodes std::cout << "we have " << nb_nodes_neigh << " impactor nodes:"; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " " << impact_nodes_val[i]; } std::cout << std::endl; UInt * master_nodes_offset_val = my_neighbor_list.master_nodes_offset.values; UInt * master_nodes_val = my_neighbor_list.master_nodes.values; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " Impactor node: " << impact_nodes_val[i] << " has master nodes:"; for(UInt mn = master_nodes_offset_val[i]; mn < master_nodes_offset_val[i+1]; ++mn) { std::cout << " " << master_nodes_val[mn]; } std::cout << std::endl; } my_contact->initSearch(); // does nothing so far std::cout << std::endl << "epsilon = " << std::numeric_limits::epsilon() << std::endl; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for(UInt s = 1; s <= max_steps; ++s) { std::cout << std::endl << "passing step " << s << "/" << max_steps << std::endl; /// apply a displacement to the slave body if(s == 2) { Real * coord = my_mesh.getNodes().values; for(UInt n = 0; n < nb_nodes; ++n) { if(coord[n*dim + 0] > 1.0) { displacement[n*dim+0] = -0.02; } } /* UInt nb_elements = my_mesh.getNbElement(element_type); UInt nb_nodes_element = my_mesh.getNbNodesPerElement(element_type); Vector element_mat = my_model.getElementMaterial(element_type); UInt * element_mat_val = element_mat.values; UInt * connectivity = my_mesh.getConnectivity(element_type).values; for(UInt el = 0; el < nb_elements; ++el) { std::cout << "element: " << el << " with mat: " << element_mat_val[el] << std::endl; if(element_mat_val[el] == impactor) { for(UInt n = 0; n < nb_nodes_element; ++n) { displacement[connectivity[el * nb_nodes_element + n]+2] = -0.2; } } }*/ } /// central difference predictor my_model.explicitPred(); /// update current position - my_model.updateCurrentPosition(); + my_model.initializeUpdateResidualData(); /// compute the penetration list std::cout << "Before solveContact" << std::endl; PenetrationList * my_penetration_list = new PenetrationList(); const_cast(my_contact->getContactSearch()).findPenetration(master, *my_penetration_list); UInt nb_nodes_pen = my_penetration_list->penetrating_nodes.getSize(); Vector pen_nodes = my_penetration_list->penetrating_nodes; UInt * pen_nodes_val = pen_nodes.values; std::cout << "we have " << nb_nodes_pen << " penetrating nodes:"; for (UInt i = 0; i < nb_nodes_pen; ++i) std::cout << " " << pen_nodes_val[i]; std::cout << std::endl; delete my_penetration_list; /// solve contact my_contact->solveContact(); /// compute the penetration list std::cout << "After solveContact" << std::endl; PenetrationList * my_penetration_list_2 = new PenetrationList(); const_cast(my_contact->getContactSearch()).findPenetration(master, *my_penetration_list_2); UInt nb_nodes_pen_2 = my_penetration_list_2->penetrating_nodes.getSize(); Vector pen_nodes_2 = my_penetration_list_2->penetrating_nodes; UInt * pen_nodes_2_val = pen_nodes_2.values; std::cout << "we have " << nb_nodes_pen_2 << " penetrating nodes:"; for (UInt i = 0; i < nb_nodes_pen_2; ++i) std::cout << " " << pen_nodes_2_val[i]; std::cout << std::endl; delete my_penetration_list_2; /// compute the residual my_model.updateResidual(false); /// compute the acceleration my_model.updateAcceleration(); /// central difference corrector my_model.explicitCorr(); } delete my_contact; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_sliding_cube_2d.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_sliding_cube_2d.cc index 9b24d2d95..b5a084cf1 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_sliding_cube_2d.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_sliding_cube_2d.cc @@ -1,311 +1,326 @@ /** * @file test_contact_rigid_sliding_cube_2d.cc * @author David Kammer * @date Tue Feb 08 15:56:42 2011 * * @brief test rigid contact for a 2d sliding cube in explicit * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_rigid.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" #include "contact_search.hh" #include "contact_search_explicit.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 2; const ElementType element_type = _triangle_3; const UInt paraview_type = TRIANGLE1; //UInt max_steps = 200000; UInt imposing_steps = 100000; Real max_displacement = -0.01; UInt damping_steps = 10000; UInt damping_interval = 100; Real damping_ratio = 0.5; Real imposed_velocity = 401.92; Real needed_time = 0.005; UInt additional_steps = 20000; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; mesh_io.read("sliding_cube_2d.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); UInt nb_nodes = my_mesh.getNbNodes(); /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors memset(my_model.getForce().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getBoundary().values, false, dim*nb_nodes*sizeof(bool)); + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); UInt nb_element = my_model.getFEM().getMesh().getNbElement(element_type); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); UInt max_steps = static_cast(needed_time*10/time_step) + imposing_steps + damping_steps + additional_steps; std::cout << "The number of time steps is found to be: " << max_steps << std::endl; Real * velocity_val = my_model.getVelocity().values; Real * acceleration_val = my_model.getAcceleration().values; my_model.assembleMassLumped(); /// contact declaration Contact * contact = Contact::newContact(my_model, _ct_rigid, _cst_expli, _cnst_regular_grid); ContactRigid * my_contact = dynamic_cast(contact); my_contact->initContact(false); Surface master = 1; my_contact->addMasterSurface(master); /*const RegularGridNeighborStructure<2> & my_rgns = dynamic_cast &>(my_contact->getContactSearch().getContactNeighborStructure(master)); const_cast&>(my_rgns).setGridSpacing(0.075, 0); const_cast&>(my_rgns).setGridSpacing(0.075, 1); */ my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); my_contact->initSearch(); // does nothing so far // boundary conditions Surface impactor = 0; Vector * top_nodes = new Vector(0, 1); Real * coordinates = my_mesh.getNodes().values; Real * displacement = my_model.getDisplacement().values; bool * boundary = my_model.getBoundary().values; UInt * surface_to_nodes_offset = my_contact->getSurfaceToNodesOffset().values; UInt * surface_to_nodes = my_contact->getSurfaceToNodes().values; // symetry boundary conditions for(UInt n = surface_to_nodes_offset[impactor]; n < surface_to_nodes_offset[impactor+1]; ++n) { UInt node = surface_to_nodes[n]; Real y_coord = coordinates[node*dim + 1]; if (y_coord > -0.00001) { boundary[node*dim + 1] = true; top_nodes->push_back(node); } } // ground boundary conditions for(UInt n = surface_to_nodes_offset[master]; n < surface_to_nodes_offset[master+1]; ++n) { UInt node = surface_to_nodes[n]; Real y_coord = coordinates[node*dim + 1]; - if (y_coord < -1.2) + if (y_coord < -1.49) boundary[node*dim] = true; boundary[node*dim + 1] = true; } UInt * top_nodes_val = top_nodes->values; #ifdef AKANTU_USE_IOHELPER /// initialize the paraview output DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(my_model.getFEM().getMesh().getNodes().values, dim, nb_nodes, "coord_sliding_cube_2d"); dumper.SetConnectivity((int *)my_model.getFEM().getMesh().getConnectivity(element_type).values, paraview_type, nb_element, C_MODE); dumper.AddNodeDataField(my_model.getDisplacement().values, dim, "displacements"); dumper.AddNodeDataField(my_model.getVelocity().values, dim, "velocity"); dumper.AddNodeDataField(my_model.getResidual().values, dim, "force"); dumper.AddNodeDataField(my_model.getForce().values, dim, "applied_force"); dumper.AddElemDataField(my_model.getMaterial(0).getStrain(element_type).values, dim*dim, "strain"); dumper.AddElemDataField(my_model.getMaterial(0).getStress(element_type).values, dim*dim, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/sliding_cube_2d/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER std::ofstream out_info; out_info.open("sliding_cube_2d.csv"); out_info << "%id,ftop,fcont,zone,stickNode,contNode" << std::endl; std::ofstream energy; energy.open("sliding_cube_2d_energy.csv"); energy << "%id,kin,pot,tot" << std::endl; Real * current_position = my_model.getCurrentPosition().values; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for(UInt s = 1; s <= max_steps; ++s) { if(s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; if(s % 200 == 0){ my_model.updateCurrentPosition(); my_contact->updateContact(); } // impose normal displacement if(s <= imposing_steps) { Real current_displacement = max_displacement/(static_cast(imposing_steps))*s; for(UInt n=0; ngetSize(); ++n) { UInt node = top_nodes_val[n]; displacement[node*dim + 1] = current_displacement; } } // damp velocity in order to find equilibrium if(s > imposing_steps && s < imposing_steps+damping_steps && s%damping_interval == 0) { for (UInt i=0; i < nb_nodes; ++i) { for (UInt j=0; j < dim; ++j) velocity_val[i*dim + j] *= damping_ratio; } } // give initial velocity if(s == imposing_steps+damping_steps) { - const Vector * act_imp_node = my_contact->getActiveImpactorNodes(); - bool * stick = my_contact->getNodeIsSticking()->values; - for(UInt i=0; i < act_imp_node->getSize(); ++i) { + Int master_index = -1; + for (UInt i=0; i < my_contact->getImpactorsInformation().size(); ++i) { + if (my_contact->getImpactorsInformation().at(i)->master_id == master) { + master_index = i; + break; + } + } + ContactRigid::ImpactorNodesInfoPerMaster * imp_info = my_contact->getImpactorsInformation().at(master_index); + bool * stick = imp_info->node_is_sticking->values; + for(UInt i=0; i < imp_info->active_impactor_nodes->getSize(); ++i) { stick[i*2] = false; stick[i*2+1] = false; } /*for (UInt i=0; i < nb_nodes; ++i) { if(coordinates[i*dim + 1] > -1) { velocity_val[i*dim] = imposed_velocity; } }*/ UInt * connectivity = my_mesh.getConnectivity(element_type).values; UInt nb_nodes_elem = Mesh::getNbNodesPerElement(element_type); for(UInt i=0; i < nb_element; ++i) { Real barycenter[dim]; Real * barycenter_p = &barycenter[0]; my_mesh.getBarycenter(i,element_type,barycenter_p); if(barycenter_p[1] > -1) { for (UInt j=0; j < nb_nodes_elem; ++j) { velocity_val[connectivity[i*nb_nodes_elem+j]*dim] = imposed_velocity; acceleration_val[connectivity[i*nb_nodes_elem +j]*dim] = 0.; } } } } my_model.explicitPred(); my_model.initializeUpdateResidualData(); my_contact->solveContact(); my_model.updateResidual(false); my_contact->avoidAdhesion(); my_contact->addFriction(); // find the total force applied at the imposed displacement surface (top) Real * residual = my_model.getResidual().values; Real top_force = 0.; for(UInt n=0; ngetSize(); ++n) { UInt node = top_nodes_val[n]; top_force += residual[node*dim + 1]; } + // find index of master surface in impactors_information + Int master_index = -1; + for (UInt i=0; i < my_contact->getImpactorsInformation().size(); ++i) { + if (my_contact->getImpactorsInformation().at(i)->master_id == master) { + master_index = i; + break; + } + } + // find the total contact force and contact area - const Vector * active_imp_nodes = my_contact->getActiveImpactorNodes(); - UInt * active_imp_nodes_val = active_imp_nodes->values; + ContactRigid::ImpactorNodesInfoPerMaster * imp_info = my_contact->getImpactorsInformation().at(master_index); + UInt * active_imp_nodes_val = imp_info->active_impactor_nodes->values; Real contact_force = 0.; Real contact_zone = 0.; - for (UInt i = 0; i < active_imp_nodes->getSize(); ++i) { + for (UInt i = 0; i < imp_info->active_impactor_nodes->getSize(); ++i) { UInt node = active_imp_nodes_val[i]; contact_force += residual[node*dim + 1]; contact_zone = std::max(contact_zone, current_position[node*dim]); } - //delete my_penetration_list; out_info << s << "," << top_force << "," << contact_force << "," << contact_zone << ","; my_model.updateAcceleration(); my_contact->addSticking(); - const Vector * sticking_nodes = my_contact->getNodeIsSticking(); + const Vector * sticking_nodes = imp_info->node_is_sticking; bool * sticking_nodes_val = sticking_nodes->values; UInt nb_sticking_nodes = 0; - for (UInt i = 0; i < active_imp_nodes->getSize(); ++i) { + for (UInt i = 0; i < imp_info->active_impactor_nodes->getSize(); ++i) { if(sticking_nodes_val[i*2]) nb_sticking_nodes++; } - out_info << nb_sticking_nodes << "," << active_imp_nodes->getSize() << std::endl; + out_info << nb_sticking_nodes << "," << imp_info->active_impactor_nodes->getSize() << std::endl; my_model.explicitCorr(); Real epot = my_model.getPotentialEnergy(); Real ekin = my_model.getKineticEnergy(); energy << s << "," << ekin << "," << epot << "," << ekin+epot << std::endl; #ifdef AKANTU_USE_IOHELPER if(s % 100 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER } out_info.close(); energy.close(); delete my_contact; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_sliding_cube_2d.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_sliding_cube_3d.cc similarity index 83% copy from test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_sliding_cube_2d.cc copy to test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_sliding_cube_3d.cc index 9b24d2d95..a4726e31f 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_sliding_cube_2d.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_rigid/test_contact_rigid_sliding_cube_3d.cc @@ -1,311 +1,318 @@ /** - * @file test_contact_rigid_sliding_cube_2d.cc + * @file test_contact_rigid_sliding_cube_3d.cc * @author David Kammer - * @date Tue Feb 08 15:56:42 2011 + * @date Mon Feb 14 10:56:42 2011 * - * @brief test rigid contact for a 2d sliding cube in explicit + * @brief test rigid contact for a 3d sliding cube in explicit * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_rigid.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" #include "contact_search.hh" #include "contact_search_explicit.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { - int dim = 2; - const ElementType element_type = _triangle_3; - const UInt paraview_type = TRIANGLE1; + int dim = 3; + const ElementType element_type = _tetrahedron_4; + const UInt paraview_type = TETRA1; //UInt max_steps = 200000; UInt imposing_steps = 100000; Real max_displacement = -0.01; UInt damping_steps = 10000; UInt damping_interval = 100; Real damping_ratio = 0.5; Real imposed_velocity = 401.92; Real needed_time = 0.005; UInt additional_steps = 20000; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; - mesh_io.read("sliding_cube_2d.msh", my_mesh); + mesh_io.read("sliding_cube_3d.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); UInt nb_nodes = my_mesh.getNbNodes(); /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors memset(my_model.getForce().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getBoundary().values, false, dim*nb_nodes*sizeof(bool)); + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); UInt nb_element = my_model.getFEM().getMesh().getNbElement(element_type); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); UInt max_steps = static_cast(needed_time*10/time_step) + imposing_steps + damping_steps + additional_steps; std::cout << "The number of time steps is found to be: " << max_steps << std::endl; Real * velocity_val = my_model.getVelocity().values; Real * acceleration_val = my_model.getAcceleration().values; my_model.assembleMassLumped(); /// contact declaration Contact * contact = Contact::newContact(my_model, _ct_rigid, _cst_expli, _cnst_regular_grid); ContactRigid * my_contact = dynamic_cast(contact); my_contact->initContact(false); Surface master = 1; my_contact->addMasterSurface(master); - /*const RegularGridNeighborStructure<2> & my_rgns = dynamic_cast &>(my_contact->getContactSearch().getContactNeighborStructure(master)); - const_cast&>(my_rgns).setGridSpacing(0.075, 0); - const_cast&>(my_rgns).setGridSpacing(0.075, 1); - */ my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); my_contact->initSearch(); // does nothing so far // boundary conditions Surface impactor = 0; Vector * top_nodes = new Vector(0, 1); Real * coordinates = my_mesh.getNodes().values; Real * displacement = my_model.getDisplacement().values; bool * boundary = my_model.getBoundary().values; UInt * surface_to_nodes_offset = my_contact->getSurfaceToNodesOffset().values; UInt * surface_to_nodes = my_contact->getSurfaceToNodes().values; - // symetry boundary conditions + // top boundary conditions -> imposing displacement for(UInt n = surface_to_nodes_offset[impactor]; n < surface_to_nodes_offset[impactor+1]; ++n) { UInt node = surface_to_nodes[n]; Real y_coord = coordinates[node*dim + 1]; if (y_coord > -0.00001) { boundary[node*dim + 1] = true; top_nodes->push_back(node); } } // ground boundary conditions for(UInt n = surface_to_nodes_offset[master]; n < surface_to_nodes_offset[master+1]; ++n) { UInt node = surface_to_nodes[n]; Real y_coord = coordinates[node*dim + 1]; - if (y_coord < -1.2) + if (y_coord < -1.49) boundary[node*dim] = true; boundary[node*dim + 1] = true; + boundary[node*dim + 2] = true; } UInt * top_nodes_val = top_nodes->values; #ifdef AKANTU_USE_IOHELPER /// initialize the paraview output DumperParaview dumper; dumper.SetMode(TEXT); - dumper.SetPoints(my_model.getFEM().getMesh().getNodes().values, dim, nb_nodes, "coord_sliding_cube_2d"); + dumper.SetPoints(my_model.getFEM().getMesh().getNodes().values, dim, nb_nodes, "coord_sliding_cube_3d"); dumper.SetConnectivity((int *)my_model.getFEM().getMesh().getConnectivity(element_type).values, paraview_type, nb_element, C_MODE); dumper.AddNodeDataField(my_model.getDisplacement().values, dim, "displacements"); dumper.AddNodeDataField(my_model.getVelocity().values, dim, "velocity"); dumper.AddNodeDataField(my_model.getResidual().values, dim, "force"); dumper.AddNodeDataField(my_model.getForce().values, dim, "applied_force"); dumper.AddElemDataField(my_model.getMaterial(0).getStrain(element_type).values, dim*dim, "strain"); dumper.AddElemDataField(my_model.getMaterial(0).getStress(element_type).values, dim*dim, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); - dumper.SetPrefix("paraview/sliding_cube_2d/"); + dumper.SetPrefix("paraview/sliding_cube_3d/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER std::ofstream out_info; - out_info.open("sliding_cube_2d.csv"); + out_info.open("sliding_cube_3d.csv"); out_info << "%id,ftop,fcont,zone,stickNode,contNode" << std::endl; std::ofstream energy; - energy.open("sliding_cube_2d_energy.csv"); + energy.open("sliding_cube_3d_energy.csv"); energy << "%id,kin,pot,tot" << std::endl; Real * current_position = my_model.getCurrentPosition().values; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for(UInt s = 1; s <= max_steps; ++s) { if(s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; if(s % 200 == 0){ my_model.updateCurrentPosition(); my_contact->updateContact(); } // impose normal displacement if(s <= imposing_steps) { Real current_displacement = max_displacement/(static_cast(imposing_steps))*s; for(UInt n=0; ngetSize(); ++n) { UInt node = top_nodes_val[n]; displacement[node*dim + 1] = current_displacement; } } // damp velocity in order to find equilibrium if(s > imposing_steps && s < imposing_steps+damping_steps && s%damping_interval == 0) { for (UInt i=0; i < nb_nodes; ++i) { for (UInt j=0; j < dim; ++j) velocity_val[i*dim + j] *= damping_ratio; } } // give initial velocity if(s == imposing_steps+damping_steps) { - const Vector * act_imp_node = my_contact->getActiveImpactorNodes(); - bool * stick = my_contact->getNodeIsSticking()->values; - for(UInt i=0; i < act_imp_node->getSize(); ++i) { + Int master_index = -1; + for (UInt i=0; i < my_contact->getImpactorsInformation().size(); ++i) { + if (my_contact->getImpactorsInformation().at(i)->master_id == master) { + master_index = i; + break; + } + } + ContactRigid::ImpactorNodesInfoPerMaster * imp_info = my_contact->getImpactorsInformation().at(master_index); + bool * stick = imp_info->node_is_sticking->values; + for(UInt i=0; i < imp_info->active_impactor_nodes->getSize(); ++i) { stick[i*2] = false; stick[i*2+1] = false; } - /*for (UInt i=0; i < nb_nodes; ++i) { - if(coordinates[i*dim + 1] > -1) { - velocity_val[i*dim] = imposed_velocity; - } - }*/ UInt * connectivity = my_mesh.getConnectivity(element_type).values; UInt nb_nodes_elem = Mesh::getNbNodesPerElement(element_type); for(UInt i=0; i < nb_element; ++i) { Real barycenter[dim]; Real * barycenter_p = &barycenter[0]; my_mesh.getBarycenter(i,element_type,barycenter_p); if(barycenter_p[1] > -1) { for (UInt j=0; j < nb_nodes_elem; ++j) { velocity_val[connectivity[i*nb_nodes_elem+j]*dim] = imposed_velocity; acceleration_val[connectivity[i*nb_nodes_elem +j]*dim] = 0.; } } } } my_model.explicitPred(); my_model.initializeUpdateResidualData(); my_contact->solveContact(); my_model.updateResidual(false); my_contact->avoidAdhesion(); my_contact->addFriction(); // find the total force applied at the imposed displacement surface (top) Real * residual = my_model.getResidual().values; Real top_force = 0.; for(UInt n=0; ngetSize(); ++n) { UInt node = top_nodes_val[n]; top_force += residual[node*dim + 1]; } - + + // find index of master surface in impactors_information + Int master_index = -1; + for (UInt i=0; i < my_contact->getImpactorsInformation().size(); ++i) { + if (my_contact->getImpactorsInformation().at(i)->master_id == master) { + master_index = i; + break; + } + } + // find the total contact force and contact area - const Vector * active_imp_nodes = my_contact->getActiveImpactorNodes(); - UInt * active_imp_nodes_val = active_imp_nodes->values; + ContactRigid::ImpactorNodesInfoPerMaster * imp_info = my_contact->getImpactorsInformation().at(master_index); + UInt * active_imp_nodes_val = imp_info->active_impactor_nodes->values; Real contact_force = 0.; Real contact_zone = 0.; - for (UInt i = 0; i < active_imp_nodes->getSize(); ++i) { + for (UInt i = 0; i < imp_info->active_impactor_nodes->getSize(); ++i) { UInt node = active_imp_nodes_val[i]; contact_force += residual[node*dim + 1]; contact_zone = std::max(contact_zone, current_position[node*dim]); } - //delete my_penetration_list; out_info << s << "," << top_force << "," << contact_force << "," << contact_zone << ","; my_model.updateAcceleration(); my_contact->addSticking(); - const Vector * sticking_nodes = my_contact->getNodeIsSticking(); + const Vector * sticking_nodes = imp_info->node_is_sticking; bool * sticking_nodes_val = sticking_nodes->values; UInt nb_sticking_nodes = 0; - for (UInt i = 0; i < active_imp_nodes->getSize(); ++i) { + for (UInt i = 0; i < imp_info->active_impactor_nodes->getSize(); ++i) { if(sticking_nodes_val[i*2]) nb_sticking_nodes++; } - out_info << nb_sticking_nodes << "," << active_imp_nodes->getSize() << std::endl; + out_info << nb_sticking_nodes << "," << imp_info->active_impactor_nodes->getSize() << std::endl; my_model.explicitCorr(); Real epot = my_model.getPotentialEnergy(); Real ekin = my_model.getKineticEnergy(); energy << s << "," << ekin << "," << epot << "," << ekin+epot << std::endl; #ifdef AKANTU_USE_IOHELPER if(s % 100 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER } out_info.close(); energy.close(); delete my_contact; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_contact_search_explicit/test_tetrahedron_4.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_contact_search_explicit/test_tetrahedron_4.cc index e1299595b..800d5e9ef 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_contact_search_explicit/test_tetrahedron_4.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_contact_search_explicit/test_tetrahedron_4.cc @@ -1,208 +1,208 @@ /** * @file test_tetrahedron_4.cc * @author David Kammer * @date Fri Dec 03 12:11:42 2010 * * @brief test contact_search_explicit for 3d case with tetrahedron_4 elements * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" #include "contact_search.hh" #include "contact_search_explicit.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 3; const ElementType element_type = _tetrahedron_4; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; mesh_io.read("cubes.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); UInt max_steps = 2; unsigned int nb_nodes = my_mesh.getNbNodes(); /// dump facet and surface information to paraview #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(my_mesh.getNodes().values, dim, nb_nodes, "tetrahedron_4_nodes_test-surface-extraction"); dumper.SetConnectivity((int*)my_mesh.getConnectivity(_tetrahedron_4).values, TETRA1, my_mesh.getNbElement(_tetrahedron_4), C_MODE); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors memset(my_model.getForce().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, dim*nb_nodes*sizeof(Real)); Real * displacement = my_model.getDisplacement().values; - + + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); my_model.assembleMassLumped(); /// contact declaration Contact * my_contact = Contact::newContact(my_model, _ct_3d_expli, _cst_expli, _cnst_regular_grid); my_contact->initContact(false); Surface master = 0; Surface impactor = 1; my_contact->addMasterSurface(master); my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); const NodesNeighborList & my_neighbor_list = dynamic_cast(my_contact->getContactSearch().getContactNeighborStructure(master).getNeighborList()); UInt nb_nodes_neigh = my_neighbor_list.impactor_nodes.getSize(); Vector impact_nodes = my_neighbor_list.impactor_nodes; UInt * impact_nodes_val = impact_nodes.values; /// print impactor nodes std::cout << "we have " << nb_nodes_neigh << " impactor nodes:"; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " " << impact_nodes_val[i]; } std::cout << std::endl; UInt * master_nodes_offset_val = my_neighbor_list.master_nodes_offset.values; UInt * master_nodes_val = my_neighbor_list.master_nodes.values; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " Impactor node: " << impact_nodes_val[i] << " has master nodes:"; for(UInt mn = master_nodes_offset_val[i]; mn < master_nodes_offset_val[i+1]; ++mn) { std::cout << " " << master_nodes_val[mn]; } std::cout << std::endl; } my_contact->initSearch(); // does nothing so far std::cout << std::endl << "epsilon = " << std::numeric_limits::epsilon() << std::endl; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for(UInt s = 1; s <= max_steps; ++s) { std::cout << std::endl << "passing step " << s << "/" << max_steps << std::endl; /// apply a displacement to the slave body if(s == 2) { Real * coord = my_mesh.getNodes().values; for(UInt n = 0; n < nb_nodes; ++n) { if(coord[n*dim + 2] > 1.0) { displacement[n*dim+2] = -0.01; } } /* UInt nb_elements = my_mesh.getNbElement(element_type); UInt nb_nodes_element = my_mesh.getNbNodesPerElement(element_type); Vector element_mat = my_model.getElementMaterial(element_type); UInt * element_mat_val = element_mat.values; UInt * connectivity = my_mesh.getConnectivity(element_type).values; for(UInt el = 0; el < nb_elements; ++el) { std::cout << "element: " << el << " with mat: " << element_mat_val[el] << std::endl; if(element_mat_val[el] == impactor) { for(UInt n = 0; n < nb_nodes_element; ++n) { displacement[connectivity[el * nb_nodes_element + n]+2] = -0.2; } } }*/ } /// central difference predictor my_model.explicitPred(); /// update current positions - my_model.updateCurrentPosition(); + my_model.initializeUpdateResidualData(); /// compute the penetration list PenetrationList * my_penetration_list = new PenetrationList(); const_cast(my_contact->getContactSearch()).findPenetration(master, *my_penetration_list); UInt nb_nodes_pen = my_penetration_list->penetrating_nodes.getSize(); Vector pen_nodes = my_penetration_list->penetrating_nodes; UInt * pen_nodes_val = pen_nodes.values; std::cout << "we have " << nb_nodes_pen << " penetrating nodes:"; for (UInt i = 0; i < nb_nodes_pen; ++i) std::cout << " " << pen_nodes_val[i]; std::cout << std::endl; delete my_penetration_list; /// compute the residual - my_model.updateResidual(); + my_model.updateResidual(false); /// compute the acceleration my_model.updateAcceleration(); /// central difference corrector my_model.explicitCorr(); } delete my_contact; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_contact_search_explicit/test_triangle_3.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_contact_search_explicit/test_triangle_3.cc index 3a6b53906..dc424a4b9 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_contact_search_explicit/test_triangle_3.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_contact_search/test_contact_search_explicit/test_triangle_3.cc @@ -1,208 +1,208 @@ /** * @file test_triangle_3.cc * @author David Kammer * @date Mon Jan 17 11:13:42 2011 * * @brief test contact search explicit for 2d case with triangle_3 elements * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" #include "contact_search.hh" #include "contact_search_explicit.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 2; const ElementType element_type = _triangle_3; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; mesh_io.read("squares.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); UInt max_steps = 2; unsigned int nb_nodes = my_mesh.getNbNodes(); /// dump facet and surface information to paraview #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(my_mesh.getNodes().values, dim, nb_nodes, "triangle_3_nodes_test-surface-extraction"); dumper.SetConnectivity((int*)my_mesh.getConnectivity(_triangle_3).values, TETRA1, my_mesh.getNbElement(_triangle_3), C_MODE); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors memset(my_model.getForce().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, dim*nb_nodes*sizeof(Real)); Real * displacement = my_model.getDisplacement().values; - + + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); my_model.assembleMassLumped(); /// contact declaration Contact * my_contact = Contact::newContact(my_model, _ct_3d_expli, _cst_expli, _cnst_regular_grid); my_contact->initContact(false); Surface master = 0; Surface impactor = 1; my_contact->addMasterSurface(master); my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); const NodesNeighborList & my_neighbor_list = dynamic_cast(my_contact->getContactSearch().getContactNeighborStructure(master).getNeighborList()); UInt nb_nodes_neigh = my_neighbor_list.impactor_nodes.getSize(); Vector impact_nodes = my_neighbor_list.impactor_nodes; UInt * impact_nodes_val = impact_nodes.values; /// print impactor nodes std::cout << "we have " << nb_nodes_neigh << " impactor nodes:"; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " " << impact_nodes_val[i]; } std::cout << std::endl; UInt * master_nodes_offset_val = my_neighbor_list.master_nodes_offset.values; UInt * master_nodes_val = my_neighbor_list.master_nodes.values; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " Impactor node: " << impact_nodes_val[i] << " has master nodes:"; for(UInt mn = master_nodes_offset_val[i]; mn < master_nodes_offset_val[i+1]; ++mn) { std::cout << " " << master_nodes_val[mn]; } std::cout << std::endl; } my_contact->initSearch(); // does nothing so far std::cout << std::endl << "epsilon = " << std::numeric_limits::epsilon() << std::endl; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for(UInt s = 1; s <= max_steps; ++s) { std::cout << std::endl << "passing step " << s << "/" << max_steps << std::endl; /// apply a displacement to the slave body if(s == 2) { Real * coord = my_mesh.getNodes().values; for(UInt n = 0; n < nb_nodes; ++n) { if(coord[n*dim + 0] > 1.0) { displacement[n*dim+0] = -0.02; } } /* UInt nb_elements = my_mesh.getNbElement(element_type); UInt nb_nodes_element = my_mesh.getNbNodesPerElement(element_type); Vector element_mat = my_model.getElementMaterial(element_type); UInt * element_mat_val = element_mat.values; UInt * connectivity = my_mesh.getConnectivity(element_type).values; for(UInt el = 0; el < nb_elements; ++el) { std::cout << "element: " << el << " with mat: " << element_mat_val[el] << std::endl; if(element_mat_val[el] == impactor) { for(UInt n = 0; n < nb_nodes_element; ++n) { displacement[connectivity[el * nb_nodes_element + n]+2] = -0.2; } } }*/ } /// central difference predictor my_model.explicitPred(); /// update current position - my_model.updateCurrentPosition(); + my_model.initializeUpdateResidualData(); /// compute the penetration list PenetrationList * my_penetration_list = new PenetrationList(); const_cast(my_contact->getContactSearch()).findPenetration(master, *my_penetration_list); UInt nb_nodes_pen = my_penetration_list->penetrating_nodes.getSize(); Vector pen_nodes = my_penetration_list->penetrating_nodes; UInt * pen_nodes_val = pen_nodes.values; std::cout << "we have " << nb_nodes_pen << " penetrating nodes:"; for (UInt i = 0; i < nb_nodes_pen; ++i) std::cout << " " << pen_nodes_val[i]; std::cout << std::endl; delete my_penetration_list; /// compute the residual my_model.updateResidual(false); /// compute the acceleration my_model.updateAcceleration(); /// central difference corrector my_model.explicitCorr(); } delete my_contact; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/CMakeLists.txt index 6de527372..01748232d 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/CMakeLists.txt +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/CMakeLists.txt @@ -1,53 +1,59 @@ #=============================================================================== # @file CMakeLists.txt # @author David Kammer # @date Tue Oct 26 16:40:24 2010 # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== #=============================================================================== add_mesh(test_regular_grid_2d_mesh squares.geo 2 1) register_test(test_regular_grid_triangle_3 test_regular_grid_triangle_3.cc) add_dependencies(test_regular_grid_triangle_3 test_regular_grid_2d_mesh) #=============================================================================== add_mesh(test_regular_grid_3d_mesh cubes.geo 3 1) register_test(test_regular_grid_tetrahedron_4 test_regular_grid_tetrahedron_4.cc) add_dependencies(test_regular_grid_tetrahedron_4 test_regular_grid_3d_mesh) #=============================================================================== #add_mesh(test_regular_grid_3d_mesh cubes.geo 3 1) register_test(test_regular_grid_tetrahedron_4_nodes test_regular_grid_tetrahedron_4_nodes.cc) add_dependencies(test_regular_grid_tetrahedron_4_nodes test_regular_grid_3d_mesh) +#=============================================================================== +add_mesh(test_regular_grid_3d_mesh_no_neighbors cube.geo 3 1) + +register_test(test_regular_grid_tetrahedron_4_nodes_no_neighbors test_regular_grid_tetrahedron_4_nodes_no_neighbors.cc) +add_dependencies(test_regular_grid_tetrahedron_4_nodes_no_neighbors test_regular_grid_3d_mesh_no_neighbors) + #=============================================================================== configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/material.dat ${CMAKE_CURRENT_BINARY_DIR}/material.dat COPYONLY ) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/cube.geo b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/cube.geo new file mode 100644 index 000000000..3a705e308 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/cube.geo @@ -0,0 +1,39 @@ +h = 0.05; + +Point(1) = {0, 0, 0, h}; +Point(2) = {1, 0, 0, h}; +Point(3) = {1, 1, 0, h}; +Point(4) = {0, 1, 0, h}; +Point(5) = {0, 0, 1, h}; +Point(6) = {1, 0, 1, h}; +Point(7) = {1, 1, 1, h}; +Point(8) = {0, 1, 1, h}; + +Line(1) = {5, 6}; +Line(2) = {6, 7}; +Line(3) = {7, 8}; +Line(4) = {8, 5}; +Line(5) = {1, 2}; +Line(6) = {2, 3}; +Line(7) = {3, 4}; +Line(8) = {4, 1}; +Line(9) = {5, 1}; +Line(10) = {6, 2}; +Line(11) = {7, 3}; +Line(12) = {8, 4}; + +Line Loop(25) = {4, 1, 2, 3}; +Plane Surface(26) = {25}; +Line Loop(27) = {7, 8, 5, 6}; +Plane Surface(28) = {27}; +Line Loop(29) = {12, -7, -11, 3}; +Plane Surface(30) = {29}; +Line Loop(31) = {12, 8, -9, -4}; +Plane Surface(32) = {31}; +Line Loop(33) = {1, 10, -5, -9}; +Plane Surface(34) = {33}; +Line Loop(35) = {2, 11, -6, -10}; +Plane Surface(36) = {35}; + +Surface Loop(49) = {32, 30, 28, 34, 26, 36}; +Volume(50) = {49}; diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_tetrahedron_4.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_tetrahedron_4.cc index a0f14989d..ed567ea72 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_tetrahedron_4.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_tetrahedron_4.cc @@ -1,160 +1,160 @@ /** * @file test_regular_grid_tetrahedron_4.cc * @author David Kammer * @date Tue Oct 26 16:58:42 2010 * * @brief test regular grid neighbor structure for 3d case * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 3; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; mesh_io.read("cubes.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); unsigned int nb_nodes = my_mesh.getNbNodes(); /// dump facet information to paraview #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(my_mesh.getNodes().values, dim, nb_nodes, "tetrahedron_4_test-surface-extraction"); dumper.SetConnectivity((int*)my_mesh.getConnectivity(_tetrahedron_4).values, TETRA1, my_mesh.getNbElement(_tetrahedron_4), C_MODE); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors memset(my_model.getForce().values, 0, 3*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, 3*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, 3*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, 3*nb_nodes*sizeof(Real)); - + + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); my_model.assembleMassLumped(); /// contact declaration Contact * my_contact = Contact::newContact(my_model, _ct_3d_expli, _cst_2d_expli, _cnst_regular_grid); // how to use contact and contact search types for testing the reg grid with normal nl? my_contact->initContact(false); Surface master = 0; my_contact->addMasterSurface(master); my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); const NeighborList & my_neighbor_list = my_contact->getContactSearch().getContactNeighborStructure(master).getNeighborList(); UInt nb_nodes_neigh = my_neighbor_list.impactor_nodes.getSize(); Vector impact_nodes = my_neighbor_list.impactor_nodes; UInt * impact_nodes_val = impact_nodes.values; UInt * node_to_elem_offset_val = my_neighbor_list.facets_offset[_triangle_3]->values; UInt * node_to_elem_val = my_neighbor_list.facets[_triangle_3]->values; /// print impactor nodes std::cout << "we have " << nb_nodes_neigh << " impactor nodes:" << std::endl; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " node " << impact_nodes_val[i] << " : "; for (UInt j = node_to_elem_offset_val[i]; j < node_to_elem_offset_val[i+1]; ++j) std::cout << node_to_elem_val[j] << " "; std::cout << std::endl; } std::cout << std::endl; #ifdef AKANTU_USE_IOHELPER DumperParaview dumper_neighbor; dumper_neighbor.SetMode(TEXT); dumper_neighbor.SetPoints(my_mesh.getNodes().values, dim, nb_nodes, "tetrahedron_4_test-neighbor-elements"); dumper_neighbor.SetConnectivity((int *)my_mesh.getConnectivity(_triangle_3).values, TRIANGLE1, my_mesh.getNbElement(_triangle_3), C_MODE); double * neigh_elem = new double [my_mesh.getNbElement(_triangle_3)]; for (UInt i = 0; i < my_mesh.getNbElement(_triangle_3); ++i) neigh_elem[i] = 0.0; UInt visualize_node = 7; UInt n = impact_nodes_val[visualize_node]; std::cout << "plot for node: " << n << std::endl; for (UInt i = node_to_elem_offset_val[visualize_node]; i < node_to_elem_offset_val[visualize_node+1]; ++i) neigh_elem[node_to_elem_val[i]] = 1.; dumper_neighbor.AddElemDataField(neigh_elem, 1, "neighbor id"); dumper_neighbor.SetPrefix("paraview/"); dumper_neighbor.Init(); dumper_neighbor.Dump(); delete [] neigh_elem; #endif //AKANTU_USE_IOHELPER delete my_contact; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_tetrahedron_4_nodes.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_tetrahedron_4_nodes.cc index 63862fa21..30b4750cc 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_tetrahedron_4_nodes.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_tetrahedron_4_nodes.cc @@ -1,141 +1,141 @@ /** * @file test_regular_grid_tetrahedron_4_nodes.cc * @author David Kammer * @date Tue Oct 26 16:58:42 2010 * * @brief test regular grid neighbor structure for 3d case * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 3; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; mesh_io.read("cubes.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); unsigned int nb_nodes = my_mesh.getNbNodes(); /// dump facet and surface information to paraview #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(my_mesh.getNodes().values, dim, nb_nodes, "tetrahedron_4_nodes_test-surface-extraction"); dumper.SetConnectivity((int*)my_mesh.getConnectivity(_tetrahedron_4).values, TETRA1, my_mesh.getNbElement(_tetrahedron_4), C_MODE); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors memset(my_model.getForce().values, 0, 3*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, 3*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, 3*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, 3*nb_nodes*sizeof(Real)); - + + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); my_model.assembleMassLumped(); /// contact declaration Contact * my_contact = Contact::newContact(my_model, _ct_3d_expli, _cst_expli, _cnst_regular_grid); my_contact->initContact(false); Surface master = 0; my_contact->addMasterSurface(master); my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); const NodesNeighborList & my_neighbor_list = dynamic_cast(my_contact->getContactSearch().getContactNeighborStructure(master).getNeighborList()); UInt nb_nodes_neigh = my_neighbor_list.impactor_nodes.getSize(); Vector impact_nodes = my_neighbor_list.impactor_nodes; UInt * impact_nodes_val = impact_nodes.values; /// print impactor nodes std::cout << "we have " << nb_nodes_neigh << " impactor nodes:"; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " " << impact_nodes_val[i]; } std::cout << std::endl; UInt * master_nodes_offset_val = my_neighbor_list.master_nodes_offset.values; UInt * master_nodes_val = my_neighbor_list.master_nodes.values; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " Impactor node: " << impact_nodes_val[i] << " has master nodes:"; for(UInt mn = master_nodes_offset_val[i]; mn < master_nodes_offset_val[i+1]; ++mn) { std::cout << " " << master_nodes_val[mn]; } std::cout << std::endl; } delete my_contact; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_tetrahedron_4_nodes.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_tetrahedron_4_nodes_no_neighbors.cc similarity index 93% copy from test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_tetrahedron_4_nodes.cc copy to test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_tetrahedron_4_nodes_no_neighbors.cc index 63862fa21..059f6421e 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_tetrahedron_4_nodes.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_tetrahedron_4_nodes_no_neighbors.cc @@ -1,141 +1,142 @@ /** - * @file test_regular_grid_tetrahedron_4_nodes.cc + * @file test_regular_grid_tetrahedron_4_nodes_no_neighbors.cc * @author David Kammer * @date Tue Oct 26 16:58:42 2010 * - * @brief test regular grid neighbor structure for 3d case + * @brief test regular grid neighbor structure for 3d case if there are no + * neighbors * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_neighbor_structure.hh" #include "regular_grid_neighbor_structure.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 3; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; - mesh_io.read("cubes.msh", my_mesh); + mesh_io.read("cube.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); unsigned int nb_nodes = my_mesh.getNbNodes(); /// dump facet and surface information to paraview #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(TEXT); - dumper.SetPoints(my_mesh.getNodes().values, dim, nb_nodes, "tetrahedron_4_nodes_test-surface-extraction"); + dumper.SetPoints(my_mesh.getNodes().values, dim, nb_nodes, "tetrahedron_4_nodes_no_neighbors_test-surface-extraction"); dumper.SetConnectivity((int*)my_mesh.getConnectivity(_tetrahedron_4).values, TETRA1, my_mesh.getNbElement(_tetrahedron_4), C_MODE); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors memset(my_model.getForce().values, 0, 3*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, 3*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, 3*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, 3*nb_nodes*sizeof(Real)); - + + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); my_model.assembleMassLumped(); /// contact declaration Contact * my_contact = Contact::newContact(my_model, _ct_3d_expli, _cst_expli, _cnst_regular_grid); my_contact->initContact(false); Surface master = 0; my_contact->addMasterSurface(master); my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); const NodesNeighborList & my_neighbor_list = dynamic_cast(my_contact->getContactSearch().getContactNeighborStructure(master).getNeighborList()); UInt nb_nodes_neigh = my_neighbor_list.impactor_nodes.getSize(); Vector impact_nodes = my_neighbor_list.impactor_nodes; UInt * impact_nodes_val = impact_nodes.values; /// print impactor nodes std::cout << "we have " << nb_nodes_neigh << " impactor nodes:"; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " " << impact_nodes_val[i]; } std::cout << std::endl; UInt * master_nodes_offset_val = my_neighbor_list.master_nodes_offset.values; UInt * master_nodes_val = my_neighbor_list.master_nodes.values; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " Impactor node: " << impact_nodes_val[i] << " has master nodes:"; for(UInt mn = master_nodes_offset_val[i]; mn < master_nodes_offset_val[i+1]; ++mn) { std::cout << " " << master_nodes_val[mn]; } std::cout << std::endl; } delete my_contact; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_triangle_3.cc b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_triangle_3.cc index 151da8a2f..9f48e4f10 100644 --- a/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_triangle_3.cc +++ b/test/test_model/test_solid_mechanics_model/test_contact/test_neighbor_structure/test_regular_grid_neighbor_structure/test_regular_grid_triangle_3.cc @@ -1,156 +1,156 @@ /** * @file test_regular_grid_triangle_3.cc * @author David Kammer * @date Tue Oct 26 16:58:42 2010 * * @brief test regular grid neighbor structure for 3d case * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "contact.hh" #include "contact_neighbor_structure.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 2; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; mesh_io.read("squares.msh", my_mesh); /// build facet connectivity and surface id MeshUtils::buildFacets(my_mesh,1,0); MeshUtils::buildSurfaceID(my_mesh); unsigned int nb_nodes = my_mesh.getNbNodes(); /// dump surface information to paraview #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(my_mesh.getNodes().values, dim, nb_nodes, "triangle_3_test-surface-extraction"); dumper.SetConnectivity((int*)my_mesh.getConnectivity(_triangle_3).values, TRIANGLE1, my_mesh.getNbElement(_triangle_3), C_MODE); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); /// initialize the vectors memset(my_model.getForce().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getVelocity().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getAcceleration().values, 0, dim*nb_nodes*sizeof(Real)); memset(my_model.getDisplacement().values, 0, dim*nb_nodes*sizeof(Real)); + my_model.initModel(); my_model.readMaterials("material.dat"); my_model.initMaterials(); - my_model.initModel(); Real time_step = my_model.getStableTimeStep(); my_model.setTimeStep(time_step/10.); my_model.assembleMassLumped(); /// contact declaration Contact * my_contact = Contact::newContact(my_model, _ct_2d_expli, _cst_2d_expli, _cnst_regular_grid); my_contact->initContact(false); Surface master = 0; my_contact->addMasterSurface(master); my_model.updateCurrentPosition(); // neighbor structure uses current position for init my_contact->initNeighborStructure(master); const NeighborList & my_neighbor_list = my_contact->getContactSearch().getContactNeighborStructure(master).getNeighborList(); UInt nb_nodes_neigh = my_neighbor_list.impactor_nodes.getSize(); Vector impact_nodes = my_neighbor_list.impactor_nodes; UInt * impact_nodes_val = impact_nodes.values; UInt * node_to_elem_offset_val = my_neighbor_list.facets_offset[_segment_2]->values; UInt * node_to_elem_val = my_neighbor_list.facets[_segment_2]->values; /// print impactor nodes std::cout << "we have " << nb_nodes_neigh << " impactor nodes:" << std::endl; for (UInt i = 0; i < nb_nodes_neigh; ++i) { std::cout << " node " << impact_nodes_val[i] << " : "; for (UInt j = node_to_elem_offset_val[i]; j < node_to_elem_offset_val[i+1]; ++j) std::cout << node_to_elem_val[j] << " "; std::cout << std::endl; } std::cout << std::endl; #ifdef AKANTU_USE_IOHELPER DumperParaview dumper_neighbor; dumper_neighbor.SetMode(TEXT); dumper_neighbor.SetPoints(my_mesh.getNodes().values, dim, nb_nodes, "triangle_3_test-neighbor-elements"); dumper_neighbor.SetConnectivity((int *)my_mesh.getConnectivity(_segment_2).values, LINE1, my_mesh.getNbElement(_segment_2), C_MODE); double * neigh_elem = new double [my_mesh.getNbElement(_segment_2)]; for (UInt i = 0; i < my_mesh.getNbElement(_segment_2); ++i) neigh_elem[i] = 0.0; UInt visualize_node = 1; UInt n = impact_nodes_val[visualize_node]; std::cout << "plot for node: " << n << std::endl; for (UInt i = node_to_elem_offset_val[visualize_node]; i < node_to_elem_offset_val[visualize_node+1]; ++i) neigh_elem[node_to_elem_val[i]] = 1.; dumper_neighbor.AddElemDataField(neigh_elem, 1, "neighbor id"); dumper_neighbor.SetPrefix("paraview/"); dumper_neighbor.Init(); dumper_neighbor.Dump(); delete [] neigh_elem; #endif //AKANTU_USE_IOHELPER finalize(); return EXIT_SUCCESS; }