diff --git a/src/materials/laminate_homogenisation.cc b/src/materials/laminate_homogenisation.cc
index 786ce7f..3d79aac 100644
--- a/src/materials/laminate_homogenisation.cc
+++ b/src/materials/laminate_homogenisation.cc
@@ -1,571 +1,559 @@
/**
* @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,
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);
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 6 elements obtained from
// A matrices averaging routine:
Stiffness_t ret_K = Stiffness_t::Zero();
ret_K = MaterialAnisotropic<twoD>::c_maker(c_maker_inp);
return ret_K ;
}
/*------------------------------------------------------------------*/
template<>
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,2,2,2,2), cf[0]*get(K,2,2,1,1), cf[1]*get(K,2,2,1,2),
cf[0]*get(K,2,2,1,1), cf[0]*get(K,1,1,1,1), cf[1]*get(K,1,1,1,2),
cf[1]*get(K,2,2,1,2), cf[1]*get(K,1,1,1,2), cf[2]*get(K,2,1,2,1);
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),
A12(0,2)/A12c(0,2), A11(1,0)/A11c(1,0), A11(0,2)/A11c(0,2),
A22(1,1)/A22c(1,1), A22(1,0)/A22c(1,0), A22(2,1)/A22c(2,1),
A12(1,0)/A12c(1,0), A12(2,0)/A12c(2,0),
A22(0,0)/A22c(0,0), A22(2,0)/A22c(2,0),
A12(1,1)/A12c(1,1), A12(2,1)/A12c(2,1),
A22(2,2)/A22c(2,2), A12(1,2)/A12c(1,2), A12(2,2)/A12c(2,2),
A11(1,1)/A11c(1,1), A11(2,1)/A11c(2,1),
A11(2,2)/A11c(2,2)};
// now the resultant stiffness is calculated fro 21 elements obtained from
// A matrices averaging routine :
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};
-
+ Strain_t F_1{F_coord},F_2{F_coord};
+ Parallel_strain_t F_para{}; F_para = get_parallel_strain(F_coord);
+ Serial_strain_t F_1_new_seri{get_serial_strain(F_coord)};
+ Serial_strain_t F_1_seri, F_2_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:
+ StressMap_t P_1{}, P_2{}, ret_P{};
+ Stiffness_t K_1{}, K_2{}, ret_K{};
+
do {
// updating variables:
F_1_seri = F_1_new_seri;
- F_2_seri = linear_eqs_seri(ratio, F_lam, F_1_seri);
+ F_2_seri = linear_eqs_seri(ratio, F_coord, 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()),
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);
+ 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);
+ 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(ret_P,ret_K);
- /*rotator.rotate_back(ret_P),
- rotator.rotate_back(ret_K));*/
+ return std::make_tuple(rotator.rotate_back(ret_P),
+ rotator.rotate_back(ret_K));
+ //(ret_P,ret_K);
}
/* ---------------------------------------------------------------------- */
template struct LamHomogen<twoD>;
template struct LamHomogen<threeD>;
} //muSpectre
+