diff --git a/examples/cohesive_element/cohesive_intrinsic/cohesive_intrinsic.cc b/examples/cohesive_element/cohesive_intrinsic/cohesive_intrinsic.cc index 841c35d0f..ec276ee41 100644 --- a/examples/cohesive_element/cohesive_intrinsic/cohesive_intrinsic.cc +++ b/examples/cohesive_element/cohesive_intrinsic/cohesive_intrinsic.cc @@ -1,151 +1,136 @@ /** * @file cohesive_intrinsic.cc * * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Mon Jan 18 2016 * * @brief Test for cohesive elements * * @section LICENSE * * Copyright (©) 2015 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 "element_group.hh" #include "mesh_iterators.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; -static void updateDisplacement(SolidMechanicsModelCohesive &, Array &, - ElementType, Real); +static void updateDisplacement(SolidMechanicsModelCohesive &, + const ElementGroup &, Real); /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize("material.dat", argc, argv); const UInt spatial_dimension = 2; const UInt max_steps = 350; - const ElementType type = _triangle_6; - Mesh mesh(spatial_dimension); mesh.read("triangle.msh"); SolidMechanicsModelCohesive model(mesh); model.getElementInserter().setLimit(_x, -0.26, -0.24); /// model initialization - model.initFull(); + model.initFull(_analysis_method = _explicit_lumped_mass, + _is_extrinsic = false); Real time_step = model.getStableTimeStep() * 0.8; model.setTimeStep(time_step); std::cout << "Time step: " << time_step << std::endl; Array & boundary = model.getBlockedDOFs(); UInt nb_nodes = mesh.getNbNodes(); /// boundary conditions for (UInt dim = 0; dim < spatial_dimension; ++dim) { for (UInt n = 0; n < nb_nodes; ++n) { boundary(n, dim) = true; } } model.setBaseName("intrinsic"); model.addDumpFieldVector("displacement"); model.addDumpField("velocity"); model.addDumpField("acceleration"); model.addDumpField("stress"); model.addDumpField("grad_u"); model.addDumpField("external_force"); model.addDumpField("internal_force"); model.dump(); /// update displacement - Array elements; + auto && elements = mesh.createElementGroup("diplacement"); Vector barycenter(spatial_dimension); - for_each_element(mesh, [&](auto && el) { - mesh.getBarycenter(el, barycenter); - if (barycenter(_x) > -0.25) - elements.push_back(el.element); - }); + + for_each_element(mesh, + [&](auto && el) { + mesh.getBarycenter(el, barycenter); + if (barycenter(_x) > -0.25) + elements.add(el, true); + }, + _element_kind = _ek_regular); Real increment = 0.01; - updateDisplacement(model, elements, type, increment); + updateDisplacement(model, elements, increment); /// Main loop for (UInt s = 1; s <= max_steps; ++s) { model.solveStep(); - updateDisplacement(model, elements, type, increment); - + updateDisplacement(model, elements, increment); if (s % 1 == 0) { model.dump(); std::cout << "passing step " << s << "/" << max_steps << std::endl; } } Real Ed = model.getEnergy("dissipated"); - Real Edt = 2 * sqrt(2); std::cout << Ed << " " << Edt << std::endl; if (Ed < Edt * 0.999 || Ed > Edt * 1.001 || std::isnan(Ed)) { std::cout << "The dissipated energy is incorrect" << std::endl; return EXIT_FAILURE; } finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ static void updateDisplacement(SolidMechanicsModelCohesive & model, - Array & elements, ElementType type, - Real increment) { - Mesh & mesh = model.getMesh(); - UInt nb_element = elements.size(); - UInt nb_nodes = mesh.getNbNodes(); - UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); - - const Array & connectivity = mesh.getConnectivity(type); + const ElementGroup & group, Real increment) { Array & displacement = model.getDisplacement(); - Array update(nb_nodes); - update.clear(); - - for (UInt el = 0; el < nb_element; ++el) { - for (UInt n = 0; n < nb_nodes_per_element; ++n) { - UInt node = connectivity(elements(el), n); - if (!update(node)) { - displacement(node, 0) -= increment; - // displacement(node, 1) += increment; - update(node) = true; - } - } + + for (auto && node : group.getNodeGroup().getNodes()) { + displacement(node, 0) += increment; } } diff --git a/examples/cohesive_element/cohesive_intrinsic/triangle.geo b/examples/cohesive_element/cohesive_intrinsic/triangle.geo index cbe3dd7a7..9b4ba1410 100644 --- a/examples/cohesive_element/cohesive_intrinsic/triangle.geo +++ b/examples/cohesive_element/cohesive_intrinsic/triangle.geo @@ -1,15 +1,15 @@ -h = 1.; +h = .1; Point(1) = { 1, 1, 0, h}; Point(2) = {-1, 1, 0, h}; Point(3) = {-1,-1, 0, h}; Point(4) = { 1,-1, 0, h}; Line(1) = {1, 2}; Line(2) = {2, 3}; Line(3) = {3, 4}; Line(4) = {4, 1}; Line Loop(5) = {1, 2, 3, 4}; Plane Surface(6) = {5}; Physical Surface(7) = {6}; Transfinite Line {1, 2, 3, 4} = 3; \ No newline at end of file diff --git a/examples/explicit/explicit_dynamic.cc b/examples/explicit/explicit_dynamic.cc index dd0f29a9e..76b6d19e4 100644 --- a/examples/explicit/explicit_dynamic.cc +++ b/examples/explicit/explicit_dynamic.cc @@ -1,113 +1,107 @@ /** * @file explicit_dynamic.cc * * @author Seyedeh Mohadeseh Taheri Mousavi * * @date creation: Sun Jul 12 2015 * @date last modification: Mon Jan 18 2016 * * @brief This code refers to the explicit dynamic example from the user manual * * @section LICENSE * * Copyright (©) 2015 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.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { initialize("material.dat", argc, argv); const UInt spatial_dimension = 3; const Real pulse_width = 2.; const Real A = 0.01; Real time_step; Real time_factor = 0.8; UInt max_steps = 1000; Mesh mesh(spatial_dimension); if (Communicator::getStaticCommunicator().whoAmI() == 0) mesh.read("bar.msh"); - mesh.distribute(); - - mesh.makePeriodic(_x); - mesh.makePeriodic(_y); - mesh.makePeriodic(_z); - SolidMechanicsModel model(mesh); /// model initialization model.initFull(_analysis_method = _explicit_lumped_mass); time_step = model.getStableTimeStep(); std::cout << "Time Step = " << time_step * time_factor << "s (" << time_step << "s)" << std::endl; model.setTimeStep(time_step * time_factor); /// boundary and initial conditions Array & displacement = model.getDisplacement(); const Array & nodes = mesh.getNodes(); for (UInt n = 0; n < mesh.getNbNodes(); ++n) { Real x = nodes(n) - 2; // Sinus * Gaussian Real L = pulse_width; Real k = 0.1 * 2 * M_PI * 3 / L; displacement(n) = A * sin(k * x) * exp(-(k * x) * (k * x) / (L * L)); } std::ofstream energy; energy.open("energy.csv"); energy << "id,rtime,epot,ekin,tot" << std::endl; model.setBaseName("explicit_dynamic"); model.addDumpField("displacement"); model.addDumpField("velocity"); model.addDumpField("acceleration"); model.addDumpField("stress"); model.dump(); for (UInt s = 1; s <= max_steps; ++s) { model.solveStep(); Real epot = model.getEnergy("potential"); Real ekin = model.getEnergy("kinetic"); energy << s << "," << s * time_step << "," << epot << "," << ekin << "," << epot + ekin << "," << std::endl; if (s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; model.dump(); } energy.close(); finalize(); return EXIT_SUCCESS; } diff --git a/examples/io/dumper/CMakeLists.txt b/examples/io/dumper/CMakeLists.txt index 16f1bdeaa..b8c2a4e1c 100644 --- a/examples/io/dumper/CMakeLists.txt +++ b/examples/io/dumper/CMakeLists.txt @@ -1,61 +1,63 @@ #=============================================================================== # @file CMakeLists.txt # # @author Fabian Barras # # @date creation: Fri Sep 03 2010 # @date last modification: Wed Jan 20 2016 # # @brief CMakeLists for DumperIOHelper examples # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 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 . # # @section DESCRIPTION # #=============================================================================== add_mesh(swiss_train_mesh swiss_train.geo 2 1) add_library(locomotive_tools locomotive_tools.cc locomotive_tools.hh ) package_get_include_dir(BOOST _boost_include_dir) target_link_libraries(locomotive_tools PRIVATE akantu) target_include_directories(locomotive_tools PRIVATE ${AKANTU_INCLUDE_DIRS} ${_boost_include_dir}) if(AKANTU_EXTRA_CXX_FLAGS) set_target_properties(locomotive_tools PROPERTIES COMPILE_FLAGS ${AKANTU_EXTRA_CXX_FLAGS}) endif() register_example(dumper_low_level SOURCES dumper_low_level.cc USE_PACKAGES IOHelper DEPENDS swiss_train_mesh locomotive_tools + DIRECTORIES_TO_CREATE paraview ) register_example(dumpable_interface SOURCES dumpable_interface.cc USE_PACKAGES IOHelper DEPENDS swiss_train_mesh locomotive_tools + DIRECTORIES_TO_CREATE paraview ) diff --git a/examples/new_material/local_material_damage.cc b/examples/new_material/local_material_damage.cc index bf28fd3a7..17807fd3c 100644 --- a/examples/new_material/local_material_damage.cc +++ b/examples/new_material/local_material_damage.cc @@ -1,112 +1,109 @@ /** * @file local_material_damage.cc * * @author Guillaume Anciaux * @author Marion Estelle Chambart * @author Nicolas Richart * * @date creation: Mon Jan 18 2016 * * @brief Specialization of the material class for the damage material * * @section LICENSE * * Copyright (©) 2015 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 "local_material_damage.hh" #include "solid_mechanics_model.hh" namespace akantu { /* -------------------------------------------------------------------------- */ LocalMaterialDamage::LocalMaterialDamage(SolidMechanicsModel & model, const ID & id) : Material(model, id), damage("damage", *this) { AKANTU_DEBUG_IN(); this->registerParam("E", E, 0., _pat_parsable, "Young's modulus"); this->registerParam("nu", nu, 0.5, _pat_parsable, "Poisson's ratio"); this->registerParam("lambda", lambda, _pat_readable, "First Lamé coefficient"); this->registerParam("mu", mu, _pat_readable, "Second Lamé coefficient"); this->registerParam("kapa", kpa, _pat_readable, "Bulk coefficient"); this->registerParam("Yd", Yd, 50., _pat_parsmod); this->registerParam("Sd", Sd, 5000., _pat_parsmod); damage.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void LocalMaterialDamage::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); lambda = nu * E / ((1 + nu) * (1 - 2 * nu)); mu = E / (2 * (1 + nu)); kpa = lambda + 2. / 3. * mu; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void LocalMaterialDamage::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = damage(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); computeStressOnQuad(grad_u, sigma, *dam); ++dam; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void LocalMaterialDamage::computePotentialEnergy(ElementType el_type, - GhostType ghost_type) { +void LocalMaterialDamage::computePotentialEnergy(ElementType el_type) { AKANTU_DEBUG_IN(); - if (ghost_type != _not_ghost) - return; Real * epot = potential_energy(el_type).storage(); - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, _not_ghost); computePotentialEnergyOnQuad(grad_u, sigma, *epot); epot++; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } static bool material_is_alocated_local_damage [[gnu::unused]] = MaterialFactory::getInstance().registerAllocator( "local_damage", [](UInt, const ID &, SolidMechanicsModel & model, const ID & id) -> std::unique_ptr { return std::make_unique(model, id); }); } // namespace akantu diff --git a/examples/new_material/local_material_damage.hh b/examples/new_material/local_material_damage.hh index 1534e3e5e..f797529a3 100644 --- a/examples/new_material/local_material_damage.hh +++ b/examples/new_material/local_material_damage.hh @@ -1,125 +1,118 @@ /** * @file local_material_damage.hh * * @author Guillaume Anciaux * @author Marion Estelle Chambart * @author Nicolas Richart * * @date creation: Mon Aug 10 2015 * @date last modification: Mon Jan 18 2016 * * @brief Material isotropic elastic * * @section LICENSE * * Copyright (©) 2015 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 "aka_common.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_LOCAL_MATERIAL_DAMAGE_HH__ #define __AKANTU_LOCAL_MATERIAL_DAMAGE_HH__ namespace akantu { class LocalMaterialDamage : public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: LocalMaterialDamage(SolidMechanicsModel & model, const ID & id = ""); virtual ~LocalMaterialDamage(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - void initMaterial(); + void initMaterial() override; /// constitutive law for all element of a type - void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); + void computeStress(ElementType el_type, + GhostType ghost_type = _not_ghost) override; + /// compute the potential energy for all elements + void computePotentialEnergy(ElementType el_type) override; + +protected: /// constitutive law for a given quadrature point inline void computeStressOnQuad(Matrix & grad_u, Matrix & sigma, Real & damage); - /// compute tangent stiffness - virtual void computeTangentStiffness(__attribute__((unused)) - const ElementType & el_type, - __attribute__((unused)) - Array & tangent_matrix, - __attribute__((unused)) - GhostType ghost_type = _not_ghost){}; - - /// compute the potential energy for all elements - void computePotentialEnergy(ElementType el_type, - GhostType ghost_type = _not_ghost); - /// compute the potential energy for on element inline void computePotentialEnergyOnQuad(Matrix & grad_u, Matrix & sigma, Real & epot); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// compute the celerity of the fastest wave in the material - inline Real getCelerity(const Element & element) const; + inline Real getCelerity(const Element & element) const override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ - +public: AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); private: /// the young modulus Real E; /// Poisson coefficient Real nu; /// First Lamé coefficient Real lambda; /// Second Lamé coefficient (shear modulus) Real mu; /// resistance to damage Real Yd; /// damage threshold Real Sd; /// Bulk modulus Real kpa; /// damage internal variable InternalField damage; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "local_material_damage_inline_impl.hh" #endif /* __AKANTU_LOCAL_MATERIAL_DAMAGE_HH__ */ diff --git a/examples/python/eigen_modes/eigen_modes.py b/examples/python/eigen_modes/eigen_modes.py index f71492dbe..d306f2cec 100644 --- a/examples/python/eigen_modes/eigen_modes.py +++ b/examples/python/eigen_modes/eigen_modes.py @@ -1,248 +1,248 @@ import subprocess import argparse import akantu as aka import numpy as np from image_saver import ImageSaver import matplotlib.pyplot as plt from scipy.sparse.linalg import eigsh from scipy.sparse import csr_matrix # ----------------------------------------------------------------------------- # parser # ----------------------------------------------------------------------------- parser = argparse.ArgumentParser(description='Eigen mode exo') parser.add_argument('-m', '--mode_number', type=int, required=True, help='precise the mode to study') parser.add_argument('-wL', '--wave_width', type=float, help='precise the width of the wave for ' 'the initial displacement') parser.add_argument('-L', '--Lbar', type=float, help='precise the length of the bar', default=10) parser.add_argument('-t', '--time_step', type=float, help='precise the timestep', default=None) parser.add_argument('-n', '--max_steps', type=int, help='precise the number of timesteps', default=500) parser.add_argument('-mh', '--mesh_h', type=float, help='characteristic mesh size', default=.2) args = parser.parse_args() mode = args.mode_number wave_width = args.wave_width time_step = args.time_step max_steps = args.max_steps mesh_h = args.mesh_h Lbar = args.Lbar # ----------------------------------------------------------------------------- # Mesh Generation # ----------------------------------------------------------------------------- geo_content = """ // Mesh size h = {0}; """.format(mesh_h) geo_content += """ h1 = h; h2 = h; // Dimensions of the bar Lx = 10; Ly = 1; // ------------------------------------------ // Geometry // ------------------------------------------ Point(101) = { 0.0, -Ly/2, 0.0, h1}; Point(102) = { Lx, -Ly/2, 0.0, h2}; Point(103) = { Lx, 0., 0.0, h2}; Point(104) = { Lx, Ly/2., 0.0, h2}; Point(105) = { 0.0, Ly/2., 0.0, h1}; Point(106) = { 0.0, 0., 0.0, h1}; Line(101) = {101,102}; Line(102) = {102,103}; Line(103) = {103,104}; Line(104) = {104,105}; Line(105) = {105,106}; Line(106) = {106,101}; Line(107) = {106,103}; Line Loop(108) = {101, 102, -107, 106}; Plane Surface(109) = {108}; Line Loop(110) = {103, 104, 105, 107}; Plane Surface(111) = {110}; Physical Surface(112) = {109, 111}; Transfinite Surface "*"; Recombine Surface "*"; Physical Surface(113) = {111, 109}; Physical Line("XBlocked") = {103, 102}; Physical Line("ImposedVelocity") = {105, 106}; Physical Line("YBlocked") = {104, 101}; """ mesh_file = 'bar' with open(mesh_file + '.geo', 'w') as f: f.write(geo_content) subprocess.call(['gmsh', '-2', mesh_file + '.geo']) mesh_file = mesh_file + '.msh' # ----------------------------------------------------------------------------- # Initialization # ----------------------------------------------------------------------------- spatial_dimension = 2 aka.parseInput('material.dat') mesh = aka.Mesh(spatial_dimension) mesh.read(mesh_file) model = aka.SolidMechanicsModel(mesh) model.initFull(aka._implicit_dynamic) model.setBaseName("waves-{0}".format(mode)) model.addDumpFieldVector("displacement") model.addDumpFieldVector("acceleration") model.addDumpFieldVector("velocity") model.addDumpField("blocked_dofs") # ----------------------------------------------------------------------------- # Boundary conditions # ----------------------------------------------------------------------------- internal_force = model.getInternalForce() displacement = model.getDisplacement() acceleration = model.getAcceleration() velocity = model.getVelocity() blocked_dofs = model.getBlockedDOFs() nbNodes = mesh.getNbNodes() position = mesh.getNodes() model.applyBC(aka.FixedValue(0.0, aka._x), "XBlocked") model.applyBC(aka.FixedValue(0.0, aka._y), "YBlocked") # ------------------------------------------------------------------------ # timestep value computation # ------------------------------------------------------------------------ time_factor = 0.8 stable_time_step = model.getStableTimeStep() * time_factor if time_step: print("Required Time Step = {0}".format(time_step)) if stable_time_step * time_factor < time_step: print("Stable Time Step = {0}".format(stable_time_step)) raise RuntimeError("required time_step too large") print("Required Time Step = {0}".format(time_step)) else: print("Stable Time Step = {0}".format(stable_time_step)) time_step = stable_time_step * time_factor model.setTimeStep(time_step) disp_sav = ImageSaver(mesh, displacement, 0, Lbar) velo_sav = ImageSaver(mesh, velocity, 0, Lbar) # ------------------------------------------------------------------------ # compute the eigen modes # ------------------------------------------------------------------------ model.assembleStiffnessMatrix() model.assembleMass() stiff = model.getDOFManager().getMatrix('K') stiff = aka.AkantuSparseMatrix(stiff).toarray() mass = model.getDOFManager().getMatrix('M') mass = aka.AkantuSparseMatrix(mass).toarray() # select the non blocked DOFs by index in the mask -mask = np.equal(blocked_dofs.flatten(), False) +mask = np.equal(blocked_dofs.flatten(), False)m Mass_star = mass[mask, :] Mass_star = csr_matrix(Mass_star[:, mask].copy()) K_star = stiff[mask, :] K_star = csr_matrix(K_star[:, mask].copy()) print('getting the eigen values') vals, vects = eigsh(K_star, M=Mass_star, which='SM', k=20) # ----------------------------------------------------------------------------- # import the initial conditions in displacement # ----------------------------------------------------------------------------- displacement.reshape(nbNodes*2)[mask] = vects[:, mode] with open('modes.txt', 'a') as f: f.write('{0} {1}\n'.format(mode, vals[mode])) model.dump() # ----------------------------------------------------------------------------- # prepare the storage of the dynamical evolution # ----------------------------------------------------------------------------- e_p = np.zeros(max_steps + 1) e_k = np.zeros(max_steps + 1) e_t = np.zeros(max_steps + 1) time = np.zeros(max_steps + 1) norm = np.zeros(max_steps + 1) epot = model.getEnergy('potential') ekin = model.getEnergy('kinetic') e_p[0] = epot e_k[0] = ekin e_t[0] = epot + ekin time[0] = 0 # ----------------------------------------------------------------------------- # loop for evolution of motion dynamics # ----------------------------------------------------------------------------- for step in range(1, max_steps + 1): model.solveStep() # outputs epot = model.getEnergy('potential') ekin = model.getEnergy('kinetic') print(step, '/', max_steps, epot, ekin, epot + ekin) e_p[step] = epot e_k[step] = ekin e_t[step] = epot + ekin time[step] = (step + 1) * time_step disp_sav.storeStep() velo_sav.storeStep() if step % 10 == 0: model.dump() # ----------------------------------------------------------------------------- # plot figures for global evolution # ----------------------------------------------------------------------------- # energy norms plt.figure(1) plt.plot(time, e_t, 'r', time, e_p, 'b', time, e_k, 'g') # space-time diagram for diplacements plt.figure(2) plt.imshow(disp_sav.getImage(), extent=(0, Lbar, max_steps * time_step, 0)) plt.xlabel("Space ") plt.ylabel("Time ") # space-time diagram for velocities plt.figure(3) plt.imshow(velo_sav.getImage(), extent=(0, Lbar, max_steps * time_step, 0)) plt.xlabel("Velocity") plt.ylabel("Time") plt.show() diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index dc253b310..afb50ac4b 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -1,80 +1,89 @@ #=============================================================================== # @file CMakeLists.txt # # @author Nicolas Richart # # @date creation: Fri Dec 12 2014 # @date last modification: Mon Jan 18 2016 # # @brief CMake file for the python wrapping of akantu # # @section LICENSE # # Copyright (©) 2015 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 . # #=============================================================================== set(PYAKANTU_SRCS py_aka_common.cc py_aka_error.cc py_akantu.cc py_boundary_conditions.cc py_fe_engine.cc py_group_manager.cc py_mesh.cc py_model.cc py_parser.cc ) package_is_activated(iohelper _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_dumpable.cc ) endif() package_is_activated(solid_mechanics _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_solid_mechanics_model.cc py_material.cc ) endif() package_is_activated(cohesive_element _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_solid_mechanics_model_cohesive.cc ) endif() package_is_activated(heat_transfer _is_activated) if (_is_activated) list(APPEND PYAKANTU_SRCS py_heat_transfer_model.cc ) endif() -add_library(pyakantu OBJECT ${PYAKANTU_SRCS}) +if(CMAKE_VERSION VERSION_LESS 3.12) + add_library(pyakantu STATIC ${PYAKANTU_SRCS}) +else() + add_library(pyakantu OBJECT ${PYAKANTU_SRCS}) +endif() + target_link_libraries(pyakantu PUBLIC akantu pybind11) target_include_directories(pyakantu INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}) set_target_properties(pyakantu PROPERTIES POSITION_INDEPENDENT_CODE TRUE) -pybind11_add_module(py11_akantu $) +if(CMAKE_VERSION VERSION_LESS 3.12) + pybind11_add_module(py11_akantu py11_akantu.cc) +else() + pybind11_add_module(py11_akantu $) +endif() target_link_libraries(py11_akantu PRIVATE pyakantu) set_target_properties(py11_akantu PROPERTIES DEBUG_POSTFIX "")