diff --git a/test/test_model/test_solid_mechanics_model/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/CMakeLists.txt
index da7f1059a..90943a0a3 100644
--- a/test/test_model/test_solid_mechanics_model/CMakeLists.txt
+++ b/test/test_model/test_solid_mechanics_model/CMakeLists.txt
@@ -1,107 +1,108 @@
#===============================================================================
# @file CMakeLists.txt
#
# @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
#
# @date creation: Fri Sep 03 2010
# @date last modification: Tue Jan 19 2016
#
# @brief configuratio for SolidMechanicsModel tests
#
# @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 <http://www.gnu.org/licenses/>.
#
# @section DESCRIPTION
#
#===============================================================================
add_akantu_test(test_materials "test_materials")
add_akantu_test(patch_tests "patch_tests")
add_akantu_test(test_cohesive "cohesive_test")
add_akantu_test(test_embedded_interface "test_embedded_interface")
add_akantu_test(test_energies "test energies")
#===============================================================================
#===============================================================================
add_mesh(test_cube3d_two_mat_mesh cube_two_materials.geo 3 1)
register_test(test_solid_mechanics_model_reassign_material
SOURCES test_solid_mechanics_model_reassign_material.cc
DEPENDS test_cube3d_two_mat_mesh
FILES_TO_COPY two_materials.dat
PACKAGE parallel implicit
PARALLEL
)
#===============================================================================
register_test(test_solid_mechanics_model_material_eigenstrain
SOURCES test_solid_mechanics_model_material_eigenstrain.cc
FILES_TO_COPY cube_3d_tet_4.msh; material_elastic_plane_strain.dat
PACKAGE implicit
)
#===============================================================================
register_test(test_material_selector
SOURCES test_material_selector.cc
FILES_TO_COPY material_selector.dat material_selector.msh
PACKAGE core
)
add_mesh(bar_segment_2 patch_tests/data/bar_segment.geo
DIM 1 ORDER 1 OUTPUT bar_segment_2.msh)
add_mesh(bar_segment_3 patch_tests/data/bar_segment.geo
DIM 1 ORDER 2 OUTPUT bar_segment_3.msh)
add_mesh(bar_triangle_3 patch_tests/data/bar_triangle.geo
DIM 2 ORDER 1 OUTPUT bar_triangle_3.msh)
add_mesh(bar_triangle_6 patch_tests/data/bar_triangle.geo
DIM 2 ORDER 2 OUTPUT bar_triangle_6.msh)
add_mesh(bar_quadrangle_4 patch_tests/data/bar_quadrangle.geo
DIM 2 ORDER 1 OUTPUT bar_quadrangle_4.msh)
add_mesh(bar_quadrangle_8 patch_tests/data/bar_quadrangle.geo
DIM 2 ORDER 2 OUTPUT bar_quadrangle_8.msh)
add_mesh(bar_tetrahedron_4 patch_tests/data/bar_tetrahedron.geo
DIM 3 ORDER 1 OUTPUT bar_tetrahedron_4.msh)
add_mesh(bar_tetrahedron_10 patch_tests/data/bar_tetrahedron.geo
DIM 3 ORDER 2 OUTPUT bar_tetrahedron_10.msh)
add_mesh(bar_hexahedron_8 patch_tests/data/bar_hexahedron.geo
DIM 3 ORDER 1 OUTPUT bar_hexahedron_8.msh)
add_mesh(bar_hexahedron_20 patch_tests/data/bar_hexahedron.geo
DIM 3 ORDER 2 OUTPUT bar_hexahedron_20.msh)
add_mesh(bar_pentahedron_6 patch_tests/data/bar_pentahedron.geo
DIM 3 ORDER 1 OUTPUT bar_pentahedron_6.msh)
add_mesh(bar_pentahedron_15 patch_tests/data/bar_pentahedron.geo
DIM 3 ORDER 2 OUTPUT bar_pentahedron_15.msh)
#===============================================================================
# dynamics tests
#===============================================================================
register_gtest_sources(
SOURCES test_solid_mechanics_model_dynamics.cc
FILES_TO_COPY test_solid_mechanics_model_dynamics_material.dat
PACKAGE core
+ PARALLEL
)
register_gtest_test(test_solid_mechanics_model
DEPENDS bar_segment_2 bar_segment_3
bar_triangle_3 bar_triangle_6
bar_quadrangle_4 bar_quadrangle_8
bar_tetrahedron_4 bar_tetrahedron_10
bar_hexahedron_8 bar_hexahedron_20
bar_pentahedron_6 bar_pentahedron_15
)
diff --git a/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_dynamics.cc b/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_dynamics.cc
index 72690cadb..11fa0cfd8 100644
--- a/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_dynamics.cc
+++ b/test/test_model/test_solid_mechanics_model/test_energies/test_solid_mechanics_model_work_dynamics.cc
@@ -1,158 +1,158 @@
/**
* @file test_solid_mechanics_model_work_dynamics.cc
*
* @author Tobias Brink <tobias.brink@epfl.ch>
*
* @date creation: Tue Dec 15 2017
* @date last modification: Dec Nov 15 2017
*
* @brief test work in dynamic simulations
*
* @section description
*
* Assuming that the kinetic energy and the potential energy of the
* linear elastic material are bug free, the work in a dynamic
* simulation must equal the change in internal energy (first law of
* thermodynamics). Work in dynamics is an infinitesimal work Fds,
* thus we need to integrate it and compare at the end. In this test,
* we use one Dirichlet boundary condition (with u = 0.0, 0.01, and
* -0.01) and one Neumann boundary condition for F on the opposite
* side. Then we do a few steps to get reference energies for work and
* internal energy. After more steps, the change in both work and
* internal energy must be equal.
*
* @section LICENSE
*
* Copyright (©) 2017 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "../test_solid_mechanics_model_fixture.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
namespace {
template <typename type_>
class TestSMMFixtureWorkDynamic : public TestSMMFixture<type_> {
public:
void SetUp() override {
this->mesh_file = "../bar" + aka::to_string(this->type) + ".msh";
TestSMMFixture<type_>::SetUp();
getStaticParser().parse("test_solid_mechanics_model_"
"work_material.dat");
/// model initialization
this->model->initFull();
/// Create a node group for Neumann BCs.
auto & apply_force_grp = this->mesh->createNodeGroup("apply_force");
auto & fixed_grp = this->mesh->createNodeGroup("fixed");
const auto & pos = this->mesh->getNodes();
const auto & lower = this->mesh->getLowerBounds();
const auto & upper = this->mesh->getUpperBounds();
UInt i = 0;
for (auto && posv : make_view(pos, this->spatial_dimension)) {
if (posv(_x) > upper(_x) - 1e-6) {
apply_force_grp.add(i);
} else if (posv(_x) < lower(_x) + 1e-6) {
fixed_grp.add(i);
}
++i;
}
this->mesh->createElementGroupFromNodeGroup("el_apply_force", "apply_force",
this->spatial_dimension - 1);
this->mesh->createElementGroupFromNodeGroup("el_fixed", "fixed",
this->spatial_dimension - 1);
Vector<Real> surface_traction(this->spatial_dimension);
surface_traction(_x) = 0.5;
if (this->spatial_dimension == 1) {
// TODO: this is a hack to work
// around non-implemented
// BC::Neumann::FromTraction for 1D
auto & force = this->model->getForce();
for (auto && pair : zip(make_view(pos, this->spatial_dimension),
make_view(force, this->spatial_dimension))) {
auto & posv = std::get<0>(pair);
auto & forcev = std::get<1>(pair);
if (posv(_x) > upper(_x) - 1e-6) {
forcev(_x) = surface_traction(_x);
}
}
} else {
this->model->applyBC(BC::Neumann::FromTraction(surface_traction),
"el_apply_force");
}
/// set up timestep
auto time_step = this->model->getStableTimeStep() * 0.1;
this->model->setTimeStep(time_step);
}
};
-TYPED_TEST_CASE(TestSMMFixtureWorkDynamic, types);
+TYPED_TEST_CASE(TestSMMFixtureWorkDynamic, gtest_element_types);
/* TODO: this is currently disabled for terrible results and performance
TYPED_TEST(TestSMMFixtureBar, WorkImplicit) {
test_body(*(this->model), *(this->mesh), _implicit_dynamic, 500);
}
*/
// model.assembleMassLumped();
TYPED_TEST(TestSMMFixtureWorkDynamic, WorkExplicit) {
/// Do the sim
std::vector<Real> displacements{0.00, 0.01, -0.01};
for (auto && u : displacements) {
this->model->applyBC(BC::Dirichlet::FixedValue(u, _x), "el_fixed");
// First, "equilibrate" a bit to get a reference state of total
// energy and work. This is needed when we have a Dirichlet with
// finite displacement on one side.
for (UInt i = 0; i < 25; ++i) {
this->model->solveStep();
}
// Again, work reported by Akantu is infinitesimal (dW) and we
// need to integrate a while to get a decent value.
double Etot0 =
this->model->getEnergy("potential") + this->model->getEnergy("kinetic");
double W = 0.0;
for (UInt i = 0; i < 200; ++i) {
/// Solve.
this->model->solveStep();
const auto dW = this->model->getEnergy("external work");
W += dW;
}
// Finally check.
const auto Epot = this->model->getEnergy("potential");
const auto Ekin = this->model->getEnergy("kinetic");
EXPECT_NEAR(W, Ekin + Epot - Etot0, 5e-2);
// Sadly not very exact for such a coarse mesh.
}
}
}
diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc
index 4d8b3a158..c049452ed 100644
--- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc
+++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_dynamics.cc
@@ -1,281 +1,314 @@
/**
* @file test_solid_mechanics_model_cube3d.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date creation: Wed Aug 04 2010
* @date last modification: Thu Aug 06 2015
*
* @brief test of the class SolidMechanicsModel on the 3d cube
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "boundary_condition_functor.hh"
#include "test_solid_mechanics_model_fixture.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
namespace {
constexpr Real A = 1e-1;
constexpr Real E = 1.;
constexpr Real poisson = .3;
constexpr Real lambda = E * poisson / (1. + poisson) / (1. - 2. * poisson);
constexpr Real mu = E / 2 / (1. + poisson);
constexpr Real rho = 1.;
const Real cp = std::sqrt((lambda + 2 * mu) / rho);
const Real cs = std::sqrt(mu / rho);
const Real c = std::sqrt(E / rho);
const Vector<Real> k = {.5, 0., 0.};
const Vector<Real> psi1 = {0., 0., 1.};
const Vector<Real> psi2 = {0., 1., 0.};
const Real knorm = k.norm();
-using verif_func = std::function<void(
- Vector<Real> & disp, const Vector<Real> & coord, Real current_time)>;
-
-template <UInt spatial_dimension> verif_func solution_disp;
-template <UInt spatial_dimension> verif_func solution_vel;
-
-template <>
-verif_func solution_disp<1> =
- [](Vector<Real> & disp, const Vector<Real> & coord, Real current_time) {
- const auto & x = coord(_x);
- const Real omega = c * k[0];
- disp(_x) = A * cos(k[0] * x - omega * current_time);
- };
-
-template <>
-verif_func solution_vel<1> =
- [](Vector<Real> & vel, const Vector<Real> & coord, Real current_time) {
- const auto & x = coord(_x);
- const Real omega = c * k[0];
- vel(_x) = omega * A * sin(k[0] * x - omega * current_time);
- };
-
-template <>
-verif_func solution_disp<2> =
- [](Vector<Real> & disp, const Vector<Real> & X, Real current_time) {
- Vector<Real> kshear = {k[1], k[0]};
- Vector<Real> kpush = {k[0], k[1]};
-
- const Real omega_p = knorm * cp;
- const Real omega_s = knorm * cs;
-
- Real phase_p = X.dot(kpush) - omega_p * current_time;
- Real phase_s = X.dot(kpush) - omega_s * current_time;
-
- disp = A * (kpush * cos(phase_p) + kshear * cos(phase_s));
- };
-
-template <>
-verif_func solution_vel<2> = [](Vector<Real> & vel, const Vector<Real> & X,
- Real current_time) {
- Vector<Real> kshear = {k[1], k[0]};
- Vector<Real> kpush = {k[0], k[1]};
-
- const Real omega_p = knorm * cp;
- const Real omega_s = knorm * cs;
-
- Real phase_p = X.dot(kpush) - omega_p * current_time;
- Real phase_s = X.dot(kpush) - omega_s * current_time;
-
- vel = A * (kpush * omega_p * cos(phase_p) + kshear * omega_s * cos(phase_s));
+/* -------------------------------------------------------------------------- */
+template <UInt dim> struct Verification {};
+
+/* -------------------------------------------------------------------------- */
+template <> struct Verification<1> {
+ void displacement(Vector<Real> & disp, const Vector<Real> & coord,
+ Real current_time) {
+ const auto & x = coord(_x);
+ const Real omega = c * k[0];
+ disp(_x) = A * cos(k[0] * x - omega * current_time);
+ }
+
+ void velocity(Vector<Real> & vel, const Vector<Real> & coord,
+ Real current_time) {
+ const auto & x = coord(_x);
+ const Real omega = c * k[0];
+ vel(_x) = omega * A * sin(k[0] * x - omega * current_time);
+ }
};
-template <>
-verif_func solution_disp<3> =
- [](Vector<Real> & disp, const Vector<Real> & coord, Real current_time) {
- const auto & X = coord;
- Vector<Real> kpush = k;
- Vector<Real> kshear1(3);
- Vector<Real> kshear2(3);
- kshear1.crossProduct(k, psi1);
- kshear2.crossProduct(k, psi2);
+/* -------------------------------------------------------------------------- */
+template <> struct Verification<2> {
+ void displacement(Vector<Real> & disp, const Vector<Real> & X,
+ Real current_time) {
+ Vector<Real> kshear = {k[1], k[0]};
+ Vector<Real> kpush = {k[0], k[1]};
- const Real omega_p = knorm * cp;
- const Real omega_s = knorm * cs;
+ const Real omega_p = knorm * cp;
+ const Real omega_s = knorm * cs;
- Real phase_p = X.dot(kpush) - omega_p * current_time;
- Real phase_s = X.dot(kpush) - omega_s * current_time;
+ Real phase_p = X.dot(kpush) - omega_p * current_time;
+ Real phase_s = X.dot(kpush) - omega_s * current_time;
- disp = A * (kpush * cos(phase_p) + kshear1 * cos(phase_s) +
- kshear2 * cos(phase_s));
+ disp = A * (kpush * cos(phase_p) + kshear * cos(phase_s));
+ }
+
+ void velocity(Vector<Real> & vel, const Vector<Real> & X, Real current_time) {
+ Vector<Real> kshear = {k[1], k[0]};
+ Vector<Real> kpush = {k[0], k[1]};
- };
+ const Real omega_p = knorm * cp;
+ const Real omega_s = knorm * cs;
-template <>
-verif_func solution_vel<3> = [](Vector<Real> & vel, const Vector<Real> & coord,
- Real current_time) {
- const auto & X = coord;
- Vector<Real> kpush = k;
- Vector<Real> kshear1(3);
- Vector<Real> kshear2(3);
- kshear1.crossProduct(k, psi1);
- kshear2.crossProduct(k, psi2);
+ Real phase_p = X.dot(kpush) - omega_p * current_time;
+ Real phase_s = X.dot(kpush) - omega_s * current_time;
- const Real omega_p = knorm * cp;
- const Real omega_s = knorm * cs;
+ vel =
+ A * (kpush * omega_p * cos(phase_p) + kshear * omega_s * cos(phase_s));
+ }
+};
- Real phase_p = X.dot(kpush) - omega_p * current_time;
- Real phase_s = X.dot(kpush) - omega_s * current_time;
+/* -------------------------------------------------------------------------- */
+template <> struct Verification<3> {
+ void displacement(Vector<Real> & disp, const Vector<Real> & coord,
+ Real current_time) {
+ const auto & X = coord;
+ Vector<Real> kpush = k;
+ Vector<Real> kshear1(3);
+ Vector<Real> kshear2(3);
+ kshear1.crossProduct(k, psi1);
+ kshear2.crossProduct(k, psi2);
+
+ const Real omega_p = knorm * cp;
+ const Real omega_s = knorm * cs;
+
+ Real phase_p = X.dot(kpush) - omega_p * current_time;
+ Real phase_s = X.dot(kpush) - omega_s * current_time;
+
+ disp = A * (kpush * cos(phase_p) + kshear1 * cos(phase_s) +
+ kshear2 * cos(phase_s));
+ }
- vel = A * (kpush * omega_p * cos(phase_p) + kshear1 * omega_s * cos(phase_s) +
+ void velocity(Vector<Real> & vel, const Vector<Real> & coord,
+ Real current_time) {
+ const auto & X = coord;
+ Vector<Real> kpush = k;
+ Vector<Real> kshear1(3);
+ Vector<Real> kshear2(3);
+ kshear1.crossProduct(k, psi1);
+ kshear2.crossProduct(k, psi2);
+
+ const Real omega_p = knorm * cp;
+ const Real omega_s = knorm * cs;
+
+ Real phase_p = X.dot(kpush) - omega_p * current_time;
+ Real phase_s = X.dot(kpush) - omega_s * current_time;
+
+ vel =
+ A * (kpush * omega_p * cos(phase_p) + kshear1 * omega_s * cos(phase_s) +
kshear2 * omega_s * cos(phase_s));
+ }
};
+/* -------------------------------------------------------------------------- */
template <ElementType _type>
class SolutionFunctor : public BC::Dirichlet::DirichletFunctor {
public:
SolutionFunctor(Real current_time, SolidMechanicsModel & model)
: current_time(current_time), model(model) {}
public:
- // static constexpr UInt dim = DimensionHelper<_type>::dim;
static constexpr UInt dim = ElementClass<_type>::getSpatialDimension();
inline void operator()(UInt node, Vector<bool> & flags, Vector<Real> & primal,
const Vector<Real> & coord) const {
flags(0) = true;
auto & vel = model.getVelocity();
- auto itVel = vel.begin(model.getSpatialDimension());
- Vector<Real> v = itVel[node];
- solution_disp<dim>(primal, coord, current_time);
- solution_vel<dim>(v, coord, current_time);
+ auto it = vel.begin(model.getSpatialDimension());
+ Vector<Real> v = it[node];
+
+ Verification<dim> verif;
+ verif.displacement(primal, coord, current_time);
+ verif.velocity(v, coord, current_time);
}
private:
Real current_time;
SolidMechanicsModel & model;
};
-template <ElementType _type, typename AM>
-void test_body(SolidMechanicsModel & model, AM analysis_method) {
- // constexpr UInt dim = DimensionHelper<_type>::dim;
- static constexpr UInt dim = ElementClass<_type>::getSpatialDimension();
+/* -------------------------------------------------------------------------- */
+// This fixture uses somewhat finer meshes representing bars.
+template <typename type_>
+class TestSMMFixtureBar
+ : public TestSMMFixture<std::tuple_element_t<0, type_>> {
+ using parent = TestSMMFixture<std::tuple_element_t<0, type_>>;
- getStaticParser().parse("test_solid_mechanics_model_"
- "dynamics_material.dat");
-
- model.initFull(SolidMechanicsModelOptions(analysis_method));
-
- bool dump_paraview = false;
-
- if (dump_paraview) {
- std::stringstream base_name;
- base_name << "bar" << analysis_method << _type;
- model.setBaseName(base_name.str());
- model.addDumpFieldVector("displacement");
- model.addDumpField("mass");
- model.addDumpField("velocity");
- model.addDumpField("acceleration");
- model.addDumpFieldVector("external_force");
- model.addDumpFieldVector("internal_force");
- model.addDumpField("stress");
- model.addDumpField("strain");
- }
+public:
+ void SetUp() override {
+ this->mesh_file = "bar" + aka::to_string(this->type) + ".msh";
+ parent::SetUp();
+
+ getStaticParser().parse("test_solid_mechanics_model_"
+ "dynamics_material.dat");
+
+ auto analysis_method = std::tuple_element_t<1, type_>::value;
+ this->model->initFull(_analysis_method = analysis_method);
+
+ if (this->dump_paraview) {
+ std::stringstream base_name;
+ base_name << "bar" << analysis_method << this->type;
+ this->model->setBaseName(base_name.str());
+ this->model->addDumpFieldVector("displacement");
+ this->model->addDumpField("mass");
+ this->model->addDumpField("velocity");
+ this->model->addDumpField("acceleration");
+ this->model->addDumpFieldVector("external_force");
+ this->model->addDumpFieldVector("internal_force");
+ this->model->addDumpField("stress");
+ this->model->addDumpField("strain");
+ }
+
+ time_step = this->model->getStableTimeStep() / 10.;
+ this->model->setTimeStep(time_step);
+ // std::cout << "timestep: " << time_step << std::endl;
- Real time_step = model.getStableTimeStep() / 10.;
- model.setTimeStep(time_step);
- std::cout << "timestep: " << time_step << std::endl;
+ const auto & position = this->mesh->getNodes();
+ auto & displacement = this->model->getDisplacement();
+ auto & velocity = this->model->getVelocity();
- UInt nb_nodes = model.getMesh().getNbNodes();
- UInt spatial_dimension = model.getSpatialDimension();
+ constexpr auto dim = parent::spatial_dimension;
- auto & nodes = model.getMesh().getNodes();
- auto & disp = model.getDisplacement();
- auto & vel = model.getVelocity();
+ Verification<dim> verif;
- Array<Real> disp_solution(nb_nodes, spatial_dimension);
+ for (auto && tuple :
+ zip(make_view(position, dim), make_view(displacement, dim),
+ make_view(velocity, dim))) {
+ verif.displacement(std::get<1>(tuple), std::get<0>(tuple), 0.);
+ verif.velocity(std::get<2>(tuple), std::get<0>(tuple), 0.);
+ }
+
+ if (dump_paraview)
+ this->model->dump();
+
+ /// boundary conditions
+ this->model->applyBC(SolutionFunctor<parent::type>(0., *this->model), "BC");
- Real current_time = 0;
+ wave_velocity = 1.; // sqrt(E/rho) = sqrt(1/1) = 1
+ simulation_time = 5 / wave_velocity;
- auto itNodes = nodes.begin(spatial_dimension);
- auto itDisp = disp.begin(spatial_dimension);
- auto itVel = vel.begin(spatial_dimension);
- for (UInt n = 0; n < nb_nodes; ++n, ++itNodes, ++itDisp, ++itVel) {
- solution_disp<dim>(*itDisp, *itNodes, current_time);
- solution_vel<dim>(*itVel, *itNodes, current_time);
+ max_steps = simulation_time / time_step; // 100
+ auto ndump = 50;
+ dump_freq = max_steps / ndump;
}
- if (dump_paraview)
- model.dump();
+protected:
+ Real time_step;
+ Real wave_velocity;
+ Real simulation_time;
+ UInt max_steps;
+ UInt dump_freq;
+ bool dump_paraview{false};
+};
- /// boundary conditions
- model.applyBC(SolutionFunctor<_type>(current_time, model), "BC");
+template <AnalysisMethod t>
+using analysis_method_t = std::integral_constant<AnalysisMethod, t>;
- Real max_error = 0.;
- Real wave_velocity = 1.; // sqrt(E/rho) = sqrt(1/1) = 1
- Real simulation_time = 5 / wave_velocity;
+#ifdef AKANTU_IMPLICIT
+using TestAnalysisTypes =
+ std::tuple<analysis_method_t<_implicit_dynamic, _explicit_lumped_mass>>;
+#else
+using TestAnalysisTypes = std::tuple<analysis_method_t<_explicit_lumped_mass>>;
+#endif
- UInt max_steps = simulation_time / time_step; // 100
- UInt ndump = 50;
- UInt dump_freq = max_steps / ndump;
+using TestTypes =
+ gtest_list_t<cross_product_t<TestElementTypes, TestAnalysisTypes>>;
- std::cout << "max_steps: " << max_steps << std::endl;
+TYPED_TEST_CASE(TestSMMFixtureBar, TestTypes);
- for (UInt s = 0; s < max_steps; ++s, current_time += time_step) {
+/* -------------------------------------------------------------------------- */
+
+TYPED_TEST(TestSMMFixtureBar, DynamicsExplicit) {
+ constexpr auto dim = TestFixture::spatial_dimension;
+ Real current_time = 0.;
+ const auto & position = this->mesh->getNodes();
+ const auto & displacement = this->model->getDisplacement();
+ UInt nb_nodes = this->mesh->getNbNodes();
+ UInt nb_global_nodes = this->mesh->getNbGlobalNodes();
+ Real max_error{0.};
- if (s % dump_freq == 0 && dump_paraview)
- model.dump();
+ Array<Real> displacement_solution(nb_nodes, dim);
+
+ Verification<dim> verif;
+
+ for (UInt s = 0; s < this->max_steps; ++s, current_time += this->time_step) {
+ if (s % this->dump_freq == 0 && this->dump_paraview)
+ this->model->dump();
/// boundary conditions
- model.applyBC(SolutionFunctor<_type>(current_time, model), "BC");
+ this->model->applyBC(
+ SolutionFunctor<TestFixture::type>(current_time, *this->model), "BC");
// compute the disp solution
- auto itDispSolution = disp_solution.begin(spatial_dimension);
- itNodes = nodes.begin(spatial_dimension);
- for (UInt n = 0; n < nb_nodes; ++n, ++itNodes, ++itDispSolution) {
- solution_disp<dim>(*itDispSolution, *itNodes, current_time);
+ for (auto && tuple :
+ zip(make_view(position, dim), make_view(displacement_solution, dim))) {
+ verif.displacement(std::get<1>(tuple), std::get<0>(tuple), current_time);
}
+
// compute the error solution
- itDispSolution = disp_solution.begin(spatial_dimension);
- itDisp = disp.begin(spatial_dimension);
Real disp_error = 0.;
- for (UInt n = 0; n < nb_nodes; ++n, ++itDispSolution, ++itDisp) {
- auto diff = *itDispSolution - *itDisp;
-
- for (UInt i = 0; i < spatial_dimension; ++i) {
- disp_error += diff(i) * diff(i);
- }
+ for (auto && tuple : zip(make_view(displacement, dim),
+ make_view(displacement_solution, dim))) {
+ auto diff = std::get<1>(tuple) - std::get<0>(tuple);
+ disp_error += diff.dot(diff);
}
- disp_error = sqrt(disp_error) / nb_nodes;
+
+ this->mesh->getCommunicator().allReduce(disp_error,
+ SynchronizerOperation::_sum);
+
+ disp_error = sqrt(disp_error) / nb_global_nodes;
max_error = std::max(disp_error, max_error);
- ASSERT_NEAR(disp_error, 0., 1e-3);
- model.solveStep();
- }
- std::cout << "max error: " << max_error << std::endl;
-}
-#ifdef AKANTU_IMPLICIT
-TYPED_TEST(TestSMMFixtureBar, DynamicsImplicit) {
- test_body<TestFixture::type>(*(this->model), _implicit_dynamic);
-}
-#endif
+ ASSERT_NEAR(disp_error, 0., 1e-3);
-TYPED_TEST(TestSMMFixtureBar, DynamicsExplicit) {
- test_body<TestFixture::type>(*(this->model), _explicit_lumped_mass);
+ this->model->solveStep();
+ }
+ // std::cout << "max error: " << max_error << std::endl;
}
}
diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_fixture.hh b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_fixture.hh
index beef6ad66..431a03c94 100644
--- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_fixture.hh
+++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_fixture.hh
@@ -1,63 +1,61 @@
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model.hh"
#include "test_gtest_utils.hh"
+#include "communicator.hh"
/* -------------------------------------------------------------------------- */
#include <gtest/gtest.h>
#include <vector>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_TEST_SOLID_MECHANICS_MODEL_FIXTURE_HH__
#define __AKANTU_TEST_SOLID_MECHANICS_MODEL_FIXTURE_HH__
using namespace akantu;
// This fixture uses very small meshes with a volume of 1.
template <typename type_> class TestSMMFixture : public ::testing::Test {
public:
static constexpr ElementType type = type_::value;
static constexpr size_t spatial_dimension =
ElementClass<type>::getSpatialDimension();
void SetUp() override {
mesh = std::make_unique<Mesh>(this->spatial_dimension);
- mesh->read(this->mesh_file);
+#if defined(AKANTU_PARALLEL)
+ if(Communicator::getStaticCommunicator().whoAmI() == 0) {
+#endif
+
+ mesh->read(this->mesh_file);
+
+#if defined(AKANTU_PARALLEL)
+ }
+
+ mesh->distribute();
+#endif
SCOPED_TRACE(aka::to_string(this->type).c_str());
model = std::make_unique<SolidMechanicsModel>(*mesh, _all_dimensions,
aka::to_string(this->type));
}
void TearDown() override {
model.reset(nullptr);
mesh.reset(nullptr);
}
protected:
std::string mesh_file{aka::to_string(this->type) + ".msh"};
std::unique_ptr<Mesh> mesh;
std::unique_ptr<SolidMechanicsModel> model;
};
template <typename type_> constexpr ElementType TestSMMFixture<type_>::type;
template <typename type_> constexpr size_t TestSMMFixture<type_>::spatial_dimension;
-using types = gtest_list_t<TestElementTypes>;
-
-TYPED_TEST_CASE(TestSMMFixture, types);
-
-
-// This fixture uses somewhat finer meshes representing bars.
-template <typename type_>
-class TestSMMFixtureBar : public TestSMMFixture<type_> {
-public:
- void SetUp() override {
- this->mesh_file = "bar" + aka::to_string(this->type) + ".msh";
- TestSMMFixture<type_>::SetUp();
- }
-};
+using gtest_element_types = gtest_list_t<TestElementTypes>;
-TYPED_TEST_CASE(TestSMMFixtureBar, types);
+TYPED_TEST_CASE(TestSMMFixture, gtest_element_types);
#endif /* __AKANTU_TEST_SOLID_MECHANICS_MODEL_FIXTURE_HH__ */