Page MenuHomec4science
Files F91213948

diagonal_compression.cc
No OneTemporary
Actions

File Metadata

Created
Sat, Nov 9, 01:13

diagonal_compression.cc
View Options

#include <chrono>
#include <coupler_solid_cohesive_contact.hh>
#include <material_selector.hh>
#include <material_selector_cohesive.hh>
using clk = std::chrono::high_resolution_clock;
using seconds = std::chrono::duration<double>;
using milliseconds = std::chrono::duration<double, std::milli>;
using hours = std::chrono::duration<double, std::ratio<3600>>;
using namespace akantu;
int main(int argc, char * argv[]) {
initialize("wall.dat", argc, argv);
const auto & user_param = getUserParser();
Vector<Real> wall_sizes = user_param.getParameter("wall_sizes");
Vector<Real> clamp_sizes = user_param.getParameter("clamp_sizes");
int settle_steps = user_param.getParameter("settle_steps");
int max_steps = user_param.getParameter("max_steps");
double displacement_increment =
user_param.getParameter("displacement_increment");
int damping_interval = user_param.getParameter("damping_interval");
double damping_ratio = user_param.getParameter("damping_ratio");
int dumping_interval = user_param.getParameter("dumping_interval");
int spatial_dimension = 3;
auto surface_element_type = _triangle_3;
Mesh mesh(spatial_dimension);
std::string mesh_file = user_param.getParameter("mesh");
mesh.read(mesh_file);
auto & mesh_facets = mesh.initMeshFacets();
auto & contact_group =
mesh_facets.createElementGroup("contact", spatial_dimension);
std::vector<std::string> group_list{
"clamp_top_inside",
"clamp_bottom_inside",
"wall_top",
"wall_bottom",
};
for (auto && gid : group_list) {
contact_group.append(mesh_facets.getElementGroup(gid));
}
CouplerSolidCohesiveContact coupler(mesh);
auto & solid = coupler.getSolidMechanicsModelCohesive();
auto & contact = coupler.getContactMechanicsModel();
MaterialCohesiveRules rules{
{{"stones", "stones"}, "stone_stone"},
{{"stones", "mortar"}, "stone_mortar"},
{{"mortar", "mortar"}, "mortar_mortar"},
{{"clamp_top", "clamp_top"}, "clamp_clamp"},
{{"clamp_bottom", "clamp_bottom"}, "clamp_clamp"}};
auto material_selector =
std::make_shared<MaterialCohesiveRulesSelector>(solid, rules);
solid.setMaterialSelector(material_selector);
coupler.initFull(_is_extrinsic = true);
auto surface_selector = std::make_shared<AllSurfaceSelector>(mesh);
contact.getContactDetector().setSurfaceSelector(surface_selector);
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "clamp_top");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "clamp_top");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _z), "clamp_top");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "clamp_bottom");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "clamp_bottom");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _z), "clamp_bottom");
auto time_step = solid.getStableTimeStep();
time_step *= 0.1;
solid.assembleMassLumped();
auto & force = solid.getExternalForce();
auto & mass = solid.getMass();
Vector<Real> direction_down{0., -1., -1.};
direction_down.normalize();
for (auto && data : zip(make_view(force, spatial_dimension),
make_view(mass, spatial_dimension))) {
auto & force = std::get<0>(data);
auto & mass = std::get<1>(data);
force = 9.81 * mass(_x) * direction_down;
}
coupler.setTimeStep(time_step);
solid.setBaseName("wall");
solid.addDumpField("displacement");
solid.addDumpField("velocity");
solid.addDumpField("mass");
solid.addDumpField("external_force");
solid.addDumpField("internal_force");
solid.addDumpField("blocked_dofs");
solid.addDumpField("grad_u");
solid.addDumpField("material_index");
solid.addDumpField("stress");
solid.addDumpFieldVectorToDumper("cohesive elements", "displacement");
solid.addDumpFieldToDumper("cohesive elements", "damage");
solid.addDumpFieldToDumper("cohesive elements", "tractions");
solid.addDumpFieldToDumper("cohesive elements", "opening");
contact.addDumpField("contact_force");
contact.addDumpField("gaps");
contact.addDumpField("areas");
solid.dump();
contact.dump();
solid.dump("cohesive elements");
std::cout << "Materials: \n"
<< " - clamps: "
<< solid.getMaterial("clamp_top").getElementFilter().size() << " + "
<< solid.getMaterial("clamp_bottom").getElementFilter().size()
<< "\n"
<< " - mortar: "
<< solid.getMaterial("mortar").getElementFilter().size() << "\n"
<< " - stones: "
<< solid.getMaterial("stones").getElementFilter().size() << "\n"
<< std::endl;
auto & velocity = solid.getVelocity();
auto increment_vect = displacement_increment * direction_down;
std::cout << increment_vect << std::endl;
std::cout << "settle step " << 0 << "/" << settle_steps << "----- ms/step"
<< "\t\t \r";
auto start_time = clk::now();
for (auto s : arange(settle_steps)) {
// coupler.checkCohesiveStress();
coupler.solveStep();
if (s % 100 == 0) {
velocity *= damping_ratio;
}
if (s % 100 == 0) {
milliseconds loop_time = clk::now() - start_time;
std::cout << "settle step " << s << "/" << settle_steps << " - "
<< loop_time.count() / s << "ms/step"
<< "\t\t \r";
std::cout << std::flush;
}
}
solid.dump();
contact.dump();
solid.dump("cohesive elements");
auto nb_element = mesh.getNbElement(spatial_dimension);
std::cout << "passing step " << 0 << "/" << max_steps
<< " - nb_element: " << nb_element << " - nb_cohesive_element: "
<< mesh.getNbElement(spatial_dimension, _not_ghost, _ek_cohesive)
<< " - nb_dofs: " << mesh.getNbGlobalNodes() * spatial_dimension
<< " - "
<< "----- ms/step"
<< " - remaining: inf h"
<< "\t\t \r";
start_time = clk::now();
for (auto s : arange(max_steps)) {
coupler.applyBC(BC::Dirichlet::IncrementValue(increment_vect(_x), _x),
"clamp_top");
coupler.applyBC(BC::Dirichlet::IncrementValue(increment_vect(_y), _y),
"clamp_top");
coupler.applyBC(BC::Dirichlet::IncrementValue(increment_vect(_z), _z),
"clamp_top");
coupler.checkCohesiveStress();
coupler.solveStep();
if (s % damping_interval == 0) {
velocity *= damping_ratio;
}
if (s % 100 == 0) {
auto loop_time = clk::now() - start_time;
auto remaining = std::chrono::duration_cast<hours>(loop_time / (s + 1) *
(max_steps - s - 1));
std::cout << "passing step " << s << "/" << max_steps
<< " - physical_time: " << s * time_step
<< " - nb_element: " << nb_element << " - nb_cohesive_element: "
<< mesh.getNbElement(spatial_dimension, _not_ghost,
_ek_cohesive)
<< " - nb_dofs: " << mesh.getNbGlobalNodes() * spatial_dimension
<< " - "
<< std::chrono::duration_cast<milliseconds>(loop_time).count() /
(s + 1)
<< "ms/step"
<< " - remaining: " << remaining.count() << " h"
<< "\t\t \r";
std::cout << std::flush;
}
if (s % dumping_interval == 0) {
solid.dump();
contact.dump();
solid.dump("cohesive elements");
}
}
auto loop_time = clk::now() - start_time;
std::cout << "passing step " << max_steps << "/" << max_steps
<< " - physical_time: " << max_steps * time_step
<< " - nb_element: " << nb_element << " - nb_cohesive_element: "
<< mesh.getNbElement(spatial_dimension, _not_ghost, _ek_cohesive)
<< " - nb_dofs: " << mesh.getNbGlobalNodes() * spatial_dimension
<< " - "
<< std::chrono::duration_cast<milliseconds>(loop_time).count() /
max_steps
<< "ms/step"
<< "\t\t \n";
std::cout << std::endl;
for (auto s : arange(settle_steps)) {
coupler.applyBC(BC::Dirichlet::IncrementValue(-increment_vect(_x), _x),
"clamp_top");
coupler.applyBC(BC::Dirichlet::IncrementValue(-increment_vect(_y), _y),
"clamp_top");
coupler.applyBC(BC::Dirichlet::IncrementValue(-increment_vect(_z), _z),
"clamp_top");
// coupler.checkCohesiveStress();
coupler.solveStep();
if (s % 10 == 0) {
velocity *= damping_ratio;
}
}
solid.dump();
contact.dump();
solid.dump("cohesive elements");
}

Event Timeline

c4science · Help