#include #include #include #include #include "vs.h" #include "fe.h" #include "omega_h_n.h" Matrix_Representation_Couple::assembly_switch Matrix_Representation_Couple::Assembly_Switch = Matrix_Representation_Couple::ALL; #include "global_discretization_gamma_h_n.h" #include "contact_searching.h" //============================================================================== // Step 1: Global_Discretizations //============================================================================== static const double E_[2] = {1.e5, 1.e10}; static const double v_ = (0.5-1.e-12); // nearly imcompressibility static const double K_[2] = {E_[0]/3.0/(1-2.0*v_), E_[1]/3.0/(1-2.0*v_)}; static const double lambda_[2] = {v_*E_[0]/((1+v_)*(1-2*v_)), v_*E_[1]/((1+v_)*(1-2*v_))}; static const double mu_[2] = {E_[0]/(2*(1+v_)), E_[1]/(2*(1+v_))}; static const double lambda_bar[2] = // plane stress modification {2*lambda_[0]*mu_[0]/(lambda_[0]+2*mu_[0]), 2*lambda_[1]*mu_[1]/(lambda_[1]+2*mu_[1])}; ofstream ofs("temp.cpp", ios::out | ios::trunc); // DEFINITION of Global_Discretization Omega_h_i::Omega_h_i(int i) : Omega_h(0){ // discritized global domain the_index = i; if(i == 0) { // elastic foundation // define nodes double v[2]; Node *node; // first row v[0] = -30.0; v[1] = -20.0; node = new Node(0, 2, v); node_array().add(node); v[0] = -16.0; node = new Node(1, 2, v); node_array().add(node); v[0] = -8.0; node = new Node(2, 2, v); node_array().add(node); v[0] = 0.0; node = new Node(3, 2, v); node_array().add(node); v[0] = 8.0; node = new Node(4, 2, v); node_array().add(node); v[0] = 16.0; node = new Node(5, 2, v); node_array().add(node); v[0] = 30.0; node = new Node(6, 2, v); node_array().add(node); // second row v[0] = -30.0; v[1] = -8.0; node = new Node(7, 2, v); node_array().add(node); v[0] = -16.0; node = new Node(8, 2, v); node_array().add(node); v[0] = -8.0; node = new Node(9, 2, v); node_array().add(node); v[0] = 0.0; node = new Node(10, 2, v); node_array().add(node); v[0] = 8.0; node = new Node(11, 2, v); node_array().add(node); v[0] = 16.0; node = new Node(12, 2, v); node_array().add(node); v[0] = 30.0; node = new Node(13, 2, v); node_array().add(node); // thrid row v[0] = -30.0; v[1] = 0.0; node = new Node(14, 2, v); node_array().add(node); v[0] = -16.0; node = new Node(15, 2, v); node_array().add(node); v[0] = -8.0; node = new Node(16, 2, v); node_array().add(node); v[0] = 0.0; node = new Node(17, 2, v); node_array().add(node); v[0] = 8.0; node = new Node(18, 2, v); node_array().add(node); v[0] = 16.0; node = new Node(19, 2, v); node_array().add(node); v[0] = 30.0; node = new Node(20, 2, v); node_array().add(node); // define elements; int ena[4]; // element node number array Omega_eh *elem; ena[0] = 0; ena[1] = 1; ena[2] = 8; ena[3] = 7; elem = new Omega_eh(0, 0, 1, 4, ena); omega_eh_array().add(elem); ena[0] = 1; ena[1] = 2; ena[2] = 9; ena[3] = 8; elem = new Omega_eh(1, 0, 1, 4, ena); omega_eh_array().add(elem); ena[0] = 2; ena[1] = 3; ena[2] = 10; ena[3] = 9; elem = new Omega_eh(2, 0, 1, 4, ena); omega_eh_array().add(elem); ena[0] = 3; ena[1] = 4; ena[2] = 11; ena[3] = 10; elem = new Omega_eh(3, 0, 1, 4, ena); omega_eh_array().add(elem); ena[0] = 4; ena[1] = 5; ena[2] = 12; ena[3] = 11; elem = new Omega_eh(4, 0, 1, 4, ena); omega_eh_array().add(elem); ena[0] = 5; ena[1] = 6; ena[2] = 13; ena[3] = 12; elem = new Omega_eh(5, 0, 1, 4, ena); omega_eh_array().add(elem); ena[0] = 7; ena[1] = 8; ena[2] = 15; ena[3] = 14; elem = new Omega_eh(6, 0, 1, 4, ena); omega_eh_array().add(elem); ena[0] = 8; ena[1] = 9; ena[2] = 16; ena[3] = 15; elem = new Omega_eh(7, 0, 1, 4, ena); omega_eh_array().add(elem); ena[0] = 9; ena[1] = 10; ena[2] = 17; ena[3] = 16; elem = new Omega_eh(8, 0, 1, 4, ena); omega_eh_array().add(elem); ena[0] = 10; ena[1] = 11; ena[2] = 18; ena[3] = 17; elem = new Omega_eh(9, 0, 1, 4, ena); omega_eh_array().add(elem); ena[0] = 11; ena[1] = 12; ena[2] = 19; ena[3] = 18; elem = new Omega_eh(10, 0, 1, 4, ena); omega_eh_array().add(elem); ena[0] = 12; ena[1] = 13; ena[2] = 20; ena[3] = 19; elem = new Omega_eh(11, 0, 1, 4, ena); omega_eh_array().add(elem); } else if(i == 1) { // rigid sled // define nodes double v[2]; Node *node; // first row v[0] = -10.0; v[1] = 0.0; node = new Node(0, 2, v); node_array().add(node); v[0] = -4.0; node = new Node(1, 2, v); node_array().add(node); v[0] = 4.0; node = new Node(2, 2, v); node_array().add(node); v[0] = 10.0; node = new Node(3, 2, v); node_array().add(node); // second row v[0] = -10.0; v[1] = 2.0; node = new Node(4, 2, v); node_array().add(node); v[0] = -4.0; node = new Node(5, 2, v); node_array().add(node); v[0] = 4.0; node = new Node(6, 2, v); node_array().add(node); v[0] = 10.0; node = new Node(7, 2, v); node_array().add(node); // thrid row v[0] = -10.0; v[1] = 4.0; node = new Node(8, 2, v); node_array().add(node); v[0] = -4.0; node = new Node(9, 2, v); node_array().add(node); v[0] = 4.0; node = new Node(10, 2, v); node_array().add(node); v[0] = 10.0; node = new Node(11, 2, v); node_array().add(node); // define elements; int ena[4]; // element node number array Omega_eh *elem; ena[0] = 0; ena[1] = 1; ena[2] = 5; ena[3] = 4; elem = new Omega_eh(0, 0, 0, 4, ena); omega_eh_array().add(elem); ena[0] = 1; ena[1] = 2; ena[2] = 6; ena[3] = 5; elem = new Omega_eh(1, 0, 0, 4, ena); omega_eh_array().add(elem); ena[0] = 2; ena[1] = 3; ena[2] = 7; ena[3] = 6; elem = new Omega_eh(2, 0, 0, 4, ena); omega_eh_array().add(elem); ena[0] = 4; ena[1] = 5; ena[2] = 9; ena[3] = 8; elem = new Omega_eh(3, 0, 0, 4, ena); omega_eh_array().add(elem); ena[0] = 5; ena[1] = 6; ena[2] = 10; ena[3] = 9; elem = new Omega_eh(4, 0, 0, 4, ena); omega_eh_array().add(elem); ena[0] = 6; ena[1] = 7; ena[2] = 11; ena[3] = 10; elem = new Omega_eh(5, 0, 0, 4, ena); omega_eh_array().add(elem); } } Contact_Omega_h_i::Contact_Omega_h_i(int i) : Omega_h_i(i) { // do nothing } Gamma_h_i::Gamma_h_i(int i, Omega_h &oh) : Omega_h(oh), the_index(i) { int ena[2]; Omega_eh *elem; if(i == 0) { // elastic foundation upper boundary as potential contact surface ena[0] = 15; ena[1] = 16; elem = new Omega_eh(0, 0, 0, 2, ena); the_omega_eh_array.add(elem); ena[0] = 16; ena[1] = 17; elem = new Omega_eh(1, 0, 0, 2, ena); the_omega_eh_array.add(elem); ena[0] = 17; ena[1] = 18; elem = new Omega_eh(2, 0, 0, 2, ena); the_omega_eh_array.add(elem); ena[0] = 18; ena[1] = 19; elem = new Omega_eh(3, 0, 0, 2, ena); the_omega_eh_array.add(elem); ena[0] = 19; ena[1] = 20; elem = new Omega_eh(4, 0, 0, 2, ena); the_omega_eh_array.add(elem); } else if(i == 1) { // rigid sled lower bounary as potential contact surface ena[0] = 0; ena[1] = 1; elem = new Omega_eh(0, 0, 0, 2, ena); the_omega_eh_array.add(elem); ena[0] = 1; ena[1] = 2; elem = new Omega_eh(1, 0, 0, 2, ena); the_omega_eh_array.add(elem); ena[0] = 2; ena[1] = 3; elem = new Omega_eh(2, 0, 0, 2, ena); the_omega_eh_array.add(elem); } } gh_on_Gamma_h_i::gh_on_Gamma_h_i(int i, int df, Omega_h& omega_h) : gh_on_Gamma_h() { gh_on_Gamma_h::__initialization(df, omega_h); if(i == 0) { // elastic foundation B.C. the_gh_array[node_order(0)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(0)][0] = 0.0; 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(5)](1) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(5)][1] = 0.0; the_gh_array[node_order(6)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(6)][0] = 0.0; the_gh_array[node_order(6)](1) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(6)][1] = 0.0; the_gh_array[node_order(7)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(7)][0] = 0.0; the_gh_array[node_order(13)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(13)][0] = 0.0; the_gh_array[node_order(14)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(14)][0] = 0.0; the_gh_array[node_order(20)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(20)][0] = 0.0; } else if(i == 1) { // rigid sled pushes down kinematically the_gh_array[node_order(8)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(8)][0] = 0.0; the_gh_array[node_order(8)](1) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(8)][1] = -0.2; the_gh_array[node_order(9)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(9)][0] = 0.0; the_gh_array[node_order(9)](1) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(9)][1] = -0.2; the_gh_array[node_order(10)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(10)][0] = 0.0; the_gh_array[node_order(10)](1) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(10)][1] = -0.2; the_gh_array[node_order(11)](0) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(11)][0] = 0.0; the_gh_array[node_order(11)](1) = gh_on_Gamma_h::Dirichlet; the_gh_array[node_order(11)][1] = -0.2; } } // INSTANTIATION of Global_Discretization(s) static int ndf = 2; //-elastic foundation----------------------------------------------------------- static Omega_h_i elastic_foundation_oh(0); static gh_on_Gamma_h_i elastic_foundation_gh(0, ndf, elastic_foundation_oh); static U_h u_0_h(ndf, elastic_foundation_oh); static Global_Discretization *elastic_foundation_type = new Global_Discretization(); static Global_Discretization elastic_foundation_gd(elastic_foundation_oh, elastic_foundation_gh, u_0_h, elastic_foundation_type); // potential elastic foundation contact segment static Gamma_h_i elastic_foundation_gamma_h(0, elastic_foundation_oh); static Global_Discretization_Gamma_h_i elastic_foundation_gamma_h_gd(0, elastic_foundation_gd, elastic_foundation_type, elastic_foundation_gamma_h); //-end of elastic foundation---------------------------------------------------- //-rigid sled------------------------------------------------------------------ static Omega_h_i rigid_sled_oh(1); static gh_on_Gamma_h_i rigid_sled_gh(1, ndf, rigid_sled_oh); static U_h u_1_h(ndf, rigid_sled_oh); static Global_Discretization *rigid_sled_type = new Global_Discretization(); static Global_Discretization rigid_sled_gd(rigid_sled_oh, rigid_sled_gh, u_1_h, rigid_sled_type); // potentail rigid sled contact segment static Gamma_h_i rigid_sled_gamma_h(1, rigid_sled_oh); static Global_Discretization_Gamma_h_i rigid_sled_gamma_h_gd(1, rigid_sled_gd, rigid_sled_type, rigid_sled_gamma_h); //-end of rigid sled----------------------------------------------------------- //-interface between the elastic foundation and the rigid sled----------------- static Contact_Omega_h_i interface_oh(2); static gh_on_Gamma_h_i interface_gh(2, ndf, interface_oh); // get initialized static U_h lambda_h(ndf, interface_oh); static Global_Discretization interface_gd(interface_oh, interface_gh, lambda_h); //-end of interface------------------------------------------------------------- //-Global_Discretization Couples------------------------------------------------ static Global_Discretization_Couple *interface_elastic_foundation_gamma_h_type = new Global_Discretization_Couple(); static Global_Discretization_Couple interface_elastic_foundation_gdc(interface_gd, elastic_foundation_gamma_h_gd, interface_elastic_foundation_gamma_h_type); static Global_Discretization_Couple *interface_rigid_sled_gamma_h_type = new Global_Discretization_Couple(); static Global_Discretization_Couple interface_rigid_sled_gdc(interface_gd, rigid_sled_gamma_h_gd, interface_rigid_sled_gamma_h_type); //-end of global_discretization couples----------------------------------------- //============================================================================== // Step 2. Element_Formulation-- // Q4 quadri-lateral plane elasticity //============================================================================== // DEFINITION of Element_Formulation class Elastic_Contact_Q4 : public Element_Formulation_Couple { public: Elastic_Contact_Q4(Element_Type_Register a) : Element_Formulation_Couple(a) {} // diagonal block formulation; standard Q4 seletive reduced integration Element_Formulation *make(int, Global_Discretization&); Elastic_Contact_Q4(int, Global_Discretization&); // off-diagonal block formulation; contact element Element_Formulation_Couple *make(int, Global_Discretization_Couple&); Elastic_Contact_Q4(int, Global_Discretization_Couple&); }; // for diagonal block formulation Element_Formulation* Elastic_Contact_Q4::make(int en, Global_Discretization& gd) { return new Elastic_Contact_Q4(en,gd); } Elastic_Contact_Q4::Elastic_Contact_Q4(int en, Global_Discretization& gd) : Element_Formulation_Couple(en, gd) { int index; if(gd.type() == elastic_foundation_type) index = 0; else if(gd.type() == rigid_sled_type) index = 1; Quadrature qp2(2, 4); // 2-dimension, 2x2 integration points H1 Z(2, (double*)0, qp2), // Natrual Coordinates N = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE( // Shape Functions "int, int, Quadrature", 4/*nen*/, 2/*nsd*/, qp2), 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 // 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 dv2(d(X).det()); // the Jacobian #ifdef __TEST_BASIS_FORMULATION //============================================================================== // Selective Reduced Integration //============================================================================== // Bubnov-Galerkin Formulation--Shape == Weighting // see Hughes(1987) "Implementation 3" on p.154 C0 e = BASIS("int", ndf), // basis vectors E = BASIS("int", nen), u = e*E, U = (e%e)*(E%E); // deviatoric; 2x2 integration H0 W_x = INTEGRABLE_SUBMATRIX("int, int, H0&", 1, nsd, Nx), Wx, Wy; Wx &= W_x[0][0]; Wy &= W_x[0][1]; C0 stiff_dev = mu_[index]* ( +( ((2.0*Wx*~Wx)+(Wy*~Wy))*((e[0]%e[0])*(E%E)) + (Wy*~Wx) *((e[0]%e[1])*(E%E)) +(Wx*~Wy) *((e[1]%e[0])*(E%E)) + ((2.0*Wy*~Wy)+(Wx*~Wx))*((e[1]%e[1])*(E%E)) ) | dv2); // volumetric; 1-point integration Quadrature qp1(2, 1); // 2-dimension, 1-point integration points H1 z(2, (double*)NULL, qp1), // Natrual Coordinates n = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE( // Shape Functions "int, int, Quadrature", 4/*nen*/, 2/*nsd*/, qp1), zai, eta; zai &= z[0]; eta &= z[1]; n[0] = (1-zai)*(1-eta)/4; n[1] = (1+zai)*(1-eta)/4; n[2] = (1+zai)*(1+eta)/4; n[3] = (1-zai)*(1+eta)/4; H1 x = n*xl; H0 nx = d(n) * d(x).inverse(); // Weighting derivative J dv1(d(x).det()); // the Jacobian // volumetric; 1-point integration H0 w_x = INTEGRABLE_SUBMATRIX("int, int, H0&", 1, nsd, nx), wx, wy; wx &= w_x[0][0]; wy &= w_x[0][1]; // lambda = (K -2/3 mu); for near incompressiblity , K/mu ->infinity, => lambda ~= K (= -B) // lambda = 10^4 mu -> 10^8 mu, for penalty method C0 stiff_vol = lambda_bar[index]*( +( wx*~wx*U[0][0]+wx*~wy*U[0][1]+ wy*~wx*U[1][0]+wy*~wy*U[1][1] ) | dv1 ); stiff &= stiff_vol + stiff_dev; #else 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]); 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(qp2); N_vp = 1.0; H0 mN_vp = N_vp & N_vp & zero; // 3x1 C0 C = ((~mN_vp)*B) | dv2, // 1x3 x 3x2 = 1x2 E = (N_vp*N_vp) | dv2, // 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; if(Matrix_Representation::Assembly_Switch == Matrix_Representation::NODAL_SCALAR) { // compute J2 on guass points H0 e = INTEGRABLE_VECTOR("int, Quadrature", 3, qp2); // deviatoric strain rate e = 0.0; for(int i = 0; i < nen; i++) e += B_dev(i*ndf)*ul[i*ndf] + B_dev(i*ndf+1)*ul[i*ndf+1]; H0 s = 2.0*mu_[index]*e, // s = 2 mu e; J2 = s*s/2.0; // s = sqrt(2*J2(s)); second invariant of deviatoric stress // Project J2 from guass points to nodes using method often used for "Stress smoothing" // 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 the_element_nodal_value &= C0(nen, (double*)0); // nodal J2; returned for global assembly H0 unit(qp2); unit = 1.0; C0 volume = unit | dv2; 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)) | dv2; sum_of_lumped_mass += lumped_mass[i]; } C0 normalized_factor = volume/sum_of_lumped_mass; for(i = 0; i < nen; i++) { lumped_mass[i] *= normalized_factor; the_element_nodal_value[i] = ( (((H0)N[i])*J2) | dv2 ) / lumped_mass[i]; } } else {// plane stress D matrix double a_ = E_[index] / (1.0-pow(v_,2)); double Dv[3][3] = { {a_, a_*v_, 0.0 }, {a_*v_, a_, 0.0 }, {0.0, 0.0, a_*(1.0-v_)/2.0} }; // plane strain D matrix //double Dv[3][3] = { // {mu_[index]*(2.0-2.0/3.0)+K_[index], mu_[index]*-2.0/3.0+K_[index], 0.0 }, // {mu_[index]*-2.0/3.0+K_[index], mu_[index]*(2.0-2.0/3.0)+K_[index], 0.0 }, // {0.0, 0.0, mu_[index]} // }; C0 D = MATRIX("int, int, const double*", 3, 3, Dv[0]); stiff &= ((~B_bar)*D*B_bar)|dv2; } #endif } // off-diagonal block formulation; contact element Element_Formulation_Couple* Elastic_Contact_Q4::make(int en, Global_Discretization_Couple& gdc) { return new Elastic_Contact_Q4(en,gdc); } Elastic_Contact_Q4::Elastic_Contact_Q4(int en, Global_Discretization_Couple& gdc) : Element_Formulation_Couple(en, gdc) { Omega_eh *elem = &(gdc.principal_discretization().omega_h()(en)); double *pp = ((Contact_Omega_eh*)elem)->projection_parameter(), *nv = ((Contact_Omega_eh*)elem)->normal_vector(); C0 n_0(2,nv), n_1(2, nv+2); // unit normal // check lambda_bar < 0 is only useful for two pass treatment, where at the second pass // no new contact search is performed; therefore, ul will contain the values of the last // calculation // lagrangian parameter greater than zero, the constraint is dropped out of active set // surface traction * normal vector = lambda; the lagrangian parameter C0 t(2, 2, (double*)0); // surface traction t[0][0] = ul[0]; t[0][1] = ul[1]; t[1][0] = ul[2]; t[1][1] = ul[3]; double lambda_0 = (double)(t[0]*n_0), lambda_1 = (double)(t[1]*n_1), lambda_bar = (lambda_0+lambda_1)/2.0; // averaged lagrangian parameter static const double EPSILON = 1.e-12; if(lambda_bar > EPSILON) { // lambda > 0, dropped out of the active set stiff &= C0(4, 4, (double*)0); // zero matrix } else { // lambda <= 0, add to the active set by forming stiffness matrix if(gdc.type() == interface_rigid_sled_gamma_h_type) pp += 2; C0 alpha(2,pp); Quadrature qp(1, 2); // 1-D, 2-points guass integration H1 Z(qp), // Natrual Coordinates N = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE( // Shape Functions "int, int, Quadrature", 2/*nen*/, 1/*nsd*/, qp); N[0] = (1.0-Z)/2.0; N[1] = (1.0+Z)/2.0; // coordinate transformation H1 X = N*xl; // Physical Coordinates; dim v = (2 x 2x2} = 2, dim dv = { 1x2 x 2x2 } = 1x2 J d_l(norm(d(X)(0))); // the normalized length C0 zero(0.0); H0 N_lambda = ((~(H0)N) || zero ) & // 2x4 (zero || (~(H0)N)); //------------------------------------------------------------------------------ // | u_0 | | (1-alpha[0]) n_0 alpha[0] n_0 | | u_0_bar | // | | = | | | | // | u_1 | | alpha[1] n_1 (1-alpha[1] n_1 | | u_1_bar | // // u_0 and u_1 are nodes on contact elements // u_0_bar and u_1_bar are nodes on two bodies // N_u_bar = N_u * "alpha-n" matrix in the above //------------------------------------------------------------------------------ H0 N_bar_0 = ( ((H0)N[0])*((1.0-alpha[0])*n_0) + ((H0)N[1])*(alpha[1]*n_1) ), N_bar_1 = ( ((H0)N[0])*(alpha[0]*n_0) + ((H0)N[1])*((1.0-alpha[1])*n_1) ), N_u_bar = ( N_bar_0[0] | zero | N_bar_1[0] | zero ) & ( zero | N_bar_0[1] | zero | N_bar_1[1] ); if(gdc.type() == interface_elastic_foundation_gamma_h_type) stiff &= -(((~N_lambda)*N_u_bar)|d_l); else stiff &= (((~N_lambda)*N_u_bar)|d_l); } } // REGISTRATION of Element(s) Element_Formulation* Element_Formulation::type_list = 0; Element_Type_Register element_type_register_instance; static Elastic_Contact_Q4 elastic_contact_q4_instance(element_type_register_instance); // element type # 0 /* remove the C-style comment in the following code segment or include as a separate file: "\vs\ex\fe\contact_frictionless_sled_on_elastic_foundation\omega_h_n.cpp" // overwirte definition of Element_Formulation_Couple::Element_Formulation_Couple( // int en, Global_Discretization_Couple& gdc) in file omega_h_n.cpp Element_Formulation_Couple::Element_Formulation_Couple(int en, Global_Discretization_Couple& gdc) : Element_Formulation(en, gdc.principal_discretization()), the_subordinate_discretization((Global_Discretization)(gdc.subordinate_discretization())) { int index = 0; if(gdc.type() == interface_rigid_sled_gamma_h_type) index = 1; Omega_eh *elem = &(gdc.principal_discretization().omega_h()(en)); s_element_no = ((Contact_Omega_eh*)elem)->associate_element()[index]; Omega_eh& s_element = the_subordinate_discretization.omega_h()(s_element_no); s_nen = s_element.element_node_no(), s_nsd = the_subordinate_discretization.omega_h()[ s_element[0] ].no_of_dim(); s_ndf = the_subordinate_discretization.u_h()[s_element[0]].no_of_dim(); s_material_type_no = s_element.material_type_no(); s_xl &= the_subordinate_discretization.element_coordinate(s_element_no); s_gl &= the_subordinate_discretization.element_fixed_variable(s_element_no); s_ul &= the_subordinate_discretization.element_free_variable(s_element_no); } void Matrix_Representation_Couple::assembly(int element_no1,int element_no2, int element_incr, int nodal_load_flag) { // add nodal force (natural boundary conditions) from the principal rhs int total_node_no = (Matrix_Representation::global_discretization()).omega_h().total_node_no(), // no distinct number ndf = (Matrix_Representation::global_discretization()).u_h()[0].no_of_dim(); if(nodal_load_flag) { for(int i = 0; i < total_node_no; i++) { Nodal_Constraint& gh = (Matrix_Representation::global_discretization()).gh_on_gamma_h()[i]; for(int j = 0; j < ndf; j++) { int equation_no = eqn[i][j]; if(gh(j) == gh_on_Gamma_h::Neumann && gh[j] != 0.0) the_rhs->global_tensor()[equation_no] += gh[j]; } } } // add stiffness, reaction (essential boundary conditions) and distributed load if(element_no2 == -1) element_no2 = element_no1; // assemble only one element if(element_no1 == -1) { // assemble all distinct elements int total_element_no = principal_discretization().omega_h().omega_eh_array().length(); for(int i = 0; i < total_element_no; i++) { Omega_eh& elem = principal_discretization().omega_h()(i); int etn = elem.element_type_no(); Element_Formulation *element_type = Element_Formulation::type_list; if(etn > 0) for(int j = 0; j < etn; j++) element_type = element_type->next; Element_Formulation ef; ef.rep() = ((Element_Formulation_Couple*)element_type)->make(i, // form element Global_Discretization_Couple(principal_discretization(), the_subordinate_discretization, gdc_type_id) ); if(Assembly_Switch == LHS || Assembly_Switch == LHS_AND_RHS || Assembly_Switch == ALL) { // a. add element stiffness matrix Element_Tensor element_lhs_matrix(i, ((Element_Formulation_Couple*)(ef.rep()))->lhs(), *this); (*the_lhs) += element_lhs_matrix; } if(Assembly_Switch == RHS || Assembly_Switch == LHS_AND_RHS || Assembly_Switch == ALL) { // b. add reaction and distributed load Element_Tensor element_rhs_vector(i, ((Element_Formulation_Couple*)(ef.rep()))->rhs(), *this); (*the_rhs) += element_rhs_vector; } if(Assembly_Switch == SUBORDINATE_RHS || Assembly_Switch == ALL) { // c. add reaction and distributed load Element_Tensor element_subordinate_rhs_vector(i, ((Element_Formulation_Couple*)(ef.rep()))->subordinate_rhs(), *this); (*the_subordinate_rhs) += element_subordinate_rhs_vector; } } } else { // assemble specified range of elements for(int i = element_no1; i <= element_no2; i+=element_incr) { Omega_eh& elem = principal_discretization().omega_h()(i); int etn = elem.element_type_no(); Element_Formulation *element_type = Element_Formulation::type_list; if(etn > 0) for(int j = 0; j < etn; j++) element_type = element_type->next; Element_Formulation ef; ef.rep() = ((Element_Formulation_Couple*)element_type)->make(i, // form element Global_Discretization_Couple(principal_discretization(), the_subordinate_discretization, gdc_type_id) ); if(Assembly_Switch == LHS || Assembly_Switch == LHS_AND_RHS || Assembly_Switch == ALL) { // a. add element stiffness matrix Element_Tensor element_lhs_matrix(i, ((Element_Formulation_Couple*)(ef.rep()))->lhs(), *this); (*the_lhs) += element_lhs_matrix; } if(Assembly_Switch == RHS || Assembly_Switch == LHS_AND_RHS || Assembly_Switch == ALL) { // b. add reaction and distributed load Element_Tensor element_rhs_vector(i, ((Element_Formulation_Couple*)(ef.rep()))->rhs(), *this); (*the_rhs) += element_rhs_vector; } if(Assembly_Switch == SUBORDINATE_RHS || Assembly_Switch == ALL) { // c. add reaction and distributed load Element_Tensor element_subordinate_rhs_vector(i, ((Element_Formulation_Couple*)(ef.rep()))->subordinate_rhs(), *this); (*the_subordinate_rhs) += element_subordinate_rhs_vector; } } } } // overwirte definition of Coupled_Global_Tensor::operator+=( // Element_Tensor& element_tensor) in file omega_h_n.cpp C0& Coupled_Global_Tensor::operator+=(Element_Tensor& element_tensor) { int en = element_tensor.element_no(); int index = 0; if(the_matrix_representation_couple.subordinate_discretization().type() == rigid_sled_type) index = 1; Omega_eh *elem = &(the_matrix_representation_couple.principal_discretization().omega_h()(en)); int s_element_no = ((Contact_Omega_eh*)elem)->associate_element()[index]; //Omega_eh elem = element_tensor.matrix_representation().global_discretization().omega_h()(en); int nen = elem->element_node_no(), first_node = (*elem)[0], ndf = the_matrix_representation_couple.global_discretization().u_h()[first_node].no_of_dim(); if(element_tensor.element_tensor().type() == vector_type) { for(int i = 0; i < nen; i++) { int node_no = (*elem)[i]; for(int j = 0; j < ndf; j++) { int eq_no = the_matrix_representation_couple.eqn[node_no][j]; if(eq_no == Matrix_Representation::Tied_Node) { int tie_node_no = the_matrix_representation_couple.global_discretization().u_h().u_h_array()[node_no].tie_node_no(); eq_no = the_matrix_representation_couple.eqn[tie_node_no][j]; } if(eq_no >=0 ) global_tensor()[eq_no] += element_tensor.element_tensor()[i*ndf+j]; } } } else if(element_tensor.element_tensor().type() == matrix_type) { Omega_eh *s_elem = &( the_matrix_representation_couple.subordinate_discretization().omega_h()(s_element_no)); int s_nen = s_elem->element_node_no(), s_first_node = (*s_elem)[0], s_ndf = the_matrix_representation_couple.subordinate_discretization().u_h()[s_first_node].no_of_dim(); for(int i = 0; i < nen; i++) for(int j = 0; j < s_nen; j++) { int node_no1 = (*elem)[i], node_no2 = (*s_elem)[j]; for(int k = 0; k < ndf; k++) for(int l = 0; l < s_ndf; l++) { int eqn1 = the_matrix_representation_couple.eqn[node_no1][k], eqn2 = the_matrix_representation_couple.var[node_no2][l]; // converts tied node if(eqn1 == Matrix_Representation::Tied_Node) { int tie_node_no = the_matrix_representation_couple.global_discretization().u_h().u_h_array()[node_no1].tie_node_no(); eqn1 = the_matrix_representation_couple.eqn[tie_node_no][k]; } if(eqn2 == Matrix_Representation::Tied_Node) { int tie_node_no = the_matrix_representation_couple.subordinate_discretization().u_h().u_h_array()[node_no2].tie_node_no(); eqn2 = the_matrix_representation_couple.var[tie_node_no][l]; } if(eqn1 >= 0 && eqn2 >= 0) // if none of the two dof is a fixed degree of freedom global_tensor()[eqn1][eqn2] += element_tensor.element_tensor()[i*ndf+k][j*s_ndf+l]; } } } return global_tensor(); } */ #include "\vs\ex\fe\contact_frictionless_sled_on_elastic_foundation\omega_h_n.cpp" //============================================================================== // Step 3: Matrix Representation and Solution Phase //============================================================================== // | K_0 Q_0^t 0 | | u_0 | | f_0 | // | Q_0 Q_1 | | lambda | = | f_i | // | 0 Q_1^t K_1 | | u_1 | | f_1 | // w/ rigid body mode constraints K_0 and K_1 are non-singular; u_0 and u_1 can // both be eliminated. //============================================================================== #include void output(Global_Discretization&, Global_Discretization&, ofstream&); void femap_output_model(Global_Discretization&, Global_Discretization&, ofstream&); void femap_output_data(int, Global_Discretization&, Global_Discretization&, ofstream&, C0&, C0&); int main() { C0 K_0, Q_0, K_1, Q_1, f_0, f_i, f_1, lambda, u_0, u_1; time_t time0, time1; time0 = time(0); // initial contact searching interface_oh.contact_searching_algorithm(interface_oh, elastic_foundation_gamma_h_gd, rigid_sled_gamma_h_gd, lambda_h, interface_gh, 0.05); // tolerance gap // form initial matrix representations Matrix_Representation mr_K_0(elastic_foundation_gd); Matrix_Representation mr_K_1(rigid_sled_gd); Matrix_Representation_Couple mrc_Q_0(interface_elastic_foundation_gdc, 0, 0, &(mr_K_0.rhs()) ); Matrix_Representation_Couple mrc_Q_1(interface_rigid_sled_gdc, 0, &(mrc_Q_0.rhs()), &(mr_K_1.rhs()) ); ofstream ofs1("D:\\FEMAP601\\contact.neu", ios::out | ios::trunc); femap_output_model(elastic_foundation_gd, rigid_sled_gd, ofs1); // loop of incremental horizaontal sliding for(int i = 0; i < 9; i++) { // assemble matrices mr_K_0.assembly(); K_0 &= (C0)(mr_K_0.lhs()); ofs << "K_0:" << endl << K_0 << endl; time1 = time(0); ofs << "elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; mr_K_1.assembly(); K_1 &= (C0)(mr_K_1.lhs()); ofs << "K_1:" << endl << K_1 << endl; time1 = time(0); ofs << "elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; // assemble off-diagonal matrices mrc_Q_0.assembly(); Q_0 &= (C0)(mrc_Q_0.lhs()); ofs << "Q_0:" << endl << Q_0 << endl; time1 = time(0); ofs << "elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; mrc_Q_1.assembly(); Q_1 &= (C0)(mrc_Q_1.lhs()); f_0 &= (C0)(mr_K_0.rhs()); f_1 &= (C0)(mr_K_1.rhs()); f_i &= (C0)(mrc_Q_0.rhs()); ofs << "Q_1:" << endl << Q_1 << endl; ofs << "f_0:" << endl << f_0 << endl; ofs << "f_1:" << endl << f_1 << endl; ofs << "f_i:" << endl << f_i << endl; time1 = time(0); ofs << "elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; // static condensation ofs << endl << "solution phase begins... " << endl; Cholesky dK0(K_0); C0 K0_inv = dK0.inverse(); C0 QK_inv_0 = Q_0*K0_inv; C0 QKQ_0 = QK_inv_0*(~Q_0); Cholesky dK1(K_1); C0 K1_inv = dK1.inverse(); C0 QK_inv_1 = Q_1*K1_inv; C0 QKQ_1 = QK_inv_1*(~Q_1); C0 QKQ = QKQ_0+QKQ_1; ofs << QKQ << endl; Cholesky mdQKQ(QKQ, 1.e-8); double integer; modf(log10(mdQKQ.cond()), &integer); ofs << "Warnning: lost of "; if(integer > 4.0) ofs << int(integer); else ofs << "4"; ofs << " significant digits in modified Cholesky decomposition." << endl; lambda &= mdQKQ * (QK_inv_0*f_0+QK_inv_1*f_1-f_i); // static condensation u_0 &= dK0*(f_0-(~Q_0)*lambda); u_1 &= dK1*(f_1-(~Q_1)*lambda); time1 = time(0); ofs << "elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; // update solution lambda_h = lambda; lambda_h = interface_gd.gh_on_gamma_h(); u_0_h = u_0; u_0_h = elastic_foundation_gd.gh_on_gamma_h(); u_1_h = u_1; u_1_h = rigid_sled_gd.gh_on_gamma_h(); // compute J2 projection on nodes Matrix_Representation::Assembly_Switch = Matrix_Representation::NODAL_SCALAR; mr_K_0.assembly(FALSE); mr_K_1.assembly(FALSE); C0 nodal_J2_0 = mr_K_0.global_nodal_value(), nodal_J2_1 = mr_K_1.global_nodal_value(); femap_output_data(i, elastic_foundation_gd, rigid_sled_gd, ofs1, nodal_J2_0, nodal_J2_1); Matrix_Representation::Assembly_Switch = Matrix_Representation::ALL; // set back to default // output solution ofs << "contact forces:" << endl << lambda_h; ofs << "elastic foundation displacement:" << endl << u_0_h; ofs << "rigid sled displacement:" << endl << u_1_h; if(i == 2) { ofstream ofs1("temp1.cpp", ios::out | ios::trunc); output(elastic_foundation_gd, rigid_sled_gd, ofs1); } else if(i == 5) { ofstream ofs2("temp2.cpp", ios::out | ios::trunc); output(elastic_foundation_gd, rigid_sled_gd, ofs2); } else if(i == 8) { ofstream ofs3("temp3.cpp", ios::out | ios::trunc); output(elastic_foundation_gd, rigid_sled_gd, ofs3); } if(i < 8) { // not the last step for(int j = 8; j < 12; j++) { // incremental horizaontal motion by update B.C. on the rigid sled int node_order = rigid_sled_gd.gh_on_gamma_h().node_order(j); rigid_sled_gd.gh_on_gamma_h().gh_array()[node_order][0] = ((double)(i+1)); } ofs << "step " << i << " :" << endl; rigid_sled_gd.u_h() = rigid_sled_gd.gh_on_gamma_h(); // contact searching with lagrangian parameter and matrix represenation redefined for each step interface_oh.contact_searching_algorithm(interface_oh, elastic_foundation_gamma_h_gd, rigid_sled_gamma_h_gd, lambda_h, interface_gh, mr_K_0, mr_K_1, mrc_Q_0, mrc_Q_1, 0.05); // tolerance gap } } time1 = time(0); ofs << "total elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; ofs.close(); return 0; } #include "output_contact.cpp" //AutoCAD r13 ADS interface program void output(Global_Discretization& elastic_foundation_gd, Global_Discretization& rigid_sled_gd, ofstream& ofs1) { // input for AutoCAD ADS program // include header Omega_h &elastic_foundation_oh = elastic_foundation_gd.omega_h(), &rigid_sled_oh = rigid_sled_gd.omega_h(); U_h &elastic_foundation_uh = elastic_foundation_gd.u_h(), &rigid_sled_uh = rigid_sled_gd.u_h(); gh_on_Gamma_h &elastic_foundation_gh = elastic_foundation_gd.gh_on_gamma_h(), &rigid_sled_gh = rigid_sled_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 = " << (elastic_foundation_oh.total_node_no()+rigid_sled_oh.total_node_no()) << ";" << endl << " define_node(" << (elastic_foundation_oh.total_node_no()+rigid_sled_oh.total_node_no()) << ");" << endl << " window_range(center_x, center_y);" << endl << " domain.total_element_no = " << (elastic_foundation_oh.total_element_no()+rigid_sled_oh.total_element_no()) << ";" << endl << " define_element(" << (elastic_foundation_oh.total_element_no()+rigid_sled_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 < elastic_foundation_oh.total_node_no(); i++) { ofs1 << " na = &(domain.node_array[" << i << "]);" << endl << " na->node_no = " << elastic_foundation_oh.node_array()[i].node_no() << "; na->ndf = " << elastic_foundation_uh.u_h_array()[i].no_of_dim() << ";" << "; na->nsd = " << elastic_foundation_oh.node_array()[i].no_of_dim() << ";" << endl << " Spoint(na->coord, " << elastic_foundation_oh.node_array()[i][0] << ", " << elastic_foundation_oh.node_array()[i][1] << ", 0.0);" << endl << " na->u = (ads_real*)calloc(" << elastic_foundation_uh.u_h_array()[i].no_of_dim() << ", sizeof(ads_real));" << endl; for(int j = 0; j < elastic_foundation_uh.u_h_array()[i].no_of_dim(); j++) ofs1 << " na->u[" << j << "] = " << elastic_foundation_uh.u_h_array()[i][j] << "; "; ofs1 << endl; ofs1 << " na->constraint_flag = (int*)calloc(" << elastic_foundation_gh.gh_array()[i].no_of_dim() << ", sizeof(int));" << endl; for(j = 0; j < elastic_foundation_gh.gh_array()[i].no_of_dim(); j++) ofs1 << " na->constraint_flag[" << j << "] = " << elastic_foundation_gh.gh_array()[i](j) << "; "; ofs1 << endl; //ofs1 << " na->f = (ads_real*)calloc(" << elastic_foundation_gh.gh_array()[i].no_of_dim() << ", sizeof(ads_real));" << endl; //for(int j = 0; j < elastic_foundation_gh.gh_array()[i].no_of_dim(); j++) // ofs1 << " na->f[" << j << "] = " << elastic_foundation_gh.gh_array()[i][j] << "; "; //ofs1 << endl; } for(i = 0; i < rigid_sled_oh.total_node_no(); i++) { ofs1 << " na = &(domain.node_array[" << (elastic_foundation_oh.total_node_no()+i) << "]);" << endl << " na->node_no = " << (elastic_foundation_oh.total_node_no()+rigid_sled_oh.node_array()[i].node_no()) << "; na->ndf = " << rigid_sled_uh.u_h_array()[i].no_of_dim() << ";" << "; na->nsd = " << rigid_sled_oh.node_array()[i].no_of_dim() << ";" << endl << " Spoint(na->coord, " << rigid_sled_oh.node_array()[i][0] << ", " << rigid_sled_oh.node_array()[i][1] << ", 0.0);" << endl << " na->u = (ads_real*)calloc(" << rigid_sled_uh.u_h_array()[i].no_of_dim() << ", sizeof(ads_real));" << endl; for(int j = 0; j < rigid_sled_uh.u_h_array()[i].no_of_dim(); j++) ofs1 << " na->u[" << j << "] = " << rigid_sled_uh.u_h_array()[i][j] << "; "; ofs1 << endl; ofs1 << " na->constraint_flag = (int*)calloc(" << rigid_sled_gh.gh_array()[i].no_of_dim() << ", sizeof(int));" << endl; for(j = 0; j < rigid_sled_gh.gh_array()[i].no_of_dim(); j++) ofs1 << " na->constraint_flag[" << j << "] = " << rigid_sled_gh.gh_array()[i](j) << "; "; ofs1 << endl; //ofs1 << " na->f = (ads_real*)calloc(" << rigid_sled_gh.gh_array()[i].no_of_dim() << ", sizeof(ads_real));" << endl; //for(int j = 0; j < rigid_sled_gh.gh_array()[i].no_of_dim(); j++) // ofs1 << " na->f[" << j << "] = " << rigid_sled_gh.gh_array()[i][j] << "; "; //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(" << (elastic_foundation_oh.total_element_no()+rigid_sled_oh.total_element_no()) << ", sizeof(Omega_eh));" << endl; for(i = 0; i < elastic_foundation_oh.total_element_no(); i++) { ofs1 << " ea = &(domain.omega_eh_array[" << i << "]);" << endl << " ea->element_no = " << elastic_foundation_oh.omega_eh_array()[i].element_no() << "; ea->element_node_no = " << elastic_foundation_oh.omega_eh_array()[i].element_node_no() << "; " << endl << " ea->element_type_no = " << elastic_foundation_oh.omega_eh_array()[i].element_type_no() << "; ea->material_type_no = " << elastic_foundation_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 < elastic_foundation_oh.omega_eh_array()[i].element_node_no(); j++) ofs1 << "ea->node_no_array[" << j << "] = " << elastic_foundation_oh.omega_eh_array()[i][j] << "; "; ofs1 << endl; } int node_no = elastic_foundation_oh.total_node_no(); for(i = 0; i < rigid_sled_oh.total_element_no(); i++) { ofs1 << " ea = &(domain.omega_eh_array[" << (elastic_foundation_oh.total_element_no()+i) << "]);" << endl << " ea->element_no = " << (elastic_foundation_oh.total_element_no()+rigid_sled_oh.omega_eh_array()[i].element_no()) << "; ea->element_node_no = " << rigid_sled_oh.omega_eh_array()[i].element_node_no() << "; " << endl << " ea->element_type_no = " << rigid_sled_oh.omega_eh_array()[i].element_type_no() << "; ea->material_type_no = " << rigid_sled_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 < rigid_sled_oh.omega_eh_array()[i].element_node_no(); j++) ofs1 << "ea->node_no_array[" << j << "] = " << (node_no+rigid_sled_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(); }