diff --git a/cohesive_extrinsic.cc b/cohesive_extrinsic.cc index f8a8fe7..a2a37e7 100644 --- a/cohesive_extrinsic.cc +++ b/cohesive_extrinsic.cc @@ -1,271 +1,275 @@ /** * @file cohesive_extrinsic.cc * * @author Zineb Fouad * @author Nicolas Richart * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Wed Feb 06 2019 * * @brief Cohesive element examples in extrinsic * * * @section LICENSE * * Copyright (©) 2015-2021 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 #include #include -#include /* -------------------------------------------------------------------------- */ using clk = std::chrono::high_resolution_clock; using std::chrono::duration_cast; -using std::chrono::microseconds ; +using std::chrono::microseconds; //#define AKANTU_VERSION_MAJOR 2 class Chrono { public: inline void start() { _start = clk::now(); }; - inline void store_time(const std::string & type) { + inline void store_time(const std::string &type) { clk::time_point _end = clk::now(); if (measures.find(type) == measures.end()) { measures[type] = duration_cast(_end - _start); nb_measures[type] = 1; } else { measures[type] += duration_cast(_end - _start); ++nb_measures[type]; } _start = clk::now(); } - virtual void printself(std::ostream & stream, int indent = 0) const { + virtual void printself(std::ostream &stream, int indent = 0) const { std::string space(indent, AKANTU_INDENT); stream << space << "Chrono [" << std::endl; std::map::const_iterator it; for (it = measures.begin(); it != measures.end(); ++it) { - const std::pair & p = *it; - const unsigned int & nb_measure = nb_measures.find(p.first)->second; + const std::pair &p = *it; + const unsigned int &nb_measure = nb_measures.find(p.first)->second; stream << space << " + " << p.first << "\t: " << std::setw(7) << std::fixed << std::setprecision(0) << p.second.count() / double(nb_measure) << "us" << std::endl; } stream << space << "]" << std::endl; } private: clk::time_point _start; std::map measures; std::map nb_measures; }; -inline std::ostream & operator<<(std::ostream & stream, const Chrono & _this) { +inline std::ostream &operator<<(std::ostream &stream, const Chrono &_this) { _this.printself(stream); return stream; } using namespace akantu; int main(int argc, char *argv[]) { initialize("material.dat", argc, argv); const UInt spatial_dimension = 3; Chrono chrono; chrono.start(); clk::time_point start_time = clk::now(); Mesh mesh(spatial_dimension); auto whoami = Communicator::getStaticCommunicator().whoAmI(); if (whoami == 0) { mesh.read("cube.msh"); } mesh.distribute(); chrono.store_time("initialize_mesh"); SolidMechanicsModelCohesive model(mesh); /// model initialization model.initFull(_analysis_method = _static, _is_extrinsic = true); chrono.store_time("init_full"); auto &blocked_dofs = model.getBlockedDOFs(); auto &force = model.getExternalForce(); /// boundary conditions model.applyBC(BC::Dirichlet::FixedValue(0.0, _z), "bottom"); // face model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "right line"); // line blocked_dofs(3, _y) = true; // point const auto &usersect = getUserParser(); const Real c = usersect.getParameter("compression"); const Real s = usersect.getParameter("shear"); const Real inc_s = usersect.getParameter("inc_shear"); const bool output_energy = usersect.getParameter("output_energy", true); const bool output_paraview = usersect.getParameter("output_energy", true); const UInt max_steps = usersect.getParameter("max_steps"); Matrix compression{{c, 0., 0.}, {0., c, 0.}, {0., 0., c}}; Matrix shear{{0., 0., s}, {0., 0., 0.}, {s, 0., 0.}}; force.zero(); model.applyBC(BC::Neumann::FromHigherDim(compression), "top"); model.applyBC(BC::Neumann::FromHigherDim(compression), "bottom"); model.applyBC(BC::Neumann::FromHigherDim(shear), "top"); model.applyBC(BC::Neumann::FromHigherDim(shear), "bottom"); model.applyBC(BC::Neumann::FromHigherDim(shear), "side left"); model.applyBC(BC::Neumann::FromHigherDim(shear), "side right"); - if(output_paraview) { + if (output_paraview) { model.setBaseName("extrinsic"); model.addDumpFieldVector("displacement"); model.addDumpField("internal_force"); model.addDumpField("external_force"); model.addDumpField("stress"); model.addDumpField("blocked_dofs"); model.addDumpField("grad_u"); model.addDumpFieldToDumper("cohesive elements", "displacement"); model.addDumpFieldToDumper("cohesive elements", "tractions"); model.dump(); model.dump("cohesive elements"); } chrono.store_time("initial_condition"); model.solveStep("static"); chrono.store_time("static_solve"); std::ofstream fout; fout.open("energies.csv", std::ofstream::out | std::ofstream::trunc); fout << "step, ed, ep, ek, ew, et" << std::endl; Real Ed{0}, Ep{0}, Ek{0}, Ew{0}; if (output_energy) { Ep = model.getEnergy("potential"); - //Ek = model.getEnergy("kinetic"); + // Ek = model.getEnergy("kinetic"); Ew += model.getEnergy("external work"); } auto Et = Ed + Ep + Ek - Ew; if (output_energy and whoami == 0) { fout << 0 << ", " << Ed << ", " << Ep << ", " << Ek << ", " << Ew << ", " << Et << std::endl; } model.initNewSolver(_explicit_lumped_mass); Real time_step = model.getStableTimeStep() * 0.05; model.setTimeStep(time_step); if (whoami == 0) { std::cout << "Time step: " << time_step << std::endl; } if (output_paraview) { model.addDumpField("velocity"); model.addDumpField("acceleration"); model.addDumpFieldToDumper("cohesive elements", "velocity"); model.dump(); model.dump("cohesive elements"); } Matrix new_shear{{0., 0., inc_s}, {0., 0., 0.}, {inc_s, 0., 0.}}; auto init_time = duration_cast(clk::now() - start_time); chrono.store_time("before_step"); start_time = clk::now(); /// Main loop for (auto s : arange(1, max_steps)) { if (s % 100 == 0 and s < 10000) { model.applyBC(BC::Neumann::FromHigherDim(new_shear), "top"); model.applyBC(BC::Neumann::FromHigherDim(new_shear), "bottom"); model.applyBC(BC::Neumann::FromHigherDim(new_shear), "side left"); model.applyBC(BC::Neumann::FromHigherDim(new_shear), "side right"); } model.checkCohesiveStress(); chrono.store_time("check_cohesive_stress"); model.solveStep("explicit_lumped"); chrono.store_time("solve_step"); - if (output_energy and whoami == 0) { + if (output_energy) { Ed = model.getEnergy("dissipated"); Ep = model.getEnergy("potential"); Ek = model.getEnergy("kinetic"); Ew += model.getEnergy("external work"); Et = Ed + Ep + Ek - Ew; - fout << s << ", " << Ed << ", " << Ep << ", " << Ek << ", " << Ew << ", " - << Et << std::endl; + if (whoami == 0) { + fout << s << ", " << Ed << ", " << Ep << ", " << Ek << ", " << Ew + << ", " << Et << std::endl; + } chrono.store_time("energies"); } if (output_paraview and (s % 100 == 0)) { model.dump(); model.dump("cohesive elements"); if (whoami == 0) { std::cout << "passing step " << s << "/" << max_steps << std::endl; } chrono.store_time("dumpers"); } } auto loop_time = duration_cast(clk::now() - start_time); if (whoami == 0) { std::cout << "The dissipated energy is " << Ed << std::endl; } if (whoami == 0) { std::cout << chrono << std::endl; - std::cout << "Nb proc: " << Communicator::getStaticCommunicator().getNbProc() << std::endl; + std::cout << "Nb proc: " + << Communicator::getStaticCommunicator().getNbProc() << std::endl; std::cout << "Full time: " << (init_time + loop_time).count() << std::endl; std::cout << "Init time: " << init_time.count() << std::endl; - std::cout << "Step time: " << loop_time.count() << " - nb_steps: " - << (max_steps -1) << " - time_per_step: " - << (loop_time.count()/(max_steps -1)) << std::endl; + std::cout << "Step time: " << loop_time.count() + << " - nb_steps: " << (max_steps - 1) + << " - time_per_step: " << (loop_time.count() / (max_steps - 1)) + << std::endl; } fout.close(); finalize(); return EXIT_SUCCESS; } diff --git a/cube.geo b/cube.geo index 4c495b2..9979109 100644 --- a/cube.geo +++ b/cube.geo @@ -1,64 +1,64 @@ -h = .3; +h = .1; Point(1) = { 1, 1, -1, h}; Point(2) = {-1, 1, -1, h}; Point(3) = {-1, -1, -1, h}; Point(4) = { 1, -1, -1, h}; Point(5) = {0, -1, -1, h}; Point(6) = {0, 1, -1, h}; Line(1) = {1, 6}; Line(2) = {6, 2}; Line(3) = {2, 3}; Line(4) = {3, 5}; Line(5) = {5, 4}; Line(6) = {4, 1}; Line(7) = {5, 6}; Line Loop(5) = {1, -7, 5, 6}; Line Loop(6) = {2, 3, 4, 7}; Plane Surface(5) = {-5}; Plane Surface(6) = {-6}; Point(11) = { 1, 1, 1, h}; Point(12) = {-1, 1, 1, h}; Point(13) = {-1, -1, 1, h}; Point(14) = { 1, -1, 1, h}; Line(11) = {11, 12}; Line(12) = {12, 13}; Line(13) = {13, 14}; Line(14) = {14, 11}; Line Loop(15) = {11, 12, 13, 14}; Plane Surface(16) = {15}; Line(15) = {2, 12}; Line(16) = {11, 1}; Line(17) = {3, 13}; Line(18) = {14, 4}; Curve Loop(16) = {17, -12, -15, 3}; Plane Surface(17) = {16}; Curve Loop(17) = {4, 5, -18, -13, -17}; Plane Surface(18) = {17}; Curve Loop(18) = {18, 6, -16, -14}; Plane Surface(19) = {18}; Curve Loop(19) = {1, 2, 15, -11, 16}; Plane Surface(20) = {19}; Surface Loop(1) = {5, 6, 20, 17, 18, 19, 16}; Volume(1) = {1}; Physical Volume("Body") = {1}; Physical Surface("top") = {16}; Physical Surface("bottom") = {5, 6}; Physical Curve("mid line") = {7}; Physical Curve("right line") = {6}; //+ Physical Surface("side left") = {17}; //+ Physical Surface("side right") = {19};