diff --git a/cohesive_extrinsic.cc b/cohesive_extrinsic.cc index 2fd46b7..0c5e8b0 100644 --- a/cohesive_extrinsic.cc +++ b/cohesive_extrinsic.cc @@ -1,293 +1,296 @@ /** * @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 #include /* -------------------------------------------------------------------------- */ using clk = std::chrono::high_resolution_clock; using seconds = std::chrono::duration; using milliseconds = std::chrono::duration; //#define AKANTU_VERSION_MAJOR 2 class Chrono { public: inline void start() { _start = clk::now(); }; inline void store_time(const std::string &type) { clk::time_point _end = clk::now(); if (measures.find(type) == measures.end()) { measures[type] = _end - _start; nb_measures[type] = 1; } else { measures[type] += _end - _start; ++nb_measures[type]; } _start = clk::now(); } virtual void printself(std::ostream &stream, int indent = 0) const { std::string space(indent, AKANTU_INDENT); stream << space << "Chrono [" << std::endl; - for (auto it = measures.begin(); it != measures.end(); ++it) { - const auto &p = *it; - const unsigned int &nb_measure = nb_measures.find(p.first)->second; - stream << space << " + " << p.first << "\t: " << std::setw(25) + for (auto && measure : measures) { + const unsigned int &nb_measure = nb_measures.find(measure.first)->second; + stream << space << " + " << measure.first << "\t: " << std::setw(25) << std::fixed << std::setprecision(16) - << p.second.count() / double(nb_measure) + << measure.second.count() / double(nb_measure) << "us - nb_repetition: " << nb_measure << 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) { _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; + + 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_paraview", true); + const UInt max_steps = usersect.getParameter("max_steps"); + const std::string mesh_filename = usersect.getParameter("mesh"); + 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.read(mesh_filename); + chrono.store_time("read_mesh"); } mesh.distribute(); - chrono.store_time("init_mesh"); + chrono.store_time("dist_mesh"); - // SolidMechanicsModelCohesive model(mesh); + //SolidMechanicsModelCohesive model(mesh); SolidMechanicsModel model(mesh); /// model initialization - model.initFull(_analysis_method = _static, _is_extrinsic = true); + model.initFull(_analysis_method = _explicit_lumped_mass, + _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) { 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("init_cond"); - - // model.solveStep("static"); + chrono.store_time("initial_conditons"); - chrono.store_time("static_solve"); + //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); + //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.}}; seconds init_time = clk::now() - start_time; chrono.store_time("before_step"); auto nb_dofs_start = mesh.getNbGlobalNodes() * spatial_dimension; start_time = clk::now(); /// Main loop for (auto s : arange(1, max_steps + 1)) { 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"); + chrono.store_time("boundary_conditions"); } - // model.checkCohesiveStress(); - // chrono.store_time("check_cohesive_stress"); + //model.checkCohesiveStress(); + //chrono.store_time("check_cohesive_stress"); model.solveStep("explicit_lumped"); chrono.store_time("solve_step"); if (output_energy) { - // Ed = model.getEnergy("dissipated"); + //Ed = model.getEnergy("dissipated"); Ep = model.getEnergy("potential"); Ek = model.getEnergy("kinetic"); Ew += model.getEnergy("external work"); Et = Ed + Ep + Ek - Ew; 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) { milliseconds loop_time = clk::now() - start_time; std::cout << "passing step " << s << "/" << max_steps - << " - nb_cohesive_element: " - << mesh.getNbElement(spatial_dimension, _not_ghost, - _ek_cohesive) + // << " - nb_cohesive_element: " + // << mesh.getNbElement(spatial_dimension, _not_ghost, + // _ek_cohesive) << " - nb_dofs: " << mesh.getNbGlobalNodes() * spatial_dimension << " - " << loop_time.count() / s << "\t\t \r"; std::cout.flush(); } chrono.store_time("dumpers"); } } auto nb_dofs_end = mesh.getNbGlobalNodes() * spatial_dimension; seconds loop_time = clk::now() - start_time; + std::cout << std::endl; + if (whoami == 0) { - std::cout << "The dissipated energy is " << Ed << std::endl; - } - - if (whoami == 0) { + //std::cout << "The dissipated energy is " << Ed << std::endl; std::cout << chrono << std::endl; std::cout << "Nb proc: " << Communicator::getStaticCommunicator().getNbProc() - << " - nb_dofs:" << mesh.getNbGlobalNodes() * spatial_dimension + << " - nb_dofs: " << mesh.getNbGlobalNodes() * spatial_dimension << 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) << std::endl; - // std::cout << "Ns DOFs - start: " << nb_dofs_start - // << " - end: " << nb_dofs_end << std::endl; + std::cout << "Ns DOFs - start: " << nb_dofs_start + << " - end: " << nb_dofs_end << std::endl; } fout.close(); finalize(); return EXIT_SUCCESS; } diff --git a/cube.geo b/cube.geo index 9979109..55aa22b 100644 --- a/cube.geo +++ b/cube.geo @@ -1,64 +1,64 @@ -h = .1; +h = .002; 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}; diff --git a/material.dat b/material.dat index 2f7dd84..f70dd97 100644 --- a/material.dat +++ b/material.dat @@ -1,37 +1,38 @@ user section [ compression = -1e5 shear = -6e7 inc_shear = -1e6 - max_steps = 10000 + max_steps = 1000 output_energy = true output_paraview = true + mesh = cube.msh ] model solid_mechanics_model [ material elastic [ name = body rho = 2500 # density E = 70e9 # young's modulus nu = 0.22 # poisson's ratio finite_deformation = true ] # material cohesive_linear [ # name = insertion # beta = 3 # G_c = 7.6 # sigma_c = 70e6 # delta_c = 0.217e-6 # ] # cohesive_inserter [ # bounding_box = [[-1., 1.], [-1., 1.], [-1., 1.]] # ] - time_step_solver static [ - name = static - non_linear_solver linear [ - type = linear - ] - ] +# time_step_solver static [ +# name = static +# non_linear_solver linear [ +# type = linear +# ] +# ] ]