Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/tests/test_model.cpp
23774fdb7372features/netcdf-mpi

Last Change
Raw File

test_model.cpp
View Options

/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "isotropic_hardening.hh"
#include "model_factory.hh"
#include "test.hh"
#include <random>
/* -------------------------------------------------------------------------- */
using namespace tamaas;
TEST(TestModel, applyElasticity) {
auto model =
ModelFactory::createModel(model_type::volume_2d, {1., 1., 1.}, {1, 1, 1});
Grid<Real, 3> gradient({1, 1, 1}, 9), stress({1, 1, 1}, 9);
// Random data objects
std::random_device rnd;
std::mt19937 mt(rnd());
std::normal_distribution<> dis(0, 1);
// Filling gradient with random data (sequential)
for (auto& du : gradient)
du = dis(mt);
std::uniform_real_distribution<> unif(0, 0.5);
model->setElasticity(std::abs(dis(mt)) * unif(mt), unif(mt));
// Computing stresses with seperate array
model->applyElasticity(stress, gradient);
// Checking correct isotropic elasticity
auto mu = model->getShearModulus(), nu = model->getPoissonRatio();
auto lambda = 2 * mu * nu / (1 - 2 * nu);
MatrixProxy<Real, 3, 3> grad(gradient(0));
auto trace = grad.trace();
Matrix<Real, 3, 3> sigma;
for (UInt i = 0; i < 3; ++i)
for (UInt j = 0; j < 3; ++j)
sigma(i, j) = (i == j) * lambda * trace + mu * (grad(i, j) + grad(j, i));
EXPECT_TRUE(compare(stress, sigma, AreFloatEqual())) << "Elasticity fail";
// Computing stress with same array
model->applyElasticity(gradient, gradient);
EXPECT_TRUE(compare(gradient, stress, AreFloatEqual()))
<< "Applying elasticity in-place fail";
}
TEST(TestModel, applyElasticitySym) {
auto model =
ModelFactory::createModel(model_type::volume_2d, {1., 1., 1.}, {1, 1, 1});
Grid<Real, 3> gradient({1, 1, 1}, 6), stress({1, 1, 1}, 6);
// Random data objects
std::random_device rnd;
std::mt19937 mt(rnd());
std::normal_distribution<> dis(0, 1);
// Filling gradient with random data (sequential)
for (auto& du : gradient)
du = dis(mt);
std::uniform_real_distribution<> unif(0, 0.5);
model->setElasticity(std::abs(dis(mt)) * unif(mt), unif(mt));
// Computing stresses with seperate array
model->applyElasticity(stress, gradient);
// Checking correct isotropic elasticity
auto mu = model->getShearModulus(), nu = model->getPoissonRatio();
auto lambda = 2 * mu * nu / (1 - 2 * nu);
SymMatrixProxy<Real, 3> sym_grad(gradient(0));
auto grad = dense(sym_grad);
auto trace = grad.trace();
Matrix<Real, 3, 3> sigma;
for (UInt i = 0; i < 3; ++i)
for (UInt j = 0; j < 3; ++j)
sigma(i, j) = (i == j) * lambda * trace + mu * (grad(i, j) + grad(j, i));
EXPECT_TRUE(compare(stress, symmetrize(sigma), AreFloatEqual()))
<< "Elasticity fail";
// Computing stress with same array
model->applyElasticity(gradient, gradient);
EXPECT_TRUE(compare(gradient, stress, AreFloatEqual()))
<< "Applying elasticity in-place fail";
}
TEST(TestIsotropicHardening, computePlasticIncrement) {
auto model =
ModelFactory::createModel(model_type::volume_2d, {1., 1., 1.}, {1, 1, 1});
auto sigma_0 = 1., h = 0.1;
IsotropicHardening<model_type::volume_2d> hardening(model.get(), sigma_0, h);
Grid<Real, 3> strain({1, 1, 1}, 6), strain_increment({1, 1, 1}, 6),
plastic_strain_increment({1, 1, 1}, 6), solution({1, 1, 1}, 6);
auto E = model->getYoungModulus();
auto nu = model->getPoissonRatio();
auto mu = model->getShearModulus();
// uniform tension state
Real sigma = 0.9;
strain_increment(0, 0, 0, 0) = sigma / E;
strain_increment(0, 0, 0, 1) = -nu * sigma / E;
strain_increment(0, 0, 0, 2) = -nu * sigma / E;
hardening.computePlasticIncrement<false>(plastic_strain_increment, strain,
strain_increment);
// just checking everything is zero
EXPECT_TRUE(compare(solution, plastic_strain_increment, AreFloatEqual()))
<< "Elastic radial return fail (plastic increment)";
// uniform tension state
sigma = 1.1;
strain_increment(0, 0, 0, 0) = sigma / E;
strain_increment(0, 0, 0, 1) = -nu * sigma / E;
strain_increment(0, 0, 0, 2) = -nu * sigma / E;
hardening.computePlasticIncrement<false>(plastic_strain_increment, strain,
strain_increment);
Real a = (sigma - sigma_0) / (3 * mu + h);
solution(0, 0, 0, 0) = a;
solution(0, 0, 0, 1) = -a / 2;
solution(0, 0, 0, 2) = -a / 2;
EXPECT_TRUE(compare(solution, plastic_strain_increment, AreFloatEqual()))
<< "Plastic radial return fail (plastic increment)";
}
c4science · Help