diff --git a/src/materials/laminate_homogenisation.cc b/src/materials/laminate_homogenisation.cc
index 82e76d2..48808f2 100644
--- a/src/materials/laminate_homogenisation.cc
+++ b/src/materials/laminate_homogenisation.cc
@@ -1,555 +1,569 @@
/**
* @file lamminate_hemogenization.hh
*
* @author Ali Falsafi <ali.falsafi@epfl.ch>
*
* @date 28 Sep 2018
*
* @brief
*
* Copyright © 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre 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
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "common/geometry.hh"
#include "common/common.hh"
#include "laminate_homogenisation.hh"
#include "material_anisotropic.hh"
#include <tuple>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template<>
constexpr auto LamHomogen<threeD>::get_serial_indexes()
->std::array<Dim_t, 2*threeD-1>{
constexpr Dim_t Dim {3};
return std::array<Dim_t, 2*Dim-1>{0, 1, 2, 3, 6};
}
template<>
constexpr auto LamHomogen<twoD>::get_serial_indexes()
->std::array<Dim_t, 2*twoD-1>{
constexpr Dim_t Dim {2};
return std::array<Dim_t, 2*Dim-1>{0, 1, 2};
}
/* ---------------------------------------------------------------------- */
template<>
constexpr auto LamHomogen<threeD>::get_parallel_indexes()
->std::array<Dim_t, (threeD-1)*(threeD-1)>{
constexpr Dim_t Dim {3};
return std::array<Dim_t, (Dim-1)*(Dim-1)>{4, 5, 7, 8};;
}
template<>
constexpr auto LamHomogen<twoD>::get_parallel_indexes()
->std::array<Dim_t, (twoD-1)*(twoD-1)>{
constexpr Dim_t Dim {2};
return std::array<Dim_t, (Dim-1)*(Dim-1)>{3};
}
/* ---------------------------------------------------------------------- */
template<Dim_t Dim>
auto LamHomogen<Dim>::
get_serial_strain(const Eigen::Ref<Strain_t>& Fs)->Serial_strain_t{
Serial_strain_t ret_F;
auto serial_indexes{get_serial_indexes()};
for (int i{0}; i < ret_F.size(); ++i){
ret_F(i) = Fs(serial_indexes[i]);
}
return ret_F;
}
/* ---------------------------------------------------------------------- */
template<Dim_t Dim>
auto LamHomogen<Dim>::
get_parallel_strain(const Eigen::Ref<Strain_t>& Fs)->Parallel_strain_t{
Parallel_strain_t ret_F;
auto parallel_indexes(get_parallel_indexes());
for (int i{0}; i < ret_F.size(); ++i){
ret_F(i) = Fs(parallel_indexes[i]);
}
return ret_F;
}
/* ---------------------------------------------------------------------- */
template<Dim_t Dim>
auto LamHomogen<Dim>::
get_total_strain(const Eigen::Ref<Serial_strain_t>& F_seri,
const Eigen::Ref<Parallel_strain_t>& F_para)
->Strain_t{
Strain_t ret_F;
auto serial_indexes{get_serial_indexes()};
auto parallel_indexes(get_parallel_indexes());
for (unsigned int i{0}; i < serial_indexes.size(); ++i){
ret_F(serial_indexes[i]) = F_seri(i);
}
for (unsigned int i{0}; i < parallel_indexes.size(); ++i){
ret_F(parallel_indexes[i]) = F_para(i);
}
return ret_F;
}
/* ---------------------------------------------------------------------- */
template<Dim_t Dim>
auto LamHomogen<Dim>::
get_serial_stiffness(const Eigen::Ref<Stiffness_t>& K)
->Serial_stiffness_t{
Serial_stiffness_t ret_K{};
auto serial_indexes{get_serial_indexes()};
for (int i{0}; i < (2*Dim-1); ++i){
for (int j{0}; j < (2*Dim-1); ++j){
ret_K(i,j) = K(serial_indexes[i], serial_indexes[j]);
}
}
return ret_K;
}
/* ---------------------------------------------------------------------- */
template<Dim_t Dim>
auto LamHomogen<Dim>::
get_parallel_stiffness(const Eigen::Ref<Stiffness_t>& K)
->Parallel_stiffness_t{
Parallel_stiffness_t ret_K{};
auto parallel_indexes{get_parallel_indexes()};
for (int i{0}; i < (Dim-1)*(Dim-1); ++i){
for (int j{0}; j < (Dim-1)*(Dim-1); ++j){
ret_K(i,j) = K(parallel_indexes[i], parallel_indexes[j]);
}
}
return ret_K;
}
/* ---------------------------------------------------------------------- */
template<Dim_t Dim>
auto LamHomogen<Dim>::
get_total_stiffness(const Eigen::Ref<Serial_stiffness_t>& K_seri,
const Eigen::Ref<Parallel_stiffness_t>& K_para)
->Stiffness_t{
Stiffness_t ret_K{};
auto serial_indexes{get_serial_indexes()};
auto parallel_indexes{get_parallel_indexes()};
for (int i{0}; i < (Dim-1)*(Dim-1); ++i){
for (int j{0}; j < (Dim-1)*(Dim-1); ++j){
ret_K(parallel_indexes[i],parallel_indexes[j]) = K_para(i, j);
}
}
for (int i{0}; i < (Dim-1)*(Dim-1); ++i){
for (int j{0}; j < (Dim-1)*(Dim-1); ++j){
ret_K(serial_indexes[i],serial_indexes[j]) = K_seri(i, j);
}
}
return ret_K;
}
/* ---------------------------------------------------------------------- */
// It does not to be especialised because all functions used init are
// especialised. this two functions compute the strai in the second
// layer of laminate material from strain in the first layer and total strain
template<Dim_t Dim>
auto LamHomogen<Dim>::linear_eqs (Real ratio,
const Eigen::Ref<Strain_t>& F_lam,
const Eigen::Ref<Strain_t>& F_1)
->Strain_t{
Strain_t ret_F;
Serial_strain_t F_lam_seri = get_serial_strain(F_lam);
Serial_strain_t F_1_seri = get_serial_strain(F_1);
Serial_strain_t F_2_seri{};
for (int i{0}; i < 2 * Dim - 1; i++){
F_2_seri(i) = (F_lam_seri(i) - ratio * F_1_seri(i)) / (1-ratio);
}
Parallel_strain_t F_lam_para = get_parallel_strain(F_lam);
Parallel_strain_t F_1_para = get_parallel_strain(F_1);
Parallel_strain_t F_2_para{};
for (int i{0}; i < 2 * (Dim-1)*(Dim-1); i++){
F_2_para(i) = (F_lam_para(i) - ratio * F_1_para(i)) / (1-ratio);
}
Strain_t F_2 = get_total_strain(F_2_seri, F_2_para);
return F_2;
}
template<Dim_t Dim>
auto LamHomogen<Dim>::
linear_eqs_seri (Real ratio,
const Eigen::Ref<Strain_t> & F_lam,
const Eigen::Ref<Serial_strain_t> & F_1)
->Serial_strain_t{
Serial_strain_t F_2_seri;
Serial_strain_t F_lam_seri = get_serial_strain(F_lam);
for (int i{0}; i < 2 * Dim - 1; i++){
F_2_seri(i) = (F_lam_seri(i) - ratio * F_1(i)) / (1-ratio);
}
return F_2_seri;
}
/* ---------------------------------------------------------------------- */
template<Dim_t Dim>
auto LamHomogen<Dim>::
delta_stress_eval (const Function_t &mat_1_stress_eval,
- const Function_t &mat_2_stress_eval,
- const Eigen::Ref<StrainMap_t>& F_1,
- const Eigen::Ref<StrainMap_t>& F_2)
+ const Function_t &mat_2_stress_eval,
+ const Eigen::Ref<StrainMap_t>& F_1,
+ const Eigen::Ref<StrainMap_t>& F_2,
+ RotatorNormal<Dim> rotator)
->std::tuple<Serial_stress_t, Serial_stiffness_t>{
auto P_K_1 = mat_1_stress_eval(F_1);
auto P_K_2 = mat_2_stress_eval(F_2);
StressMap_t del_P_map; del_P_map = std::get<0>(P_K_1) - std::get<0>(P_K_2);
-
- auto del_P_seri = get_serial_strain(Eigen::Map<Stress_t>(del_P_map.data()));
+ del_P_map = rotator.rotate(del_P_map);
+ auto del_P_seri = get_serial_strain(Eigen::Map<Stress_t>
+ (del_P_map.data())
+ );
Stiffness_t del_K; del_K = std::get<1>(P_K_1) - std::get<1>(P_K_2);
+ del_K = rotator.rotate(del_K);
auto del_K_seri = get_serial_stiffness(del_K);
return std::make_tuple(del_P_seri, del_K_seri);
}
/* ---------------------------------------------------------------------- */
template<Dim_t Dim>
inline auto LamHomogen<Dim>::
del_energy_eval (const Real del_F_norm,
const Real delta_P_norm)
->Real{
return del_F_norm * delta_P_norm;
}
-
+ /*------------------------------------------------------------------------*/
/* ---------------------------------------------------------------------- */
template<>
auto LamHomogen<twoD>::
lam_K_combine(const Eigen::Ref<Stiffness_t>& K_1,
const Eigen::Ref<Stiffness_t>& K_2,
const Real ratio)
->Stiffness_t{
using Mat_A_t = Eigen::Matrix<Real, twoD, twoD>;
using Vec_A_t = Eigen::Matrix<Real, twoD, 1>;
using Vec_AT_t = Eigen::Matrix<Real, 1, twoD>;
std::array<double, 3> cf = {1.0, sqrt(2.0), 2.0};
auto get_A11 = [cf](const Eigen::Ref<Stiffness_t>& K){
Mat_A_t A11 = Mat_A_t::Zero() ;
A11<<
cf[0]*get(K,0,0,0,0), cf[1]*get(K,0,0,0,1),
cf[1]*get(K,0,0,0,1), cf[2]*get(K,0,1,0,1);
return A11;
};
Mat_A_t A11c;
A11c<<
cf[0], cf[1],
cf[1], cf[2];
auto get_A12 = [cf](const Eigen::Ref<Stiffness_t>& K){
Vec_A_t A12 = Vec_A_t::Zero() ;
A12<<
cf[0]*get(K,0,0,1,1),
cf[1]*get(K,1,1,0,1);
return A12;
};
Vec_A_t A12c = Vec_A_t::Zero();
A12c <<
cf[0],
cf[1];
Mat_A_t A11_1 {get_A11(K_1)};
Mat_A_t A11_2 {get_A11(K_2)};
Vec_A_t A12_1 = {get_A12(K_1)}; Vec_AT_t A21_1 = A12_1.transpose();
Vec_A_t A12_2 = {get_A12(K_2)}; Vec_AT_t A21_2 = A12_2.transpose();
Real A22_1 = get(K_1,1,1,1,1);
Real A22_2 = get(K_2,1,1,1,1);
auto get_inverse_average =
[&ratio](const Eigen::Ref<Mat_A_t> &matrix_1,
const Eigen::Ref<Mat_A_t> &matrix_2){
return (( ratio * matrix_1.inverse() +
(1-ratio)* matrix_2.inverse()).inverse());
};
auto get_average =
[&ratio](Real A_1,
Real A_2){
return ratio*A_1 + (1-ratio)*A_2;
};
auto get_average_vec =
[&ratio](Vec_A_t A_1,
Vec_A_t A_2){
return ratio*A_1 + (1-ratio)*A_2;};
auto get_average_vecT =
[&ratio](Vec_AT_t A_1,
Vec_AT_t A_2){
return ratio*A_1 + (1-ratio)*A_2;};
// calculating average of A matrices of the materials
Mat_A_t A11 = get_inverse_average(A11_1, A11_2);
Vec_A_t A12 = A11 * get_average_vec(A11_1.inverse()*A12_1,
A11_2.inverse()*A12_2);
Real A22 = get_average(A22_1 - A21_1 * A11_1.inverse() * A12_1,
A22_2 - A21_2 * A11_2.inverse() * A12_2) +
get_average_vecT(A21_1 * A11_1.inverse(),
A21_2 * A11_2.inverse()) *
A11 * get_average_vec(A11_1.inverse() * A12_1,
A11_2.inverse() * A12_2);
std::vector<Real> c_maker_inp =
{A11(0,0)/A11c(0,0),A12(0,0)/A12c(0,0),A11(0,1)/A11c(0,1),
A22 ,A12(1,0)/A12c(1,0),
A11(1,1)/A11c(1,1)};
// now the resultant stiffness is calculated fro 21 elements obtained from
// A matrices :
Stiffness_t ret_K = Stiffness_t::Zero();
ret_K = MaterialAnisotropic<twoD>::c_maker(c_maker_inp);
return ret_K ;
- return K_1 + K_2 * ratio;
}
+
+ /*------------------------------------------------------------------*/
template<>
- auto LamHomogen<threeD>::
- lam_K_combine(const Eigen::Ref<Stiffness_t>& K_1,
- const Eigen::Ref<Stiffness_t>& K_2,
- const Real ratio)
+ auto LamHomogen<threeD>::
+ lam_K_combine(const Eigen::Ref<Stiffness_t>& K_1,
+ const Eigen::Ref<Stiffness_t>& K_2,
+ const Real ratio)
->Stiffness_t{
// the combination method is obtained form P. 163 of
// "Theory of Composites"
// Author : Milton_G_W
//constructing "A" matrices( A11, A12, A21, A22) according to the
//procedure from the book:
// this type of matrix will be used in calculating the combinatio of the
// Stiffness matrixes
using Mat_A_t = Eigen::Matrix<Real, threeD, threeD>;
// these coeffs are used in constructing A matrices from k and vice versa
std::array<double, 3> cf = {1.0, sqrt(2.0), 2.0};
auto get_A11 = [cf](const Eigen::Ref<Stiffness_t>& K){
Mat_A_t A11 = Mat_A_t::Zero() ;
A11<<
cf[0]*get(K,0,0,0,0), cf[1]*get(K,0,0,0,2), cf[1]*get(K,0,0,0,1),
cf[1]*get(K,0,0,0,2), cf[2]*get(K,0,2,0,2), cf[2]*get(K,0,2,0,1),
cf[1]*get(K,0,0,0,1), cf[2]*get(K,0,2,0,1), cf[2]*get(K,0,1,0,1);
return A11;
};
Mat_A_t A11c;
A11c<<
cf[0], cf[1], cf[1],
cf[1], cf[2], cf[2],
cf[1], cf[2], cf[2];
auto get_A12 = [cf](const Eigen::Ref<Stiffness_t>& K){
Mat_A_t A12 = Mat_A_t::Zero();
A12 <<
cf[0]*get(K,0,0,1,1), cf[0]*get(K,0,0,2,2), cf[1]*get(K,0,0,1,2),
cf[1]*get(K,1,1,0,2), cf[1]*get(K,2,2,0,2), cf[2]*get(K,1,2,0,2),
cf[1]*get(K,1,1,0,1), cf[1]*get(K,2,2,0,1), cf[2]*get(K,1,2,0,1);
return A12;
};
Mat_A_t A12c;
A12c<<
cf[0], cf[0], cf[1],
cf[1], cf[1], cf[2],
cf[1], cf[1], cf[2];
auto get_A22 = [cf](const Eigen::Ref<Stiffness_t>& K){
Mat_A_t A22 = Mat_A_t::Zero();
A22 <<
cf[0]*get(K,1,1,1,1), cf[0]*get(K,1,1,2,2), cf[1]*get(K,1,1,1,2),
cf[0]*get(K,1,1,2,2), cf[0]*get(K,2,2,2,2), cf[1]*get(K,2,2,1,2),
cf[1]*get(K,1,1,1,2), cf[1]*get(K,2,2,1,2), cf[2]*get(K,1,2,1,2);
return A22;
};
Mat_A_t A22c;
A22c <<
cf[0], cf[0], cf[1],
cf[0], cf[0], cf[1],
cf[1], cf[1], cf[2];
Mat_A_t A11_1 {get_A11(K_1)};
Mat_A_t A11_2 {get_A11(K_2)};
Mat_A_t A12_1 {get_A12(K_1)}; Mat_A_t A21_1 {A12_1.transpose()};
Mat_A_t A12_2 {get_A12(K_2)}; Mat_A_t A21_2 {A12_2.transpose()};
Mat_A_t A22_1 {get_A22(K_1)};
Mat_A_t A22_2 {get_A22(K_2)};
// these two functions are routines to compute average of "A" matrices
auto get_inverse_average =
[&ratio](const Eigen::Ref<Mat_A_t> &matrix_1,
const Eigen::Ref<Mat_A_t> &matrix_2){
return (( ratio * matrix_1.inverse() +
(1-ratio)* matrix_2.inverse()).inverse());
};
auto get_average =
[&ratio](const Mat_A_t &matrix_1,
const Mat_A_t &matrix_2){
return (ratio * matrix_1 + (1-ratio) * matrix_2);
};
// calculating average of A matrices of the materials
Mat_A_t A11 = get_inverse_average(A11_1, A11_2);
Mat_A_t A12 = A11 * get_average(A11_1.inverse()*A12_1,
A11_2.inverse()*A12_2);
Mat_A_t A22 = get_average(A22_1 - A21_1 * A11_1.inverse() * A12_1,
A22_2 - A21_2 * A11_2.inverse() * A12_2) +
get_average(A21_1 * A11_1.inverse(),
A21_2 * A11_2.inverse()) *
A11 * get_average(A11_1.inverse() * A12_1,
A11_2.inverse() * A12_2);
std::vector<Real> c_maker_inp =
{A11(0,0)/A11c(0,0), A12(0,0)/A12c(0,0), A12(0,1)/A12c(0,1), A11(0,2)/A11c(0,2)
, A11(1,0)/A11c(1,0), A12(0,2)/A12c(0,2),
A22(0,0)/A22c(0,0), A22(0,1)/A22c(0,1), A12(2,0)/A12c(2,0), A12(1,0)/A12c(1,0)
, A22(0,2)/A22c(0,2),
A22(1,1)/A22c(1,1), A12(2,1)/A12c(2,1), A12(1,1)/A12c(1,1), A22(1,2)/A22c(1,2),
A11(2,2)/A11c(2,2), A11(1,2)/A11c(1,2), A12(2,2)/A12c(2,2),
A11(1,1)/A11c(1,1), A12(1,2)/A11c(1,2),
A22(2,2)/A22c(2,2)};
// now the resultant stiffness is calculated fro 21 elements obtained from
// A matrices :
Stiffness_t ret_K = Stiffness_t::Zero();
ret_K = MaterialAnisotropic<threeD>::c_maker(c_maker_inp);
return ret_K ;
}
/* ---------------------------------------------------------------------- */
template<Dim_t Dim>
auto LamHomogen<Dim>::
lam_P_combine(Eigen::Ref<StressMap_t> P_1,
Eigen::Ref<StressMap_t> P_2,
const Real ratio)
->StressMap_t{
StressMap_t P_tmp{};
//return P_tmp ;
Serial_stress_t P_1_seri =
get_serial_strain(Eigen::Map<Stress_t>(P_1.data()));
Parallel_stress_t P_1_para =
get_parallel_strain(Eigen::Map<Stress_t>(P_1.data()));
+
Parallel_stress_t P_2_para =
get_parallel_strain(Eigen::Map<Stress_t>(P_2.data()));
Parallel_stress_t P_para = ratio * P_1_para + (1-ratio)* P_2_para ;
auto ret_P = get_total_strain(P_1_seri, P_para);
return Eigen::Map<StressMap_t> (ret_P.data());
}
/* ---------------------------------------------------------------------- */
template<Dim_t Dim>
auto LamHomogen<Dim>::
laminate_solver(Eigen::Ref<Strain_t> F_coord,
Function_t mat_1_stress_eval,
Function_t mat_2_stress_eval,
Real ratio,
const Eigen::Ref<Vec_t>& normal_vec,
Real tol,
Dim_t max_iter)
->std::tuple<StressMap_t, Stiffness_t>{
/**
* here we rotate the strain such that the laminate intersection normal
* would align with the x-axis F_lam is the total strain in the new
* coordinates.
*/
Real del_energy;
Dim_t iter{0};
RotatorNormal<Dim> rotator(normal_vec);
StrainMap_t F_rot = rotator.rotate(Eigen::Map<StrainMap_t>(F_coord.data()));
Eigen::Map<Strain_t> F_lam(F_rot.data());
+ F_lam = F_coord;
Strain_t F_1{F_lam},F_2{F_lam};
Parallel_strain_t F_para{}; F_para = get_parallel_strain
(Eigen::Map<Strain_t>(F_lam.data()));
StressMap_t P_1{}, P_2{}, ret_P{};
Stiffness_t K_1{}, K_2{}, ret_K{};
Serial_strain_t F_1_new_seri ; F_1_new_seri = get_serial_strain(F_lam);
Serial_strain_t F_1_seri{F_1_new_seri}, F_2_seri{F_1_new_seri};
std::tuple<Serial_stress_t, Serial_stiffness_t> delta_P_K_seri;
Serial_stress_t delta_P_seri{};
Serial_stiffness_t delta_K_seri{};
std::tuple<StressMap_t, Stiffness_t> P_K_1, P_K_2;
// TODO: to make sure that the loop is not performed if not necessary;
// maybe change it to while loop:
do {
// updating variables:
F_1_seri = F_1_new_seri;
F_2_seri = linear_eqs_seri(ratio, F_lam, F_1_seri);
F_1 = get_total_strain(F_1_seri, F_para);
F_2 = get_total_strain(F_2_seri, F_para);
/**
* P_1 - P_2 and its Jacobian is calculated here
*/
delta_P_K_seri = delta_stress_eval(mat_1_stress_eval,
mat_2_stress_eval,
Eigen::Map<StrainMap_t> (F_1.data()),
- Eigen::Map<StrainMap_t> (F_2.data()));
+ Eigen::Map<StrainMap_t> (F_2.data()),
+ rotator);
delta_P_seri = std::get<0> (delta_P_K_seri);
delta_K_seri = std::get<1> (delta_P_K_seri);
//solving for ∇P(x_n+1 - x_n)=-P
F_1_new_seri = F_1_seri -
delta_K_seri.colPivHouseholderQr().solve(delta_P_seri);
// calculating the criterion for ending the loop
del_energy = del_energy_eval((F_1_new_seri - F_1_seri).norm(),
delta_P_seri.norm());
iter ++;
} while(del_energy > tol && iter < max_iter);
// storing the reuktant strain in each layer in F_1 and F_2
F_1 = get_total_strain(F_1_seri, F_para);
F_2 = get_total_strain(F_2_seri, F_para);
// computing stress and stiffness of each layer according to it's
// correspondent strin calculated in the loop
P_K_1 = mat_1_stress_eval(Eigen::Map<StrainMap_t> (F_1.data()));
P_K_2 = mat_2_stress_eval(Eigen::Map<StrainMap_t> (F_2.data()));
P_1 = std::get<0>(P_K_1); P_2 = std::get<0>(P_K_2);
K_1 = std::get<1>(P_K_1); K_2 = std::get<1>(P_K_2);
+
// combine the computed strains and tangents to have the resultant
// stress and tangent of the pixel
+ //P_1 = rotator.rotate(P_1);
+ //P_2 = rotator.rotate(P_2);
ret_P = lam_P_combine(P_1, P_2, ratio);
+ //K_1 = rotator.rotate(K_1);
+ //K_2 = rotator.rotate(K_2);
ret_K = lam_K_combine(K_1, K_2, ratio);
//ret_P = Eigen::Map<Stress_t> (P_rot.data());
//ret_K = rotator.rotate_back(ret_K);
- return std::make_tuple(rotator.rotate_back(ret_P),
- rotator.rotate_back(ret_K));
+ return std::make_tuple(ret_P,ret_K);
+ /*rotator.rotate_back(ret_P),
+ rotator.rotate_back(ret_K));*/
}
/* ---------------------------------------------------------------------- */
- template struct LamHomogen<twoD>;
+//template struct LamHomogen<twoD>;
template struct LamHomogen<threeD>;
} //muSpectre
diff --git a/src/materials/laminate_homogenisation.hh b/src/materials/laminate_homogenisation.hh
index 4231286..d1173ad 100644
--- a/src/materials/laminate_homogenisation.hh
+++ b/src/materials/laminate_homogenisation.hh
@@ -1,192 +1,193 @@
/**
* @file lamminate_hemogenization.hh
*
* @author Ali Falsafi <ali.falsafi@epfl.ch>
*
* @date 28 Sep 2018
*
* @brief Laminatehomogenisation enables one to obtain the resulting stress
* and stiffness tensors of a laminate pixel that is consisted of two
* materialswith a certain normal vector of their interface plane.
* note that it is supposed to be used in static way. so it does note
* any data member. It is merely a collection of functions used to
* calculate effective stress and stiffness.
*
*
* Copyright © 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre 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
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef LAMINATE_HOMOGENISATION_H
#define LAMINATE_HOMOGENISATION_H
#include "common/geometry.hh"
#include "common/common.hh"
#include "common/field_map.hh"
namespace muSpectre{
// TODO: to template with the formulation as well
template<Dim_t Dim>
struct LamHomogen{
//! typedefs for data handled by this interface
using Vec_t = Eigen::Matrix<Real, Dim, 1>;
using Stiffness_t = T4Mat<Real, Dim>;
using Serial_stiffness_t = Eigen::Matrix<Real, 2*Dim-1, 2*Dim-1>;
using Parallel_stiffness_t = T4Mat<Real, (Dim-1)>;
using Strain_t = Eigen::Matrix<Real, Dim*Dim, 1>;
using Serial_strain_t = Eigen::Matrix<Real, 2*Dim-1, 1>;
using Parallel_strain_t = Eigen::Matrix<Real, (Dim-1)*(Dim-1), 1>;
using Stress_t = Strain_t;
using Serial_stress_t = Serial_strain_t;
using Parallel_stress_t = Parallel_strain_t;
using StrainMap_t = Eigen::Matrix<Real, Dim, Dim>;
using StressMap_t = StrainMap_t;
using Function_t = std::function<std::tuple<StressMap_t,
Stiffness_t>
(const Eigen::Ref<StrainMap_t> &)>;
// these two functions return the indexes in column major strain and stress
// tensors that the behavior is either serial or parallel concerning
// combining the laminate layers
static constexpr auto get_serial_indexes()->std::array<Dim_t, 2*Dim-1>;
static constexpr auto get_parallel_indexes()->std::array<Dim_t, (Dim-1)*(Dim-1)>;
static constexpr auto get_serial_elimination_indexes()
->std::array<Dim_t, (Dim-1)*(Dim-1)>;
static constexpr auto get_parallel_elimination_indexes()
->std::array<Dim_t, 2*Dim-1>;
static constexpr auto get_if_serial()->std::array<bool, Dim*Dim>;
// these functions return the parts of a stress or strain tensor that
// behave either serial or parallel in a laminate structure
// they are used to obtain a certian part of tensor which is needed by the
// solver used in this struct to calculate the effective stress and stiffness
// or compose the complete tensor from its sub sets.
// obtain the serial part of stress or strain tensor
static auto get_serial_strain(const Eigen::Ref<Strain_t>& F) -> Serial_strain_t;
// obtain the parallel part of stress or strain tensor
static auto get_parallel_strain(const Eigen::Ref<Strain_t>& F)
-> Parallel_strain_t;
// compose the complete strain or stress tensor from
//its serial and parallel parts
static auto get_total_strain(const Eigen::Ref<Serial_strain_t>& F_seri,
const Eigen::Ref<Parallel_strain_t>& F_para)->Strain_t;
// obtain the serial part of stifness tensor
static auto get_serial_stiffness(const Eigen::Ref<Stiffness_t>& K)
->Serial_stiffness_t;
// obtain the parallel part of stifness tensor
static auto get_parallel_stiffness(const Eigen::Ref<Stiffness_t>& K)
->Parallel_stiffness_t;
// compose the complete stiffness tensor from
//its serial and parallel parts
static auto get_total_stiffness(const Eigen::Ref<Serial_stiffness_t>& K_seri,
const Eigen::Ref<Parallel_stiffness_t>& K_para)
->Stiffness_t;
/**
* this functions calculate the strain in the second materials considering
* the fact that:
* 1. In the directions that the laminate layers are in series the weighted
* average of the strains should be equal to the total strain of the laminate
* 2. In the directions that the laminate layers are in parallel the strain
* is constant and equal in all of the layers.
*/
static auto linear_eqs(Real ratio,
const Eigen::Ref<Strain_t>& F_lam,
const Eigen::Ref<Strain_t> & F_1)->Strain_t;
static auto linear_eqs_seri(Real ratio,
const Eigen::Ref<Strain_t> & F_lam,
const Eigen::Ref<Serial_strain_t> & F_1)->Serial_strain_t;
/**
* the objective in hemgenisation of a single laminate pixel is equating the
* stress in the serial directions so the difference of stress between their
* layers should tend to zero. this function return the stress difference and
* the difference of Stiffness matrices which is used as the Jacobian in the
* solution process
*/
static auto delta_stress_eval(const Function_t &mat_1_stress_eval,
const Function_t &mat_2_stress_eval,
const Eigen::Ref<StrainMap_t>& F_1,
- const Eigen::Ref<StrainMap_t>& F_2)
+ const Eigen::Ref<StrainMap_t>& F_2,
+ RotatorNormal<Dim> rotator)
->std::tuple<Serial_stress_t, Serial_stiffness_t>;
/**
* the following function s claculates the energy dimension error of the
* solution. it will beused in each step of the solution to detremine the
* relevant difference that implementation of that step has had on
* convergence to the solution.
*/
inline static auto del_energy_eval(const Real del_F_norm,
const Real delta_P_norm)
->Real;
/**
* This functions calculate the resultant stress and tangent matrices
* according to the computed F_1 and F_2 from the solver.
*
*/
static auto lam_P_combine(Eigen::Ref<StressMap_t> P_1,
Eigen::Ref<StressMap_t> P_2,
const Real ratio)
->StressMap_t;
static auto lam_K_combine(const Eigen::Ref<Stiffness_t>& K_1,
const Eigen::Ref<Stiffness_t>& K_2,
const Real ratio)
->Stiffness_t;
/**
* This is the main solver function that might be called staically from
* an external file. this will return the resultant stress and stiffness
* tensor according to interanl "equilibrium" of the lamiante.
* The inputs are :
* 1- global Strain
* 2- stress calculation function of the layer 1
* 3- stress calculation function of the layer 2
* 4- the ratio of the first material in the laminate sturucture of the pixel
* 5- the normal vector of the interface of two layers
* 6- the tolerance error for the internal solution of the laminate pixel
* 7- the maximum iterations for the internal solution of the laminate pixel
*/
static auto laminate_solver(Eigen::Ref<Strain_t> F_coord,
Function_t mat_1_stress_eval,
Function_t mat_2_stress_eval,
Real ratio,
const Eigen::Ref<Vec_t>& normal_vec,
Real tol = 1e-6,
Dim_t max_iter = 1000)
->std::tuple<StressMap_t, Stiffness_t>;
}; // LamHomogen
} //muSpectre
#endif /* LAMINATE_HOMOGENISATION_H */
diff --git a/tests/split_test_laminate_solver.cc b/tests/split_test_laminate_solver.cc
index f8a79c0..b08236d 100644
--- a/tests/split_test_laminate_solver.cc
+++ b/tests/split_test_laminate_solver.cc
@@ -1,145 +1,175 @@
/**
* @file test_laminate_solver.cc
*
* @author Till Junge <ali.falsafi@epfl.ch>
*
* @date 19 Oct 2018
*
* @brief Tests for the large-strain, laminate homogenisation algorithm
*
* Copyright © 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre 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
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <type_traits>
#include <boost/mpl/list.hpp>
#include <boost/range/combine.hpp>
#include "materials/material_linear_elastic1.hh"
#include "materials/material_linear_elastic2.hh"
+#include "materials/material_orthotropic.hh"
#include "materials/laminate_homogenisation.hh"
#include "tests.hh"
#include "tests/test_goodies.hh"
#include "common/field_collection.hh"
#include "common/iterators.hh"
namespace muSpectre{
BOOST_AUTO_TEST_SUITE(laminate_homogenisation);
template <class Mat1_t, class Mat2_t, Dim_t Dim>
struct MaterialsFixture
{
using Vec_t = Eigen::Matrix<Real, Dim, 1>;
using Stiffness_t = T4Mat<Real, Dim>;
using Serial_stiffness_t = Eigen::Matrix<Real, 2*Dim-1, 2*Dim-1>;
using Parallel_stiffness_t = T4Mat<Real, (Dim-1)>;
using Strain_t = Eigen::Matrix<Real, Dim*Dim, 1>;
using Serial_strain_t = Eigen::Matrix<Real, 2*Dim-1, 1>;
using Parallel_strain_t = Eigen::Matrix<Real, (Dim-1)*(Dim-1), 1>;
using Stress_t = Strain_t;
using Serial_stress_t = Serial_strain_t;
using Paralle_stress_t = Parallel_strain_t;
using StrainMap_t = Eigen::Matrix<Real, Dim, Dim>;
using StressMap_t = StrainMap_t;
using T4MatMap_t = Eigen::Map<T4Mat<Real, Dim>>;
using Function_t = std::function<std::tuple<StressMap_t,
Stiffness_t>
(const Eigen::Ref<StrainMap_t> &)>;
constexpr static double contrast = 1.0 ;
constexpr static Real lambda1{2}, mu1{1.5};
constexpr static Real lambda2{contrast*lambda1}, mu2{contrast*mu1};
constexpr static Real young1{mu1*(3*lambda1 + 2*mu1)/(lambda1 + mu1)};
constexpr static Real young2{mu2*(3*lambda2 + 2*mu2)/(lambda2 + mu2)};
constexpr static Real poisson1{lambda1/(2*(lambda1 + mu1))};
constexpr static Real poisson2{lambda2/(2*(lambda2 + mu2))};
+
+
// constructor :
- MaterialsFixture():mat1("Name1", young1, poisson1),
- mat2("Name2", young2, poisson2)
- {};
+ MaterialsFixture();
// constructor with th e normal vector:
MaterialsFixture(Vec_t normal_vec_inp):mat1("Name1", young1, poisson1),
mat2("Name2", young2, poisson2),
normal_vec{normal_vec_inp}
{};
Mat1_t mat1;
Mat2_t mat2;
Vec_t normal_vec{Vec_t::Random()};
static constexpr Dim_t fix_dim{Dim};
};
+ template<>
+ MaterialsFixture<MaterialLinearElastic1<threeD, threeD>,
+ MaterialLinearElastic1<threeD, threeD>, threeD>:: MaterialsFixture():mat1("Name1", young1, poisson1),
+ mat2("Name2", young2, poisson2){}
+ template<>
+ MaterialsFixture<MaterialLinearElastic1<twoD, twoD>,
+ MaterialLinearElastic1<twoD, twoD>, twoD>:: MaterialsFixture():mat1("Name1", young1, poisson1),
+ mat2("Name2", young2, poisson2){}
+ Real contrast{1.0};
+ std::array<Real,9> ortho_inp1 {1.1,2.1,3.1,4.1,5.1,6.1,7.1,8.1,9.1};
+ std::array<Real,9> ortho_inp2 {contrast*1.1,contrast*2.1,contrast*3.1,contrast*4.1,
+ contrast*5.1,contrast*6.1,contrast*7.1,contrast*8.1,contrast*9.1};
+ template<>
+ MaterialsFixture<MaterialOrthotropic<threeD, threeD>,
+ MaterialOrthotropic<threeD, threeD>, threeD>::
+ MaterialsFixture():mat1("Name1", std::vector<Real> {ortho_inp1.begin(), ortho_inp1.end()}),
+ mat2("Name2", std::vector<Real> {ortho_inp2.begin(), ortho_inp2.end()}){}
+
+
using fix_list = boost::mpl::list
- <MaterialsFixture<MaterialLinearElastic1<threeD, threeD>,
+ <MaterialsFixture<MaterialOrthotropic<threeD, threeD>,
+ MaterialOrthotropic<threeD, threeD>, threeD>>;//,
+ /*MaterialsFixture<MaterialLinearElastic1<threeD, threeD>,
MaterialLinearElastic1<threeD, threeD>, threeD>,
MaterialsFixture<MaterialLinearElastic1<twoD, twoD>,
- MaterialLinearElastic1<twoD, twoD>, twoD>>;
+ MaterialLinearElastic1<twoD, twoD>, twoD>>;*/
BOOST_FIXTURE_TEST_CASE_TEMPLATE(identical_material, Fix, fix_list, Fix) {
auto & mat1{Fix::mat1};
auto & mat2{Fix::mat2};
using Strain_t = typename Fix::Strain_t;
using StrainMap_t = typename Fix::StrainMap_t;
StrainMap_t F;
F = StrainMap_t::Random();
+ /*F<< 1,0,0,
+ 0,0,0,
+ 0,0,0;*/
using Fucntion_t = typename Fix::Function_t;
Fucntion_t mat1_evaluate_stress_func =
[&mat1](const Eigen::Ref<StrainMap_t> & strain){
return mat1.evaluate_stress_tangent(std::move(strain));
};
Fucntion_t mat2_evaluate_stress_func =
[&mat2](const Eigen::Ref<StrainMap_t> & strain){
return mat2.evaluate_stress_tangent(std::move(strain));
};
auto P_K_lam =
LamHomogen<Fix::fix_dim>::laminate_solver(Eigen::Map<Strain_t>(F.data()),
mat1_evaluate_stress_func,
mat2_evaluate_stress_func,
0.5,
Fix::normal_vec,
1e-8,
1e5);
auto P_lam = std::get<0> (P_K_lam);
auto K_lam = std::get<1> (P_K_lam);
auto P_K_ref = mat1_evaluate_stress_func(F);
auto P_ref = std::get<0> (P_K_ref);
auto K_ref = std::get<1> (P_K_ref);
+ //std::cout<< "p_lam:"<< std::endl << P_lam <<std::endl;
+ std::cout<< "k_lam:"<< std::endl << K_lam <<std::endl;
+
+ //std::cout<< "p_ref:"<< std::endl << P_ref <<std::endl;
+ std::cout<< "k_ref:"<< std::endl << K_ref <<std::endl;
auto err_P{(P_lam - P_ref).norm()};
auto err_K{(K_lam - K_ref).norm()};
BOOST_CHECK_LT(err_P, tol);
BOOST_CHECK_LT(err_K, tol);
};
BOOST_AUTO_TEST_SUITE_END();
} //muSpectre