diff --git a/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_hexahedron.geo b/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_hexahedron.geo
index cd7c6a709..6ca6ce290 100644
--- a/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_hexahedron.geo
+++ b/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_hexahedron.geo
@@ -1,79 +1,76 @@
h = .6;
Point(1) = {0, 0, 0, h};
Point(2) = {10, 0, 0, h};
Point(3) = {0, 1., 0, h};
Point(4) = {10, 1., 0, h};
Point(5) = {0, 0, 1., h};
Point(6) = {10, 0, 1., h};
Point(7) = {0, 1., 1., h};
Point(8) = {10, 1., 1., h};
Mesh.SecondOrderIncomplete = 1;//+
Line(1) = {7, 8};
//+
Line(2) = {8, 6};
//+
Line(3) = {6, 5};
//+
Line(4) = {5, 7};
//+
Line(5) = {3, 7};
//+
Line(6) = {3, 1};
//+
Line(7) = {1, 5};
//+
Line(8) = {4, 2};
//+
Line(9) = {4, 8};
//+
Line(10) = {2, 6};
//+
Line(11) = {4, 3};
//+
Line(12) = {1, 2};
//+
Line Loop(1) = {2, -10, -8, 9};
//+
Plane Surface(1) = {1};
//+
Line Loop(2) = {7, 4, -5, 6};
//+
Plane Surface(2) = {2};
//+
Line Loop(3) = {12, 10, 3, -7};
//+
Plane Surface(3) = {3};
//+
Line Loop(4) = {3, 4, 1, 2};
//+
Plane Surface(4) = {4};
//+
Line Loop(5) = {1, -9, 11, 5};
//+
Plane Surface(5) = {5};
//+
Line Loop(6) = {11, 6, 12, -8};
//+
Plane Surface(6) = {6};
//+
Surface Loop(1) = {3, 6, 5, 4, 2, 1};
//+
Volume(1) = {1};
//+
-Physical Surface("Left") = {2};
-//+
-Physical Surface("Right") = {1};
-//+
+Physical Surface("BC") = {1, 2, 3, 4, 5, 6};
Physical Volume("bulk") = {1};
Transfinite Surface "*";
Transfinite Volume "*";
Recombine Surface "*";
Mesh.SecondOrderIncomplete = 1;
\ No newline at end of file
diff --git a/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_quadrangle.geo b/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_quadrangle.geo
index 4826d400d..f47d77171 100644
--- a/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_quadrangle.geo
+++ b/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_quadrangle.geo
@@ -1,27 +1,26 @@
-h = .4;
+h = .6;
Point(1) = {0, 0, 0, h};
Point(2) = {10, 0, 0, h};
Point(3) = {0, 1., 0, h};
Point(4) = {10, 1., 0, h};
Line(1) = {1, 2};
//+
Line(2) = {4, 3};
//+
Line(3) = {2, 4};
//+
Line(4) = {3, 1};
//+
Line Loop(1) = {1, 3, 2, 4};
//+
Plane Surface(1) = {1};
-Physical Line("Left") = {4};
-Physical Line("Right") = {3};
-//+
+Physical Line("BC") = {4, 3, 1, 2};
Physical Surface("bulk") = {1};
-//+
+
+Transfinite Surface "*";
Recombine Surface {1};
Mesh.SecondOrderIncomplete = 1;
\ No newline at end of file
diff --git a/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_segment.geo b/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_segment.geo
index 7aa27ae0a..ab49b3850 100644
--- a/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_segment.geo
+++ b/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_segment.geo
@@ -1,8 +1,7 @@
h = .6;
Point(1) = {0, 0, 0, h};
Point(2) = {10, 0, 0, h};
Line(1) = {1, 2};
-Physical Point("Left") = {1};
-Physical Point("Right") = {2};
-//+
+
+Physical Point("BC") = {1, 2};
Physical Line("bulk") = {1};
diff --git a/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_tetrahedron.geo b/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_tetrahedron.geo
index 1df0b9749..fbcfb6e33 100644
--- a/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_tetrahedron.geo
+++ b/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_tetrahedron.geo
@@ -1,73 +1,70 @@
h = .6;
Point(1) = {0, 0, 0, h};
Point(2) = {10, 0, 0, h};
Point(3) = {0, 1., 0, h};
Point(4) = {10, 1., 0, h};
Point(5) = {0, 0, 1., h};
Point(6) = {10, 0, 1., h};
Point(7) = {0, 1., 1., h};
Point(8) = {10, 1., 1., h};
Mesh.SecondOrderIncomplete = 1;//+
Line(1) = {7, 8};
//+
Line(2) = {8, 6};
//+
Line(3) = {6, 5};
//+
Line(4) = {5, 7};
//+
Line(5) = {3, 7};
//+
Line(6) = {3, 1};
//+
Line(7) = {1, 5};
//+
Line(8) = {4, 2};
//+
Line(9) = {4, 8};
//+
Line(10) = {2, 6};
//+
Line(11) = {4, 3};
//+
Line(12) = {1, 2};
//+
Line Loop(1) = {2, -10, -8, 9};
//+
Plane Surface(1) = {1};
//+
Line Loop(2) = {7, 4, -5, 6};
//+
Plane Surface(2) = {2};
//+
Line Loop(3) = {12, 10, 3, -7};
//+
Plane Surface(3) = {3};
//+
Line Loop(4) = {3, 4, 1, 2};
//+
Plane Surface(4) = {4};
//+
Line Loop(5) = {1, -9, 11, 5};
//+
Plane Surface(5) = {5};
//+
Line Loop(6) = {11, 6, 12, -8};
//+
Plane Surface(6) = {6};
//+
Surface Loop(1) = {3, 6, 5, 4, 2, 1};
//+
Volume(1) = {1};
//+
-Physical Surface("Left") = {2};
-//+
-Physical Surface("Right") = {1};
-//+
+Physical Surface("BC") = {1, 2, 3, 4, 5, 6};
Physical Volume("bulk") = {1};
diff --git a/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_triangle.geo b/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_triangle.geo
index 240c54287..6fa5be9e8 100644
--- a/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_triangle.geo
+++ b/test/test_model/test_solid_mechanics_model/patch_tests/data/bar_triangle.geo
@@ -1,24 +1,23 @@
h = .6;
Point(1) = {0, 0, 0, h};
Point(2) = {10, 0, 0, h};
Point(3) = {0, 1., 0, h};
Point(4) = {10, 1., 0, h};
Line(1) = {1, 2};
//+
Line(2) = {4, 3};
//+
Line(3) = {2, 4};
//+
Line(4) = {3, 1};
//+
Line Loop(1) = {1, 3, 2, 4};
//+
Plane Surface(1) = {1};
-Physical Line("Left") = {4};
-Physical Line("Right") = {3};
-//+
+Physical Line("BC") = {4, 3, 1 , 2};
+
Physical Surface("bulk") = {1};
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 9744e0d8e..a099bbd09 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,243 +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 {
template <typename type_>
class TestSMMFixtureDynamics : public TestSMMFixture<type_> {
public:
static constexpr ElementType type = type_::value;
void SetUp() override {
if (this->type == _pentahedron_6 || this->type == _pentahedron_15)
return;
std::cout << "testing type " << this->type << std::endl;
TestSMMFixture<type_>::SetUp();
}
std::string makeMeshName() override {
std::stringstream element_type;
element_type << type;
SCOPED_TRACE(element_type.str().c_str());
return std::string("bar") + element_type.str() + ".msh";
}
};
-template <UInt spatial_dimension>
-auto solution_disp =
+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();
+
+template <UInt spatial_dimension> auto solution_disp = 1;
+template <UInt spatial_dimension> auto solution_vel = 1;
+
+template <>
+auto solution_disp<1> =
[](Vector<Real> & disp, const Vector<Real> & coord, Real current_time) {
const auto & x = coord(_x);
- constexpr Real k = .5;
- constexpr Real omega = k;
- disp(_x) = cos(k * x - omega * current_time);
+ const Real omega = c * k[0];
+ disp(_x) = A * cos(k[0] * x - omega * current_time);
};
-template <UInt spatial_dimension>
-auto solution_vel =
+template <>
+auto solution_vel<1> =
[](Vector<Real> & vel, const Vector<Real> & coord, Real current_time) {
-
const auto & x = coord(_x);
- constexpr Real k = .5;
- constexpr Real omega = k;
- vel(_x) = omega * sin(k * x - omega * current_time);
+ const Real omega = c * k[0];
+ vel(_x) = omega * A * sin(k[0] * x - omega * current_time);
};
-template <ElementType _type> struct DimensionHelper {
- static constexpr int dim = -1;
-};
+template <>
+auto solution_disp<2> =
+ [](Vector<Real> & disp, const Vector<Real> & coord, Real current_time) {
+
+ const auto & X = coord;
+ 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 <>
+auto solution_vel<2> = [](Vector<Real> & vel, const Vector<Real> & coord,
+ Real current_time) {
+
+ const auto & X = coord;
+ 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 <> struct DimensionHelper<_segment_2> {
- static constexpr UInt dim = 1;
-};
-template <> struct DimensionHelper<_segment_3> {
- static constexpr UInt dim = 1;
-};
-template <> struct DimensionHelper<_triangle_3> {
- static constexpr UInt dim = 2;
-};
-template <> struct DimensionHelper<_triangle_6> {
- static constexpr UInt dim = 2;
-};
-template <> struct DimensionHelper<_quadrangle_4> {
- static constexpr UInt dim = 2;
};
-template <> struct DimensionHelper<_quadrangle_8> {
- static constexpr UInt dim = 2;
+
+template <>
+auto 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);
+
+ 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));
+
+ };
+
+template <>
+auto 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);
+
+ 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);
}
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();
getStaticParser().parse("test_solid_mechanics_model_"
"dynamics_material.dat");
model.initFull(SolidMechanicsModelOptions(analysis_method));
- bool dump_paraview = false;
+ bool dump_paraview = true;
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");
}
Real time_step = model.getStableTimeStep() / 10.;
model.setTimeStep(time_step);
std::cout << "timestep: " << time_step << std::endl;
UInt nb_nodes = model.getMesh().getNbNodes();
UInt spatial_dimension = model.getSpatialDimension();
auto & nodes = model.getMesh().getNodes();
auto & disp = model.getDisplacement();
auto & vel = model.getVelocity();
Array<Real> disp_solution(nb_nodes, spatial_dimension);
Real current_time = 0;
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);
}
if (dump_paraview)
model.dump();
/// boundary conditions
- model.applyBC(SolutionFunctor<_type>(current_time, model), "Left");
- model.applyBC(SolutionFunctor<_type>(current_time, model), "Right");
+ model.applyBC(SolutionFunctor<_type>(current_time, model), "BC");
Real max_error = 0.;
Real wave_velocity = 1.; // sqrt(E/rho) = sqrt(1/1) = 1
Real simulation_time = 5 / wave_velocity;
UInt max_steps = simulation_time / time_step; // 100
+ UInt ndump = 50;
+ UInt dump_freq = max_steps / ndump;
+
std::cout << "max_steps: " << max_steps << std::endl;
for (UInt s = 0; s < max_steps; ++s, current_time += time_step) {
- if (dump_paraview)
+ if (s % dump_freq == 0 && dump_paraview)
model.dump();
/// boundary conditions
- model.applyBC(SolutionFunctor<_type>(current_time, model), "Left");
- model.applyBC(SolutionFunctor<_type>(current_time, model), "Right");
+ model.applyBC(SolutionFunctor<_type>(current_time, 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);
}
// 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);
}
}
disp_error = sqrt(disp_error) / nb_nodes;
max_error = std::max(disp_error, max_error);
- ASSERT_NEAR(disp_error, 0., 1e-1);
+ ASSERT_NEAR(disp_error, 0., 1e-3);
model.solveStep();
}
std::cout << "max error: " << max_error << std::endl;
}
TYPED_TEST_CASE(TestSMMFixtureDynamics, types);
#ifdef AKANTU_IMPLICIT
TYPED_TEST(TestSMMFixture, DynamicsImplicit) {
if (this->type != _pentahedron_6 && this->type != _pentahedron_15)
test_body<this->type>(*(this->model), _implicit_dynamic);
}
#endif
TYPED_TEST(TestSMMFixtureDynamics, DynamicsExplicit) {
if (this->type != _pentahedron_6 && this->type != _pentahedron_15)
test_body<this->type>(*(this->model), _explicit_lumped_mass);
}
}