Homec4science
Diffusion diagonal_compression (master) Browse

R12270/diagonal_compression_static.cc
5fc0a4afc319master

Last Change
Raw File

diagonal_compression_static.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;
//#define explicit_
int main(int argc, char * argv[]) {
initialize("wall_static.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);
const auto & low = mesh.getBBox().getLowerBounds();
const auto & up = mesh.getBBox().getUpperBounds();
auto & top = mesh.createElementGroup("top", spatial_dimension);
auto & bottom = mesh.createElementGroup("bottom", spatial_dimension);
for (auto && data :
enumerate(make_view(mesh.getNodes(), spatial_dimension))) {
auto && node = std::get<0>(data);
auto && pos = std::get<1>(data);
if ((std::abs(pos(_z) - up(_z)) < 1e-10) and
(pos(_y) > up(_y) - clamp_sizes(_y))) {
top.addNode(node);
}
if ((std::abs(pos(_y) - up(_y)) < 1e-10) and
(pos(_z) > up(_z) - clamp_sizes(_z))) {
top.addNode(node);
}
if ((std::abs(pos(_z) - low(_z)) < 1e-10) and
(pos(_y) < low(_y) + clamp_sizes(_y))) {
bottom.addNode(node);
}
if ((std::abs(pos(_y) - low(_y)) < 1e-10) and
(pos(_z) < low(_z) + clamp_sizes(_z))) {
bottom.addNode(node);
}
}
SolidMechanicsModel solid(mesh);
auto material_selector =
std::make_shared<MeshDataMaterialSelector<std::string>>("physical_names",
solid);
solid.setMaterialSelector(material_selector);
#ifdef explicit_
solid.initFull(_analysis_method = _explicit_lumped_mass);
#else
solid.initFull(_analysis_method = _static);
#endif
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "top");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "top");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _z), "top");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "bottom");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "bottom");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _z), "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;
}
solid.setTimeStep(time_step);
solid.setBaseName("wall");
solid.addDumpField("displacement");
solid.addDumpField("external_force");
solid.addDumpField("internal_force");
solid.addDumpField("blocked_dofs");
solid.addDumpField("grad_u");
solid.addDumpField("material_index");
solid.addDumpField("stress");
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 increment_vect = displacement_increment * direction_down;
std::cout << increment_vect << std::endl;
solid.dump();
#ifdef explicit_
auto & velocity = solid.getVelocity();
#endif
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";
auto start_time = clk::now();
#ifdef explicit_
for (auto s : arange(max_steps)) {
solid.applyBC(BC::Dirichlet::IncrementValue(increment_vect(_x), _x),
"top");
solid.applyBC(BC::Dirichlet::IncrementValue(increment_vect(_y), _y),
"top");
solid.applyBC(BC::Dirichlet::IncrementValue(increment_vect(_z), _z),
"top");
solid.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();
}
}
#else
solid.solveStep();
auto loop_time = clk::now() - start_time;
std::cout << std::chrono::duration_cast<milliseconds>(loop_time).count() << std::endl;
#endif
solid.dump();
}
c4science · Help