-----------------------------------------------------------------------------------------------------------------------------------------

  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

Reply via email to