test_elastic_materials.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Thu, Jun 13, 08:22

test_elastic_materials.cc
View Options

	/* -------------------------------------------------------------------------- */
	#include "material_elastic.hh"
	#include "material_elastic_orthotropic.hh"
	#include "solid_mechanics_model.hh"
	#include "test_material_fixtures.hh"
	/* -------------------------------------------------------------------------- */
	#include <gtest/gtest.h>
	#include <type_traits>
	/* -------------------------------------------------------------------------- */

	using namespace akantu;

	using types =
	::testing::Types<Traits<MaterialElastic, 1>, Traits<MaterialElastic, 2>,
	Traits<MaterialElastic, 3>,

	Traits<MaterialElasticOrthotropic, 2>,
	Traits<MaterialElasticOrthotropic, 3>,

	Traits<MaterialElasticLinearAnisotropic, 2>,
	Traits<MaterialElasticLinearAnisotropic, 3>>;

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElastic<1>>::testComputeStress() {
	Real E = 3.;
	setParam("E", E);

	Matrix<Real> eps = {{2}};
	Matrix<Real> sigma(1, 1);
	Real sigma_th = 2;
	this->computeStressOnQuad(eps, sigma, sigma_th);

	auto solution = E * eps(0, 0) + sigma_th;
	ASSERT_NEAR(sigma(0, 0), solution, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElastic<1>>::testEnergyDensity() {
	Real eps = 2, sigma = 2;
	Real epot = 0;
	this->computePotentialEnergyOnQuad({{eps}}, {{sigma}}, epot);
	Real solution = 2;
	ASSERT_NEAR(epot, solution, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <>
	void FriendMaterial<MaterialElastic<1>>::testComputeTangentModuli() {
	Real E = 2;
	setParam("E", E);
	Matrix<Real> tangent(1, 1);
	this->computeTangentModuliOnQuad(tangent);
	ASSERT_NEAR(tangent(0, 0), E, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElastic<1>>::testPushWaveSpeed() {
	Real E = 3., rho = 2.;
	setParam("E", E);
	setParam("rho", rho);

	auto wave_speed = this->getPushWaveSpeed(Element());
	auto solution = std::sqrt(E / rho);
	ASSERT_NEAR(wave_speed, solution, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElastic<1>>::testShearWaveSpeed() {
	ASSERT_THROW(this->getShearWaveSpeed(Element()), debug::Exception);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElastic<2>>::testComputeStress() {
	Real E = 1.;
	Real nu = .3;
	Real sigma_th = 0.3; // thermal stress
	setParam("E", E);
	setParam("nu", nu);

	Matrix<Real> rotation_matrix = getRandomRotation2d();

	auto grad_u = this->getDeviatoricStrain(1.).block(0, 0, 2, 2);

	auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix);

	Matrix<Real> sigma_rot(2, 2);
	this->computeStressOnQuad(grad_u_rot, sigma_rot, sigma_th);

	auto sigma = this->reverseRotation(sigma_rot, rotation_matrix);

	Matrix<Real> identity(2, 2);
	identity.eye();

	Matrix<Real> sigma_expected =
	0.5 * E / (1 + nu) * (grad_u + grad_u.transpose()) + sigma_th * identity;

	auto diff = sigma - sigma_expected;
	Real stress_error = diff.norm<L_inf>();
	ASSERT_NEAR(stress_error, 0., 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElastic<2>>::testEnergyDensity() {
	Matrix<Real> sigma = {{1, 2}, {2, 4}};
	Matrix<Real> eps = {{1, 0}, {0, 1}};
	Real epot = 0;
	Real solution = 2.5;
	this->computePotentialEnergyOnQuad(eps, sigma, epot);
	ASSERT_NEAR(epot, solution, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <>
	void FriendMaterial<MaterialElastic<2>>::testComputeTangentModuli() {
	Real E = 1.;
	Real nu = .3;
	setParam("E", E);
	setParam("nu", nu);
	Matrix<Real> tangent(3, 3);

	/* Plane Strain */
	// clang-format off
	Matrix<Real> solution = {
	{1 - nu, nu, 0},
	{nu, 1 - nu, 0},
	{0, 0, (1 - 2 * nu) / 2},
	};
	// clang-format on
	solution = E / ((1 + nu) (1 - 2 * nu));

	this->computeTangentModuliOnQuad(tangent);
	Real tangent_error = (tangent - solution).norm<L_2>();
	ASSERT_NEAR(tangent_error, 0, 1e-14);

	/* Plane Stress */
	this->plane_stress = true;
	this->updateInternalParameters();
	// clang-format off
	solution = {
	{1, nu, 0},
	{nu, 1, 0},
	{0, 0, (1 - nu) / 2},
	};
	// clang-format on
	solution = E / (1 - nu nu);

	this->computeTangentModuliOnQuad(tangent);
	tangent_error = (tangent - solution).norm<L_2>();
	ASSERT_NEAR(tangent_error, 0, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElastic<2>>::testPushWaveSpeed() {
	Real E = 1.;
	Real nu = .3;
	Real rho = 2;
	setParam("E", E);
	setParam("nu", nu);
	setParam("rho", rho);
	auto wave_speed = this->getPushWaveSpeed(Element());

	Real K = E / (3 * (1 - 2 * nu));
	Real mu = E / (2 * (1 + nu));
	Real sol = std::sqrt((K + 4. / 3 * mu) / rho);

	ASSERT_NEAR(wave_speed, sol, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElastic<2>>::testShearWaveSpeed() {
	Real E = 1.;
	Real nu = .3;
	Real rho = 2;
	setParam("E", E);
	setParam("nu", nu);
	setParam("rho", rho);
	auto wave_speed = this->getShearWaveSpeed(Element());

	Real mu = E / (2 * (1 + nu));
	Real sol = std::sqrt(mu / rho);

	ASSERT_NEAR(wave_speed, sol, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElastic<3>>::testComputeStress() {
	Real E = 1.;
	Real nu = .3;
	Real sigma_th = 0.3; // thermal stress
	setParam("E", E);
	setParam("nu", nu);

	Matrix<Real> rotation_matrix = getRandomRotation3d();

	auto grad_u = this->getDeviatoricStrain(1.);

	auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix);

	Matrix<Real> sigma_rot(3, 3);
	this->computeStressOnQuad(grad_u_rot, sigma_rot, sigma_th);

	auto sigma = this->reverseRotation(sigma_rot, rotation_matrix);

	Matrix<Real> identity(3, 3);
	identity.eye();

	Matrix<Real> sigma_expected =
	0.5 * E / (1 + nu) * (grad_u + grad_u.transpose()) + sigma_th * identity;

	auto diff = sigma - sigma_expected;
	Real stress_error = diff.norm<L_inf>();
	ASSERT_NEAR(stress_error, 0., 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElastic<3>>::testEnergyDensity() {
	Matrix<Real> sigma = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}};
	Matrix<Real> eps = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
	Real epot = 0;
	Real solution = 5.5;
	this->computePotentialEnergyOnQuad(eps, sigma, epot);
	ASSERT_NEAR(epot, solution, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <>
	void FriendMaterial<MaterialElastic<3>>::testComputeTangentModuli() {
	Real E = 1.;
	Real nu = .3;
	setParam("E", E);
	setParam("nu", nu);
	Matrix<Real> tangent(6, 6);

	// clang-format off
	Matrix<Real> solution = {
	{1 - nu, nu, nu, 0, 0, 0},
	{nu, 1 - nu, nu, 0, 0, 0},
	{nu, nu, 1 - nu, 0, 0, 0},
	{0, 0, 0, (1 - 2 * nu) / 2, 0, 0},
	{0, 0, 0, 0, (1 - 2 * nu) / 2, 0},
	{0, 0, 0, 0, 0, (1 - 2 * nu) / 2},
	};
	// clang-format on
	solution = E / ((1 + nu) (1 - 2 * nu));

	this->computeTangentModuliOnQuad(tangent);
	Real tangent_error = (tangent - solution).norm<L_2>();
	ASSERT_NEAR(tangent_error, 0, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElastic<3>>::testPushWaveSpeed() {
	Real E = 1.;
	Real nu = .3;
	Real rho = 2;
	setParam("E", E);
	setParam("nu", nu);
	setParam("rho", rho);
	auto wave_speed = this->getPushWaveSpeed(Element());

	Real K = E / (3 * (1 - 2 * nu));
	Real mu = E / (2 * (1 + nu));
	Real sol = std::sqrt((K + 4. / 3 * mu) / rho);

	ASSERT_NEAR(wave_speed, sol, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElastic<3>>::testShearWaveSpeed() {
	Real E = 1.;
	Real nu = .3;
	Real rho = 2;
	setParam("E", E);
	setParam("nu", nu);
	setParam("rho", rho);
	auto wave_speed = this->getShearWaveSpeed(Element());

	Real mu = E / (2 * (1 + nu));
	Real sol = std::sqrt(mu / rho);

	ASSERT_NEAR(wave_speed, sol, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticOrthotropic<2>>::testComputeStress() {
	UInt Dim = 2;
	Real E1 = 1.;
	Real E2 = 2.;
	Real nu12 = 0.1;
	Real G12 = 2.;

	setParamNoUpdate("E1", E1);
	setParamNoUpdate("E2", E2);
	setParamNoUpdate("nu12", nu12);
	setParamNoUpdate("G12", G12);

	// material frame of reference is rotate by rotation_matrix starting from
	// canonical basis
	Matrix<Real> rotation_matrix = getRandomRotation2d();

	// canonical basis as expressed in the material frame of reference, as
	// required by MaterialElasticOrthotropic class (it is simply given by the
	// columns of the rotation_matrix; the lines give the material basis expressed
	// in the canonical frame of reference)
	Vector<Real> v1(Dim);
	Vector<Real> v2(Dim);
	for (UInt i = 0; i < Dim; ++i) {
	v1[i] = rotation_matrix(i, 0);
	v2[i] = rotation_matrix(i, 1);
	}

	setParamNoUpdate("n1", v1.normalize());
	setParamNoUpdate("n2", v2.normalize());

	// set internal Cijkl matrix expressed in the canonical frame of reference
	this->updateInternalParameters();

	// gradient in material frame of reference
	auto grad_u = ( this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.) ).block(0, 0, 2, 2);

	// gradient in canonical basis (we need to rotate back to the canonical
	// basis)
	auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix);

	// stress in the canonical basis
	Matrix<Real> sigma_rot(2, 2);
	this->computeStressOnQuad(grad_u_rot, sigma_rot);

	// stress in the material reference (we need to apply the rotation)
	auto sigma = this->applyRotation(sigma_rot, rotation_matrix);

	// construction of Cijkl engineering tensor in the material frame of reference
	// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
	Real nu21 = nu12*E2/E1;
	Real gamma = 1 / ( 1 - nu12*nu21 );

	Matrix<Real> C_expected(2Dim,2Dim,0);
	C_expected(0,0) = gamma * E1 ;
	C_expected(1,1) = gamma * E2 ;
	C_expected(2,2) = G12;

	C_expected(1,0) = C_expected(0,1) = gamma * E1 * nu21;

	// epsilon is computed directly in the material frame of reference
	Matrix<Real> epsilon = 0.5 * ( grad_u + grad_u.transpose() );

	// sigma_expected is computed directly in the material frame of reference
	Matrix<Real> sigma_expected(Dim, Dim);
	for ( UInt i = 0; i < Dim; ++i ) {
	for ( UInt j = 0; j < Dim; ++j ) {
	sigma_expected(i,i) += C_expected(i,j)*epsilon(j,j);
	}
	}

	sigma_expected(0,1) = sigma_expected(1,0) = C_expected(2,2) * 2 * epsilon(0,1);

	// sigmas are checked in the material frame of reference
	auto diff = sigma - sigma_expected;
	Real stress_error = diff.norm<L_inf>();
	ASSERT_NEAR(stress_error, 0., 1e-13);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticOrthotropic<2>>::testEnergyDensity() {
	Matrix<Real> sigma = {{1, 2}, {2, 4}};
	Matrix<Real> eps = {{1, 0}, {0, 1}};
	Real epot = 0;
	Real solution = 2.5;
	this->computePotentialEnergyOnQuad(eps, sigma, epot);
	ASSERT_NEAR(epot, solution, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticOrthotropic<2>>::testComputeTangentModuli() {

	// Note: for this test material and canonical basis coincide
	Vector<Real> n1 = {1, 0};
	Vector<Real> n2 = {0, 1};
	Real E1 = 1.;
	Real E2 = 2.;
	Real nu12 = 0.1;
	Real G12 = 2.;

	setParamNoUpdate("n1", n1);
	setParamNoUpdate("n2", n2);
	setParamNoUpdate("E1", E1);
	setParamNoUpdate("E2", E2);
	setParamNoUpdate("nu12", nu12);
	setParamNoUpdate("G12", G12);

	// set internal Cijkl matrix expressed in the canonical frame of reference
	this->updateInternalParameters();

	// construction of Cijkl engineering tensor in the material frame of reference
	// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
	Real nu21 = nu12*E2/E1;
	Real gamma = 1 / ( 1 - nu12*nu21 );

	Matrix<Real> C_expected(3,3);
	C_expected(0,0) = gamma * E1 ;
	C_expected(1,1) = gamma * E2 ;
	C_expected(2,2) = G12;

	C_expected(1,0) = C_expected(0,1) = gamma * E1 * nu21;

	Matrix<Real> tangent(3,3);
	this->computeTangentModuliOnQuad(tangent);

	Real tangent_error = (tangent - C_expected).norm<L_2>();
	ASSERT_NEAR(tangent_error, 0, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticOrthotropic<2>>::testPushWaveSpeed() {

	// Note: for this test material and canonical basis coincide
	Vector<Real> n1 = {1, 0};
	Vector<Real> n2 = {0, 1};
	Real E1 = 1.;
	Real E2 = 2.;
	Real nu12 = 0.1;
	Real G12 = 2.;
	Real rho = 2.5;

	setParamNoUpdate("n1", n1);
	setParamNoUpdate("n2", n2);
	setParamNoUpdate("E1", E1);
	setParamNoUpdate("E2", E2);
	setParamNoUpdate("nu12", nu12);
	setParamNoUpdate("G12", G12);
	setParamNoUpdate("rho", rho);

	// set internal Cijkl matrix expressed in the canonical frame of reference
	this->updateInternalParameters();

	// construction of Cijkl engineering tensor in the material frame of reference
	// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
	Real nu21 = nu12*E2/E1;
	Real gamma = 1 / ( 1 - nu12*nu21 );

	Matrix<Real> C_expected(3,3);
	C_expected(0,0) = gamma * E1 ;
	C_expected(1,1) = gamma * E2 ;
	C_expected(2,2) = G12;

	C_expected(1,0) = C_expected(0,1) = gamma * E1 * nu21;

	Vector<Real> eig_expected(3);
	C_expected.eig(eig_expected);

	auto celerity_expected = std::sqrt(eig_expected(0)/rho);

	auto celerity = this->getCelerity(Element());

	ASSERT_NEAR(celerity_expected, celerity, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticOrthotropic<2>>::testShearWaveSpeed() {

	AKANTU_DEBUG_WARNING("WARNING: Shear wave speed for 2D orthotropic materials is not present in akantu."
	"Celerity is used instead to get the stable time step (see testPushWaveSpeed).");
	ASSERT_NEAR(0, 0, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticOrthotropic<3>>::testComputeStress() {
	UInt Dim = 3;
	Real E1 = 1.;
	Real E2 = 2.;
	Real E3 = 3.;
	Real nu12 = 0.1;
	Real nu13 = 0.2;
	Real nu23 = 0.3;
	Real G12 = 2.;
	Real G13 = 3.;
	Real G23 = 1.;

	setParamNoUpdate("E1", E1);
	setParamNoUpdate("E2", E2);
	setParamNoUpdate("E3", E3);
	setParamNoUpdate("nu12", nu12);
	setParamNoUpdate("nu13", nu13);
	setParamNoUpdate("nu23", nu23);
	setParamNoUpdate("G12", G12);
	setParamNoUpdate("G13", G13);
	setParamNoUpdate("G23", G23);

	// material frame of reference is rotate by rotation_matrix starting from
	// canonical basis
	Matrix<Real> rotation_matrix = getRandomRotation3d();

	// canonical basis as expressed in the material frame of reference, as
	// required by MaterialElasticOrthotropic class (it is simply given by the
	// columns of the rotation_matrix; the lines give the material basis expressed
	// in the canonical frame of reference)
	Vector<Real> v1(Dim);
	Vector<Real> v2(Dim);
	Vector<Real> v3(Dim);
	for (UInt i = 0; i < Dim; ++i) {
	v1[i] = rotation_matrix(i, 0);
	v2[i] = rotation_matrix(i, 1);
	v3[i] = rotation_matrix(i, 2);
	}

	setParamNoUpdate("n1", v1.normalize());
	setParamNoUpdate("n2", v2.normalize());
	setParamNoUpdate("n3", v3.normalize());

	// set internal Cijkl matrix expressed in the canonical frame of reference
	this->updateInternalParameters();

	// gradient in material frame of reference
	auto grad_u = this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.);

	// gradient in canonical basis (we need to rotate back to the canonical
	// basis)
	auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix);

	// stress in the canonical basis
	Matrix<Real> sigma_rot(3, 3);
	this->computeStressOnQuad(grad_u_rot, sigma_rot);

	// stress in the material reference (we need to apply the rotation)
	auto sigma = this->applyRotation(sigma_rot, rotation_matrix);

	// construction of Cijkl engineering tensor in the material frame of reference
	// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
	Real nu21 = nu12*E2/E1;
	Real nu31 = nu13*E3/E1;
	Real nu32 = nu23*E3/E2;
	Real gamma = 1 / ( 1 - nu12nu21 - nu23nu32 - nu31nu13 - 2nu21nu32nu13 );

	Matrix<Real> C_expected(6,6);
	C_expected(0,0) = gamma * E1 * ( 1-nu23*nu32 );
	C_expected(1,1) = gamma * E2 * ( 1-nu13*nu31 );
	C_expected(2,2) = gamma * E3 * ( 1-nu12*nu21 );

	C_expected(1,0) = C_expected(0,1) = gamma * E1 * ( nu21 + nu31*nu23 );
	C_expected(2,0) = C_expected(0,2) = gamma * E1 * ( nu31 + nu21*nu32 );
	C_expected(2,1) = C_expected(1,2) = gamma * E2 * ( nu32 + nu12*nu31 );

	C_expected(3,3) = G23;
	C_expected(4,4) = G13;
	C_expected(5,5) = G12;

	// epsilon is computed directly in the material frame of reference
	Matrix<Real> epsilon = 0.5 * ( grad_u + grad_u.transpose() );

	// sigma_expected is computed directly in the material frame of reference
	Matrix<Real> sigma_expected(Dim, Dim);
	for ( UInt i = 0; i < Dim; ++i ) {
	for ( UInt j = 0; j < Dim; ++j ) {
	sigma_expected(i,i) += C_expected(i,j)*epsilon(j,j);
	}
	}

	sigma_expected(0,1) = C_expected(5,5) * 2 * epsilon(0,1);
	sigma_expected(0,2) = C_expected(4,4) * 2 * epsilon(0,2);
	sigma_expected(1,2) = C_expected(3,3) * 2 * epsilon(1,2);
	sigma_expected(1,0) = sigma_expected(0,1);
	sigma_expected(2,0) = sigma_expected(0,2);
	sigma_expected(2,1) = sigma_expected(1,2);

	// sigmas are checked in the material frame of reference
	auto diff = sigma - sigma_expected;
	Real stress_error = diff.norm<L_inf>();
	ASSERT_NEAR(stress_error, 0., 1e-13);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticOrthotropic<3>>::testEnergyDensity() {
	Matrix<Real> sigma = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}};
	Matrix<Real> eps = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
	Real epot = 0;
	Real solution = 5.5;
	this->computePotentialEnergyOnQuad(eps, sigma, epot);
	ASSERT_NEAR(epot, solution, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticOrthotropic<3>>::testComputeTangentModuli() {

	// Note: for this test material and canonical basis coincide
	UInt Dim = 3;
	Vector<Real> n1 = {1, 0, 0};
	Vector<Real> n2 = {0, 1, 0};
	Vector<Real> n3 = {0, 0, 1};
	Real E1 = 1.;
	Real E2 = 2.;
	Real E3 = 3.;
	Real nu12 = 0.1;
	Real nu13 = 0.2;
	Real nu23 = 0.3;
	Real G12 = 2.;
	Real G13 = 3.;
	Real G23 = 1.;

	setParamNoUpdate("n1", n1);
	setParamNoUpdate("n2", n2);
	setParamNoUpdate("n3", n3);
	setParamNoUpdate("E1", E1);
	setParamNoUpdate("E2", E2);
	setParamNoUpdate("E3", E3);
	setParamNoUpdate("nu12", nu12);
	setParamNoUpdate("nu13", nu13);
	setParamNoUpdate("nu23", nu23);
	setParamNoUpdate("G12", G12);
	setParamNoUpdate("G13", G13);
	setParamNoUpdate("G23", G23);

	// set internal Cijkl matrix expressed in the canonical frame of reference
	this->updateInternalParameters();

	// construction of Cijkl engineering tensor in the material frame of reference
	// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
	Real nu21 = nu12*E2/E1;
	Real nu31 = nu13*E3/E1;
	Real nu32 = nu23*E3/E2;
	Real gamma = 1 / ( 1 - nu12nu21 - nu23nu32 - nu31nu13 - 2nu21nu32nu13 );

	Matrix<Real> C_expected(2Dim,2Dim,0);
	C_expected(0,0) = gamma * E1 * ( 1-nu23*nu32 );
	C_expected(1,1) = gamma * E2 * ( 1-nu13*nu31 );
	C_expected(2,2) = gamma * E3 * ( 1-nu12*nu21 );

	C_expected(1,0) = C_expected(0,1) = gamma * E1 * ( nu21 + nu31*nu23 );
	C_expected(2,0) = C_expected(0,2) = gamma * E1 * ( nu31 + nu21*nu32 );
	C_expected(2,1) = C_expected(1,2) = gamma * E2 * ( nu32 + nu12*nu31 );

	C_expected(3,3) = G23;
	C_expected(4,4) = G13;
	C_expected(5,5) = G12;

	Matrix<Real> tangent(6,6);
	this->computeTangentModuliOnQuad(tangent);

	Real tangent_error = (tangent - C_expected).norm<L_2>();
	ASSERT_NEAR(tangent_error, 0, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticOrthotropic<3>>::testPushWaveSpeed() {

	// Note: for this test material and canonical basis coincide
	UInt Dim = 3;
	Vector<Real> n1 = {1, 0, 0};
	Vector<Real> n2 = {0, 1, 0};
	Vector<Real> n3 = {0, 0, 1};
	Real E1 = 1.;
	Real E2 = 2.;
	Real E3 = 3.;
	Real nu12 = 0.1;
	Real nu13 = 0.2;
	Real nu23 = 0.3;
	Real G12 = 2.;
	Real G13 = 3.;
	Real G23 = 1.;
	Real rho = 2.3;

	setParamNoUpdate("n1", n1);
	setParamNoUpdate("n2", n2);
	setParamNoUpdate("n3", n3);
	setParamNoUpdate("E1", E1);
	setParamNoUpdate("E2", E2);
	setParamNoUpdate("E3", E3);
	setParamNoUpdate("nu12", nu12);
	setParamNoUpdate("nu13", nu13);
	setParamNoUpdate("nu23", nu23);
	setParamNoUpdate("G12", G12);
	setParamNoUpdate("G13", G13);
	setParamNoUpdate("G23", G23);
	setParamNoUpdate("rho", rho);

	// set internal Cijkl matrix expressed in the canonical frame of reference
	this->updateInternalParameters();

	// construction of Cijkl engineering tensor in the material frame of reference
	// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
	Real nu21 = nu12*E2/E1;
	Real nu31 = nu13*E3/E1;
	Real nu32 = nu23*E3/E2;
	Real gamma = 1 / ( 1 - nu12nu21 - nu23nu32 - nu31nu13 - 2nu21nu32nu13 );

	Matrix<Real> C_expected(2Dim,2Dim,0);
	C_expected(0,0) = gamma * E1 * ( 1-nu23*nu32 );
	C_expected(1,1) = gamma * E2 * ( 1-nu13*nu31 );
	C_expected(2,2) = gamma * E3 * ( 1-nu12*nu21 );

	C_expected(1,0) = C_expected(0,1) = gamma * E1 * ( nu21 + nu31*nu23 );
	C_expected(2,0) = C_expected(0,2) = gamma * E1 * ( nu31 + nu21*nu32 );
	C_expected(2,1) = C_expected(1,2) = gamma * E2 * ( nu32 + nu12*nu31 );

	C_expected(3,3) = G23;
	C_expected(4,4) = G13;
	C_expected(5,5) = G12;

	Vector<Real> eig_expected(6);
	C_expected.eig(eig_expected);

	auto celerity_expected = std::sqrt(eig_expected(0)/rho);

	auto celerity = this->getCelerity(Element());

	ASSERT_NEAR(celerity_expected, celerity, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticOrthotropic<3>>::testShearWaveSpeed() {

	AKANTU_DEBUG_WARNING("WARNING: Shear wave speed for 3D orthotropic materials is not present in akantu."
	"Celerity is used instead to get the stable time step (see testPushWaveSpeed).");
	ASSERT_NEAR(0, 0, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticLinearAnisotropic<2>>::testComputeStress() {
	UInt Dim = 2;
	Matrix<Real> C = {
	{1.0, 0.3, 0.4},
	{0.3, 2.0, 0.1},
	{0.4, 0.1, 1.5},
	};

	setParamNoUpdate("C11", C(0,0));
	setParamNoUpdate("C12", C(0,1));
	setParamNoUpdate("C13", C(0,2));
	setParamNoUpdate("C22", C(1,1));
	setParamNoUpdate("C23", C(1,2));
	setParamNoUpdate("C33", C(2,2));

	// material frame of reference is rotate by rotation_matrix starting from
	// canonical basis
	Matrix<Real> rotation_matrix = getRandomRotation2d();

	// canonical basis as expressed in the material frame of reference, as
	// required by MaterialElasticLinearAnisotropic class (it is simply given by
	// the columns of the rotation_matrix; the lines give the material basis
	// expressed in the canonical frame of reference)
	Vector<Real> v1(Dim);
	Vector<Real> v2(Dim);
	for (UInt i = 0; i < Dim; ++i) {
	v1[i] = rotation_matrix(i, 0);
	v2[i] = rotation_matrix(i, 1);
	}

	setParamNoUpdate("n1", v1.normalize());
	setParamNoUpdate("n2", v2.normalize());

	// set internal Cijkl matrix expressed in the canonical frame of reference
	this->updateInternalParameters();

	// gradient in material frame of reference
	auto grad_u = (this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.)).block(0, 0, 2, 2);

	// gradient in canonical basis (we need to rotate back to the canonical
	// basis)
	auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix);

	// stress in the canonical basis
	Matrix<Real> sigma_rot(2, 2);
	this->computeStressOnQuad(grad_u_rot, sigma_rot);

	// stress in the material reference (we need to apply the rotation)
	auto sigma = this->applyRotation(sigma_rot, rotation_matrix);

	// epsilon is computed directly in the material frame of reference
	Matrix<Real> epsilon = 0.5 * ( grad_u + grad_u.transpose() );
	Vector<Real> epsilon_voigt(3);
	epsilon_voigt(0) = epsilon(0,0);
	epsilon_voigt(1) = epsilon(1,1);
	epsilon_voigt(2) = 2 * epsilon(0,1);

	// sigma_expected is computed directly in the material frame of reference
	Vector<Real> sigma_voigt = C * epsilon_voigt;
	Matrix<Real> sigma_expected(Dim, Dim);
	sigma_expected(0,0) = sigma_voigt(0);
	sigma_expected(1,1) = sigma_voigt(1);
	sigma_expected(0,1) = sigma_expected(1,0) = sigma_voigt(2);

	// sigmas are checked in the material frame of reference
	auto diff = sigma - sigma_expected;
	Real stress_error = diff.norm<L_inf>();
	ASSERT_NEAR(stress_error, 0., 1e-13);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticLinearAnisotropic<2>>::testEnergyDensity() {
	Matrix<Real> sigma = {{1, 2}, {2, 4}};
	Matrix<Real> eps = {{1, 0}, {0, 1}};
	Real epot = 0;
	Real solution = 2.5;
	this->computePotentialEnergyOnQuad(eps, sigma, epot);
	ASSERT_NEAR(epot, solution, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticLinearAnisotropic<2>>::testComputeTangentModuli() {

	// Note: for this test material and canonical basis coincide
	Matrix<Real> C = {
	{1.0, 0.3, 0.4},
	{0.3, 2.0, 0.1},
	{0.4, 0.1, 1.5},
	};

	setParamNoUpdate("C11", C(0,0));
	setParamNoUpdate("C12", C(0,1));
	setParamNoUpdate("C13", C(0,2));
	setParamNoUpdate("C22", C(1,1));
	setParamNoUpdate("C23", C(1,2));
	setParamNoUpdate("C33", C(2,2));

	// set internal Cijkl matrix expressed in the canonical frame of reference
	this->updateInternalParameters();

	Matrix<Real> tangent(3,3);
	this->computeTangentModuliOnQuad(tangent);

	Real tangent_error = (tangent - C).norm<L_2>();
	ASSERT_NEAR(tangent_error, 0, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticLinearAnisotropic<2>>::testPushWaveSpeed() {

	// Note: for this test material and canonical basis coincide
	Matrix<Real> C = {
	{1.0, 0.3, 0.4},
	{0.3, 2.0, 0.1},
	{0.4, 0.1, 1.5},
	};

	Real rho = 2.7;

	setParamNoUpdate("C11", C(0,0));
	setParamNoUpdate("C12", C(0,1));
	setParamNoUpdate("C13", C(0,2));
	setParamNoUpdate("C22", C(1,1));
	setParamNoUpdate("C23", C(1,2));
	setParamNoUpdate("C33", C(2,2));
	setParamNoUpdate("rho", rho);

	// set internal Cijkl matrix expressed in the canonical frame of reference
	this->updateInternalParameters();

	Vector<Real> eig_expected(3);
	C.eig(eig_expected);

	auto celerity_expected = std::sqrt(eig_expected(0)/rho);

	auto celerity = this->getCelerity(Element());

	ASSERT_NEAR(celerity_expected, celerity, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticLinearAnisotropic<2>>::testShearWaveSpeed() {

	AKANTU_DEBUG_WARNING("WARNING: Shear wave speed for 2D anisotropic materials is not present in akantu."
	"Celerity is used instead to get the stable time step (see testPushWaveSpeed).");
	ASSERT_NEAR(0, 0, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticLinearAnisotropic<3>>::testComputeStress() {
	UInt Dim = 3;
	Matrix<Real> C = {
	{1.0, 0.3, 0.4, 0.3, 0.2, 0.1},
	{0.3, 2.0, 0.1, 0.2, 0.3, 0.2},
	{0.4, 0.1, 1.5, 0.1, 0.4, 0.3},
	{0.3, 0.2, 0.1, 2.4, 0.1, 0.4},
	{0.2, 0.3, 0.4, 0.1, 0.9, 0.1},
	{0.1, 0.2, 0.3, 0.4, 0.1, 1.2},
	};

	setParamNoUpdate("C11", C(0,0));
	setParamNoUpdate("C12", C(0,1));
	setParamNoUpdate("C13", C(0,2));
	setParamNoUpdate("C14", C(0,3));
	setParamNoUpdate("C15", C(0,4));
	setParamNoUpdate("C16", C(0,5));
	setParamNoUpdate("C22", C(1,1));
	setParamNoUpdate("C23", C(1,2));
	setParamNoUpdate("C24", C(1,3));
	setParamNoUpdate("C25", C(1,4));
	setParamNoUpdate("C26", C(1,5));
	setParamNoUpdate("C33", C(2,2));
	setParamNoUpdate("C34", C(2,3));
	setParamNoUpdate("C35", C(2,4));
	setParamNoUpdate("C36", C(2,5));
	setParamNoUpdate("C44", C(3,3));
	setParamNoUpdate("C45", C(3,4));
	setParamNoUpdate("C46", C(3,5));
	setParamNoUpdate("C55", C(4,4));
	setParamNoUpdate("C56", C(4,5));
	setParamNoUpdate("C66", C(5,5));

	// material frame of reference is rotate by rotation_matrix starting from
	// canonical basis
	Matrix<Real> rotation_matrix = getRandomRotation3d();

	// canonical basis as expressed in the material frame of reference, as
	// required by MaterialElasticLinearAnisotropic class (it is simply given by
	// the columns of the rotation_matrix; the lines give the material basis
	// expressed in the canonical frame of reference)
	Vector<Real> v1(Dim);
	Vector<Real> v2(Dim);
	Vector<Real> v3(Dim);
	for (UInt i = 0; i < Dim; ++i) {
	v1[i] = rotation_matrix(i, 0);
	v2[i] = rotation_matrix(i, 1);
	v3[i] = rotation_matrix(i, 2);
	}

	setParamNoUpdate("n1", v1.normalize());
	setParamNoUpdate("n2", v2.normalize());
	setParamNoUpdate("n3", v3.normalize());

	// set internal Cijkl matrix expressed in the canonical frame of reference
	this->updateInternalParameters();

	// gradient in material frame of reference
	auto grad_u = this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.);

	// gradient in canonical basis (we need to rotate back to the canonical
	// basis)
	auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix);

	// stress in the canonical basis
	Matrix<Real> sigma_rot(3, 3);
	this->computeStressOnQuad(grad_u_rot, sigma_rot);

	// stress in the material reference (we need to apply the rotation)
	auto sigma = this->applyRotation(sigma_rot, rotation_matrix);

	// epsilon is computed directly in the material frame of reference
	Matrix<Real> epsilon = 0.5 * ( grad_u + grad_u.transpose() );
	Vector<Real> epsilon_voigt(6);
	epsilon_voigt(0) = epsilon(0,0);
	epsilon_voigt(1) = epsilon(1,1);
	epsilon_voigt(2) = epsilon(2,2);
	epsilon_voigt(3) = 2 * epsilon(1,2);
	epsilon_voigt(4) = 2 * epsilon(0,2);
	epsilon_voigt(5) = 2 * epsilon(0,1);

	// sigma_expected is computed directly in the material frame of reference
	Vector<Real> sigma_voigt = C * epsilon_voigt;
	Matrix<Real> sigma_expected(Dim, Dim);
	sigma_expected(0,0) = sigma_voigt(0);
	sigma_expected(1,1) = sigma_voigt(1);
	sigma_expected(2,2) = sigma_voigt(2);
	sigma_expected(1,2) = sigma_expected(2,1) = sigma_voigt(3);
	sigma_expected(0,2) = sigma_expected(2,0) = sigma_voigt(4);
	sigma_expected(0,1) = sigma_expected(1,0) = sigma_voigt(5);

	// sigmas are checked in the material frame of reference
	auto diff = sigma - sigma_expected;
	Real stress_error = diff.norm<L_inf>();
	ASSERT_NEAR(stress_error, 0., 1e-13);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticLinearAnisotropic<3>>::testEnergyDensity() {
	Matrix<Real> sigma = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}};
	Matrix<Real> eps = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
	Real epot = 0;
	Real solution = 5.5;
	this->computePotentialEnergyOnQuad(eps, sigma, epot);
	ASSERT_NEAR(epot, solution, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticLinearAnisotropic<3>>::testComputeTangentModuli() {

	// Note: for this test material and canonical basis coincide
	Matrix<Real> C = {
	{1.0, 0.3, 0.4, 0.3, 0.2, 0.1},
	{0.3, 2.0, 0.1, 0.2, 0.3, 0.2},
	{0.4, 0.1, 1.5, 0.1, 0.4, 0.3},
	{0.3, 0.2, 0.1, 2.4, 0.1, 0.4},
	{0.2, 0.3, 0.4, 0.1, 0.9, 0.1},
	{0.1, 0.2, 0.3, 0.4, 0.1, 1.2},
	};

	setParamNoUpdate("C11", C(0,0));
	setParamNoUpdate("C12", C(0,1));
	setParamNoUpdate("C13", C(0,2));
	setParamNoUpdate("C14", C(0,3));
	setParamNoUpdate("C15", C(0,4));
	setParamNoUpdate("C16", C(0,5));
	setParamNoUpdate("C22", C(1,1));
	setParamNoUpdate("C23", C(1,2));
	setParamNoUpdate("C24", C(1,3));
	setParamNoUpdate("C25", C(1,4));
	setParamNoUpdate("C26", C(1,5));
	setParamNoUpdate("C33", C(2,2));
	setParamNoUpdate("C34", C(2,3));
	setParamNoUpdate("C35", C(2,4));
	setParamNoUpdate("C36", C(2,5));
	setParamNoUpdate("C44", C(3,3));
	setParamNoUpdate("C45", C(3,4));
	setParamNoUpdate("C46", C(3,5));
	setParamNoUpdate("C55", C(4,4));
	setParamNoUpdate("C56", C(4,5));
	setParamNoUpdate("C66", C(5,5));

	// set internal Cijkl matrix expressed in the canonical frame of reference
	this->updateInternalParameters();

	Matrix<Real> tangent(6,6);
	this->computeTangentModuliOnQuad(tangent);

	Real tangent_error = (tangent - C).norm<L_2>();
	ASSERT_NEAR(tangent_error, 0, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticLinearAnisotropic<3>>::testPushWaveSpeed() {

	// Note: for this test material and canonical basis coincide
	Matrix<Real> C = {
	{1.0, 0.3, 0.4, 0.3, 0.2, 0.1},
	{0.3, 2.0, 0.1, 0.2, 0.3, 0.2},
	{0.4, 0.1, 1.5, 0.1, 0.4, 0.3},
	{0.3, 0.2, 0.1, 2.4, 0.1, 0.4},
	{0.2, 0.3, 0.4, 0.1, 0.9, 0.1},
	{0.1, 0.2, 0.3, 0.4, 0.1, 1.2},
	};
	Real rho = 2.9;

	setParamNoUpdate("C11", C(0,0));
	setParamNoUpdate("C12", C(0,1));
	setParamNoUpdate("C13", C(0,2));
	setParamNoUpdate("C14", C(0,3));
	setParamNoUpdate("C15", C(0,4));
	setParamNoUpdate("C16", C(0,5));
	setParamNoUpdate("C22", C(1,1));
	setParamNoUpdate("C23", C(1,2));
	setParamNoUpdate("C24", C(1,3));
	setParamNoUpdate("C25", C(1,4));
	setParamNoUpdate("C26", C(1,5));
	setParamNoUpdate("C33", C(2,2));
	setParamNoUpdate("C34", C(2,3));
	setParamNoUpdate("C35", C(2,4));
	setParamNoUpdate("C36", C(2,5));
	setParamNoUpdate("C44", C(3,3));
	setParamNoUpdate("C45", C(3,4));
	setParamNoUpdate("C46", C(3,5));
	setParamNoUpdate("C55", C(4,4));
	setParamNoUpdate("C56", C(4,5));
	setParamNoUpdate("C66", C(5,5));
	setParamNoUpdate("rho", rho);

	// set internal Cijkl matrix expressed in the canonical frame of reference
	this->updateInternalParameters();

	Vector<Real> eig_expected(6);
	C.eig(eig_expected);

	auto celerity_expected = std::sqrt(eig_expected(0)/rho);

	auto celerity = this->getCelerity(Element());

	ASSERT_NEAR(celerity_expected, celerity, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	template <> void FriendMaterial<MaterialElasticLinearAnisotropic<3>>::testShearWaveSpeed() {

	AKANTU_DEBUG_WARNING("WARNING: Shear wave speed for 3D anisotropic materials is not present in akantu."
	"Celerity is used instead to get the stable time step (see testPushWaveSpeed).");
	ASSERT_NEAR(0, 0, 1e-14);
	}

	/* -------------------------------------------------------------------------- */
	namespace {

	template <typename T>
	class TestElasticMaterialFixture : public ::TestMaterialFixture<T> {};

	TYPED_TEST_CASE(TestElasticMaterialFixture, types);

	TYPED_TEST(TestElasticMaterialFixture, ElasticComputeStress) {
	this->material->testComputeStress();
	}

	TYPED_TEST(TestElasticMaterialFixture, ElasticEnergyDensity) {
	this->material->testEnergyDensity();
	}

	TYPED_TEST(TestElasticMaterialFixture, ElasticComputeTangentModuli) {
	this->material->testComputeTangentModuli();
	}

	TYPED_TEST(TestElasticMaterialFixture, ElasticComputePushWaveSpeed) {
	this->material->testPushWaveSpeed();
	}

	TYPED_TEST(TestElasticMaterialFixture, ElasticComputeShearWaveSpeed) {
	this->material->testShearWaveSpeed();
	}
	} // namespace

test_elastic_materials.ccNo OneTemporaryActions

File Metadata

test_elastic_materials.ccView Options

Event Timeline

test_elastic_materials.cc
No OneTemporary
Actions

test_elastic_materials.cc
View Options