diff --git a/examples/new_material/CMakeLists.txt b/examples/new_material/CMakeLists.txt index ed959fe8f..98a76ef50 100644 --- a/examples/new_material/CMakeLists.txt +++ b/examples/new_material/CMakeLists.txt @@ -1,47 +1,49 @@ #=============================================================================== # @file CMakeLists.txt # # @author Seyedeh Mohadeseh Taheri Mousavi # # @date creation: Mon Jan 18 2016 # @date last modification: Tue Jan 19 2016 # # @brief CMakeFile for new material example # # @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(new_local_material_barre_trou_mesh barre_trou.geo 2 2) register_example(new_local_material SOURCES local_material_damage.cc local_material_damage.hh local_material_damage_inline_impl.hh new_local_material.cc DEPENDS new_local_material_barre_trou_mesh FILES_TO_COPY material.dat ) + +add_example(viscoelastic_maxwell "Example on how to instantiate and obtaine energies from the viscoelastic maxwell material" PACKAGE core) diff --git a/examples/new_material/CMakeLists.txt b/examples/new_material/viscoelastic_maxwell/CMakeLists.txt similarity index 69% copy from examples/new_material/CMakeLists.txt copy to examples/new_material/viscoelastic_maxwell/CMakeLists.txt index ed959fe8f..5ac1bd788 100644 --- a/examples/new_material/CMakeLists.txt +++ b/examples/new_material/viscoelastic_maxwell/CMakeLists.txt @@ -1,47 +1,46 @@ #=============================================================================== # @file CMakeLists.txt # -# @author Seyedeh Mohadeseh Taheri Mousavi +# @author Emil Gallyamov # -# @date creation: Mon Jan 18 2016 -# @date last modification: Tue Jan 19 2016 +# @date creation: Tue Nov 20 2018 # -# @brief CMakeFile for new material example +# @brief CMakeFile for viscoelastic material example # # @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(new_local_material_barre_trou_mesh barre_trou.geo 2 2) -register_example(new_local_material - SOURCES - local_material_damage.cc - local_material_damage.hh - local_material_damage_inline_impl.hh - new_local_material.cc - DEPENDS - new_local_material_barre_trou_mesh - FILES_TO_COPY - material.dat + + + +add_mesh(material_viscoelastic_maxwell_mesh material_viscoelastic_maxwell_mesh.geo 2 1) +register_example(material_viscoelastic_maxwell_energies + SOURCES material_viscoelastic_maxwell_energies.cc + DEPENDS material_viscoelastic_maxwell_mesh + USE_PACKAGES iohelper + FILES_TO_COPY material_viscoelastic_maxwell.dat + DIRECTORIES_TO_CREATE paraview ) + diff --git a/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell.dat b/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell.dat new file mode 100644 index 000000000..1ccd6ae26 --- /dev/null +++ b/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell.dat @@ -0,0 +1,9 @@ +material viscoelastic_maxwell [ + name = anything + rho = 7800 # density + Einf = 15e6 # young's modulus + nu = 0.0 # poisson's ratio + Eta = [40e6] # viscosity + Ev =[10e6] # viscous stiffness + Plane_Stress = 0 +] diff --git a/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_energies.cc b/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_energies.cc new file mode 100644 index 000000000..85849fe62 --- /dev/null +++ b/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_energies.cc @@ -0,0 +1,181 @@ +/** + * @file material_viscoelastic_maxwell_energies.cc + * + * @author Emil Gallyamov + * + * @date creation: Tue Nov 20 2018 + * @date last modification: + * + * @brief Example of using viscoelastic material and computing energies + * + * @section LICENSE + * + * Copyright (©) 2010-2018 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 "material_viscoelastic_maxwell.hh" +#include "non_linear_solver.hh" +#include "solid_mechanics_model.hh" +#include "sparse_matrix.hh" + +using namespace akantu; + +/* -------------------------------------------------------------------------- */ +/* Main */ +/* -------------------------------------------------------------------------- */ + +int main(int argc, char *argv[]) { + akantu::initialize("material_viscoelastic_maxwell.dat", argc, argv); + + // sim data + Real eps = 0.1; + + const UInt dim = 2; + Real sim_time = 100.; + Real T = 10.; + Mesh mesh(dim); + mesh.read("material_viscoelastic_maxwell_mesh.msh"); + + SolidMechanicsModel model(mesh); + + /* ------------------------------------------------------------------------ */ + /* Initialization */ + /* ------------------------------------------------------------------------ */ + model.initFull(_analysis_method = _static); + std::cout << model.getMaterial(0) << std::endl; + + std::stringstream filename_sstr; + filename_sstr << "material_viscoelastic_maxwell_output.out"; + std::ofstream output_data; + output_data.open(filename_sstr.str().c_str()); + + Material &mat = model.getMaterial(0); + + Real time_step = 0.1; + + UInt nb_nodes = mesh.getNbNodes(); + const Array &coordinate = mesh.getNodes(); + Array &displacement = model.getDisplacement(); + Array &blocked = model.getBlockedDOFs(); + + /// Setting time step + + model.setTimeStep(time_step); + + model.setBaseName("dynamic"); + model.addDumpFieldVector("displacement"); + model.addDumpField("blocked_dofs"); + model.addDumpField("external_force"); + model.addDumpField("internal_force"); + model.addDumpField("grad_u"); + model.addDumpField("stress"); + model.addDumpField("strain"); + + + UInt max_steps = sim_time / time_step + 1; + Real time = 0.; + + auto &solver = model.getNonLinearSolver(); + solver.set("max_iterations", 10); + solver.set("threshold", 1e-7); + solver.set("convergence_type", _scc_residual); + + /* ------------------------------------------------------------------------ */ + /* Main loop */ + /* ------------------------------------------------------------------------ */ + for (UInt s = 0; s <= max_steps; ++s) { + + std::cout << "Time Step = " << time_step << "s" << std::endl; + std::cout << "Time = " << time << std::endl; + + // impose displacement + Real epsilon = 0; + if (time < T) { + epsilon = eps * time / T; + } else { + epsilon = eps; + } + + for (UInt n = 0; n < nb_nodes; ++n) { + if (Math::are_float_equal(coordinate(n, 0), 0.0)) { + displacement(n, 0) = 0; + blocked(n, 0) = true; + displacement(n, 1) = epsilon * coordinate(n, 1); + blocked(n, 1) = true; + } else if (Math::are_float_equal(coordinate(n, 1), 0.0)) { + displacement(n, 0) = epsilon * coordinate(n, 0); + blocked(n, 0) = true; + displacement(n, 1) = 0; + blocked(n, 1) = true; + } else if (Math::are_float_equal(coordinate(n, 0), 0.001)) { + displacement(n, 0) = epsilon * coordinate(n, 0); + blocked(n, 0) = true; + displacement(n, 1) = epsilon * coordinate(n, 1); + blocked(n, 1) = true; + } else if (Math::are_float_equal(coordinate(n, 1), 0.001)) { + displacement(n, 0) = epsilon * coordinate(n, 0); + blocked(n, 0) = true; + displacement(n, 1) = epsilon * coordinate(n, 1); + blocked(n, 1) = true; + } + } + + try { + model.solveStep(); + } catch (debug::Exception &e) { + } + + // for debugging + // auto int_force = model.getInternalForce(); + // auto &K = model.getDOFManager().getMatrix("K"); + // K.saveMatrix("K.mtx"); + + Int nb_iter = solver.get("nb_iterations"); + Real error = solver.get("error"); + bool converged = solver.get("converged"); + + if (converged) { + std::cout << "Converged in " << nb_iter << " iterations" << std::endl; + } else { + std::cout << "Didn't converge after " << nb_iter + << " iterations. Error is " << error << std::endl; + return EXIT_FAILURE; + } + + model.dump(); + + Real epot = mat.getEnergy("potential"); + Real edis = mat.getEnergy("dissipated"); + Real work = mat.getEnergy("work"); + + // data output + output_data << s * time_step << " " << epsilon << + << " " << epot << " " << edis << " " << + work << std::endl; + time += time_step; + + } + output_data.close(); + finalize(); +} diff --git a/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_mesh.geo b/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_mesh.geo new file mode 100644 index 000000000..173c1f300 --- /dev/null +++ b/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_mesh.geo @@ -0,0 +1,33 @@ +// Mesh size +h = 1; // Top cube + +// Dimensions of top cube +Lx = 0.001; +Ly = 0.001; + + +// ------------------------------------------ +// Geometry +// ------------------------------------------ + +// Base Cube +Point(101) = { 0.0, 0.0, 0.0, h}; // Bottom Face +Point(102) = { Lx, 0.0, 0.0, h}; // Bottom Face +Point(103) = { Lx, Ly, 0.0, h}; // Bottom Face +Point(104) = { 0.0, Ly, 0.0, h}; // Bottom Face + +// Base Cube +Line(101) = {101,102}; // Bottom Face +Line(102) = {102,103}; // Bottom Face +Line(103) = {103,104}; // Bottom Face +Line(104) = {104,101}; // Bottom Face + +// Base Cube +Line Loop(101) = {101:104}; + +Plane Surface(101) = {101}; + +Physical Surface(1) = {101}; + +Transfinite Surface "*"; +Recombine Surface "*"; diff --git a/test/test_model/test_solid_mechanics_model/test_materials/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_materials/CMakeLists.txt index eaee42f83..7fcecdb2d 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/CMakeLists.txt +++ b/test/test_model/test_solid_mechanics_model/test_materials/CMakeLists.txt @@ -1,87 +1,88 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # # @date creation: Fri Oct 22 2010 # @date last modification: Mon Jan 29 2018 # # @brief configuration for materials tests # # @section LICENSE # # Copyright (©) 2010-2018 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(test_local_material_barre_trou_mesh barre_trou.geo 2 2) add_mesh(test_local_material_barre_trou_mesh mesh_section_gap.geo 2 2) register_test(test_local_material SOURCES test_local_material.cc local_material_damage.cc EXTRA_FILES local_material_damage.hh local_material_damage_inline_impl.cc DEPENDS test_local_material_barre_trou_mesh FILES_TO_COPY material.dat DIRECTORIES_TO_CREATE paraview PACKAGE core ) # ============================================================================== add_mesh(test_interpolate_stress_mesh interpolation.geo 3 2) register_test(test_interpolate_stress test_interpolate_stress.cc FILES_TO_COPY material_interpolate.dat DEPENDS test_interpolate_stress_mesh DIRECTORIES_TO_CREATE paraview PACKAGE lapack core ) #=============================================================================== add_mesh(test_material_orthotropic_square_mesh square.geo 2 1) register_test(test_material_orthotropic SOURCES test_material_orthotropic.cc DEPENDS test_material_orthotropic_square_mesh FILES_TO_COPY orthotropic.dat DIRECTORIES_TO_CREATE paraview PACKAGE core lapack ) #=============================================================================== register_test(test_material_mazars SOURCES test_material_mazars.cc FILES_TO_COPY material_mazars.dat DIRECTORIES_TO_CREATE paraview PACKAGE core lapack ) # ============================================================================== add_akantu_test(test_material_viscoelastic "test the visco elastic materials") add_akantu_test(test_material_non_local "test the non-local materials") add_akantu_test(test_material_elasto_plastic_linear_isotropic_hardening "test the elasto plastic with linear isotropic hardening materials") +add_akantu_test(test_material_viscoelastic_maxwell "test the viscoelastic maxwell material") # ============================================================================== add_mesh(test_multi_material_elastic_mesh test_multi_material_elastic.geo 2 1) register_test(test_multi_material_elastic SOURCES test_multi_material_elastic.cc FILES_TO_COPY test_multi_material_elastic.dat DEPENDS test_multi_material_elastic_mesh PACKAGE implicit) # ============================================================================== # Material unit tests # ============================================================================== register_gtest_sources(SOURCES test_elastic_materials.cc PACKAGE core) register_gtest_sources(SOURCES test_finite_def_materials.cc PACKAGE core) register_gtest_sources(SOURCES test_damage_materials.cc PACKAGE core) register_gtest_sources(SOURCES test_plastic_materials.cc PACKAGE core) register_gtest_sources(SOURCES test_material_thermal.cc PACKAGE core) register_gtest_test(test_material) diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic_maxwell/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic_maxwell/CMakeLists.txt new file mode 100644 index 000000000..d90b92a4f --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic_maxwell/CMakeLists.txt @@ -0,0 +1,31 @@ +#=============================================================================== +# @file CMakeLists.txt +# +# @author Nicolas Richart +# +# @date creation: Thu Aug 09 2012 +# @date last modification: Wed Feb 03 2016 +# +# @brief configuration for viscoelastic materials tests +# +# @section LICENSE +# +# Copyright (©) 2010-2018 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 . +# +#=============================================================================== + +add_mesh(test_material_viscoelastic_maxwell_mesh + test_material_viscoelastic_maxwell.geo 2 1) + +register_test(test_material_viscoelasti_maxwell_relaxation + SOURCES test_material_viscoelastic_maxwell_relaxation.cc + DEPENDS test_material_viscoelastic_maxwell_mesh + FILES_TO_COPY material_viscoelastic_maxwell.dat + PACKAGE core + ) diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic_maxwell/material_viscoelastic_maxwell.dat b/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic_maxwell/material_viscoelastic_maxwell.dat new file mode 100644 index 000000000..1ccd6ae26 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic_maxwell/material_viscoelastic_maxwell.dat @@ -0,0 +1,9 @@ +material viscoelastic_maxwell [ + name = anything + rho = 7800 # density + Einf = 15e6 # young's modulus + nu = 0.0 # poisson's ratio + Eta = [40e6] # viscosity + Ev =[10e6] # viscous stiffness + Plane_Stress = 0 +] diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic_maxwell/test_material_viscoelastic_maxwell.geo b/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic_maxwell/test_material_viscoelastic_maxwell.geo new file mode 100644 index 000000000..173c1f300 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic_maxwell/test_material_viscoelastic_maxwell.geo @@ -0,0 +1,33 @@ +// Mesh size +h = 1; // Top cube + +// Dimensions of top cube +Lx = 0.001; +Ly = 0.001; + + +// ------------------------------------------ +// Geometry +// ------------------------------------------ + +// Base Cube +Point(101) = { 0.0, 0.0, 0.0, h}; // Bottom Face +Point(102) = { Lx, 0.0, 0.0, h}; // Bottom Face +Point(103) = { Lx, Ly, 0.0, h}; // Bottom Face +Point(104) = { 0.0, Ly, 0.0, h}; // Bottom Face + +// Base Cube +Line(101) = {101,102}; // Bottom Face +Line(102) = {102,103}; // Bottom Face +Line(103) = {103,104}; // Bottom Face +Line(104) = {104,101}; // Bottom Face + +// Base Cube +Line Loop(101) = {101:104}; + +Plane Surface(101) = {101}; + +Physical Surface(1) = {101}; + +Transfinite Surface "*"; +Recombine Surface "*"; diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic_maxwell/test_material_viscoelastic_maxwell_relaxation.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic_maxwell/test_material_viscoelastic_maxwell_relaxation.cc new file mode 100644 index 000000000..072314cc6 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic_maxwell/test_material_viscoelastic_maxwell_relaxation.cc @@ -0,0 +1,196 @@ +/** + * @file test_material_viscoelastic_maxwell_relaxation.cc + * + * @author Emil Gallyamov + * + * @date creation: Thu May 17 2018 + * @date last modification: + * + * @brief test of the viscoelastic material: relaxation + * + * @section LICENSE + * + * Copyright (©) 2010-2018 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 "material_viscoelastic_maxwell.hh" +#include "non_linear_solver.hh" +#include "solid_mechanics_model.hh" +#include "sparse_matrix.hh" + +using namespace akantu; + +/* -------------------------------------------------------------------------- */ +/* Main */ +/* -------------------------------------------------------------------------- */ + +int main(int argc, char *argv[]) { + akantu::initialize("material_viscoelastic_maxwell.dat", argc, argv); + + // sim data + Real eps = 0.1; + + const UInt dim = 2; + Real sim_time = 100.; + Real T = 10.; + Real tolerance = 1e-6; + Mesh mesh(dim); + mesh.read("test_material_viscoelastic_maxwell.msh"); + + const ElementType element_type = _quadrangle_4; + SolidMechanicsModel model(mesh); + + /* ------------------------------------------------------------------------ */ + /* Initialization */ + /* ------------------------------------------------------------------------ */ + model.initFull(_analysis_method = _static); + std::cout << model.getMaterial(0) << std::endl; + + std::stringstream filename_sstr; + filename_sstr << "test_material_viscoelastic_maxwell.out"; + std::ofstream output_data; + output_data.open(filename_sstr.str().c_str()); + + Material &mat = model.getMaterial(0); + + const Array &stress = mat.getStress(element_type); + + Vector Eta = mat.get("Eta"); + Vector Ev = mat.get("Ev"); + Real Einf = mat.get("Einf"); + Real nu = mat.get("nu"); + Real lambda = Eta(0) / Ev(0); + Real pre_mult = 1 / (1 + nu) / (1 - 2 * nu); + Real time_step = 0.1; + + UInt nb_nodes = mesh.getNbNodes(); + const Array &coordinate = mesh.getNodes(); + Array &displacement = model.getDisplacement(); + Array &blocked = model.getBlockedDOFs(); + + /// Setting time step + + model.setTimeStep(time_step); + + UInt max_steps = sim_time / time_step + 1; + Real time = 0.; + + auto &solver = model.getNonLinearSolver(); + solver.set("max_iterations", 200); + solver.set("threshold", 1e-7); + solver.set("convergence_type", _scc_residual); + + /* ------------------------------------------------------------------------ */ + /* Main loop */ + /* ------------------------------------------------------------------------ */ + for (UInt s = 0; s <= max_steps; ++s) { + + std::cout << "Time Step = " << time_step << "s" << std::endl; + std::cout << "Time = " << time << std::endl; + + // impose displacement + Real epsilon = 0; + if (time < T) { + epsilon = eps * time / T; + } else { + epsilon = eps; + } + + for (UInt n = 0; n < nb_nodes; ++n) { + if (Math::are_float_equal(coordinate(n, 0), 0.0)) { + displacement(n, 0) = 0; + blocked(n, 0) = true; + displacement(n, 1) = epsilon * coordinate(n, 1); + blocked(n, 1) = true; + } else if (Math::are_float_equal(coordinate(n, 1), 0.0)) { + displacement(n, 0) = epsilon * coordinate(n, 0); + blocked(n, 0) = true; + displacement(n, 1) = 0; + blocked(n, 1) = true; + } else if (Math::are_float_equal(coordinate(n, 0), 0.001)) { + displacement(n, 0) = epsilon * coordinate(n, 0); + blocked(n, 0) = true; + displacement(n, 1) = epsilon * coordinate(n, 1); + blocked(n, 1) = true; + } else if (Math::are_float_equal(coordinate(n, 1), 0.001)) { + displacement(n, 0) = epsilon * coordinate(n, 0); + blocked(n, 0) = true; + displacement(n, 1) = epsilon * coordinate(n, 1); + blocked(n, 1) = true; + } + } + + try { + model.solveStep(); + } catch (debug::Exception &e) { + } + + Int nb_iter = solver.get("nb_iterations"); + Real error = solver.get("error"); + bool converged = solver.get("converged"); + + if (converged) { + std::cout << "Converged in " << nb_iter << " iterations" << std::endl; + } else { + std::cout << "Didn't converge after " << nb_iter + << " iterations. Error is " << error << std::endl; + return EXIT_FAILURE; + } + + // analytical solution + Real solution_11 = 0.; + if (time < T) { + solution_11 = pre_mult * eps / T * + (Einf * time + lambda * Ev(0) * (1 - exp(-time / lambda))); + } else { + solution_11 = pre_mult * eps * + (Einf + + lambda * Ev(0) / T * + (exp((T - time) / lambda) - exp(-time / lambda))); + } + + // data output + output_data << s * time_step << " " << epsilon << " " << solution_11; + Array::const_matrix_iterator stress_it = stress.begin(dim, dim); + Array::const_matrix_iterator stress_end = stress.end(dim, dim); + for (; stress_it != stress_end; ++stress_it) { + output_data << " " << (*stress_it)(0, 0); + // test error + Real rel_error_11 = + std::abs(((*stress_it)(0, 0) - solution_11) / solution_11); + if (rel_error_11 > tolerance) { + std::cerr << "Relative error: " << rel_error_11 << std::endl; + return EXIT_FAILURE; + } + } + output_data << std::endl; + time += time_step; + + } + output_data.close(); + finalize(); + + std::cout << "Test successful!" << std::endl; + return EXIT_SUCCESS; +}