/*
---
written by Clemens Preisinger
29.10.14
---
*/

#include <math.h>
#include <limits>

#include "../../stdafx.h"
#include "../../Ensemble/Superimposition/CSForceIterator.h"
#include "UtilChecker_EC3.h"

#define _USE_MATH_DEFINES

namespace CroSecDesign {

	void
	UtilChecker_EC3_Elastic::copy(const UtilChecker& orig, DCMap* mapDict){
		UtilChecker_EC3_Simple_Elastic::copy(orig, mapDict);
		const UtilChecker_EC3_Elastic* _orig = dynamic_cast<const UtilChecker_EC3_Elastic*>(&orig);
		if (_orig && _orig != this) {
			gamma_m0_ = _orig->gamma_m0_;
			gamma_m1_ = _orig->gamma_m1_;
			util_ = _orig->util_;
			calc_ = _orig->calc_;
		}
	}

	void
	UtilChecker_EC3_Elastic::visit(const Beam3DCroSec& crosec, const Element& elem)
	{

		// if (cs_forces_->getElementAction(&elem).size() == 0) {
			// there are no cross section forces for the element
		//	return;
		// }

		real util_grp = 0;
#ifndef FREEVERSION

		real PI2 = SQR(M_PI);

		const Beam3DCroSec_I* beam_crosec_I = dynamic_cast<const Beam3DCroSec_I*>(&crosec);
		bool is_I_profile = beam_crosec_I != 0 && beam_crosec_I->isIProfile();

		real lk_y = 0;
		real lk_z = 0;
		real lk_LT = 0;

		if (can_buckle()) {
			lk_y = bkl_len_->bklLengthY()[elem.ind()];
			lk_z = bkl_len_->bklLengthZ()[elem.ind()];
			lk_LT = bkl_len_->bklLengthLT()[elem.ind()];
		}
		
		// multiply strength with factor to take account of deformation
		const Material& material = crosec.material();
		real fy = material.ft();
		real eps = material.eps();
		real E = material.E();
		real G = material.Gtr();

		// material with zero resistance is automatically insufficient
		if (fy <= 0 || E <= 0 || G <= 0) {
			utilization_ = std::numeric_limits<real>::max();
			return;
		}

		real A = crosec.A;
		real Cw = crosec.Cw;
		real Iy = crosec.Iyy;
		real Iz = crosec.Izz;
		real It = crosec.Ipp;
		real zg = crosec.zg;

		// determine the resisting cross section forces
		real N_Rd = A * fy;
		if (can_buckle()) N_Rd /= gamma_m1_; else N_Rd /= gamma_m0_;
		real Mtrd = crosec.Wt * fy / sqrt(3.0) / gamma_m0_;
		real Vy_Rd = crosec.Ay * fy / sqrt(3.0) / gamma_m0_;
		real Vz_Rd = crosec.Az * fy / sqrt(3.0) / gamma_m0_;

		real C1 = 1.0;
		real C2 = 1.6;
		real C3 = 1.0;

		// according to ÖNORM B 1993-1-1, 6.7.1 (8)
		real Mcr_pos;
		real Mcr_neg;
		if (is_I_profile && beam_crosec_I->zj_pos() != 0.0) {
			Mcr_pos = C1 * PI2 * E * Iz / SQR(lk_z)*
				(sqrt(SQR(lk_z / lk_LT)*Cw / Iz +
				SQR(lk_z)*G*It / (PI2*E*Iz) +
				SQR(C2*zg - C3*beam_crosec_I->zj_pos())) - (C2*zg - C3*beam_crosec_I->zj_pos()));
			Mcr_neg = C1 * PI2 * E * Iz / SQR(lk_z)*
				(sqrt(SQR(lk_z / lk_LT)*Cw / Iz +
				SQR(lk_z)*G*It / (PI2*E*Iz) +
				SQR(C2*zg - C3*beam_crosec_I->zj_neg())) - (C2*zg - C3*beam_crosec_I->zj_neg()));
		}
		else {
			Mcr_pos = C1 * PI2 * E * Iz / SQR(lk_z)*
				(sqrt(SQR(lk_z / lk_LT)*Cw / Iz +
				SQR(lk_z)*G*It / (PI2*E*Iz) +
				SQR(C2*zg)) - C2*zg);
			Mcr_neg = Mcr_pos;
		}

		// cross section with zero resistance is automatically insufficient
		if (N_Rd == 0 || Mtrd == 0 || Vy_Rd == 0 ||
			Vz_Rd == 0 || Mcr_pos == 0 || Mcr_neg == 0) {
			utilization_ = std::numeric_limits<real>::max();
			return;
		}

		real alpha_y = crosec.alpha_y;
		real alpha_z = crosec.alpha_z;
		real alpha_LT = crosec.alpha_LT;
		real Ncr_y = PI2*E*Iy / SQR(lk_y);
		real Ncr_z = PI2*E*Iz / SQR(lk_z);

		real elem_length = elem.characteristicLength();

		std::vector<real> psi_y_values(cs_forces_->numberOfLoadCases(), 1);
		if ((elem_length >= lk_y) && (elem_length >= lk_LT)) psi_y_values = psi(*cs_forces_, elem, Node::y_r);

		// determine the moment coefficients Cmy, Cmz and CmLT
		// according to EN 1993-1-1 Appendix B, procedure 2, Table B.3
		std::vector<real> Cmy_values(cs_forces_->numberOfLoadCases(), 1);
		std::vector<real> Cmz_values(cs_forces_->numberOfLoadCases(), 1);
		std::vector<real> CmLT_values(cs_forces_->numberOfLoadCases(), 1);
	
		if (elem_length >= lk_y) Cmy_values = CmValue(*cs_forces_, elem, Node::y_r);
		if (elem_length >= lk_z) Cmz_values = CmValue(*cs_forces_, elem, Node::z_r);
		if (elem_length >= lk_LT) CmLT_values = CmValue(*cs_forces_, elem, Node::y_r);

		// determine cross section values for calculating stresses;
		real const min_prop = 1E-8;
		real A_sig = A > min_prop ? A : min_prop;
		real Ay_sig = crosec.Ay > min_prop ? crosec.Ay : min_prop;
		real Az_sig = crosec.Az > min_prop ? crosec.Az : min_prop;
		real Wt_sig = crosec.Wt > min_prop ? crosec.Wt : min_prop;
		real Wy_sig = crosec.Wely > min_prop ? crosec.Wely : min_prop;
		real Wz_sig = crosec.Welz > min_prop ? crosec.Welz : min_prop;

		index_t sample_point_ind = 0;

		SuperImp::const_CSForceIterator cs_iter = cs_forces_->const_csForceIteratorBegin(elem);
		SuperImp::const_CSForceIterator cs_iter_end = cs_forces_->const_csForceIteratorEnd(elem);
		while (cs_iter != cs_iter_end) {
			const VecRealProxy cs_force = *cs_iter;

			index_t lc_ind = cs_iter.lc_ind();
			real psi_y = psi_y_values[lc_ind];
			real Cmy = Cmy_values[lc_ind];
			real Cmz = Cmz_values[lc_ind];
			real CmLT = CmLT_values[lc_ind];

			++cs_iter;

			real N_com = -cs_force[Node::x_t]; // compressive normal force is positive
			real N_com_mod = N_com > 0 ? N_com : 0;
			real Vy = abs(cs_force[Node::y_t]);
			real Vz = abs(cs_force[Node::z_t]);
			real Mt = abs(cs_force[Node::x_r]);
			real My = cs_force[Node::y_r];
			real Mz = cs_force[Node::z_r];

			UtilDetailsEC3 util_element_;

			// calculate stresses
			real tau_max_y = std::abs(Mt / Wt_sig) + std::abs(Vy / Ay_sig);
			real tau_max_z = std::abs(Mt / Wt_sig) + std::abs(Vz / Az_sig);
			util_element_.stresses.tau_max = tau_max_y > tau_max_z ? tau_max_y : tau_max_z;
			real sig_max_M = std::abs(My / Wy_sig) + std::abs(Mz / Wz_sig);
			real sig_max_N = -N_com / A_sig;
			real sig_max_1 = sig_max_N + sig_max_M;
			real sig_max_2 = sig_max_N - sig_max_M;
			util_element_.stresses.sig_xx_max = sig_max_1 > sig_max_2 ? sig_max_1 : sig_max_2;
			util_element_.stresses.sig_xx_min = sig_max_1 < sig_max_2 ? sig_max_1 : sig_max_2;

			// determine the resisting cross section forces
			index_t crosec_class = crosec.croSecClass(N_com, My, Mz, eps);
			real My_Rd = Wy(crosec_class, crosec) * fy;
			if (lk_y != 0) My_Rd /= gamma_m1_; else My_Rd /= gamma_m0_;
			real Mz_Rd = Wz(crosec_class, crosec) * fy;
			if (lk_z != 0) Mz_Rd /= gamma_m1_; else Mz_Rd /= gamma_m0_;

			// cross section with zero resistance is automatically insufficient
			if (My_Rd == 0 || Mz_Rd == 0) {
				utilization_ = std::numeric_limits<real>::max();
				return;
			}

			// account for shear and torsional moment
			// according to EN 1993-1-1, 6.2.7 (9)
			// simplification: largest tau_t_Ed is used for both directions Y and Z
			real tau_t_Ed = Mt / crosec.Wt;
			real Vy_pl_T_Rd = 0;
			real Vz_pl_T_Rd = 0;
			if (crosec.isHollow() == false) {
				// according to EN 1993-1-1, 6.2.7, equation 6.26
				util_element_.shear.fac_Mt = tau_t_Ed / (1.25*(fy / sqrt(3.0)) / gamma_m0_);
				if (util_element_.shear.fac_Mt < 1.0) {
					real h = sqrt(1 - util_element_.shear.fac_Mt);
					Vy_pl_T_Rd = h * Vy_Rd;
					Vz_pl_T_Rd = h * Vz_Rd;
				}
			}
			else {
				// according to EN 1993-1-1, 6.2.7, equation 6.28
				util_element_.shear.fac_Mt = tau_t_Ed / ((fy / sqrt(3.0)) / gamma_m0_);
				if (util_element_.shear.fac_Mt < 1.0) {
					real h = sqrt(1 - util_element_.shear.fac_Mt);
					Vy_pl_T_Rd = h * Vy_Rd;
					Vz_pl_T_Rd = h * Vz_Rd;
				}
			}

			// account for shear and bending moment
			// according to EN 1993-1-1, 6.2.8 (4)
			real rhoY = 1.0;
			real rhoZ = 1.0;
			if (Vy_pl_T_Rd != 0) {
				util_element_.shear.fac_Vy = Vy / Vy_pl_T_Rd;
				if (Vy <= 0.5*Vy_pl_T_Rd) {
					rhoY = 0.0;
				}
				else {
					rhoY = SQR(2 * Vy / Vy_pl_T_Rd - 1);
				}
			}
			else {
				util_element_.shear.fac_Vy = 1000;
			}

			if (Vz_pl_T_Rd != 0) {
				util_element_.shear.fac_Vz = Vz / Vz_pl_T_Rd;
				if (Vz <= 0.5*Vz_pl_T_Rd) {
					rhoZ = 0.0;
				}
				else {
					rhoZ = SQR(2 * Vz / Vz_pl_T_Rd - 1);
				}
			}
			else {
				util_element_.shear.fac_Vz = 1000;
			}

			// simplification: largest rho is used for both directions Y and Z
			real rho = rhoY > rhoZ ? rhoY : rhoZ;
			
			// when set to 1 rho causes N_Rd, My_Rd and Mz_Rd to be zero which leads to util = 3000 which is hard to interpret
			// rho = rho < 1.0 ? rho : 1.0;
			rho = rho < 1.0 ? rho : 0.99;

			// reduction factor according to EN 1993-1-1, 6.2.8 (3), equation 6.29
			real red_fac = 1.0 - rho;
			real fy_red = red_fac * fy;

			real N_Rd_red = N_Rd * red_fac;
			real My_Rd_red = My_Rd * red_fac;
			real Mz_Rd_red = Mz_Rd * red_fac;

			// according to EN 1993-1-1, 6.3.1.2, equation 6.49f
			real lam_y_bar = sqrt(A*fy_red / Ncr_y);
			real chi_y = 1.0;
			if (can_buckle() && lk_y != 0) {
				real phi_y = 0.5*(1.0 + alpha_y*(lam_y_bar - 0.2) + SQR(lam_y_bar));
				chi_y = 1 / (phi_y + sqrt(SQR(phi_y) - SQR(lam_y_bar)));
				chi_y = chi_y < 1 ? chi_y : 1;
			}

			// according to EN 1993-1-1, 6.3.1.2, equation 6.49f
			real lam_z_bar = sqrt(A*fy_red / Ncr_z);
			real chi_z = 1.0;
			if (can_buckle() && lk_z != 0) {
				real phi_z = 0.5*(1.0 + alpha_z*(lam_z_bar - 0.2) + SQR(lam_z_bar));
				chi_z = 1 / (phi_z + sqrt(SQR(phi_z) - SQR(lam_z_bar)));
				chi_z = chi_z < 1 ? chi_z : 1;
			}

			// ideal LT-buckling moment in case of cross sections with one axis of symmetry
			// according to ÖNORM B 1993-1-1, 6.7.1 (8)
			real Mcr;
			if (My > 0) {
				Mcr = Mcr_pos;
			}
			else {
				Mcr = Mcr_neg;
			}

			real chi_LT_mod = 1.0;
			// hollow cross sections are not susceptible to lateral torsional buckling (EN 1993-1-1, 6.3.2.1, (2))
			if (crosec.isHollow() == false) {
				if (can_buckle() && (lk_LT != 0)) {
					real lam_lt_0 = 0.4;
					real beta = 0.75;
					real lam_LT_bar = sqrt(Wy(crosec_class, crosec)*fy_red / Mcr);
					if (is_I_profile) {
						real phi_LT = 0.5*(1 + alpha_LT*(lam_LT_bar - lam_lt_0) + beta*SQR(lam_LT_bar));
						real chi_LT1 = 1 / (phi_LT + sqrt(SQR(phi_LT) - beta*SQR(lam_LT_bar)));
						real chi_LT2 = 1 / SQR(lam_LT_bar);
						real chi_LT = chi_LT1 < chi_LT2 ? chi_LT1 : chi_LT2;
						chi_LT = chi_LT < 1 ? chi_LT : 1;
						// according to EN 1993-1-1, 6.3.2.3, Table 6.6
						real kc = 1 / (1.33 - 0.33*psi_y);
						// according to EN 1993-1-1, 6.3.2.3 (2)
						real f = 1.0 - 0.5*(1.0 - kc)*(1 - 2.0*SQR(lam_LT_bar - 0.8));
						f = f < 1.0 ? f : 1.0;
						// according to EN 1993-1-1, 6.3.2.3, (6.58)
						chi_LT_mod = chi_LT / f;
					}
					else {
						real phi_LT = 0.5*(1 + alpha_LT*(lam_LT_bar - 0.2) + SQR(lam_LT_bar));
						chi_LT_mod = 1 / (phi_LT + sqrt(SQR(phi_LT) - SQR(lam_LT_bar)));
					}
					// according to EN 1993-1-1, 6.3.2.3, (6.58)
					chi_LT_mod = chi_LT_mod < 1.0 ? chi_LT_mod : 1.0;
				}
			}


			real kyy_ = kyy(crosec_class, Cmy, lam_y_bar, N_com_mod, chi_y, N_Rd_red);
			real kzz_ = kzz(crosec_class, Cmy, lam_z_bar, N_com_mod, chi_z, N_Rd_red, is_I_profile);
			real kyz_ = kyz(crosec_class, kzz_);
			real kzy_ = kzy(crosec_class, CmLT, lam_z_bar, N_com_mod, chi_z, N_Rd_red, is_I_profile, kyy_);

			if (N_Rd_red == 0) {
				util_element_.bending_Y.fac_N = 1000;
				util_element_.bending_Z.fac_N = 1000;
			}
			else {
				if (N_com < 0) {
					// no buckling in tension
					util_element_.bending_Y.fac_N = -N_com / (N_Rd_red);
					util_element_.bending_Z.fac_N = -N_com / (N_Rd_red);
				}
				else {
					util_element_.bending_Y.fac_N = N_com / (chi_y*N_Rd_red);
					util_element_.bending_Z.fac_N = N_com / (chi_z*N_Rd_red);
				}
			}

			if (My_Rd_red == 0) {
				util_element_.bending_Y.fac_My = 1000;
				util_element_.bending_Z.fac_My = 1000;
			}
			else {
				util_element_.bending_Y.fac_My = std::abs(kyy_ * My / (chi_LT_mod*My_Rd_red));
				util_element_.bending_Z.fac_My = std::abs(kzy_ * My / (chi_LT_mod*My_Rd_red));
			}

			if (Mz_Rd_red == 0) {
				util_element_.bending_Y.fac_Mz = 1000;
				util_element_.bending_Z.fac_Mz = 1000;
			}
			else {
				util_element_.bending_Y.fac_Mz = std::abs(kyz_ * Mz / Mz_Rd_red);
				util_element_.bending_Z.fac_Mz = std::abs(kzz_ * Mz / Mz_Rd_red);
			}

			real res_util_element_ = util_element_.util();

			if (util_grp < res_util_element_) {
				util_grp = res_util_element_;
				util_ = util_element_;
				calc_.sample_point_ind = sample_point_ind++;
				calc_.crosec_class = crosec_class;
				calc_.N_Rd = N_Rd_red;
				calc_.Vy_Rd = Vy_Rd;
				calc_.Vz_Rd = Vz_Rd;
				calc_.Mt_Rd = Mtrd;
				calc_.My_Rd = My_Rd_red;
				calc_.Mz_Rd = Mz_Rd_red;
				calc_.Mcr = Mcr;
				calc_.Ncr_y = Ncr_y;
				calc_.Ncr_z = Ncr_z;
				calc_.psi_y = psi_y;
				calc_.Cmy = Cmy;
				calc_.Cmz = Cmz;
				calc_.CmLT = CmLT;
				calc_.chi_y = chi_y;
				calc_.chi_z = chi_z;
				calc_.chi_LT_mod = chi_LT_mod;
				calc_.kyy = kyy_;
				calc_.kzz = kzz_;
				calc_.kyz = kyz_;
				calc_.kzy = kzy_;

				if (util_grp > utilization_) {
					utilization_ = util_grp;
					compressive_ = cs_force[Node::x_t] < 0;
				}

			}
		}
#endif
	}

	real
		UtilChecker_EC3_Elastic::kyy(index_t crosec_class, real Cmy, real lam_y_bar,
		real N, real chi_y, real Nrd) const
	{
		if (Nrd == 0) return 1.0;

		real kyy1 = Cmy*(1 + 0.6*lam_y_bar*N / (chi_y*Nrd));
		real kyy2 = Cmy*(1 + 0.6*N / (chi_y*Nrd));
		real kyy = kyy1 < kyy2 ? kyy1 : kyy2;
		return kyy;
	}

	void
		UtilChecker_EC3_Plastic::copy(const UtilChecker& orig, DCMap* mapDict) {
		UtilChecker_EC3_Elastic::copy(orig, mapDict);
	}
	
	real
		UtilChecker_EC3_Plastic::kyy(index_t crosec_class, real Cmy, real lam_y_bar,
		real N, real chi_y, real Nrd) const
	{
		if (Nrd == 0) return 1.0;

		if (crosec_class>2) return
			UtilChecker_EC3_Elastic::kyy(crosec_class, Cmy, lam_y_bar, N, chi_y, Nrd);

		real kyy1 = Cmy*(1 + (lam_y_bar - 0.2)*N / (chi_y*Nrd));
		real kyy2 = Cmy*(1 + 0.8*N / (chi_y*Nrd));
		real kyy = kyy1 < kyy2 ? kyy1 : kyy2;
		return kyy;
	}

	real
		UtilChecker_EC3_Elastic::kzz(index_t crosec_class, real Cmz, real lam_z_bar,
		real N, real chi_z, real Nrd, bool is_I_profil) const
	{
		if (Nrd == 0) return 1.0;

		real kzz1 = Cmz*(1 + 0.6*lam_z_bar*N / (chi_z*Nrd));
		real kzz2 = Cmz*(1 + 0.6*N / (chi_z*Nrd));
		real kzz = kzz1 < kzz2 ? kzz1 : kzz2;
		return kzz;
	}

	real
		UtilChecker_EC3_Plastic::kzz(index_t crosec_class, real Cmz, real lam_z_bar,
		real N, real chi_z, real Nrd, bool is_I_profil) const
	{
		if (Nrd == 0) return 1.0;

		if (crosec_class>2) return
			UtilChecker_EC3_Elastic::kzz(crosec_class, Cmz, lam_z_bar,
			N, chi_z, Nrd, is_I_profil);

		real kzz;
		if (is_I_profil) {
			real kzz1 = Cmz*(1 + (2 * lam_z_bar - 0.6)*N / (chi_z*Nrd));
			real kzz2 = Cmz*(1 + 1.4*N / (chi_z*Nrd));
			kzz = kzz1 < kzz2 ? kzz1 : kzz2;
		}
		else {
			real kzz1 = Cmz*(1 + (lam_z_bar - 0.2)*N / (chi_z*Nrd));
			real kzz2 = Cmz*(1 + 0.8*N / (chi_z*Nrd));
			kzz = kzz1 < kzz2 ? kzz1 : kzz2;
		}
		return kzz;
	}

	real
		UtilChecker_EC3_Elastic::kyz(index_t crosec_class, real kzz) const
	{
		return kzz;
	}

	real
		UtilChecker_EC3_Plastic::kyz(index_t crosec_class, real kzz) const
	{
		if (crosec_class>2) return
			UtilChecker_EC3_Elastic::kyz(crosec_class, kzz);

		return 0.6*kzz;;
	}

	real
		UtilChecker_EC3_Elastic::kzy(index_t crosec_class, real CmLT, real lam_z_bar,
		real N, real chi_z, real Nrd, bool is_I_profil, real kyy) const
	{
		if (Nrd == 0) return 1.0;

		if (is_I_profil) {
			real kzy1 = 1 - 0.05*lam_z_bar / (CmLT - 0.25)*N / (chi_z*Nrd);
			real kzy2 = 1 - 0.05 / (CmLT - 0.25)*N / (chi_z*Nrd);
			real kzy = kzy1 < kzy2 ? kzy1 : kzy2;
			return kzy;
		}
		else {
			return 0.8*kyy;
		}
	}

	real
		UtilChecker_EC3_Plastic::kzy(index_t crosec_class, real CmLT, real lam_z_bar,
		real N, real chi_z, real Nrd, bool is_I_profil, real kyy) const
	{
		if (Nrd == 0) return 1.0;

		if (crosec_class>2) return
			UtilChecker_EC3_Elastic::kzy(crosec_class, CmLT, lam_z_bar,
			N, chi_z, Nrd, is_I_profil, kyy);

		if (is_I_profil) {
			if (lam_z_bar < 0.4) {
				real kzy1 = 0.6 + lam_z_bar;
				real kzy2 = 1 - 0.1*lam_z_bar / (CmLT - 0.25)*N / (chi_z*Nrd);
				return kzy1 < kzy2 ? kzy1 : kzy2;
			}
			else {
				real kzy1 = 1 - 0.1*lam_z_bar / (CmLT - 0.25)*N / (chi_z*Nrd);
				real kzy2 = 1 - 0.1 / (CmLT - 0.25)*N / (chi_z*Nrd);
				return kzy1 < kzy2 ? kzy1 : kzy2;
			}
		}
		else {
			return 0.6*kyy;
		}
	}
		
	std::vector<real>
		UtilChecker_EC3_Elastic::psi(const SuperImp::CSForceCollection& cs_forces, const Element& elem, int dof) const {
		index_t lc_num = cs_forces.numberOfLoadCases();
		std::vector<real> res(lc_num);
		for (index_t lc_ind = 0; lc_ind < lc_num; ++lc_ind) {
			real M0 = (*cs_forces.const_csForceIteratorBegin(elem, lc_ind))[dof];
			real M2 = (*cs_forces.const_csForceIteratorLast(elem, lc_ind))[dof];
			real Mh = std::abs(M0) > std::abs(M2) ? M0 : M2;
			real Mpsi = std::abs(M0) < std::abs(M2) ? M0 : M2;
			real psi = 1.0;
			if (Mh != 0) psi = Mpsi / Mh;
			res[lc_ind] = psi;
		}
		return res;
	}

	std::vector<real>
		UtilChecker_EC3_Elastic::CmValue(const SuperImp::CSForceCollection& cs_forces, const Element& elem, int dof) const {
		index_t lc_num = cs_forces.numberOfLoadCases();
		std::vector<real> res(lc_num);
		for (index_t lc_ind = 0; lc_ind < lc_num; ++lc_ind) {
			real M0 = (*cs_forces.const_csForceIteratorBegin(elem, lc_ind))[dof];
			real M2 = (*cs_forces.const_csForceIteratorLast(elem, lc_ind))[dof];
			SuperImp::const_CSForceIterator iter_middle = cs_forces.const_csForceIteratorMiddle(elem, lc_ind);
			real Ms = (*iter_middle)[dof];
			if (cs_forces.nSamples(elem) % 2 == 0) {
				Ms = (Ms + (*(++iter_middle))[dof])*0.5;
			}

			real Mh = std::abs(M0) > std::abs(M2) ? M0 : M2;
			real Mpsi = std::abs(M0) < std::abs(M2) ? M0 : M2;

			real psi = 1.0;
			if (Mh != 0) psi = Mpsi / Mh;

			real cm;
			// a uniformly distributed load is assumed
			if (std::abs(Mh) > std::abs(Ms)) {
				real alpha_s = Ms / Mh;
				if (alpha_s >= 0) {
					cm = 0.2 + 0.8*alpha_s;
					cm = cm < 0.4 ? 0.4 : cm;
				}
				else {
					if (psi >= 0) {
						cm = 0.1 - 0.8*alpha_s;
						cm = cm < 0.4 ? 0.4 : cm;
					}
					else {
						cm = 0.1*(1 - psi) - 0.8*alpha_s;
						cm = cm < 0.4 ? 0.4 : cm;
					}
				}
			}
			else {
				real alpha_h = 1;
				if (Ms != 0) alpha_h = Mh / Ms;
				if (alpha_h >= 0) {
					cm = 0.95 + 0.05*alpha_h;
				}
				else {
					if (psi >= 0) {
						cm = 0.95 + 0.05*alpha_h;
					}
					else {
						cm = 0.95 + 0.05*alpha_h*(1.0 + 2.0*psi);
					}
				}
			}

			// for buckling with lateral displacement according to EN 1993-1-1, Table B.3
			// is too conservative for systems that do no displace laterally
			cm = cm < 0.9 ? 0.9 : cm;
			res[lc_ind] = cm;
		}
		return res;
	}

	real
		UtilChecker_EC3_Plastic::Wy(index_t crosec_class, const Beam3DCroSec& crosec) const {
		if (crosec_class <= 2) return crosec.Wply;
		return crosec.Wely;
	}

	real
		UtilChecker_EC3_Plastic::Wz(index_t crosec_class, const Beam3DCroSec& crosec) const {
		if (crosec_class <= 2) return crosec.Wplz;
		return crosec.Welz;
	}

}