-----------------------------------------------------------------------------------------------------------------------------------------
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 <[email protected]> wrote:
>
>
> On Fri, Dec 4, 2015 at 4:22 PM, Xujun Zhao <[email protected]> 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
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
