#include #include #include #include #include "vs.h" #include "fe.h" #include //============================================================================== // Step 1: Global_Discretizations //============================================================================== // Finite deformation elastoplasticity references in Comp. Meth. in Appl. Mech. Engrg.: // Finite Deformaiton: Simo, Taylor and Pister (1985, vol. 51, p.177) // Simo (1988, vol 66, p. 199, and 1988, vol. 68, p.1) // Simo, Kennedy and Taylor, (1989, vol. 74, p.177) // Elastoplasticity: Simo and Taylor (1985, vol. 48, p.114) static const double E_ = 70.0e6; static const double v_ = 0.2; static const double Y_infinity = 0.234e6; static const double Y_0 = 0.234e6; static const double iota = 0.1; static const double d_kapa = 0.01; // linear isotropic hardenning parameter static const double d_H_alpha = 0.0; // no kinematic hardenning static const double lambda_ = v_*E_/((1.0+v_)*(1.0-2.0*v_)); static const double mu_ = E_/(2.0*(1.0+v_)); static const double K_ = lambda_+2.0/3.0*mu_; static ofstream ofs("temp.cpp", ios::out | ios::trunc); // DEFINITION of Global_Discretization static Quadrature qp(2, 4); // 2-dimension, 2x2 integration points for history terms class Non_Linear_Omega_eh : public Omega_eh { // plastic history terms, (.)_n is at t_n; (.) is at t_n+1 C0 the_theta, // mean dilatation; averaged J over an element domain, theta = (J|dv) / meas(Omega_eh) the_theta_n, the_p, // average pressure the_p_n; H0 the_J, // the Jacobian = det(F) the_J_n, the_F, // deformation gradient tensor the_F_n, the_sigma, // Cauchy stresses; tau (Kirchhoff stresses) = theta * sigma the_sigma_n, // sigma use engineering vector convention {sigma_xx, sigma_yy, sigma_xy} the_F_p_inv, // inverse of plastic deformation gradient tensor the_F_p_inv_n, the_ep, // total plastic strain the_ep_n, the_alpha, // "back stresses" (for kinematic hardenning law) the_alpha_n; //void __initialization(); public: Non_Linear_Omega_eh(int en = 0, int etn = 0, int mtn = 0, int nen = 0, int* ena = NULL); C0& theta(); C0& p(); H0& J(); H0& F(); H0& sigma(); H0& F_p_inv(); H0& ep(); H0& alpha(); C0& theta_n(); C0& p_n(); H0& J_n(); H0& F_n(); H0& sigma_n(); H0& F_p_inv_n(); H0& ep_n(); H0& alpha_n(); void update_history_data(); }; /*void Non_Linear_Omega_eh::__initialization() { the_theta &= C0(1.0); the_J &= H0(qp); the_F &= INTEGRABLE_MATRIX("int, int, Quadrature", 2, 2, qp); //the_F &= H0(2, 2, (double*)0, qp); cout << the_F[0][0] << endl; the_sigma &= H0(3, (double*)0, qp); // engineering vector convention {sigma_xx, sigma_yy, sigma_xy} the_F_p_inv &= INTEGRABLE_MATRIX("int, int, Quadrature", 2, 2, qp); the_ep &= H0((double*)0, qp); the_alpha &= H0(3, (double*)0, qp); the_theta_n &= C0(1.0); the_J_n &= H0(qp); the_F_n &= INTEGRABLE_MATRIX("int, int, Quadrature", 2, 2, qp); the_sigma_n &= H0(3, (double*)0, qp); the_F_p_inv_n &= INTEGRABLE_MATRIX("int, int, Quadrature", 2, 2, qp); the_ep_n &= H0((double*)0, qp); the_alpha_n &= H0(3, (double*)0, qp); the_J = 1.0; the_J_n = 1.0; // set identity matrix the_F(0,0) = the_F(1,1) = the_F_n(0,0) = the_F_n(1,1) = the_F_p_inv(0,0) = the_F_p_inv(1,1) = the_F_p_inv_n(0,0) = the_F_p_inv_n(1,1) = 1.0; } */ static double integrable_identity[4][2][2] = { { {1.0, 0.0}, {0.0, 1.0} }, { {1.0, 0.0}, {0.0, 1.0} }, { {1.0, 0.0}, {0.0, 1.0} }, { {1.0, 0.0}, {0.0, 1.0} } }; Non_Linear_Omega_eh::Non_Linear_Omega_eh(int en, int etn, int mtn, int nen, int* ena) : Omega_eh(en, etn, mtn, nen, ena) { //__initialization(); the_theta &= C0(1.0); the_p &= C0(0.0); the_J &= H0(qp); the_F &= INTEGRABLE_MATRIX("int, int, const double*, Quadrature, int, int", 2, 2, integrable_identity[0][0], qp, 2, 2); the_sigma &= H0(4, (double*)0, qp); // vector convention {sigma_xx, sigma_yy, sigma_zz, sigma_xy} the_F_p_inv &= INTEGRABLE_MATRIX("int, int, const double*, Quadrature, int, int", 2, 2, integrable_identity[0][0], qp, 2, 2); the_ep &= H0((double*)0, qp); the_alpha &= H0(4, (double*)0, qp); the_theta_n &= C0(1.0); the_p_n &= C0(0.0);; the_J_n &= H0(qp); the_F_n &= INTEGRABLE_MATRIX("int, int, const double*, Quadrature, int, int", 2, 2, integrable_identity[0][0], qp, 2, 2); the_sigma_n &= H0(4, (double*)0, qp); the_F_p_inv_n &= INTEGRABLE_MATRIX("int, int, const double*, Quadrature, int, int", 2, 2, integrable_identity[0][0], qp, 2, 2); the_ep_n &= H0((double*)0, qp); the_alpha_n &= H0(4, (double*)0, qp); the_J = 1.0; the_J_n = 1.0; } C0& Non_Linear_Omega_eh::theta() { return the_theta; } C0& Non_Linear_Omega_eh::p() { return the_p; } H0& Non_Linear_Omega_eh::J() { return the_J; } H0& Non_Linear_Omega_eh::F() { return the_F; } H0& Non_Linear_Omega_eh::sigma() { return the_sigma; } H0& Non_Linear_Omega_eh::F_p_inv() { return the_F_p_inv; } H0& Non_Linear_Omega_eh::ep() { return the_ep; } H0& Non_Linear_Omega_eh::alpha() { return the_alpha; } C0& Non_Linear_Omega_eh::theta_n() { return the_theta_n; } C0& Non_Linear_Omega_eh::p_n() { return the_p_n; } H0& Non_Linear_Omega_eh::J_n() { return the_J_n; } H0& Non_Linear_Omega_eh::F_n() { return the_F_n; } H0& Non_Linear_Omega_eh::sigma_n() { return the_sigma_n; } H0& Non_Linear_Omega_eh::F_p_inv_n() { return the_F_p_inv_n; } H0& Non_Linear_Omega_eh::ep_n() { return the_ep_n; } H0& Non_Linear_Omega_eh::alpha_n() { return the_alpha_n; } void Non_Linear_Omega_eh::update_history_data() { the_theta_n = the_theta; the_p_n = the_p; the_J_n = the_J; the_F_n = the_F; the_sigma_n = the_sigma; the_F_p_inv_n = the_F_p_inv; the_ep_n = the_ep; the_alpha_n = the_alpha; } void define_node(Dynamic_Array& dna); void define_element(Dynamic_Array& dea); Omega_h::Omega_h() { // discritized global domain define_node(node_array()); define_element(omega_eh_array()); } void define_node(Dynamic_Array& dna) { double v[2]; Node *node; v[0] = 5.0; v[1] = 0.0; node = new Node(0, 2, v); dna.add(node); v[0] = 6.25; node = new Node(1, 2, v); dna.add(node); v[0] = 7.5; node = new Node(2, 2, v); dna.add(node); v[0] = 8.625; node = new Node(3, 2, v); dna.add(node); v[0] = 10.0; node = new Node(4, 2, v); dna.add(node); v[0] = 4.9039; v[1] = 0.9755; node = new Node(5, 2, v); dna.add(node); v[0] = 6.1299; v[1] = 1.2189; node = new Node(6, 2, v); dna.add(node); v[0] = 7.3559; v[1] = 1.4632; node = new Node(7, 2, v); dna.add(node); v[0] = 8.678; v[1] = 1.4632; node = new Node(8, 2, v); dna.add(node); v[0] = 10.0; v[1] = 1.4632; node = new Node(9, 2, v); dna.add(node); v[0] = 4.6194; v[1] = 1.9134; node = new Node(10, 2, v); dna.add(node); v[0] = 5.7742; v[1] = 2.3918; node = new Node(11, 2, v); dna.add(node); v[0] = 6.9291; v[1] = 2.8701; node = new Node(12, 2, v); dna.add(node); v[0] = 8.625; node = new Node(13, 2, v); dna.add(node); v[0] = 10.0; node = new Node(14, 2, v); dna.add(node); v[0] = 4.1573; v[1] = 2.7779; node = new Node(15, 2, v); dna.add(node); v[0] = 5.1967; v[1] = 3.4724; node = new Node(16, 2, v); dna.add(node); v[0] = 6.2360; v[1] = 4.1668; node = new Node(17, 2, v); dna.add(node); v[0] = 8.118; v[1] = 4.6760; node = new Node(18, 2, v); dna.add(node); v[0] = 10.0; v[1] = 5.1851; node = new Node(19, 2, v); dna.add(node); v[0] = 3.5355; v[1] = 3.5355; node = new Node(20, 2, v); dna.add(node); v[0] = 4.4194; v[1] = 4.4194; node = new Node(21, 2, v); dna.add(node); v[0] = 5.3033; v[1] = 5.3033; node = new Node(22, 2, v); dna.add(node); v[0] = 7.6517; v[1] = 6.4017; node = new Node(23, 2, v); dna.add(node); v[0] = 10.0; v[1] = 7.5; node = new Node(24, 2, v); dna.add(node); v[0] = 0.0; v[1] = 5.0; node = new Node(25, 2, v); dna.add(node); v[0] = 0.9755; v[1] = 4.9039; node = new Node(26, 2, v); dna.add(node); v[0] = 1.9134; v[1] = 4.6194; node = new Node(27, 2, v); dna.add(node); v[0] = 2.7779; v[1] = 4.1573; node = new Node(28, 2, v); dna.add(node); v[0] = 0.0; v[1] = 6.25; node = new Node(29, 2, v); dna.add(node); v[0] = 1.1959; v[1] = 6.0121; node = new Node(30, 2, v); dna.add(node); v[0] = 2.3918; v[1] = 5.7742; node = new Node(31, 2, v); dna.add(node); v[0] = 3.4056; v[1] = 5.0968; node = new Node(32, 2, v); dna.add(node); v[0] = 0.0; v[1] = 7.5; node = new Node(33, 2, v); dna.add(node); v[0] = 1.4632; v[1] = 7.3559; node = new Node(34, 2, v); dna.add(node); v[0] = 2.8701; v[1] = 6.9291; node = new Node(35, 2, v); dna.add(node); v[0] = 4.1668; v[1] = 6.2360; node = new Node(36, 2, v); dna.add(node); v[0] = 0.0; v[1] = 8.75; node = new Node(37, 2, v); dna.add(node); v[0] = 1,4351; node = new Node(38, 2, v); dna.add(node); v[0] = 2.8701; node = new Node(39, 2, v); dna.add(node); v[0] = 4.0867; node = new Node(40, 2, v); dna.add(node); v[0] = 5.3033; node = new Node(41, 2, v); dna.add(node); v[0] = 7.6517; node = new Node(42, 2, v); dna.add(node); v[0] = 10.0; node = new Node(43, 2, v); dna.add(node); v[0] = 0.0; v[1] = 10.0; node = new Node(44, 2, v); dna.add(node); v[0] = 1,4351; node = new Node(45, 2, v); dna.add(node); v[0] = 2.8701; node = new Node(46, 2, v); dna.add(node); v[0] = 4.0867; node = new Node(47, 2, v); dna.add(node); v[0] = 5.3033; node = new Node(48, 2, v); dna.add(node); v[0] = 7.6517; node = new Node(49, 2, v); dna.add(node); v[0] = 10.0; node = new Node(50, 2, v); dna.add(node); v[0] = 0.0; v[1] = 11.25; node = new Node(51, 2, v); dna.add(node); v[0] = 1,4351; node = new Node(52, 2, v); dna.add(node); v[0] = 2.8701; node = new Node(53, 2, v); dna.add(node); v[0] = 4.0867; node = new Node(54, 2, v); dna.add(node); v[0] = 5.3033; node = new Node(55, 2, v); dna.add(node); v[0] = 7.6517; node = new Node(56, 2, v); dna.add(node); v[0] = 10.0; node = new Node(57, 2, v); dna.add(node); v[0] = 0.0; v[1] = 12.6; node = new Node(58, 2, v); dna.add(node); v[0] = 1,4351; node = new Node(59, 2, v); dna.add(node); v[0] = 2.8701; node = new Node(60, 2, v); dna.add(node); v[0] = 4.0867; node = new Node(61, 2, v); dna.add(node); v[0] = 5.3033; node = new Node(62, 2, v); dna.add(node); v[0] = 7.6517; node = new Node(63, 2, v); dna.add(node); v[0] = 10.0; node = new Node(64, 2, v); dna.add(node); v[0] = 0.0; v[1] = 13.95; node = new Node(65, 2, v); dna.add(node); v[0] = 1,4351; node = new Node(66, 2, v); dna.add(node); v[0] = 2.8701; node = new Node(67, 2, v); dna.add(node); v[0] = 4.0867; node = new Node(68, 2, v); dna.add(node); v[0] = 5.3033; node = new Node(69, 2, v); dna.add(node); v[0] = 7.6517; node = new Node(70, 2, v); dna.add(node); v[0] = 10.0; node = new Node(71, 2, v); dna.add(node); v[0] = 0.0; v[1] = 15.3; node = new Node(72, 2, v); dna.add(node); v[0] = 1,4351; node = new Node(73, 2, v); dna.add(node); v[0] = 2.8701; node = new Node(74, 2, v); dna.add(node); v[0] = 4.0867; node = new Node(75, 2, v); dna.add(node); v[0] = 5.3033; node = new Node(76, 2, v); dna.add(node); v[0] = 7.6517; node = new Node(77, 2, v); dna.add(node); v[0] = 10.0; node = new Node(78, 2, v); dna.add(node); v[0] = 0.0; v[1] = 16.65; node = new Node(79, 2, v); dna.add(node); v[0] = 1,4351; node = new Node(80, 2, v); dna.add(node); v[0] = 2.8701; node = new Node(81, 2, v); dna.add(node); v[0] = 4.0867; node = new Node(82, 2, v); dna.add(node); v[0] = 5.3033; node = new Node(83, 2, v); dna.add(node); v[0] = 7.6517; node = new Node(84, 2, v); dna.add(node); v[0] = 10.0; node = new Node(85, 2, v); dna.add(node); v[0] = 0.0; v[1] = 18.0; node = new Node(86, 2, v); dna.add(node); v[0] = 1,4351; node = new Node(87, 2, v); dna.add(node); v[0] = 2.8701; node = new Node(88, 2, v); dna.add(node); v[0] = 4.0867; node = new Node(89, 2, v); dna.add(node); v[0] = 5.3033; node = new Node(90, 2, v); dna.add(node); v[0] = 7.6517; node = new Node(91, 2, v); dna.add(node); v[0] = 10.0; node = new Node(92, 2, v); dna.add(node); } void define_element(Dynamic_Array& dea) { // define elements; four-node int ena[4]; // element node number array Omega_eh *elem; ena[0] = 0; ena[1] = 1; ena[2] = 6; ena[3] = 5; elem = new Non_Linear_Omega_eh(0, 0, 0, 4, ena); dea.add(elem); ena[0] = 1; ena[1] = 2; ena[2] = 7; ena[3] = 6; elem = new Non_Linear_Omega_eh(1, 0, 0, 4, ena); dea.add(elem); ena[0] = 2; ena[1] = 3; ena[2] = 8; ena[3] = 7; elem = new Non_Linear_Omega_eh(2, 0, 0, 4, ena); dea.add(elem); ena[0] = 3; ena[1] = 4; ena[2] = 9; ena[3] = 8; elem = new Non_Linear_Omega_eh(3, 0, 0, 4, ena); dea.add(elem); ena[0] = 5; ena[1] = 6; ena[2] = 11; ena[3] = 10; elem = new Non_Linear_Omega_eh(4, 0, 0, 4, ena); dea.add(elem); ena[0] = 6; ena[1] = 7; ena[2] = 12; ena[3] = 11; elem = new Non_Linear_Omega_eh(5, 0, 0, 4, ena); dea.add(elem); ena[0] = 7; ena[1] = 8; ena[2] = 13; ena[3] = 12; elem = new Non_Linear_Omega_eh(6, 0, 0, 4, ena); dea.add(elem); ena[0] = 8; ena[1] = 9; ena[2] = 14; ena[3] = 13; elem = new Non_Linear_Omega_eh(7, 0, 0, 4, ena); dea.add(elem); ena[0] = 10; ena[1] = 11; ena[2] = 16; ena[3] = 15; elem = new Non_Linear_Omega_eh(8, 0, 0, 4, ena); dea.add(elem); ena[0] = 11; ena[1] = 12; ena[2] = 17; ena[3] = 16; elem = new Non_Linear_Omega_eh(9, 0, 0, 4, ena); dea.add(elem); ena[0] = 12; ena[1] = 13; ena[2] = 18; ena[3] = 17; elem = new Non_Linear_Omega_eh(10, 0, 0, 4, ena); dea.add(elem); ena[0] = 13; ena[1] = 14; ena[2] = 19; ena[3] = 18; elem = new Non_Linear_Omega_eh(11, 0, 0, 4, ena); dea.add(elem); ena[0] = 15; ena[1] = 16; ena[2] = 21; ena[3] = 20; elem = new Non_Linear_Omega_eh(12, 0, 0, 4, ena); dea.add(elem); ena[0] = 16; ena[1] = 17; ena[2] = 22; ena[3] = 21; elem = new Non_Linear_Omega_eh(13, 0, 0, 4, ena); dea.add(elem); ena[0] = 17; ena[1] = 18; ena[2] = 23; ena[3] = 22; elem = new Non_Linear_Omega_eh(14, 0, 0, 4, ena); dea.add(elem); ena[0] = 18; ena[1] = 19; ena[2] = 24; ena[3] = 23; elem = new Non_Linear_Omega_eh(15, 0, 0, 4, ena); dea.add(elem); ena[0] = 25; ena[1] = 26; ena[2] = 30; ena[3] = 29; elem = new Non_Linear_Omega_eh(16, 0, 0, 4, ena); dea.add(elem); ena[0] = 26; ena[1] = 27; ena[2] = 31; ena[3] = 30; elem = new Non_Linear_Omega_eh(17, 0, 0, 4, ena); dea.add(elem); ena[0] = 27; ena[1] = 28; ena[2] = 32; ena[3] = 31; elem = new Non_Linear_Omega_eh(18, 0, 0, 4, ena); dea.add(elem); ena[0] = 28; ena[1] = 20; ena[2] = 21; ena[3] = 32; elem = new Non_Linear_Omega_eh(19, 0, 0, 4, ena); dea.add(elem); ena[0] = 29; ena[1] = 30; ena[2] = 34; ena[3] = 33; elem = new Non_Linear_Omega_eh(20, 0, 0, 4, ena); dea.add(elem); ena[0] = 30; ena[1] = 31; ena[2] = 35; ena[3] = 34; elem = new Non_Linear_Omega_eh(21, 0, 0, 4, ena); dea.add(elem); ena[0] = 31; ena[1] = 32; ena[2] = 36; ena[3] = 35; elem = new Non_Linear_Omega_eh(22, 0, 0, 4, ena); dea.add(elem); ena[0] = 32; ena[1] = 21; ena[2] = 22; ena[3] = 36; elem = new Non_Linear_Omega_eh(23, 0, 0, 4, ena); dea.add(elem); ena[0] = 33; ena[1] = 34; ena[2] = 38; ena[3] = 37; elem = new Non_Linear_Omega_eh(24, 0, 0, 4, ena); dea.add(elem); ena[0] = 34; ena[1] = 35; ena[2] = 39; ena[3] = 38; elem = new Non_Linear_Omega_eh(25, 0, 0, 4, ena); dea.add(elem); ena[0] = 35; ena[1] = 36; ena[2] = 40; ena[3] = 39; elem = new Non_Linear_Omega_eh(26, 0, 0, 4, ena); dea.add(elem); ena[0] = 36; ena[1] = 22; ena[2] = 41; ena[3] = 40; elem = new Non_Linear_Omega_eh(27, 0, 0, 4, ena); dea.add(elem); ena[0] = 22; ena[1] = 23; ena[2] = 42; ena[3] = 41; elem = new Non_Linear_Omega_eh(28, 0, 0, 4, ena); dea.add(elem); ena[0] = 23; ena[1] = 24; ena[2] = 43; ena[3] = 42; elem = new Non_Linear_Omega_eh(29, 0, 0, 4, ena); dea.add(elem); ena[0] = 37; ena[1] = 38; ena[2] = 45; ena[3] = 44; elem = new Non_Linear_Omega_eh(30, 0, 0, 4, ena); dea.add(elem); ena[0] = 38; ena[1] = 39; ena[2] = 46; ena[3] = 45; elem = new Non_Linear_Omega_eh(31, 0, 0, 4, ena); dea.add(elem); ena[0] = 39; ena[1] = 40; ena[2] = 47; ena[3] = 46; elem = new Non_Linear_Omega_eh(32, 0, 0, 4, ena); dea.add(elem); ena[0] = 40; ena[1] = 41; ena[2] = 48; ena[3] = 47; elem = new Non_Linear_Omega_eh(33, 0, 0, 4, ena); dea.add(elem); ena[0] = 41; ena[1] = 42; ena[2] = 49; ena[3] = 48; elem = new Non_Linear_Omega_eh(34, 0, 0, 4, ena); dea.add(elem); ena[0] = 42; ena[1] = 43; ena[2] = 50; ena[3] = 49; elem = new Non_Linear_Omega_eh(35, 0, 0, 4, ena); dea.add(elem); ena[0] = 44; ena[1] = 45; ena[2] = 52; ena[3] = 51; elem = new Non_Linear_Omega_eh(36, 0, 0, 4, ena); dea.add(elem); ena[0] = 45; ena[1] = 46; ena[2] = 53; ena[3] = 52; elem = new Non_Linear_Omega_eh(37, 0, 0, 4, ena); dea.add(elem); ena[0] = 46; ena[1] = 47; ena[2] = 54; ena[3] = 53; elem = new Non_Linear_Omega_eh(38, 0, 0, 4, ena); dea.add(elem); ena[0] = 47; ena[1] = 48; ena[2] = 55; ena[3] = 54; elem = new Non_Linear_Omega_eh(39, 0, 0, 4, ena); dea.add(elem); ena[0] = 48; ena[1] = 49; ena[2] = 56; ena[3] = 55; elem = new Non_Linear_Omega_eh(40, 0, 0, 4, ena); dea.add(elem); ena[0] = 49; ena[1] = 50; ena[2] = 57; ena[3] = 56; elem = new Non_Linear_Omega_eh(41, 0, 0, 4, ena); dea.add(elem); ena[0] = 51; ena[1] = 52; ena[2] = 59; ena[3] = 58; elem = new Non_Linear_Omega_eh(42, 0, 0, 4, ena); dea.add(elem); ena[0] = 52; ena[1] = 53; ena[2] = 60; ena[3] = 59; elem = new Non_Linear_Omega_eh(43, 0, 0, 4, ena); dea.add(elem); ena[0] = 53; ena[1] = 54; ena[2] = 61; ena[3] = 60; elem = new Non_Linear_Omega_eh(44, 0, 0, 4, ena); dea.add(elem); ena[0] = 54; ena[1] = 55; ena[2] = 62; ena[3] = 61; elem = new Non_Linear_Omega_eh(45, 0, 0, 4, ena); dea.add(elem); ena[0] = 55; ena[1] = 56; ena[2] = 63; ena[3] = 62; elem = new Non_Linear_Omega_eh(46, 0, 0, 4, ena); dea.add(elem); ena[0] = 56; ena[1] = 57; ena[2] = 64; ena[3] = 63; elem = new Non_Linear_Omega_eh(47, 0, 0, 4, ena); dea.add(elem); ena[0] = 58; ena[1] = 59; ena[2] = 66; ena[3] = 65; elem = new Non_Linear_Omega_eh(48, 0, 0, 4, ena); dea.add(elem); ena[0] = 59; ena[1] = 60; ena[2] = 67; ena[3] = 66; elem = new Non_Linear_Omega_eh(49, 0, 0, 4, ena); dea.add(elem); ena[0] = 60; ena[1] = 61; ena[2] = 68; ena[3] = 67; elem = new Non_Linear_Omega_eh(50, 0, 0, 4, ena); dea.add(elem); ena[0] = 61; ena[1] = 62; ena[2] = 69; ena[3] = 68; elem = new Non_Linear_Omega_eh(51, 0, 0, 4, ena); dea.add(elem); ena[0] = 62; ena[1] = 63; ena[2] = 70; ena[3] = 69; elem = new Non_Linear_Omega_eh(52, 0, 0, 4, ena); dea.add(elem); ena[0] = 63; ena[1] = 64; ena[2] = 71; ena[3] = 70; elem = new Non_Linear_Omega_eh(53, 0, 0, 4, ena); dea.add(elem); ena[0] = 65; ena[1] = 66; ena[2] = 73; ena[3] = 72; elem = new Non_Linear_Omega_eh(54, 0, 0, 4, ena); dea.add(elem); ena[0] = 66; ena[1] = 67; ena[2] = 74; ena[3] = 73; elem = new Non_Linear_Omega_eh(55, 0, 0, 4, ena); dea.add(elem); ena[0] = 67; ena[1] = 68; ena[2] = 75; ena[3] = 74; elem = new Non_Linear_Omega_eh(56, 0, 0, 4, ena); dea.add(elem); ena[0] = 68; ena[1] = 69; ena[2] = 76; ena[3] = 75; elem = new Non_Linear_Omega_eh(57, 0, 0, 4, ena); dea.add(elem); ena[0] = 69; ena[1] = 70; ena[2] = 77; ena[3] = 76; elem = new Non_Linear_Omega_eh(58, 0, 0, 4, ena); dea.add(elem); ena[0] = 70; ena[1] = 71; ena[2] = 78; ena[3] = 77; elem = new Non_Linear_Omega_eh(59, 0, 0, 4, ena); dea.add(elem); ena[0] = 72; ena[1] = 73; ena[2] = 80; ena[3] = 79; elem = new Non_Linear_Omega_eh(60, 0, 0, 4, ena); dea.add(elem); ena[0] = 73; ena[1] = 74; ena[2] = 81; ena[3] = 80; elem = new Non_Linear_Omega_eh(61, 0, 0, 4, ena); dea.add(elem); ena[0] = 74; ena[1] = 75; ena[2] = 82; ena[3] = 81; elem = new Non_Linear_Omega_eh(62, 0, 0, 4, ena); dea.add(elem); ena[0] = 75; ena[1] = 76; ena[2] = 83; ena[3] = 82; elem = new Non_Linear_Omega_eh(63, 0, 0, 4, ena); dea.add(elem); ena[0] = 76; ena[1] = 77; ena[2] = 84; ena[3] = 83; elem = new Non_Linear_Omega_eh(64, 0, 0, 4, ena); dea.add(elem); ena[0] = 77; ena[1] = 78; ena[2] = 85; ena[3] = 84; elem = new Non_Linear_Omega_eh(65, 0, 0, 4, ena); dea.add(elem); ena[0] = 79; ena[1] = 80; ena[2] = 87; ena[3] = 86; elem = new Non_Linear_Omega_eh(66, 0, 0, 4, ena); dea.add(elem); ena[0] = 80; ena[1] = 81; ena[2] = 88; ena[3] = 87; elem = new Non_Linear_Omega_eh(67, 0, 0, 4, ena); dea.add(elem); ena[0] = 81; ena[1] = 82; ena[2] = 89; ena[3] = 88; elem = new Non_Linear_Omega_eh(68, 0, 0, 4, ena); dea.add(elem); ena[0] = 82; ena[1] = 83; ena[2] = 90; ena[3] = 89; elem = new Non_Linear_Omega_eh(69, 0, 0, 4, ena); dea.add(elem); ena[0] = 83; ena[1] = 84; ena[2] = 91; ena[3] = 90; elem = new Non_Linear_Omega_eh(70, 0, 0, 4, ena); dea.add(elem); ena[0] = 84; ena[1] = 85; ena[2] = 92; ena[3] = 91; elem = new Non_Linear_Omega_eh(71, 0, 0, 4, ena); dea.add(elem); } gh_on_Gamma_h::gh_on_Gamma_h(int df, Omega_h& omega_h) { // boundary conditions __initialization(df, omega_h); the_gh_array[node_order(0)](1) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(0)][1] = 0.0; the_gh_array[node_order(1)](1) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(1)][1] = 0.0; the_gh_array[node_order(2)](1) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(2)][1] = 0.0; the_gh_array[node_order(3)](1) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(3)][1] = 0.0; the_gh_array[node_order(4)](1) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(4)][1] = 0.0; the_gh_array[node_order(25)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(25)][0] = 0.0; the_gh_array[node_order(29)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(29)][0] = 0.0; the_gh_array[node_order(33)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(33)][0] = 0.0; the_gh_array[node_order(37)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(37)][0] = 0.0; the_gh_array[node_order(44)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(44)][0] = 0.0; the_gh_array[node_order(51)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(51)][0] = 0.0; the_gh_array[node_order(58)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(58)][0] = 0.0; the_gh_array[node_order(65)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(65)][0] = 0.0; the_gh_array[node_order(72)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(72)][0] = 0.0; the_gh_array[node_order(79)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(79)][0] = 0.0; the_gh_array[node_order(86)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(86)][0] = 0.0; // incremental loading (pull apart on top line) for(int j = 86; j <= 92; j++) { the_gh_array[node_order(j)](1) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(j)][1] = 0.0; // initial values; redefined at each time step } } // INSTANTIATION of Global_Discretization static int ndf = 2; static Omega_h oh; static gh_on_Gamma_h gh(ndf, oh); static U_h uh(ndf, oh); static Global_Discretization gd(oh, gh, uh); static U_h delta_uh(ndf, oh); // total increment of displacement of one incremental loading step static gh_on_Gamma_h delta_gh(ndf, oh); static Global_Discretization gd_delta_uh(oh, delta_gh, delta_uh); static U_h d_uh(ndf, oh); // increment of increment of displacement within one incremental loading step static gh_on_Gamma_h d_gh(ndf, oh); static Global_Discretization gd_d_uh(oh, d_gh, d_uh); static double identity_44[4][4] = { {1.0, 0.0, 0.0, 0.0}, {0.0, 1.0, 0.0, 0.0}, {0.0, 0.0, 1.0, 0.0}, {0.0, 0.0, 0.0, 1.0} }; static double identity_33[3][3] = { {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0} }; static double identity_22[2][2] = { {1.0, 0.0}, {0.0, 1.0} }; static double one[4] = {1.0, 1.0, 1.0, 0.0}; static const double i_dev[4][4] = { {(1.0-1.0/3.0), -1.0/3.0, -1.0/3.0, 0.0}, {-1.0/3.0, (1.0-1.0/3.0), -1.0/3.0, 0.0}, {-1.0/3.0, -1.0/3.0, (1.0-1.0/3.0), 0.0}, {0.0, 0.0, 0.0, 1.0} }; static const double i_mu[4][4] = { {1.0, 0.0, 0.0, 0.0}, {0.0, 1.0, 0.0, 0.0}, {0.0, 0.0, 1.0, 0.0}, {0.0, 0.0, 0.0, 0.5} }; C0 I_44 = MATRIX("int, int, const double*", 4, 4, identity_44[0]), I_33 = MATRIX("int, int, const double*", 3, 3, identity_33[0]), I_22 = MATRIX("int, int, const double*", 2, 2, identity_22[0]), One = VECTOR("int, const double*", 4, one), I_dev = MATRIX("int, int, const double*", 4, 4, i_dev[0]), I_mu = MATRIX("int, int, const double*", 4, 4, i_mu[0]); //============================================================================== // Step 2. Element_Formulation-- // quadri-lateral finite deformation plane strain elastoplastic element //============================================================================== // DEFINITION of Element_Formulation class Finite_Deformation_Elastoplastic_B_bar_Q41 : public Element_Formulation { // history data C0 theta, theta_n, // mean dilatation p, p_n; // average pressure H0 j, J_n, // the Jacobian F, F_n, // deformation gradient tensor sigma, sigma_n, // Cauchy stresses F_p_inv, F_p_inv_n, // inverse of plastic deformation gradient tensor ep, ep_n, // total plastic strain alpha, alpha_n; // back stresses (for kinematic hardenning law) void __initialization(int); public: Finite_Deformation_Elastoplastic_B_bar_Q41(Element_Type_Register a) : Element_Formulation(a) {} Element_Formulation *make(int, Global_Discretization&); Finite_Deformation_Elastoplastic_B_bar_Q41(int, Global_Discretization&); }; void Finite_Deformation_Elastoplastic_B_bar_Q41::__initialization(int en) { // ul at t_n+1 ul += gl; // the default procedure behind the scene is the reaction due to fixed B. C.; rhs -= tangent stiffness * d_gl, // this step should be done only in the first iteration of a incremental loading step, // d_gh is then set to zero after the first iteration gl = gd_d_uh.element_fixed_variable(en); // history data Omega_eh *elem = &(oh(en)); theta &= ((Non_Linear_Omega_eh*) elem)->theta(); // at t_n+1 p &= ((Non_Linear_Omega_eh*) elem)->p(); j &= ((Non_Linear_Omega_eh*) elem)->J(); F &= ((Non_Linear_Omega_eh*) elem)->F(); sigma &= ((Non_Linear_Omega_eh*) elem)->sigma(); F_p_inv &= ((Non_Linear_Omega_eh*) elem)->F_p_inv(); ep &= ((Non_Linear_Omega_eh*) elem)->ep(); alpha &= ((Non_Linear_Omega_eh*) elem)->alpha(); theta_n &= ((Non_Linear_Omega_eh*) elem)->theta_n(); // at t_n p_n &= ((Non_Linear_Omega_eh*) elem)->p_n(); J_n &= ((Non_Linear_Omega_eh*) elem)->J_n(); F_n &= ((Non_Linear_Omega_eh*) elem)->F_n(); sigma_n &= ((Non_Linear_Omega_eh*) elem)->sigma_n(); F_p_inv_n &= ((Non_Linear_Omega_eh*) elem)->F_p_inv_n(); ep_n &= ((Non_Linear_Omega_eh*) elem)->ep_n(); alpha_n &= ((Non_Linear_Omega_eh*) elem)->alpha_n(); } Element_Formulation* Finite_Deformation_Elastoplastic_B_bar_Q41::make(int en, Global_Discretization& gd) { return new Finite_Deformation_Elastoplastic_B_bar_Q41(en,gd); } // ad-hoc utility functions static H0 tr(H0& a) { return (a[0][0]+a[1][1]+a[2][2]); } static C0 tr(C0& a) { return (a[0][0]+a[1][1]+a[2][2]); } static H0 acos(H0& a) { // cos^-1 Quadrature& qp = a.quadrature_point(); H0 ret_val(qp); for(int i = 0; i < qp.no_of_quadrature_point(); i++) ret_val.quadrature_point_value(i) = acos( (double)(a.quadrature_point_value(i)) ); return ret_val; } static int new_time_flag = FALSE; Finite_Deformation_Elastoplastic_B_bar_Q41::Finite_Deformation_Elastoplastic_B_bar_Q41(int en, Global_Discretization& gd) : Element_Formulation(en, gd) { // retreive history data from last time step __initialization(en); H1 Z(2, (double*)0, qp), // Natrual Coordinates N = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE( // Shape Functions "int, int, Quadrature", 4/*nen*/, 2/*nsd*/, qp), Zai, Eta; Zai &= Z[0]; Eta &= Z[1]; N[0] = (1.0-Zai)*(1.0-Eta)/4.0; N[1] = (1.0+Zai)*(1.0-Eta)/4.0; N[2] = (1.0+Zai)*(1.0+Eta)/4.0; N[3] = (1.0-Zai)*(1.0+Eta)/4.0; //=========================================================================== // // Simo, Taylor and Pister, 1985, Comp. Meth. in Appl. Mech. and Engrg. // //=========================================================================== // 1. elastic response is exact, incremental objective integration scheme is // avoid entirely. This is achieved by multiplicative slpit of elastic and // plastic deformation gradient tensor. // 2. current configuration (phi) and its coordinates (x) are defined at // time n+1, iteration i+1, from the x-coordinate the derivatives of // shape function are evaluated // Consequences of 2: // ---In elastic predictor phase, no plastic flow is assumed, therefore, // plastic internal variables remain unchanged. // ---In plastic corrector phase, the process occurs at the "fixed current // configuration" (phi_n+1 does not change, d := d^e + d^p = 0), therefore, // the integration scheme reduces to be identical to infinitesimal case // ---the objectivity is followed directly from the objectivity of the // constitutive law at the current configuration // Pro-and-con of 2, from computational point of view: // Advantage: plastic internal variables integration scheme reduces to the // infinitesimal case, no push-forward/pull-back is necessary // Disadvantage: at the end of each time step, need to compute intermediate // configuration with (F_^p)^-1 @ (n, i+1) from history stress // data, this step is very complicated and costy // // N.B. Simo, 1988, same Journal, use phi defines at time_n, objectivity // requires push-forward/pull-back of left Cauchy-Green tensor and back // stresses with deformation gradient of incremental displacement, no // inverted stresses for computing inverse of plastic deformation // gradient tensor (Step 3 below) is necessary. // //=========================================================================== //=========================================================================== // coordinate transformation // Physical Coordinates; // x is the current configuration @(time n+1, iteration i+1) C0 u(nen, ndf, (double*)0); for(int i = 0; i < nen; i ++) for(int j = 0; j < ndf; j++) u[i][j] = ul[i*ndf+j]; H1 x = N*(xl+u); // N.B. d(N)/d(X) = d(N)/d(Z)*(d(X)/d(Z))^-1 // => {nen x nsd}*{nsd x nsd}={nen x nsd} H0 Nx = d(N) * d(x).inverse(); // Weighting derivative, J dv(d(x).det()); // the Jacobian //=========================================================================== //=========================================================================== // standard B-matrix formulation // dim ~Nx = {2x4}, dim ~wx = {{2x1}x4} H0 w_x = INTEGRABLE_SUBMATRIX("int, int, H0&", 1, 2/*nsd*/, Nx), wx, wy, B, B_bar; wx &= w_x[0][0]; wy &= w_x[0][1]; C0 zero(0.0); Matrix_Representation::axisymmetric_flag = TRUE; // flag for B-matrix expressed in axisymmetric compatible format (4x2) H0 N_r(4, (double*)0, qp); N_r = 0.0; // if plain strian //H0 N_r &= ~((H0)N)/(H0)x; // if axisymmetric B &= (~wx || zero ) & // axisymmetric compatible format (zero || ~wy ) & // dim B = {4x8}, where dim ~wx[][] = {1x4} (~N_r || zero ) & (~wy || ~wx ); //=========================================================================== //=========================================================================== // B_bar formulation--with mean dilatation approximation for Q4/1 element // pressure-volumetric strain shape function H0 N_vp(qp); N_vp = 1.0; H0 mN_vp = N_vp & N_vp & N_vp & zero; // 4x1; axisymmetric compatible format C0 C = ((~mN_vp)*B) | dv, // 1x4 x 4x2 = 1x2; or C = (~One*B) | dv E = (N_vp) | dv, // 1; area of the element W = C/E; // 1x2; average shape derivatives H0 B_vol = (1.0/3.0) * (mN_vp % W[0]), // 4x1 x 1x2; for computing pressure B_dev = I_dev * B; B_bar &= B_dev + B_vol; //=========================================================================== // Finite Deformation Kinematics============================================= // // Multiplicative split of elasto-plastic deformation: F = F^e * F^p // Multiplicative kinematic split: F_hat = J^(-1/3) * F, where J := det F // // Step 0: compute deformation gradient tensor and the Jacobian at t_n+1, // iteration i+1 (F_i+1, J_i+1) from u; Note J_i+1 = det(F_i+1) // Remark 1: F = dx/dX = d(X+u)/dX = I + du/dx dx/dX = I + (grad u) * F // post-multiply F^-1 on both sides => F F^-1 = F^-1 + (grad u)*F*F^-1 // => F^-1 = I - grad u, therefore get F^-1, then get F from F^-1 // 0.1 compute F_i+1^-1 from F_i+1^-1 = I - grad u // 0.2 compute F_i+1 and det(F_i+1) from F_i+1^-1 //=========================================================================== { // Kinematics of finite deformation // 0.1 compute F_i+1^-1 from F_i+1^-1 = I - grad u H0 F_inv = I_22 - (~u)*Nx; // F_inv @ (time n+1, iteration i+1); // 0.2 compute F_i+1 and det(F_i+1) from F_i+1^-1 F = F_inv.inverse(); // F @ (time n+1, iteration i+1) j = F.det(); // J @ (time n+1, iteration i+1) } // average dilatation, volumetric store-energy function, and mean pressure H0 unit(qp); unit = 1.0; C0 VOL = (1.0/j) | dv, // volume at material configuration vol = unit | dv, // volume at current configuration theta = vol / VOL; // mean dilatation; volume change averaged over the element // volumetric store energy function: U(J) = 1/2 K log(J)^2; C2 THETA(theta), // mean dilatation; re-declared to be twice differentiable U = K_*((log(THETA)).pow(2))/2.0; // volumetric store energy function; restricted to be a poly-convex function p = d(U); // mean pressure; p = U' = (dU/d J); N.B. U'' is used in computing tangent stiffness matrix H0 J_two_third = exp(log(j)*2.0/3.0); // J^(2/3) //=========================================================================== //=========================================================================== // // Step 3: compute (F_n^p)^-1 @ (n, i+1) // //=========================================================================== // N.B. "Step 3" is activated only at the beginning of a time step, instead // of only at the end of the time step as described in Simo et al., 1985 // since Step 3 at the last time step does not need to be computed. // // if new time step (n+1, i = 0) // Note (F_n^p)^-1 @ (n, i+1) is un-available at the beginning of t_n+1 // compute (F_n^p)^-1 @ (n, i+1) from F and sigma @ (n, i+1) by // (F_n^p)^-1 = F_n^e * F_n^-1; // (F_i^p)^-1 @ (n+1, i = 0) := (F_n^p)^-1 @ (n, i+1) // else if a new iteration (i+1) within the current time step (n+1, i+1) // skip Step 3, (F_i^p)^-1 @ (n+1, i) is retreived from the history data // // @ (n+1, i = 0) // 3.a1 compute dev b^hat_n^e from stress history data of sigma_n @ (n, i+1) // Remark 3.a1: tau_n = J sigma_n, // dev b^hat_n^e = dev tau_n / mu // where sigma_n: Cauchy stresses from the history data @ (n, i+1) // tau: Kirchhoff stresses // b: left Cauchy-Green tensor // 3.a2 compute (b^hat_n^e) from (dev b^hat_n^e), requring det(b^hat_n^e) = 1; // i.e., add trace to an uni-modular deviatoric tensor // Remark 3.a2: For those who are familiar with solid mechanics, // The solution for the above problem is the same // as the solution for "principle deviatoric stresses" // see Y.C. Fung, "Foundation of Solid Mechanics", p. 80-81, // Problem: b_hat^e = dev b_hat^e + theta * 1, where // theta = 1/3 tr b_hat^e , the volumetric part // Objective: compute theta then add to diagonal of deviator // to obtain b_hat^e, (requring det b_hat^e = 1) // Lemma: compute volumetric part giving deviatoric part and // determinant = 1 // in plain strain // | b'11 b'12 | // dev b_hat^e = | | // | b'21 b'22 | // from C.Y. Fung Eqn. 7, with dilatation (theta_0) in place of pressure (sigma_0) // theta_0^3 - J2 * theta_0 + (J3-I3) = 0 -------------(**) // by definitions // I3 := det b_hat^e = 1 as given // 3rd invariant of a tensor // J3 := (b'11*b'22 - b'12^2)*b'33 // 3rd invariant of the deviatoric part of a tensor // J2 := 1/2 (b'11^2 + b'22^2 + b'33^2) + b'12^2 // 2nd invariant of the deviatoric part of a tensor // => J3 = (b'11*b'22 - (J2-1/2 (b'11^2+b'22^2+b'33^2)) )*b'33 // = (1/2 (b'11+b'22^2+b'33)^2 -b'11*b'33-b'22*b'33- J2)*b'33 // = (b'33*b'33-J2)*b'33 // N.B. J1 := b'11 + b'22 + b'33 = 0; // 1st invaraint of the deviatoric part of a tensor // b'23 = b'13 = 0 // i.e., J3 = (b'33^2-J2)*b'33 // substitute J3 = (b'33^2-J2)*b'33 into (**) // theta_0^3 - J2*theta_0 - ((J2-b'33^2)*b'33+1) = 0 // define J3' = (J2-b'33^2)*b'33+1 // theta_0^3 - J2*theta_0 - J3' = 0 ---------------------(***) // one has a cubic algebraic equation that is similar // to the one for solving principle deviatric stresses in C.Y. Fung // // cos 3 alpha = J'3 * sqrt(2)/b'o^3, where J2 = 3/2 b'o^2 // => theta_0 = sqrt(2) b'0 cos alpha // // 3.b Polar Decomposition; compute b_n^e from b^hat_n^e, then V_n^e from b_n^e, // then set F_n^e = V_n^e, where V is the left strecth tensor // b_n^e = J^2/3 * b^hat_n^e // Remark 3.b: see e.g. Marsden and Hughes, 1983, "Mathematical // Foundations of Elasticity", p. 55 // from Cayley-Hamilton theorem in linear algebra // V^2 - Iv * V + IIv * 1 = 0 // Ib = trace(b), IIb = det(b) // V = (b+sqrt(IIb)*1)/sqrt(Ib+2sqrt(IIb)) // // 3.c compute (F_n^p)^-1 = F_n^e * F_n^-1, // where F_n(n, i+1) is from history data //=========================================================================== if(new_time_flag) { // flag to switch on Step 3 only at a new time step // Step 3: procedure to find the "intermediate" configuration // 3.a1 invert stress to get strain from constitutive law H0 dev_b_hat_n_e = INTEGRABLE_MATRIX("int, int, Quadrature", 3, 3, qp), tau_n = J_n*sigma_n, // Kirchhoff stresses compute from history "the Jacobian" and history Cauchy stresses dev_tau_n = INTEGRABLE_MATRIX("int, int, Quadrature", 3, 3, qp); // deviatoric Kirchhoff stresses H0 tr_tau_n = tau_n[0]+tau_n[1]+tau_n[2]; dev_tau_n[0][0] = tau_n[0] - (1.0/3.0)*tr_tau_n; dev_tau_n[1][1] = tau_n[1] - (1.0/3.0)*tr_tau_n; dev_tau_n[2][2] = tau_n[2] - (1.0/3.0)*tr_tau_n; dev_tau_n[0][1] = dev_tau_n[1][0] = tau_n[3]; dev_b_hat_n_e = dev_tau_n / mu_; // deviatoric part of the constitutive law // 3.a2 add one third of trace to an uni-modular deviatoric tensor H0 bp = dev_b_hat_n_e, // alias; b_prime := dev b_hat_n_e J2 = (bp[0][0].pow(2)+bp[1][1].pow(2)+bp[2][2].pow(2))/2.0 + // 2nd invariant of a deviatoric tensor bp[0][1].pow(2); H0 theta_0(qp); #if defined(__NUMERICAL_ROOT_FINDING) for(int i = 0; i < qp.no_of_quadrature_point(); i++) { C0 J2_q = J2.quadrature_point_value(i), theta_0_q = theta_0.quadrature_point_value(i), bp_q = bp.quadrature_point_value(i), J3_prime = (J2_q-bp_q[2][2].pow(2))*bp_q[2][2] + 1.0; // solution for a cubic algebraic equation: // // theta_0^3 - J2*theta_0 - J3' = 0; // // where theta_0 is "one third of trace"; // i.e., theta_0 := 1/3 tr(b_hat_n_e) int MAX_ITERATION_NO = 50, count = 0; double EPSILON = 1.e-6; C1 X(1.0), // X_0 = 1.0; initial value f; C0 d_X(0.0); // Root Finding Algorithm: Taylor's series expasion to 1st-order // f(X_0+d_X) = f(X_0)+ f'(X_0)* d_X + O(X^2); // find root for f(X) = 0, one get d_X = -f(X_0)/f'(X_0) // where X = X_0 + d_X do { f &= X.pow(3) - J2_q*X - J3_prime; d_X = -((C0)f)/d(f); ((C0)X) += d_X; } while((double)norm(d_X) > EPSILON && ++count < MAX_ITERATION_NO); if(count == MAX_ITERATION_NO) ofs << "Warning: No convergence achieved for solution of a cubic algebraic equation!" << endl; theta_0_q = ((C0)X); } #else // Closed form for(int i = 0; i < qp.no_of_quadrature_point(); i++) { C0 J2_q = J2.quadrature_point_value(i), theta_0_q = theta_0.quadrature_point_value(i); if((double)J2_q < 1.e-6) { theta_0_q = 1.0; } else { C0 bp_q = bp.quadrature_point_value(i), J3_prime = (J2_q-bp_q[2][2].pow(2))*bp_q[2][2] + 1.0, // from definition above bp_oct = sqrt((2.0/3.0)*J2_q); // equivlent to octahedral stresses in principle dev stresses soln. C0 temp = (-J2_q/3.0).pow(3)+(J3_prime/2.0).pow(2); if((double)temp > 0.0) { // if cos^2 (3 alpha) > 1.0 double arg1 = (double)(J3_prime/2.0+sqrt(temp)), arg2 = (double)(J3_prime/2.0-sqrt(temp)), a1, a2; if(fabs(arg1) < 1.e-6) a1 = 0.0; else a1 = ((arg1>=0.0) ? 1.0 : -1.0) * exp(logl(fabs(arg1))/3.0); if(fabs(arg2) < 1.e-6) a2 = 0.0; else a2 = ((arg2>=0.0) ? 1.0 : -1.0) * exp(logl(fabs(arg2))/3.0); theta_0_q = a1 + a2; } else { C0 cos_3alpha = J3_prime * sqrt(2.0) / bp_oct.pow(3), alpha = acos(cos_3alpha)/3.0; theta_0_q = sqrt(2.0) * bp_oct * cos(alpha); } } } #endif // add "one-third of trace" (theta_0) to an uni-modular deviatoric tensor (dev b_hat_n_e) { H0 b_hat_n_e = dev_b_hat_n_e + theta_0 * I_33; // 3.b compute b_n^e from b^hat_n^e, then V_n^e from b_n^e H0 b_n_e = J_two_third * b_hat_n_e; // compute V_n^e from b_n^e // Polar Decomposition; use 2-D Caley-Hamilton principle to compute elastic left stretch tensor H0 V_n_e = INTEGRABLE_MATRIX("int, int, Quadrature", 2, 2, qp); H0 I_b = b_n_e[0][0]+b_n_e[1][1], // 2x2 trace II_b =b_n_e[0][0]*b_n_e[1][1]-b_n_e[0][1]*b_n_e[1][0]; // 2x2 determinant; V_n_e[0][0] = (b_n_e[0][0]+sqrt(II_b)) / sqrt(I_b+2.0*sqrt(II_b)); V_n_e[1][1] = (b_n_e[1][1]+sqrt(II_b)) / sqrt(I_b+2.0*sqrt(II_b)); V_n_e[0][1] = b_n_e[0][1] / sqrt(I_b+2.0*sqrt(II_b)); V_n_e[1][0] = b_n_e[1][0] / sqrt(I_b+2.0*sqrt(II_b)); //V_n_e[2][2] = sqrt((double)b_n_e[2][2]); // this component does not need to be computed // elastic deformation gradient H0 F_n_e = V_n_e; // 3.c compute (F_n^p)^-1 = F_n^e * F_n^-1 F_p_inv_n = F_n_e * F_n.inverse(); // inverse of plastic deformation gradient tensor in history data } }// end if new time step //=========================================================================== // Step 1: Elastic Predictor // compute "trial" TR-F^e @ (n+1, i+1), then TR-b_hat^e, // assuming no plastic flow occurs for the trial step // => F_n+1^p := F_n^p // TR-F^e@(n+1, i+1) = F@(n+1, i+1) * (F_n+1^p)^-1 // = F@(n+1, i+1) * (F_n^p)^-1 // // 1.a1 compute TR-F^e@(n+1, i+1) = F@(n+1, i+1) * (F_n^p)^-1 // 1.a2 TR-b_hat^e@(n+1, i+1) // = J^(-2/3)@(n+1, i+1) * TR-b^e@(n+1, i+1) // = J^(-2/3)@(n+1, i+1) * (TR-F^e * (TR-F^e)^T)@(n+1, i+1) // 1.b1 compute dev TR-b_hat^e@(n+1, i+1) // 1.b2 dev TR-tau@(n+1, i+1) = mu dev TR-b_hat^e@(n+1, i+1) //=========================================================================== H0 s(4, (double*)0, qp); // use engineering convention; // s := {s_xx, s_yy, s_zz, s_xy} = "trial elastic" deviatoric Kirchhoff stresses H0 mu_bar; // parameter for computing deviatoric tangent moduli { F_p_inv = F_p_inv_n; // assume no plastic flow occurs at trial elastic step H0 TR_F_e = F * F_p_inv, // trial elastic "deformation gradient tensor" TR_b_e = INTEGRABLE_MATRIX("int, int, Quadrature", 3, 3, qp), // trial elastic "left Cauchy-Green tensor" ref_TR_b_e =INTEGRABLE_MATRIX("int, int, H0&, int, int, Quadrature", 2, 2, TR_b_e, 0, 0, qp); ref_TR_b_e = TR_F_e *(~TR_F_e); // b^e = F^e*(F^e)^t TR_b_e[2][2] = 1.0; // plain strain; the 3rd-Dimension dx/dX = 1.0 H0 TR_b_hat_e = TR_b_e / J_two_third, // trial elastic "modified" left Cauchy-Green tensor dev_TR_b_hat_e = TR_b_hat_e - (1.0/3.0) * tr(TR_b_hat_e) * I_33, // trial elastic "deviatoric" modified left Cauchy-Green tensor dev_TR_tau = mu_ * dev_TR_b_hat_e; // trial elastic "deviatoric Kirchhoff stresses" mu_bar &= (mu_/3.0) * tr(TR_b_hat_e); // convert to engineering convention s[0] = dev_TR_tau[0][0]; // s_xx s[1] = dev_TR_tau[1][1]; // s_yy s[2] = dev_TR_tau[2][2]; // s_zz s[3] = dev_TR_tau[0][1]; // s_xy } // deviatoric tangent moduli; 4x4 H0 a_dev(4, 4, (double*)0, qp), // trial deviatoric tangent moduli; 4x4 TR_a_dev = 2.0*mu_bar*(I_mu - (One%One)/3.0) - (2.0/3.0)*(s%One+One%s); //=========================================================================== // Elastoplastic material model; //=========================================================================== // material constants: Y, d_kapa, d_H_alpha; initial alpha(back-stresses) = 0, // where, kapa = Y+d_kapa*ep is the "linear isotropic // hardenning law" and, alpha_n+1 = alpha_n + // sqrt(2/3) d_H_alpha lambda is the "linear kinematic // hardenning law" // from t_n(store in memory): s_n(deviatoric stresses), e^p sub n(plastic strain), // alpha_n(back_stresses) // compute for t_n+1: yield ratio, sigma(for residual), tangent operator(for // stiffness) //=========================================================================== H0 ZAI = s - alpha_n, // (2.4);(i-2) kinematic hardenning law ZAI_norm(qp); #define qv quadrature_point_value for(i = 0; i < qp.no_of_quadrature_point(); i++) { // aliases for quadrature point value C0 ZAI_0_q, ZAI_1_q, ZAI_2_q, ZAI_3_q; ZAI_0_q &= ZAI[0].qv(i); ZAI_1_q &= ZAI[1].qv(i); ZAI_2_q &= ZAI[2].qv(i); ZAI_3_q &= ZAI[3].qv(i); ZAI_norm.qv(i) = sqrt((double)(ZAI_0_q.pow(2)+ZAI_1_q.pow(2)+ZAI_2_q.pow(2)+2.0*ZAI_3_q.pow(2))); } // isotropic hardenning law for kapa(e^p) H0 KAPA(qp); KAPA = Y_infinity - (Y_infinity - Y_0) * exp(-iota*ep); H0 radius = sqrt(2.0/3.0)*KAPA, // R = sqrt(2/3) kapa(e^p) yield_ratio = ZAI_norm/radius; H0 lambda(qp); //=========================================================================== // solve plastic flow law with "local" Newton iteration on each quadrature point //=========================================================================== for(i = 0; i < qp.no_of_quadrature_point(); i++) { // aliases for quadrature point value------------------------------------- C0 s_q, ep_q, ep_n_q, J_q, alpha_q, alpha_n_q, ZAI_q, ZAI_norm_q, sigma_q, a_dev_q, TR_a_dev_q, mu_bar_q, B_bar_q, B_vol_q, radius_q, lambda_q, yield_ratio_q; s_q &= s.qv(i); ep_q &= ep.qv(i); ep_n_q = ep_n.qv(i); J_q = j.qv(i); // history terms alpha_q &= alpha.qv(i); alpha_n_q = alpha_n.qv(i); ZAI_q &= ZAI.qv(i); ZAI_norm_q &= ZAI_norm.qv(i); // stress deviator sigma_q &= sigma.qv(i); // internal stresses and tangent moduli a_dev_q &= a_dev.qv(i); TR_a_dev_q = TR_a_dev.qv(i); mu_bar_q &= mu_bar.qv(i); B_bar_q = B_bar.qv(i); B_vol_q = B_vol.qv(i); // B matrices radius_q &= radius.qv(i); lambda_q &= lambda.qv(i); // yield function parameters yield_ratio_q &= yield_ratio.qv(i); //------------------------------------------------------------------------ if(ZAI_norm_q > radius_q) { // plastic state //=========================================================================== // Step 2: Plastic Corrector is exactly the same as in the infinitesimal case // N.B. yield funtion phi = || dev tau || - sqrt(2/3) * kapa(ep) //=========================================================================== ofs << "***Plastic state: element # " << en << ", quadrature point # " << i << ", initial yield ratio: " << yield_ratio_q << endl; const double EPSILON = 1.e-20; const int MAX_ITERATION_NO = 50; double d_lambda_norm, d_lambda_norm_0; C0 d_lambda(0.0); int count = 0; C1 phi, // yield function, phi(LAMBDA; ...) LAMBDA(0.0); // plastic consistenet parameter (gamma times delta_t) do { // isotropic hardenning law C1 EP = ep_n_q+sqrt(2.0/3.0)*LAMBDA, KAPA = Y_infinity - (Y_infinity - Y_0)*exp(-iota*EP), RADIUS = sqrt(2.0/3.0)*KAPA, // linear kinematic hardenning law delta_H_alpha = sqrt(2.0/3.0)*d_H_alpha*LAMBDA; // phi(lambda), the yield function; (3.8) phi &= -RADIUS + ZAI_norm_q - (2.0*mu_*(LAMBDA) + sqrt(2.0/3.0)*delta_H_alpha); d_lambda = - ((C0)phi)/d(phi); // Newton's formula; or Taylor's series expansion up to 1st-order ((C0) LAMBDA) += d_lambda; // update with increment d_lambda_norm = norm(d_lambda); if(count == 0) d_lambda_norm_0 = d_lambda_norm; } while(d_lambda_norm > (d_lambda_norm_0*EPSILON) && count++ < MAX_ITERATION_NO); // until converge lambda_q = (C0) LAMBDA; if(count == MAX_ITERATION_NO) { ofs << " Warning: local Newton iteration failed to converge!" << endl << " element # " << en << ", quadrature point # " << i << endl; } // update history terms================================================ ep_q = ep_n_q + sqrt(2.0/3.0)*lambda_q; // (2.4) or (3.4);(iii), integrate the total plastic strain C0 n = ZAI_q/ZAI_norm_q; // (3.2); (ii) unit vector normal to the yield surface double delta_H_alpha = sqrt(2.0/3.0)*d_H_alpha*(double)lambda_q; alpha_q = alpha_n_q + sqrt(2.0/3.0)*delta_H_alpha*n; // (3.5); (iv-1), use linear kinematic hardenning law C1 EP((double)ep_q), KAPA = Y_infinity - (Y_infinity - Y_0)*exp(-iota*EP); double factor = 1.0+1.e-8; // force the stress state to stay slightly outside the yield surface s_q = alpha_q + (sqrt(2.0/3.0)*((double)((C0)KAPA))*factor)*n; // (3.3); (iv-2), //===================================================================== //===================================================================== // consistent tangent elastoplastic moduli //===================================================================== // dev n^2; // N.B. n is expressed in engineering convention all vector space operation // that is defined in mathematical convention doesn't apply C0 n2(4, (double*)0); // n^2 n2[0] = n[0]*n[0]+n[3]*n[3]; n2[1] = n[3]*n[3]+n[1]*n[1]; n2[2] = n[2]*n[2]; n2[3] = n[3]*(n[0]+n[1]); C0 trn2 = n2[0]+n2[1]+n2[2]; // tr(n2) C0 dev_n2(4, (double*)0); dev_n2 = n2; for(int i = 0; i < 3; i++) dev_n2[i] -= trn2/3.0; // scaling factors; // Simo, 1988, p.7, without "pull-back" and "push forwrad" of back stresses C0 mu_2bar = mu_bar_q - (alpha_n_q[0]+alpha_n_q[1]+alpha_n_q[2])/3.0, f_0 = 2.0*mu_bar_q*lambda_q / ZAI_norm_q, delta_0 = (1.0+d_H_alpha/3.0/mu_+d(KAPA)/3.0/mu_2bar), f_1 = (1.0/delta_0-f_0), delta_1 = f_1*2.0*mu_2bar - (1.0/delta_0*(1.0+d_H_alpha/3.0/mu_)-1.0)*(4.0/3.0)*lambda_q, delta_2 = 2.0*ZAI_norm_q*f_1; // trial deviatoric hardenning moduli; 3x3 C0 TR_h_dev = 2.0 * mu_2bar * (I_mu - (One%One)/3.0) - (2.0/3.0)*(ZAI_q%One+One%ZAI_q); a_dev_q = TR_a_dev_q - f_0 * TR_h_dev - delta_1*(n%n) - (delta_2/2.0)*(n%dev_n2+dev_n2%n); } else { // elastic state a_dev_q = TR_a_dev_q; ofs << " Elastic state: element # " << en << ", quadrature point # " << i << ", yield ratio: " << yield_ratio_q << endl; } // end if plastic or elastic state // compute stresses at t_n+1 (v); deviatoric stresses plus elastic volumetric stresses sigma_q = s_q/J_q + p*One; // (v); deviatoric Cauchy stresses plus mean pressure } // end for each quadrature point // post-processing code for reaction, stresses, and strain if(Matrix_Representation::Assembly_Switch == Matrix_Representation::REACTION) { C0 stiff_geometrical, stiff_material; // A. geometrical part (or the so-called initial stress part, also see Zienkewicz and Taylor 1991) // grad N : (sigma grad N) { C0 e = BASIS("int", ndf), // basis vectors E = BASIS("int", nen), U = (e%e)*(E%E); H0 fact = wx*(sigma[0]*~wx+sigma[3]*~wy) + // grad N : (sigma grad N) wy*(sigma[3]*~wx+sigma[1]*~wy); stiff_geometrical &= (+(fact*U[0][0]+ fact*U[1][1])) |dv; // 8x8 } // B. material part { C0 dN_bar = (( (~wx) || (~wy) ) | dv) / vol; // 3.9b average shape function derivatives; 1x8 stiff_material &= ( ((~B) * (a_dev/j + p*((One%One) - 2.0*I_mu)) * B) | dv ) // 3.11; retain B-matrix formulation + dd(U) * theta * ((~dN_bar)*dN_bar) * vol; // J^2 U''/J } stiff &= stiff_geometrical + stiff_material; // 8x8 the_element_nodal_value &= stiff * ul; // reaction due to free and fixed degree of freedom } else if(Matrix_Representation::Assembly_Switch == Matrix_Representation::STRAIN) { H0 epsilon_hat = INTEGRABLE_VECTOR("int, Quadrature", 3, qp); epsilon_hat = 0.0; for(int i = 0; i < nen; i++) epsilon_hat += B_bar(i*ndf)*ul[i*ndf] + B_bar(i*ndf+1)*ul[i*ndf+1]; int nqp = qp.no_of_quadrature_point(); ofs.setf(ios::left,ios::adjustfield); ofs << setw(9) << " elem #, " << setw(14) << "x-coor.," << setw(14) << "y-coor.," << setw(14) << "epsilon-x," << setw(14) << "epsilon-y," << setw(14) << "epsilon-xy" << endl; for(i = 0; i < nqp; i++) { ofs << setw(9) << en << setw(14) << ((H0)x[0]).quadrature_point_value(i) << setw(14) << ((H0)x[1]).quadrature_point_value(i) << setw(14) << (epsilon_hat[0].quadrature_point_value(i)) << setw(14) << (epsilon_hat[1].quadrature_point_value(i)) << setw(14) << (epsilon_hat[2].quadrature_point_value(i)) << endl; } } else if(Matrix_Representation::Assembly_Switch == Matrix_Representation::NODAL_STRAIN) { int strain_no; if(!Matrix_Representation::axisymmetric_flag) strain_no = (ndf+1)*ndf/2; else strain_no = 4; the_element_nodal_value &= C0(nen*strain_no, (double*)0); // returned for global assembly C0 projected_nodal_strain = SUBVECTOR("int, C0&", strain_no, the_element_nodal_value); H0 epsilon_hat = INTEGRABLE_VECTOR("int, Quadrature", strain_no, qp); epsilon_hat = 0.0; for(int i = 0; i < nen; i++) epsilon_hat += B_bar(i*ndf)*ul[i*ndf] + B_bar(i*ndf+1)*ul[i*ndf+1]; // Least Sqaures Method = Bubnov-Galerkin Method for Q8 interpolation to Smooth stresses // lumped mass procedure; see Zienkiewicz and Taylor, vol.1, p.607, Eq.A8.5 and A8.6 H0 unit(qp); unit = 1.0; C0 area = unit | dv; C0 lumped_mass = VECTOR("int", nen); C0 sum_of_lumped_mass(0.0); for(i = 0; i < nen; i++) { lumped_mass[i] = (((H0)N[i]).pow(2)) | dv; sum_of_lumped_mass += lumped_mass[i]; } C0 normalized_factor = area/sum_of_lumped_mass; for(i = 0; i < nen; i++) { lumped_mass[i] *= normalized_factor; projected_nodal_strain(i) = ( (((H0)N[i])*epsilon_hat) | dv ) / lumped_mass[i]; } } else if(Matrix_Representation::Assembly_Switch == Matrix_Representation::STRESS) { /* H0 sigma_hat = INTEGRABLE_VECTOR("int, Quadrature", 3, qp); sigma_hat = 0.0; H0 Dt = a_dev/j + p*(One%One-2.0*I_44); H0 DtB_bar = Dt*B_bar; // 3x3 x 3x8 = 3x8 for(int i = 0; i < nen; i++) sigma_hat += DtB_bar(i*ndf)*ul[i*ndf] + DtB_bar(i*ndf+1)*ul[i*ndf+1]; */ int nqp = qp.no_of_quadrature_point(); ofs.setf(ios::left,ios::adjustfield); ofs << setw(9) << " elem #, " << setw(14) << "x-coor.," << setw(14) << "y-coor.," << setw(14) << "sigma-x," << setw(14) << "sigma-y," << setw(14) << "sigma-xy" << endl; for(int i = 0; i < nqp; i++) { ofs << setw(9) << en << setw(14) << ((H0)x[0]).quadrature_point_value(i) << setw(14) << ((H0)x[1]).quadrature_point_value(i) << setw(14) << (sigma[0].quadrature_point_value(i)) << setw(14) << (sigma[1].quadrature_point_value(i)) << setw(14) << (sigma[2].quadrature_point_value(i)) << endl; } } else if(Matrix_Representation::Assembly_Switch == Matrix_Representation::NODAL_STRESS) { int stress_no; if(!Matrix_Representation::axisymmetric_flag) stress_no = (ndf+1)*ndf/2; else stress_no = 4; the_element_nodal_value &= C0(nen*stress_no, (double*)0); // returned for global assembly C0 projected_nodal_stress = SUBVECTOR("int, C0&", stress_no, the_element_nodal_value); /* H0 sigma_hat = INTEGRABLE_VECTOR("int, Quadrature", stress_no, qp); sigma_hat = 0.0; H0 Dt = a_dev/j + p*((One%One)-2.0*I_44); H0 DtB_bar = Dt*B_bar; // 3x3 x 3x8 = 3x8 for(int i = 0; i < nen; i++) sigma_hat += DtB_bar(i*ndf)*ul[i*ndf] + DtB_bar(i*ndf+1)*ul[i*ndf+1]; */ // Least Sqaures Method = Bubnov-Galerkin Method for Q8 interpolation to Smooth stresses // lumped mass procedure; see Zienkiewicz and Taylor, vol.1, p.607, Eq.A8.5 and A8.6 H0 unit(qp); unit = 1.0; C0 area = unit | dv; C0 lumped_mass = VECTOR("int", nen); C0 sum_of_lumped_mass(0.0); for(int i = 0; i < nen; i++) { lumped_mass[i] = (((H0)N[i]).pow(2)) | dv; sum_of_lumped_mass += lumped_mass[i]; } C0 normalized_factor = area/sum_of_lumped_mass; for(i = 0; i < nen; i++) { lumped_mass[i] *= normalized_factor; projected_nodal_stress(i) = ( (((H0)N[i])*sigma) | dv ) / lumped_mass[i]; } } else if(Matrix_Representation::Assembly_Switch == Matrix_Representation::NODAL_SCALAR) { the_element_nodal_value &= C0(nen, (double*)0); // returned for global assembly // Least Sqaures Method = Bubnov-Galerkin Method for Q8 interpolation to Smooth stresses // lumped mass procedure; see Zienkiewicz and Taylor, vol.1, p.607, Eq.A8.5 and A8.6 H0 unit(qp); unit = 1.0; C0 area = unit | dv; C0 lumped_mass = VECTOR("int", nen); C0 sum_of_lumped_mass(0.0); for(int i = 0; i < nen; i++) { lumped_mass[i] = (((H0)N[i]).pow(2)) | dv; sum_of_lumped_mass += lumped_mass[i]; } C0 normalized_factor = area/sum_of_lumped_mass; for(i = 0; i < nen; i++) { lumped_mass[i] *= normalized_factor; the_element_nodal_value[i] = ( (((H0)N[i])*yield_ratio) | dv ) / lumped_mass[i]; } } else { // if not post-processing, compute tangent stiffness and form residual vector // non-linear mixed formulation (3-field Hu-Washizu variational principle) tangent stiffness matrix // for tangent stiffness matrix, B_bar formulation does not carry over to geometrical // non-linear case (see Simo et al. 1985); additional terms occurs that // only make the formulation even more complicated. Follow Simo(1988) formula in that // the B matrix formulation is retained. The formulation is much more simple. C0 stiff_geometrical, stiff_material; // A. geometrical part (or the so-called initial stress part, see also Zienkewicz and Taylor 1991) // grad N : (sigma grad N) { C0 e = BASIS("int", ndf), // basis vectors E = BASIS("int", nen), U = (e%e)*(E%E); H0 fact = wx*(sigma[0]*~wx+sigma[3]*~wy) + // grad N : (sigma grad N) wy*(sigma[3]*~wx+sigma[1]*~wy); stiff_geometrical &= (+(fact*U[0][0]+ fact*U[1][1])) |dv; // 8x8 } // B. material part { C0 dN_bar = (( (~wx) || (~wy) ) | dv) / vol; // 3.9b average shape function derivatives; 1x8 stiff_material &= ( ((~B) * (a_dev/j + p*(/*(One%One)*/ - 2.0*I_mu)) * B) | dv ) // 3.11; retain B-matrix formulation + dd(U) * theta * ((~dN_bar)*dN_bar) * vol; // J^2 U''/J } stiff &= stiff_geometrical + stiff_material; // 8x8 // compute the residual (non-linear problem only); internal stress divergence term force &= -((~B)*sigma)|dv; } } // REGISTRATION of Element(s) Element_Formulation* Element_Formulation::type_list = 0; Element_Type_Register element_type_register_instance; static Finite_Deformation_Elastoplastic_B_bar_Q41 finite_deformation_elastoplastic_b_bar_q41_instance(element_type_register_instance); // element type # 0 //============================================================================== // Step 3: Matrix Formulation and Solution Phase //============================================================================== void femap_output_model(Global_Discretization&, ofstream&); void femap_output_data(int, Global_Discretization&, ofstream&, C0&, C0&, C0&); int main() { time_t time0 = time(0); Matrix_Representation mr(gd); delta_uh.matrix_representation() = &mr; d_uh.matrix_representation() = &mr; time_t time1 = time(0); ofs << "elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; C0 delta_u, // total increment of current increment loading step du; // increment of one Newton iteration step; delta_u = sumation of du ofstream ofs1("D:\\FEMAP601\\finite.neu", ios::out | ios::trunc); femap_output_model(gd, ofs1); // incremental loading (pull apart on top line) for(int i = 0; i < 60; i++) { // incremental loading loop if(i != 0) new_time_flag = TRUE; // flag to turn on Step 3 of the finite deformation procedure //======================================================================== //"global" Newton-Raphson Method for elasto-plastic minimun energy solution //======================================================================== int count = 0; const int MAX_ITERATION_NO = 40; const double EPSILON = 1.e-9; double residual_norm, energy, energy_0; #ifdef __GOLDEN_SECTION_LINE_SEARCH double energy_old; #endif //=initial state of a incremental loading step ======================== if(i != 0) { delta_u = 0.0; // initial free total increment of displacement set to zero delta_uh = delta_u; } for(int j = 86; j <= 92; j++) { d_gh[j][1] = 0.1; // displacement increment on top delta_gh[j][1] = 0.0; } d_uh = d_gh; delta_uh = delta_gh; uh = gh; //===================================================================== do { ofs << (i+1) << "-time step, " << (++count) << "-iteration, " << endl; cout << (i+1) << "-time step, " << (count) << "-iteration, " << endl; #ifdef __GOLDEN_SECTION_LINE_SEARCH Matrix_Representation::Assembly_Switch = Matrix_Representation::ALL; mr.assembly(); // form tangent stiffness matrix(lhs) and residual vector(rhs) new_time_flag = FALSE; // flag to turn off Step 3 of the finite deformation procedure, after first iteration time1 = time(0); ofs << " after aseembly elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; C0 p = ((C0)(mr.rhs())) / ((C0)(mr.lhs())); // Newton's formula for search direction energy = fabs( (double)(p * ((C0)(mr.rhs())) ) ); //=begin update 1====================================================== // update fixed dof if(count == 1) { // first iteration only for(int j = 86; j <= 92; j++) { gh[j][1] += d_gh[j][1]; delta_gh[j][1] = d_gh[j][1]; d_gh[j][1] = 0.0; } uh = gh; delta_uh = delta_gh; d_uh = d_gh; } //=end update 1======================================================== if(count != 1 && energy > energy_old) { // energy norm diverges; abandant the current solution and activate line search cout << " energy: " << energy << endl; U_h uh_temp(ndf, oh); U_h delta_uh_temp(ndf, oh); // save current total displacement for(int j = 0; j < uh_temp.total_node_no(); j++) { // and total increment of displacement uh_temp[j] = uh[j]; delta_uh_temp[j] = delta_uh[j]; } double alpha = 0.0, // line search parameter left = 0.0, right = 1.0, length = right-left; // initial bracket int line_search_counter = 0; do { (C0)(mr.rhs()) = 0.0; double golden_alpha; alpha = golden_alpha = (left + 0.618 * length); // temporary update free dof // define size of delta_u; set to zero for a new time step for(int j = 0; j < uh_temp.total_node_no(); j++) { // restore total displacement and total increment of displacement uh[j] = uh_temp[j]; delta_uh[j] = delta_uh_temp[j]; } delta_uh += alpha*p; // total increment of displacement uh += alpha*p; // total displacement Matrix_Representation::Assembly_Switch = Matrix_Representation::RHS; mr.assembly(); // form residual vector (rhs) double golden_residual = norm( (C0)(mr.rhs()) ); C0 mr_rhs_temp = (C0)(mr.rhs()); (C0)(mr.rhs()) = 0.0; double left_alpha; alpha = left_alpha = (left + 0.382 * length); // temporary update free dof // define size of delta_u; set to zero for a new time step for(j = 0; j < uh_temp.total_node_no(); j++) { // restore total displacement and total increment of displacement uh[j] = uh_temp[j]; delta_uh[j] = delta_uh_temp[j]; } delta_uh += alpha*p; // total increment of displacement uh += alpha*p; // total displacement Matrix_Representation::Assembly_Switch = Matrix_Representation::RHS; mr.assembly(); // form residual vector (rhs) double left_residual = norm( (C0)(mr.rhs()) ); if(golden_residual < left_residual) { left = left + 0.382 * length; alpha = golden_alpha; (C0)(mr.rhs()) = mr_rhs_temp; } else { right = left + 0.618 * length; alpha = left_alpha; } length = right-left; //cout << left_residual << ", " << golden_residual << endl; //ofs << line_search_counter << " alpha: " << alpha << endl; cout << line_search_counter << " alpha: " << alpha << " current energy: " << fabs( (double)((alpha*p)*((C0)(mr.rhs()))) ) << endl; } while (line_search_counter++ < 10); du = alpha * p; for(j = 0; j < uh_temp.total_node_no(); j++) { // restore total displacement and total increment of displacement uh[j] = uh_temp[j]; delta_uh[j] = delta_uh_temp[j]; } } else du = p; residual_norm = norm((C0)(mr.rhs())); //energy = fabs( (double)(du * ((C0)(mr.rhs())) ) ); energy_old = energy; if(count == 1) energy_0 = energy; ofs << " residual norm: " << residual_norm << " energy: " << fabs( (double)(du*((C0)(mr.rhs()))) ) << endl; cout << " energy: " << fabs( (double)(du*((C0)(mr.rhs()))) ) << endl; (C0)(mr.lhs()) = 0.0; (C0)(mr.rhs()) = 0.0; #else mr.assembly(); // form tangent stiffness matrix(rhs) and residual vector(lhs) new_time_flag = FALSE; // flag to turn off Step 3 of the finite deformation procedure, after first iteration time1 = time(0); ofs << " after aseembly elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; du = ((C0)(mr.rhs())) / ((C0)(mr.lhs())); // Newton's formula residual_norm = norm((C0)(mr.rhs())); energy = fabs( (double)(du * ((C0)(mr.rhs())) ) ); if(count == 1) energy_0 = energy; ofs << " residual norm: " << residual_norm << " energy: " << energy << endl; cout << " energy: " << energy << endl; //ofs << uh << endl; (C0)(mr.lhs()) = 0.0; (C0)(mr.rhs()) = 0.0; //=begin update ====================================================== // update fixed dof if(count == 1) { // first iteration only for(int j = 86; j <= 92; j++) { gh[j][1] += d_gh[j][1]; delta_gh[j][1] = d_gh[j][1]; d_gh[j][1] = 0.0; } uh = gh; delta_uh = delta_gh; d_uh = d_gh; } #endif // update free dof // define size of delta_u; set to zero for a new time step if(i == 0 && count == 1) delta_u &= C0(du.length(), (double*)0); delta_uh = delta_u += du; // total increment of displacement uh += du; // total displacement //=end=update========================================================== } while(energy > (EPSILON*energy_0) && // convergence criterion (count < MAX_ITERATION_NO) ); if(count == MAX_ITERATION_NO) { ofs << " Warning: global Newton iteration failed!" << endl << " current solution: " << endl << uh << endl; return 0; } #if defined(__REACTION) // compute reactions Matrix_Representation::Assembly_Switch = Matrix_Representation::REACTION; mr.assembly(FALSE); // without loading (flag == FALSE) ofs << "reactions:" << endl << (mr.global_nodal_value()) << endl; Matrix_Representation::Assembly_Switch = Matrix_Representation::ALL; // reset #endif //======================================================================== //if(i == 59) { // compute stress projection on nodes Matrix_Representation::Assembly_Switch = Matrix_Representation::NODAL_STRESS; mr.assembly(FALSE); C0 nodal_stress = mr.global_nodal_value(); // compute strain projection on nodes Matrix_Representation::Assembly_Switch = Matrix_Representation::NODAL_STRAIN; mr.assembly(FALSE); C0 nodal_strain = mr.global_nodal_value(); // compute yield ratio projection on nodes Matrix_Representation::Assembly_Switch = Matrix_Representation::NODAL_SCALAR; mr.assembly(FALSE); C0 nodal_yield_ratio = mr.global_nodal_value(); femap_output_data(i, gd, ofs1, nodal_stress, nodal_strain, nodal_yield_ratio); // ofstream ofs1("temp59.cpp", ios::out | ios::trunc); //output(gd, ofs1, nodal_stress, nodal_strain, nodal_yield_ratio); Matrix_Representation::Assembly_Switch = Matrix_Representation::ALL; //} for(j = 0; j < oh.total_element_no(); j++) { // update history terms Omega_eh *elem = &(oh(j)); ((Non_Linear_Omega_eh*) elem)->update_history_data(); } // output solution ofs << endl << (i+1) << "-time step solution: " << endl; ofs << uh << endl; } ofs1.close(); //============================================================================== // post-processing to compute reactions, stresses and strains //============================================================================== #if defined(__GAUSS_STRAIN) // compute stresses at gauss integration point Matrix_Representation::Assembly_Switch = Matrix_Representation::STRAIN; mr.assembly(FALSE); time1 = time(0); ofs << "elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; #endif #if defined(__NODAL_STRAIN) // compute stress projection on nodes Matrix_Representation::Assembly_Switch = Matrix_Representation::NODAL_STRAIN; mr.assembly(FALSE); ofs << "nodal strains:" << endl << (mr.global_nodal_value()) << endl; time1 = time(0); ofs << "elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; #endif #if defined(__GAUSS_STRESS) // compute stresses at gauss integration point Matrix_Representation::Assembly_Switch = Matrix_Representation::STRESS; mr.assembly(FALSE); time1 = time(0); ofs << "elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; #endif #if defined(__NODAL_STRESS) // compute stress projection on nodes Matrix_Representation::Assembly_Switch = Matrix_Representation::NODAL_STRESS; mr.assembly(FALSE); ofs << "nodal stresses:" << endl << (mr.global_nodal_value()) << endl; time1 = time(0); ofs << "elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; #endif //============================================================================== time1 = time(0); ofs << "total elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; ofs.close(); return 0; } #include "output3.cpp" /*void output(Global_Discretization& gd, ofstream& ofs1, C0& sig, C0& eps, C0& yield) { // input for AutoCAD ADS program // include header Omega_h &oh = gd.omega_h(); U_h &uh = gd.u_h(); gh_on_Gamma_h &gh = gd.gh_on_gamma_h(); ofs1 << "#include \"c:\\fe\\acad\\fe2d\\test\\grafe.h\"" << endl << endl; ofs1 << "extern Omega_h domain;" << endl; // void input(center_x, center_y) ofs1 << "void input(center_x, center_y)" << endl << " ads_real *center_x, *center_y;" << endl << "{" << endl << " domain.index = 0;" << endl << " domain.total_node_no = " << oh.total_node_no() << ";" << endl << " define_node(" << oh.total_node_no() << ");" << endl << " window_range(center_x, center_y);" << endl << " domain.total_element_no = " << oh.total_element_no() << ";" << endl << " define_element(" << oh.total_element_no() << ");" << endl << "}" << endl << endl; // void define_node(total_node_no) ofs1 << "void define_node(total_node_no)" << endl << " int total_node_no;" << endl << "{" << endl << " Node *na;" << endl << " domain.node_array = (Node*) calloc(total_node_no, sizeof(Node));" << endl; for(int i = 0; i < oh.total_node_no(); i++) { ofs1 << " na = &(domain.node_array[" << i << "]);" << endl << " na->node_no = " << oh.node_array()[i].node_no() << "; na->ndf = " << uh.u_h_array()[i].no_of_dim() << ";" << "; na->nsd = " << oh.node_array()[i].no_of_dim() << ";" << endl << " Spoint(na->coord, " << oh.node_array()[i][0] << ", " << oh.node_array()[i][1] << ", 0.0);" << endl << " na->u = (ads_real*)calloc(" << uh.u_h_array()[i].no_of_dim() << ", sizeof(ads_real));" << endl; for(int j = 0; j < uh.u_h_array()[i].no_of_dim(); j++) ofs1 << " na->u[" << j << "] = " << uh.u_h_array()[i][j] << "; "; ofs1 << endl; ofs1 << " na->constraint_flag = (int*)calloc(" << gh.gh_array()[i].no_of_dim() << ", sizeof(int));" << endl; for(j = 0; j < gh.gh_array()[i].no_of_dim(); j++) ofs1 << " na->constraint_flag[" << j << "] = " << gh.gh_array()[i](j) << "; "; ofs1 << endl; ofs1 << " na->f = (ads_real*)calloc(" << gh.gh_array()[i].no_of_dim() << ", sizeof(ads_real));" << endl; for(j = 0; j < gh.gh_array()[i].no_of_dim(); j++) ofs1 << " na->f[" << j << "] = " << gh.gh_array()[i][j] << "; "; ofs1 << endl; ofs1 << " na->sig = (ads_real*)calloc(3, sizeof(ads_real));" << endl; for(j = 0; j < 3; j++) { // skip sigma_r int index = (j == 2) ? 3 : j; ofs1 << " na->sig[" << j << "] = " << sig[i][index] << "; "; } ofs1 << endl; ofs1 << " na->eps = (ads_real*)calloc(3, sizeof(ads_real));" << endl; for(j = 0; j < 3; j++) { // skip epsilon_r int index = (j == 2) ? 3 : j; ofs1 << " na->eps[" << j << "] = " << eps[i][index] << "; "; } ofs1 << endl; ofs1 << " na->yield = " << yield[i] << "; "; ofs1 << endl; } ofs1 << "}" << endl << endl; // void define_element(total_element_no) ofs1 << "void define_element(total_element_no)" << endl << " int total_element_no;" << endl << "{" << endl << " Omega_eh *ea;" << endl << " domain.omega_eh_array = (Omega_eh*) calloc(" << oh.total_element_no() << ", sizeof(Omega_eh));" << endl; for(i = 0; i < oh.total_element_no(); i++) { ofs1 << " ea = &(domain.omega_eh_array[" << i << "]);" << endl << " ea->element_no = " << oh.omega_eh_array()[i].element_no() << "; ea->element_node_no = " << oh.omega_eh_array()[i].element_node_no() << "; " << endl << " ea->element_type_no = " << oh.omega_eh_array()[i].element_type_no() << "; ea->material_type_no = " << oh.omega_eh_array()[i].material_type_no() << ";" << endl << " ea->node_no_array = (int*) calloc(ea->element_node_no, sizeof(int));" << endl << " "; for(int j = 0; j < oh.omega_eh_array()[i].element_node_no(); j++) ofs1 << "ea->node_no_array[" << j << "] = " << oh.omega_eh_array()[i][j] << "; "; ofs1 << endl; } ofs1 << "}" << endl << endl; //void clean_exit(void) ofs1 << "void clean_exit()" << endl << "{" << endl << " int i;" << endl << " for(i = 0; i < domain.total_node_no; i++) {" << endl << " free(domain.node_array[i].u);" << endl << " free(domain.node_array[i].constraint_flag);" << endl << " free(domain.node_array[i].f);" << endl << " free(domain.node_array[i].sig);" << endl << " free(domain.node_array[i].eps);" << endl << " }" << endl << " free(domain.node_array);" << endl << " for(i = 0; i < domain.total_element_no; i++) {" << endl << " free(domain.omega_eh_array[i].node_no_array);" << endl << " }" << endl << " free(domain.omega_eh_array);" << endl << "}" << endl; ofs1.close(); }*/