#include #include #include #include #include #include "\vs\include\vs.h" #include "\vs\ex\fe\include\fe.h" static ofstream ofs("temp.cpp", ios::out | ios::trunc); // define node, element, constraint, and load //============================================================================== // Step 1: Global_Discretizations //============================================================================== #include "model.cpp" //============================================================================== // Step 2. Element_Formulation--tetraheda element //============================================================================== // DEFINITION of Element_Formulation class ElasticT4 : public Element_Formulation { public: ElasticT4(Element_Type_Register a) : Element_Formulation(a) {} Element_Formulation *make(int, Global_Discretization&); ElasticT4(int, Global_Discretization&); }; Element_Formulation* ElasticT4::make(int en, Global_Discretization& gd) { return new ElasticT4(en,gd); } // 3d D matrix static const double E_ = 2.9e7, v_ = 0.32, mu_ = E_/2.0/(1.+v_); // ANSI 4130 Steel static const double a_ = E_*(1.-v_)/(1.-v_)/(1.-2.*v_); static const double Dv[6][6] = { {a_, a_* v_/(1.-v_), a_* v_/(1.-v_), 0.0, 0.0, 0.0 }, {a_* v_/(1.-v_), a_, a_* v_/(1.-v_), 0.0, 0.0, 0.0 }, {a_* v_/(1.-v_), a_* v_/(1.-v_), a_, 0.0, 0.0, 0.0 }, {0.0, 0.0, 0.0, a_* (1.-2.*v_)/2./(1.-v_), 0.0, 0.0 }, {0.0, 0.0, 0.0, 0.0, a_* (1.-2.*v_)/2./(1.-v_), 0.0 }, {0.0, 0.0, 0.0, 0.0, 0.0, a_* (1.-2.*v_)/2./(1.-v_)} }; C0 D = MATRIX("int, int, const double*", 6, 6, Dv[0]); static const double i_dev[6][6] = { {(1.0-1.0/3.0), -1.0/3.0, -1.0/3.0, 0.0, 0.0, 0.0}, {-1.0/3.0, (1.0-1.0/3.0), -1.0/3.0, 0.0, 0.0, 0.0}, {-1.0/3.0, -1.0/3.0, (1.0-1.0/3.0), 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 1.0} }; C0 I_dev = MATRIX("int, int, const double*", 6, 6, i_dev[0]); ElasticT4::ElasticT4(int en, Global_Discretization& gd) : Element_Formulation(en, gd) { // A. Declarations #if defined(__DEGENERATED_CUBE) // Fourth-order accurate, 2x2x2 quadrature rule in 3-d, on 6 face centers double s3 = 1./sqrt(3.0), w[8] = {1., 1., 1., 1., 1., 1., 1., 1.}, coord[8][3] = {{-s3, -s3, -s3}, {s3, -s3, -s3}, {s3, s3, -s3}, {-s3, s3, -s3}, {-s3, -s3, s3}, {s3, -s3, s3}, {s3, s3, s3}, {-s3, s3, s3}}; Quadrature qp(w, coord[0], 8, 3/*nsd*/); // Fourth-order accurate, six-point quadrature rule in 3-d, on 6 face centers; this // quadrature rule produces singular d(X).inverse() for tetrahedra degenerated from cube //double w[6] = {4./3., 4./3., 4./3., 4./3., 4./3., 4./3.}, // coord[6][3] = {{1., 0., 0.}, {-1., 0., 0.}, {0., 1., 0.}, // {0., -1., 0.}, {0., 0., 1.}, {0., 0., -1.}}; //Quadrature qp(w, coord[0], 6, 3/*nsd*/); H1 Z(3, (double*)0, qp), // Natrual Coordinates Zai, Eta, Zeta; Zai &= Z[0]; Eta &= Z[1], Zeta = Z[2]; // B. Definitions // B1. Shape function for trilinear hexahedral element p. 124 in Hughes[1987] double z[8][3] = { {-1., -1., -1.}, {1., -1., -1.}, {1., 1., -1.}, {-1., 1., -1.}, {-1., -1., 1.}, {1., -1., 1.}, {1., 1., 1.}, {-1., 1., 1.}}; H1 n = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 8/*nen*/, 3/*nsd*/, qp); for(int i = 0; i < 8; i++) n[i] = (1.+z[i][0]*Zai)*(1.+z[i][1]*Eta)*(1.+z[i][2]*Zeta)/8.; // Shape function of wedge-shaped element formed by degenerating the trilinear hexahedral element, p.125 in Hughes[1987] H1 np = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 6, 3, qp); np[0] = n[0]; np[1] = n[1]; np[2] = n[4]; np[3] = n[5]; np[4] = n[2] + n[3]; np[5] = n[6] + n[7]; // Shape function of tetrahedra element formed by degenerating the wedge-shaped element H1 N = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 4, 3, qp); N[0] = np[0]; N[1] = np[1]; N[2] = np[4]; N[3] = np[2] + np[3] + np[5]; #else // Shape function of tetrahedra element using tetrahedra natural coordinates, // Zienkiewicz & Taylor, 1989, p.177 double w[1] = {1.}, coord[1][3] = {{.25, .25, .25}}; Quadrature qp(w, coord[0], 1, 3); H1 L(3, (double*)0, qp), N = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 4, 3, qp); N[0] = L[0]; N[1] = L[1]; N[2] = L[2]; N[3] = 1.0 - L[0] - L[1] - L[2]; #endif // B2. coordinate transformation H1 X = N*xl; // Physical Coordinates // d(N)/d(X) = d(N)/d(Z)*(d(X)/d(Z))^-1 => {4x3}*{3x3}={4x3} 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, nsd, Nx), // dim ~Nx = {3x4}, dim wx = {{3x1}x4} Wx, Wy, Wz, b[4], B; Wx &= W_x[0][0]; Wy &= W_x[0][1], Wz &= W_x[0][2]; // dim Wx = {4x1} C0 zero = C0(0.0); for(int j = 0; j < nen; j++) { H0 wx = Wx[j][0], wy = Wy[j][0], wz = Wz[j][0]; b[j] &= (wx | zero | zero) & (zero | wy | zero) & (zero | (zero | wz ) ) & (wy | wx | zero) & (zero | wz | wy ) & (wz | zero | wx ); } B &= b[0] | b[1] | b[2] | b[3]; // dim B = {6x12} if(Matrix_Representation::Assembly_Switch == Matrix_Representation::NODAL_SCALAR) { // compute J2 on guass points H0 B_dev = I_dev * B; H0 e = INTEGRABLE_VECTOR("int, Quadrature", 6, qp); // 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_*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(qp); unit = 1.0; C0 volume = 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 = 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) | dv ) / lumped_mass[i]; } } else stiff &= ((~B) * (D * B)) | dv; // {12x6}*{6x6}*{6x12}={12x12} } // REGISTRATION of Element(s) Element_Formulation* Element_Formulation::type_list = 0; Element_Type_Register element_type_register_instance; static ElasticT4 tetrahedra_instance(element_type_register_instance); // element type # 0 //============================================================================== // Step 3: Matrix Formulation and Solution Phase //============================================================================== //void output(Global_Discretization&); #include "output.cpp" int main() { // INSTANTIATION of Global_Discretization int ndf = 3; Omega_h oh; gh_on_Gamma_h gh(ndf, oh); U_h uh(ndf, oh); Global_Discretization gd(oh, gh, uh); time_t time0 = time(0); Matrix_Representation mr(gd); ofstream ofs1("D:\\FEMAP601\\test1.neu", ios::out | ios::trunc); femap_output_model(gd, ofs1); mr.assembly(); C0 u = ((C0)(mr.rhs())) / ((C0)(mr.lhs())); gd.u_h() = u; // update free degree of freedom gd.u_h() = gd.gh_on_gamma_h(); // update fixed degree of freedom ofs << "nodal displacement:" << endl; ofs << gd.u_h(); // output solution // post-processing for compute J2 projection from guass point to nodes Matrix_Representation::Assembly_Switch = Matrix_Representation::NODAL_SCALAR; mr.assembly(FALSE); C0 nodal_J2 = mr.global_nodal_value(); ofs << "nodal J2:" << endl; for(int i = 0; i < oh.total_node_no(); i++) ofs << "{" << oh.node_array()[i].node_no() << " | " << nodal_J2[i] << " }" << endl; Matrix_Representation::Assembly_Switch = Matrix_Representation::ALL; // set back to default femap_output_data(gd, ofs1, nodal_J2); time_t time1 = time(0); ofs << "elipsed time: " << (difftime(time1, time0)/60) << " minutes." << endl; ofs.close(); //output(gd); return 0; }