-----------------------------------------------------------------------------------------------------------------------------------------
EquationSystems elasticity_es (mesh); // mesh is read from files. ElasticitySystem& elasticity_system = elasticity_es.add_system<ElasticitySystem> ("Elasticity"); unsigned int U_var = 0, V_var = 0, W_var = 0; U_var = elasticity_system.add_variable ("U", FIRST); V_var = elasticity_system.add_variable ("V", FIRST); if(dim==3) W_var = elasticity_system.add_variable ("W", FIRST); elasticity_es.init (); ---------------------------------------------------- subroutine -------------------------------------------------- // ================================================================================== void ElasticitySystem::build_nodal_force_gravity() { START_LOG("build_nodal_force_gravity()", "ElasticitySystem"); // Get a constant reference to the mesh object and mesh dimension. const MeshBase& mesh = this->get_mesh(); const unsigned int dim = mesh.mesh_dimension(); const Real force_den = 1.0; // force density of the solid this->rhs->zero(); // Numeric ids corresponding to each variable in the system const unsigned int u_var = this->variable_number ("U"); // u_var = 0 const unsigned int v_var = this->variable_number ("V"); // v_var = 1 unsigned int w_var = 0; if(dim==3) w_var = this->variable_number ("W"); // w_var = 2 // Build a FEBase object of the specified type // Note: typedef FEGenericBase<Real> FEBase in libMesh; FEType fe_vel_type = this->variable_type(u_var); AutoPtr<FEBase> fe_vel (FEBase::build(dim, fe_vel_type) ); //fe_vel->print_info(libMesh::out); // Gauss Quadrature rule and tell the finite element objects QGauss qrule (dim, fe_vel_type.default_quadrature_order()); //3^dim pts //QGauss qrule (dim, FOURTH ); // THIRD/FOURTH/FIFTH fe_vel->attach_quadrature_rule (&qrule); // The element Jacobian * quadrature weight at each integration point. const std::vector<Real>& JxW = fe_vel->get_JxW(); const std::vector<std::vector<Real> >& phi = fe_vel->get_phi(); //const std::vector<std::vector<RealGradient> >& dphi = fe_vel->get_dphi(); // A reference to the DofMap object for this system. const DofMap & dof_map = this->get_dof_map(); std::vector<dof_id_type> dof_indices; std::vector<dof_id_type> dof_indices_u, dof_indices_v, dof_indices_w; // element nodal vector DenseVector<Number> Fe; DenseSubVector<Number> Fu(Fe), Fv(Fe), Fw(Fe); // Now we will loop over all the elements in the mesh that live // on the local processor, and compute the element rhs Fe. MeshBase::const_element_iterator el = mesh.active_local_elements_begin(); const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end(); for ( ; el != end_el; ++el) { // Store a pointer to the element we are currently working on. const Elem* elem = *el; // Get the degree of freedom indices for the current element. dof_map.dof_indices (elem, dof_indices); dof_map.dof_indices (elem, dof_indices_u, u_var); dof_map.dof_indices (elem, dof_indices_v, v_var); const unsigned int n_dofs = dof_indices.size(); const unsigned int n_u_dofs = dof_indices_u.size(); const unsigned int n_v_dofs = dof_indices_v.size(); unsigned int n_w_dofs = 0; if(dim==3) { dof_map.dof_indices (elem, dof_indices_w, w_var); n_w_dofs = dof_indices_w.size(); } // NOTE: here JxW and dphi and other element quantities are not computed. // This will be done in the elem loop after fe->reinit() fe_vel->reinit (elem); //elem->print_info(libMesh::out); libMesh::out<<"\n\n"; //fe_vel->print_info(libMesh::out); libMesh::out<<"\n"; // Zero the element right-hand-side vector before summing them. // DenseSubVector.reposition () member takes the (row_offset, row_size) Fe.resize (n_dofs); // Ke.resize (n_dofs, n_dofs); Fu.reposition (u_var*n_u_dofs, n_u_dofs); Fv.reposition (v_var*n_u_dofs, n_v_dofs); if(dim==3) Fw.reposition (w_var*n_u_dofs, n_w_dofs); // Loop over each Gauss pt for (unsigned int qp=0; qp<qrule.n_points(); qp++) { // Assemble the subvector for (unsigned int i=0; i<n_u_dofs; i++) { Fu(i) += JxW[qp]*phi[i][qp]*force_den; } } // If this assembly program were to be used on an adaptive mesh, // we would have to apply any hanging node constraint equations. dof_map.constrain_element_vector (Fe, dof_indices); // Add the element matrix and rhs vector to the global system. this->rhs->add_vector (Fe, dof_indices); } // end for el-loop this->rhs->close(); // test view vector this->solution = this->rhs->clone(); // this->solution->print_global(); STOP_LOG("build_nodal_force_gravity()", "ElasticitySystem"); } On Fri, Dec 4, 2015 at 5:24 PM, John Peterson <jwpeter...@gmail.com> wrote: > > > On Fri, Dec 4, 2015 at 4:22 PM, Xujun Zhao <xzha...@gmail.com> wrote: > >> Hi all, >> >> I try to evaluate nodal force according to a given force density, say, it >> is a constant 1. When I use HEX8 elem, the results are correct. but if I >> use TET4 elem, it gives rise to wrong results. I am using Cubit to >> generate >> mesh for a solid sphere, then it is read by libMesh. I doubt if there are >> any inconsistencies between them, for example, different node # order. Or >> I >> made some mistakes that I didn't notice. Anyone has any ideas? >> >> On the other hand, can I generate TET4 elements from libMesh? so that I >> can >> make a test. >> > > Post your code. There's no way to tell what's going wrong without it... > > Either make a gist on github or share a google doc, the mailing list does > not allow attachments. > > -- > John > ------------------------------------------------------------------------------ Go from Idea to Many App Stores Faster with Intel(R) XDK Give your users amazing mobile app experiences with Intel(R) XDK. Use one codebase in this all-in-one HTML5 development environment. Design, debug & build mobile apps & 2D/3D high-impact games for multiple OSs. http://pubads.g.doubleclick.net/gampad/clk?id=254741911&iu=/4140 _______________________________________________ Libmesh-users mailing list Libmesh-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/libmesh-users