#include #include #include #include #include "vs.h" #include "fe.h" #include //============================================================================== // Step 1: Global_Discretizations //============================================================================== // From numerical example in Simo and Taylor (1985, Comput. Meths. in Appl. Mech. Engrg., 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 rheology history terms H0 // (.)_old is at t_n; (.) is at t_n+1 the_s, // deviatoric stress, s_xx, s_yy, s_xy the_s_old, the_ep, // total plastic strain the_ep_old, the_alpha, // "back stress" (for kinematic hardening law) the_alpha_old; void __initialization(); public: Non_Linear_Omega_eh(int en = 0, int etn = 0, int mtn = 0, int nen = 0, int* ena = NULL) : Omega_eh(en, etn, mtn, nen, ena) { __initialization(); } H0& s_old() { return the_s_old; } H0& ep_old() { return the_ep_old; } H0& alpha_old() { return the_alpha_old; } H0& s() { return the_s; } H0& ep() { return the_ep; } H0& alpha() { return the_alpha; } void update_history_data(); }; void Non_Linear_Omega_eh::__initialization() { the_s &= H0(3, (double*)0, qp); the_ep &= H0((double*)0, qp); the_alpha &= H0(3, (double*)0, qp); the_s_old &= H0(3, (double*)0, qp); the_ep_old &= H0((double*)0, qp); the_alpha_old &= H0(3, (double*)0, qp); } void Non_Linear_Omega_eh::update_history_data() { the_s_old = the_s; the_ep_old = the_ep; the_alpha_old = the_alpha; } Omega_h::Omega_h() { // discritized global domain // define nodes double v[2]; Node *node; v[0] = 5.0; v[1] = 0.0; node = new Node(0, 2, v); node_array().add(node); v[0] = 6.25; node = new Node(1, 2, v); node_array().add(node); v[0] = 7.5; node = new Node(2, 2, v); node_array().add(node); v[0] = 8.625; node = new Node(3, 2, v); node_array().add(node); v[0] = 10.0; node = new Node(4, 2, v); node_array().add(node); v[0] = 4.9039; v[1] = 0.9755; node = new Node(5, 2, v); node_array().add(node); v[0] = 6.1299; v[1] = 1.2189; node = new Node(6, 2, v); node_array().add(node); v[0] = 7.3559; v[1] = 1.4632; node = new Node(7, 2, v); node_array().add(node); v[0] = 8.678; v[1] = 1.4632; node = new Node(8, 2, v); node_array().add(node); v[0] = 10.0; v[1] = 1.4632; node = new Node(9, 2, v); node_array().add(node); v[0] = 4.6194; v[1] = 1.9134; node = new Node(10, 2, v); node_array().add(node); v[0] = 5.7742; v[1] = 2.3918; node = new Node(11, 2, v); node_array().add(node); v[0] = 6.9291; v[1] = 2.8701; node = new Node(12, 2, v); node_array().add(node); v[0] = 8.625; node = new Node(13, 2, v); node_array().add(node); v[0] = 10.0; node = new Node(14, 2, v); node_array().add(node); v[0] = 4.1573; v[1] = 2.7779; node = new Node(15, 2, v); node_array().add(node); v[0] = 5.1967; v[1] = 3.4724; node = new Node(16, 2, v); node_array().add(node); v[0] = 6.2360; v[1] = 4.1668; node = new Node(17, 2, v); node_array().add(node); v[0] = 8.118; v[1] = 4.6760; node = new Node(18, 2, v); node_array().add(node); v[0] = 10.0; v[1] = 5.1851; node = new Node(19, 2, v); node_array().add(node); v[0] = 3.5355; v[1] = 3.5355; node = new Node(20, 2, v); node_array().add(node); v[0] = 4.4194; v[1] = 4.4194; node = new Node(21, 2, v); node_array().add(node); v[0] = 5.3033; v[1] = 5.3033; node = new Node(22, 2, v); node_array().add(node); v[0] = 7.6517; v[1] = 6.4017; node = new Node(23, 2, v); node_array().add(node); v[0] = 10.0; v[1] = 7.5; node = new Node(24, 2, v); node_array().add(node); v[0] = 0.0; v[1] = 5.0; node = new Node(25, 2, v); node_array().add(node); v[0] = 0.9755; v[1] = 4.9039; node = new Node(26, 2, v); node_array().add(node); v[0] = 1.9134; v[1] = 4.6194; node = new Node(27, 2, v); node_array().add(node); v[0] = 2.7779; v[1] = 4.1573; node = new Node(28, 2, v); node_array().add(node); v[0] = 0.0; v[1] = 6.25; node = new Node(29, 2, v); node_array().add(node); v[0] = 1.1959; v[1] = 6.0121; node = new Node(30, 2, v); node_array().add(node); v[0] = 2.3918; v[1] = 5.7742; node = new Node(31, 2, v); node_array().add(node); v[0] = 3.4056; v[1] = 5.0968; node = new Node(32, 2, v); node_array().add(node); v[0] = 0.0; v[1] = 7.5; node = new Node(33, 2, v); node_array().add(node); v[0] = 1.4632; v[1] = 7.3559; node = new Node(34, 2, v); node_array().add(node); v[0] = 2.8701; v[1] = 6.9291; node = new Node(35, 2, v); node_array().add(node); v[0] = 4.1668; v[1] = 6.2360; node = new Node(36, 2, v); node_array().add(node); v[0] = 0.0; v[1] = 8.75; node = new Node(37, 2, v); node_array().add(node); v[0] = 1,4351; node = new Node(38, 2, v); node_array().add(node); v[0] = 2.8701; node = new Node(39, 2, v); node_array().add(node); v[0] = 4.0867; node = new Node(40, 2, v); node_array().add(node); v[0] = 5.3033; node = new Node(41, 2, v); node_array().add(node); v[0] = 7.6517; node = new Node(42, 2, v); node_array().add(node); v[0] = 10.0; node = new Node(43, 2, v); node_array().add(node); v[0] = 0.0; v[1] = 10.0; node = new Node(44, 2, v); node_array().add(node); v[0] = 1,4351; node = new Node(45, 2, v); node_array().add(node); v[0] = 2.8701; node = new Node(46, 2, v); node_array().add(node); v[0] = 4.0867; node = new Node(47, 2, v); node_array().add(node); v[0] = 5.3033; node = new Node(48, 2, v); node_array().add(node); v[0] = 7.6517; node = new Node(49, 2, v); node_array().add(node); v[0] = 10.0; node = new Node(50, 2, v); node_array().add(node); v[0] = 0.0; v[1] = 11.25; node = new Node(51, 2, v); node_array().add(node); v[0] = 1,4351; node = new Node(52, 2, v); node_array().add(node); v[0] = 2.8701; node = new Node(53, 2, v); node_array().add(node); v[0] = 4.0867; node = new Node(54, 2, v); node_array().add(node); v[0] = 5.3033; node = new Node(55, 2, v); node_array().add(node); v[0] = 7.6517; node = new Node(56, 2, v); node_array().add(node); v[0] = 10.0; node = new Node(57, 2, v); node_array().add(node); v[0] = 0.0; v[1] = 12.6; node = new Node(58, 2, v); node_array().add(node); v[0] = 1,4351; node = new Node(59, 2, v); node_array().add(node); v[0] = 2.8701; node = new Node(60, 2, v); node_array().add(node); v[0] = 4.0867; node = new Node(61, 2, v); node_array().add(node); v[0] = 5.3033; node = new Node(62, 2, v); node_array().add(node); v[0] = 7.6517; node = new Node(63, 2, v); node_array().add(node); v[0] = 10.0; node = new Node(64, 2, v); node_array().add(node); v[0] = 0.0; v[1] = 13.95; node = new Node(65, 2, v); node_array().add(node); v[0] = 1,4351; node = new Node(66, 2, v); node_array().add(node); v[0] = 2.8701; node = new Node(67, 2, v); node_array().add(node); v[0] = 4.0867; node = new Node(68, 2, v); node_array().add(node); v[0] = 5.3033; node = new Node(69, 2, v); node_array().add(node); v[0] = 7.6517; node = new Node(70, 2, v); node_array().add(node); v[0] = 10.0; node = new Node(71, 2, v); node_array().add(node); v[0] = 0.0; v[1] = 15.3; node = new Node(72, 2, v); node_array().add(node); v[0] = 1,4351; node = new Node(73, 2, v); node_array().add(node); v[0] = 2.8701; node = new Node(74, 2, v); node_array().add(node); v[0] = 4.0867; node = new Node(75, 2, v); node_array().add(node); v[0] = 5.3033; node = new Node(76, 2, v); node_array().add(node); v[0] = 7.6517; node = new Node(77, 2, v); node_array().add(node); v[0] = 10.0; node = new Node(78, 2, v); node_array().add(node); v[0] = 0.0; v[1] = 16.65; node = new Node(79, 2, v); node_array().add(node); v[0] = 1,4351; node = new Node(80, 2, v); node_array().add(node); v[0] = 2.8701; node = new Node(81, 2, v); node_array().add(node); v[0] = 4.0867; node = new Node(82, 2, v); node_array().add(node); v[0] = 5.3033; node = new Node(83, 2, v); node_array().add(node); v[0] = 7.6517; node = new Node(84, 2, v); node_array().add(node); v[0] = 10.0; node = new Node(85, 2, v); node_array().add(node); v[0] = 0.0; v[1] = 18.0; node = new Node(86, 2, v); node_array().add(node); v[0] = 1,4351; node = new Node(87, 2, v); node_array().add(node); v[0] = 2.8701; node = new Node(88, 2, v); node_array().add(node); v[0] = 4.0867; node = new Node(89, 2, v); node_array().add(node); v[0] = 5.3033; node = new Node(90, 2, v); node_array().add(node); v[0] = 7.6517; node = new Node(91, 2, v); node_array().add(node); v[0] = 10.0; node = new Node(92, 2, v); node_array().add(node); // define elements; eight-node serendipity element 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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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); omega_eh_array().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 int ndf = 2; Omega_h oh; gh_on_Gamma_h gh(ndf, oh); U_h uh(ndf, oh); Global_Discretization gd(oh, gh, uh); U_h delta_uh(ndf, oh); // total increment of displacement of one incremental loading step gh_on_Gamma_h delta_gh(ndf, oh); Global_Discretization gd_delta_uh(oh, delta_gh, delta_uh); U_h d_uh(ndf, oh); // increment of increment of displacement within one incremental loading step gh_on_Gamma_h d_gh(ndf, oh); Global_Discretization gd_d_uh(oh, d_gh, d_uh); //============================================================================== // Step 2. Element_Formulation-- // quadri-lateral plane strain elastoplastic element //============================================================================== // DEFINITION of Element_Formulation class Elastoplastic_B_bar_Q41 : public Element_Formulation { C0 delta_ul; // total increment of displacement of one incremental loading step (fixed and free) H0 s, s_old, ep, ep_old, alpha, alpha_old; // history terms void __initialization(int); public: Elastoplastic_B_bar_Q41(Element_Type_Register a) : Element_Formulation(a) {} Element_Formulation *make(int, Global_Discretization&); Elastoplastic_B_bar_Q41(int, Global_Discretization&); }; void Elastoplastic_B_bar_Q41::__initialization(int en) { // delta_ul: free and fixed increment of displacement at t_n+1; to compute current incremental strain rate // delta_gh = 0.0 at first iteration, delta_gh(after first iteration) = d_gh(at first iteration; d_gh set to zero after first iteration) delta_ul &= gd_delta_uh.element_free_variable(en) + gd_delta_uh.element_fixed_variable(en); // ul at t_n+1 to compute internal stress divergence term 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 set to zero after first iteration gl = gd_d_uh.element_fixed_variable(en); // history terms Omega_eh *elem = &(oh(en)); s &= ((Non_Linear_Omega_eh*) elem)->s(); // at t_n+1 ep &= ((Non_Linear_Omega_eh*) elem)->ep(); alpha &= ((Non_Linear_Omega_eh*) elem)->alpha(); s_old &= ((Non_Linear_Omega_eh*) elem)->s_old(); // at t_n ep_old &= ((Non_Linear_Omega_eh*) elem)->ep_old(); alpha_old &= ((Non_Linear_Omega_eh*) elem)->alpha_old(); } Element_Formulation* Elastoplastic_B_bar_Q41::make(int en, Global_Discretization& gd) { return new Elastoplastic_B_bar_Q41(en,gd); } static const double i_dev[3][3] = { {(1.0-1.0/3.0), -1.0/3.0, 0.0}, {-1.0/3.0, (1.0-1.0/3.0),0.0}, {0.0, 0.0, 1.0} }; C0 I_dev = MATRIX("int, int, const double*", 3, 3, i_dev[0]); static const double i_dev2[3][3] = { {(1.0-1.0/3.0), -1.0/3.0, 0.0}, {-1.0/3.0, (1.0-1.0/3.0),0.0}, {0.0, 0.0, 0.5} }; C0 I_dev2 = MATRIX("int, int, const double*", 3, 3, i_dev2[0]); //static double identity[3][3] = { // {1.0, 0.0, 0.0}, // {0.0, 1.0, 0.0}, // {0.0, 0.0, 1.0} // }; static double one[3] = {1.0, 1.0, 0.0}; C0 //I = MATRIX("int, int, const double*", 3, 3, identity[0]), One = VECTOR("int, const double*", 3, one); static double identity_mu[3][3] = { {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 0.5} }; C0 I_mu = MATRIX("int, int, const double*", 3, 3, identity_mu[0]); Elastoplastic_B_bar_Q41::Elastoplastic_B_bar_Q41(int en, Global_Discretization& gd) : Element_Formulation(en, gd) { Quadrature qp(2, 4); // 2-dimension, 2x2 integration points 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; // coordinate transformation H1 X = N*xl; // Physical Coordinates // N.B. d(N)/d(X) = d(N)/d(Z)*(d(X)/d(Z))^-1 => {8x2}*{2x2}={8x2} H0 Nx = d(N) * d(X).inverse(); // Weighting derivative J dv(d(X).det()); // the Jacobian H0 w_x = INTEGRABLE_SUBMATRIX("int, int, H0&", 1, 2/*nsd*/, Nx), // dim ~Nx = {2x4}, dim ~wx = {{2x1}x4} wx, wy, B, B_bar; wx &= w_x[0][0]; wy &= w_x[0][1]; C0 zero(0.0); B &= (~wx || zero) & // dim B = {3x8}, where dim ~wx[][] = {1x4} (zero || ~wy ) & (~wy || ~wx ); // pressure-volumetric strain shape function H0 N_vp(qp); //H0 zai, eta; // zai &= ((H0)Z[0]); eta &= ((H0)Z[1]); N_vp = 1.0; H0 mN_vp = N_vp & N_vp & zero; // 3x1 C0 C = ((~mN_vp)*B) | dv, // 1x3 x 3x2 = 1x2 E = (N_vp*N_vp) | dv, // 1 W = C/E; // 1x2 H0 B_vol = (1.0/3.0) * (mN_vp % W[0]), // 3x1 x 1x2 B_dev = I_dev * B; B_bar &= B_dev + B_vol; //=========================================================================== // Elastoplastic material model; Radial Return Algorithm //=========================================================================== // material constants: Y, d_kapa, d_H_alpha; initial alpha(back-stress) = 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 stress), e^p sub n(plastic strain), alpha_n(back_stress) // compute for t_n+1: yield ratio, sigma(for residual), tangent operator(for stiffness) //=========================================================================== __initialization(en); // retreive incerment of displacement and history terms from last time step //=========================================================================== // trial elastic stress--may shoot out of the yeild surface //=========================================================================== H0 delta_e = INTEGRABLE_VECTOR("int, Quadrature", 3, qp); // deviatoric strain rate delta_e = 0.0; for(int i = 0; i < nen; i++) // delta_ul: element total increment of displacement (fixed and free) delta_e += B_dev(i*ndf)*delta_ul[i*ndf] + B_dev(i*ndf+1)*delta_ul[i*ndf+1]; s = s_old + delta_e*(2.0*mu_*I_mu); // "TRIAL elastic" deviatoric stress; (i-1) H0 ZAI = s - alpha_old, // (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_0_q &= ZAI[0].qv(i); ZAI_1_q &= ZAI[1].qv(i); ZAI_2_q &= ZAI[2].qv(i); ZAI_norm.qv(i) = sqrt((double)(ZAI_0_q.pow(2)+ZAI_1_q.pow(2)+2.0*ZAI_2_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 Dt(3, 3, (double*)0, qp), // tangent moduli // sigma is for computing a residual(i.e., the rhs internal stress divergence term) sigma = INTEGRABLE_VECTOR("int, Quadrature", 3, qp), 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_old_q, alpha_q, alpha_old_q, ZAI_q, ZAI_norm_q, sigma_q, Dt_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_old_q = ep_old.qv(i); // history terms alpha_q &= alpha.qv(i); alpha_old_q = alpha_old.qv(i); ZAI_q &= ZAI.qv(i); ZAI_norm_q &= ZAI_norm.qv(i); // stress deviator sigma_q &= sigma.qv(i); Dt_q &= Dt.qv(i); // internal stress and tangent moduli 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 ofs << "***Plastic state: element # " << en << ", quadrature point # " << i << ", 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_old_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_old_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_old_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), //===================================================================== double gamma = 1.0/(1.0+((((double)d(KAPA))+d_H_alpha)/3.0/mu_)); #if defined(__TEST_CONSISTENT_TANGENT_OPERATOR) // consistent tangent elastoplastic moduli double beta = sqrt(2.0/3.0)*((double)((C0)KAPA) + delta_H_alpha)/(double)ZAI_norm_q, gamma_bar = gamma - (1.0-beta); Dt_q = lambda_*(One%One) + // lambda_ = K_-2/3 mu_ 2.0 * mu_ * (beta * I_mu - gamma_bar * (n%n)); #else // continuum tangent elastoplastic moduli Dt_q = lambda_*(One%One) + 2.0 * mu_ * (I_mu - gamma * (n%n)); #endif } else { // elastic state Dt_q = lambda_*(One%One) + 2.0*mu_*I_mu; ofs << " Elastic state: element # " << en << ", quadrature point # " << i << ", yield ratio: " << yield_ratio_q << endl; } // end if plastic or elastic state // compute stress at t_n+1 (v); deviatoric stress plus elastic volumetric stress C0 eps_v(3, (double*)0); for(int i = 0; i < nen; i++) eps_v += B_vol_q(i*ndf)*ul[i*ndf] + B_vol_q(i*ndf+1)*ul[i*ndf+1]; eps_v *= 3.0; sigma_q = s_q + K_*eps_v; // (v); deviatoric stress plus elastic volumetric pressure } // end for each quadrature point // post-processing code for reaction, stress, and strain if(Matrix_Representation::Assembly_Switch == Matrix_Representation::REACTION) { stiff &= ((~B_bar)*Dt*B_bar)|dv; // 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 = (ndf+1)*ndf/2; 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", 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]; // 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 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-xx," << setw(14) << "sigma-yy," << 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)) // replace sigma_hat with sigma that is already available << 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 = (ndf+1)*ndf/2; 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", 3, qp); sigma_hat = 0.0; 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]; // replace sigma_hat with sigma that is already available } } 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 stiff &= ((~B_bar)*Dt*B_bar)|dv; // 8x8 force &= -((~B_bar)*sigma)|dv; // compute the residual (non-linear problem only); internal stress divergence term } } // REGISTRATION of Element(s) Element_Formulation* Element_Formulation::type_list = 0; Element_Type_Register element_type_register_instance; static Elastoplastic_B_bar_Q41 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.0) << " 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\\test.neu", ios::out | ios::trunc); femap_output_model(gd, ofs1); // incremental loading (pull apart on top line) for(int i = 0; i < 10; i++) { // incremental loading loop //======================================================================== //"global" Newton-Raphson Method for elasto-plastic minimun energy solution //======================================================================== int count = 0; #if defined(__TEST_CONJUGATE_GRADIENT_METHOD) const int MAX_ITERATION_NO = 20; const double EPSILON = 1.e-9; //=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.01; // displacement increment on top delta_gh[j][1] = 0.0; } d_uh = d_gh; delta_uh = delta_gh; uh = gh; //===================================================================== double norm_p_0; C0 p; do { ofs << (i+1) << "-time step, " << (++count) << "-iteration, " << endl; cout << (i+1) << "-time step, " << (count) << "-iteration, " << endl; //=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======================================================== C0 g, Q; for(int k = 0; k < MAX_ITERATION_NO; k++) { if(k == 0) { Matrix_Representation::Assembly_Switch = Matrix_Representation::ALL; mr.assembly(); // form tangent stiffness matrix(lhs) and residual vector(rhs) time1 = time(0); ofs << " after aseembly elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; Q = ((C0)(mr.lhs())); p = ((C0)(mr.rhs())); g = -p; // negative gradient as search direction; p = -g (C0)(mr.rhs()) = 0.0; (C0)(mr.lhs()) = 0.0; } double alpha = (-g*p)/(p*(Q*p)); du = alpha*p; // x_k+1 = x_k + alpha p if(i == 0 && count == 1 && k == 0) delta_u &= C0(p.length(), (double*)0); delta_uh = delta_u += du; // total increment of displacement uh += du; // total displacement // compute new search direction; p_i+1 if(k != 29) { Matrix_Representation::Assembly_Switch = Matrix_Representation::ALL; mr.assembly(); // form residual vector (rhs) C0 g1 = -(C0)(mr.rhs()); // next gradient Q = ((C0)(mr.lhs())); (C0)(mr.rhs()) = 0.0; (C0)(mr.lhs()) = 0.0; double beta = g1*(Q*p)/(p*(Q*p)); // Luenberger //double beta = (g1*g1)/g.pow(2);//(g1*p)/p.pow(2); // Fletcher-Reeves //double beta = ((g1-g)*g1)/g.pow(2); // Polak-Ribiere p = -g1 + beta * p; // new search direction g = g1; // update gradient } if(count == 1 && k == 0) norm_p_0 = (double)norm(p); if(((double)norm(p)) < (norm_p_0*EPSILON)) break; cout << k << " p norm: " << ((double)norm(p)) << endl; } // end of partial conjugate gradient for loop } while(((double)norm(p)) > (EPSILON*norm_p_0) && // convergence criterion (count < MAX_ITERATION_NO) ); #elif defined(__QUASI_NEWTON_BFGS_METHOD) const int MAX_ITERATION_NO = 10; const double EPSILON = 1.e-9; //=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.01; // displacement increment on top delta_gh[j][1] = 0.0; } d_uh = d_gh; delta_uh = delta_gh; uh = gh; //===================================================================== double energy, energy_0; do { ofs << (i+1) << "-time step, " << (++count) << "-partial quasi-newton loop, " << endl; cout << (i+1) << "-time step, " << (count) << "-partial quasi-newton loop, " << endl; //=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======================================================== C0 g, Q, H, d_; for(int k = 0; k < MAX_ITERATION_NO; k++) { // partial quasi-newton for loop if(k == 0) { Matrix_Representation::Assembly_Switch = Matrix_Representation::ALL; mr.assembly(); // form tangent stiffness matrix(lhs) and residual vector(rhs) Q = ((C0)(mr.lhs())); // Hessian g = -((C0)(mr.rhs())); // gradient Cholesky dQ(Q, 1.e-12); // modified cholesky method for semi-positive definite symmetric matrix H = dQ.inverse(); // inverse of Hessian d_ = -H*g; // search direction (C0)(mr.rhs()) = 0.0; (C0)(mr.lhs()) = 0.0; } //if(count == 1 && k == 0) energy_0 = norm(d_*g); //ofs << k << " bfgs iteration energy: " << norm(d_*g) << " residual: " << norm(g) << endl; //cout << k << " bfgs iteration energy: " << norm(d_*g) << " residual: " << norm(g) << endl; #ifdef __GOLDEN_SECTION_LINE_SEARCH // gloden section line search for alpha begins====================== 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 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*d_; // total increment of displacement uh += alpha*d_; // 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(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*d_; // total increment of displacement uh += alpha*d_; // 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 << " (" << line_search_counter << ", " << alpha << ", " << norm( (C0)(mr.rhs()) ) << ") "; } while (line_search_counter++ < 10); //cout << endl; du = alpha * d_; // gloden section line search for alpha ends======================== double residual_norm = norm((C0)(mr.rhs())); energy = fabs( (double)(du * ((C0)(mr.rhs())) ) ); if(count == 1 && k == 0) energy_0 = energy; ofs << " residual norm: " << residual_norm << " energy: " << energy << endl; ofs << k << " bfgs iteration energy: " << energy << " residual: " << residual_norm << endl; cout << k << " bfgs iteration energy: " << energy << " residual: " << residual_norm << endl; (C0)(mr.rhs()) = 0.0; 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]; } #else du = d_; double residual_norm = norm(g); energy = fabs( (double)(du * g) ); if(count == 1 && k == 0) energy_0 = energy; ofs << " residual norm: " << residual_norm << " energy: " << energy << endl; ofs << k << " bfgs iteration energy: " << energy << " residual: " << residual_norm << endl; cout << k << " bfgs iteration energy: " << energy << " residual: " << residual_norm << endl; #endif if(i == 0 && count == 1 && k == 0) delta_u &= C0(d_.length(), (double*)0); delta_uh = delta_u += du; // total increment of displacement uh += du; // total displacement if(k != (MAX_ITERATION_NO-1)) { Matrix_Representation::Assembly_Switch = Matrix_Representation::RHS; mr.assembly(); // form residual vector (rhs) C0 g1 = -(C0)(mr.rhs()); // next gradient (C0)(mr.rhs()) = 0.0; C0 p, q, Hq; p &= du; q = g1 -g; Hq = H*q; // BFGS update H += (1.0+(q*Hq)/(q*p))*((p%p)/(p*q)) - (p%Hq+Hq%p)/(q*p); // update inverse of Hessian g = g1; // update gradient d_ = -H*g; // update search direction } if(energy < (energy_0*EPSILON)) break; } // end partial quasi-newton for loop } while(energy > (energy_0*EPSILON) && count < MAX_ITERATION_NO); #else double residual_norm, energy, energy_0; const int MAX_ITERATION_NO = 50; const double EPSILON = 1.e-9; //=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.01; // 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; mr.assembly(); // form tangent stiffness matrix(rhs) and residual vector(lhs) time1 = time(0); ofs << " after aseembly elipsed time: " << (difftime(time1, time0)/60.0) << " 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; (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; } // 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) ); #endif if(count == MAX_ITERATION_NO) { ofs << " Warning: global Newton iteration failed!" << endl << " current solution: " << endl << uh; //return 0; } //======================================================================== //if(i == 9) { // 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); Matrix_Representation::Assembly_Switch = Matrix_Representation::ALL; // set back to default //} 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 reaction, stresses and strains //============================================================================== #if defined(__REACTION) // compute reaction Matrix_Representation::Assembly_Switch = Matrix_Representation::REACTION; mr.assembly(FALSE); // without loading (flag = FALSE) ofs << "reaction:" << endl << (mr.global_nodal_value()) << endl; time1 = time(0); ofs << "elipsed time: " << (difftime(time1, time0)/60.0) << " minutes." << endl; #endif #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.0) << " 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.0) << " 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.0) << " 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.0) << " minutes." << endl; #endif //============================================================================== time1 = time(0); ofs << "total elipsed time: " << (difftime(time1, time0)/60.0) << " minutes." << endl; ofs.close(); return 0; } #include "output2.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++) ofs1 << " na->sig[" << j << "] = " << sig[i][j] << "; "; ofs1 << endl; ofs1 << " na->eps = (ads_real*)calloc(3, sizeof(ads_real));" << endl; for(j = 0; j < 3; j++) ofs1 << " na->eps[" << j << "] = " << eps[i][j] << "; "; 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(); }*/