diff --git a/cohesive_extrinsic.cc b/cohesive_extrinsic.cc index ad11a50..dd8b9f1 100644 --- a/cohesive_extrinsic.cc +++ b/cohesive_extrinsic.cc @@ -1,180 +1,179 @@ /** * @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 "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char *argv[]) { initialize("material.dat", argc, argv); const UInt spatial_dimension = 3; const UInt max_steps = 10000; Mesh mesh(spatial_dimension); auto whoami = Communicator::getStaticCommunicator().whoAmI(); if (whoami == 0) { mesh.read("cube.msh"); } mesh.distribute(); SolidMechanicsModelCohesive model(mesh); /// model initialization model.initFull(_analysis_method = _static, _is_extrinsic = true); - model.initNewSolver(_explicit_lumped_mass); - - Real time_step = model.getStableTimeStep() * 0.05; - model.setTimeStep(time_step); - std::cout << "Time step: " << time_step << std::endl; - - auto &inserter = model.getElementInserter(); - inserter.setLimit(_y, 0.30, 0.20); - model.updateAutomaticInsertion(); - const auto &position = mesh.getNodes(); auto &velocity = model.getVelocity(); auto &displacement = model.getDisplacement(); 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 Real c = -1e6; - Real s = -3e8; + Real s = -2e8; 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"); model.setBaseName("extrinsic"); model.addDumpFieldVector("displacement"); model.addDumpField("velocity"); model.addDumpField("acceleration"); 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", "velocity"); model.addDumpFieldToDumper("cohesive elements", "tractions"); auto lower = mesh.getLowerBounds(); auto upper = mesh.getUpperBounds(); - /// initial conditions - Real loading_rate = 0.5; - Real disp_update = loading_rate * time_step; - model.dump(); model.dump("cohesive elements"); model.solveStep("static"); + 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}; + Ep = model.getEnergy("potential"); + Ew += model.getEnergy("external work"); + + auto Et = Ed + Ep + Ek - Ew; + + fout << s << ", " << 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; + } + model.dump(); model.dump("cohesive elements"); s = -1e6; Matrix new_shear{{0., 0., s}, {0., 0., 0.}, {s, 0., 0.}}; - Real Ed, Ep, Ek, Ew{0}; - - Ew += model.getEnergy("external work"); - - std::ofstream fout; - fout.open("energies.csv", std::ofstream::out | std::ofstream::trunc); - fout << "step, ed, ep, ek, ew, et" << std::endl; - /// 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(); model.solveStep("explicit_lumped"); Ed = model.getEnergy("dissipated"); Ep = model.getEnergy("potential"); Ek = model.getEnergy("kinetic"); Ew += model.getEnergy("external work"); - auto Et = Ed + Ep + Ek - Ew; + Et = Ed + Ep + Ek - Ew; fout << s << ", " << Ed << ", " << Ep << ", " << Ek << ", " << Ew << ", " << Et << std::endl; if (s % 100 == 0) { model.dump(); model.dump("cohesive elements"); if (whoami == 0) { std::cout << "passing step " << s << "/" << max_steps << std::endl; } } } if (whoami == 0) { std::cout << "The dissipated energy is " << Ed << std::endl; } fout.close(); finalize(); return EXIT_SUCCESS; }