The source file systems_of_equations_ex4.C with comments:

Systems Example 4 - Linear Elastic Cantilever

By David Knezevic

In this example we model a homogeneous isotropic cantilever using the equations of linear elasticity. We set the Poisson ratio to \nu = 0.3 and clamp the left boundary and apply a vertical load at the right boundary.



C++ include files that we need
        #include <iostream>
        #include <algorithm>
        #include <math.h>
        
libMesh includes
        #include "libmesh/libmesh.h"
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_generation.h"
        #include "libmesh/exodusII_io.h"
        #include "libmesh/gnuplot_io.h"
        #include "libmesh/linear_implicit_system.h"
        #include "libmesh/equation_systems.h"
        #include "libmesh/fe.h"
        #include "libmesh/quadrature_gauss.h"
        #include "libmesh/dof_map.h"
        #include "libmesh/sparse_matrix.h"
        #include "libmesh/numeric_vector.h"
        #include "libmesh/dense_matrix.h"
        #include "libmesh/dense_submatrix.h"
        #include "libmesh/dense_vector.h"
        #include "libmesh/dense_subvector.h"
        #include "libmesh/perf_log.h"
        #include "libmesh/elem.h"
        #include "libmesh/boundary_info.h"
        #include "libmesh/zero_function.h"
        #include "libmesh/dirichlet_boundaries.h"
        #include "libmesh/string_to_enum.h"
        #include "libmesh/getpot.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
Matrix and right-hand side assemble
        void assemble_elasticity(EquationSystems& es,
                                 const std::string& system_name);
        
Define the elasticity tensor, which is a fourth-order tensor i.e. it has four indices i,j,k,l
        Real eval_elasticity_tensor(unsigned int i,
                                    unsigned int j,
                                    unsigned int k,
                                    unsigned int l);
        
Begin the main program.
        int main (int argc, char** argv)
        {
Initialize libMesh and any dependent libaries
          LibMeshInit init (argc, argv);
        
Initialize the cantilever mesh
          const unsigned int dim = 2;
        
Skip this 2D example if libMesh was compiled as 1D-only.
          libmesh_example_assert(dim <= LIBMESH_DIM, "2D support");
        
Create a 2D mesh distributed across the default MPI communicator.
          Mesh mesh(init.comm(), dim);
          MeshTools::Generation::build_square (mesh,
                                               50, 10,
                                               0., 1.,
                                               0., 0.2,
                                               QUAD9);
        
        
Print information about the mesh to the screen.
          mesh.print_info();
        
        
Create an equation systems object.
          EquationSystems equation_systems (mesh);
        
Declare the system and its variables. Create a system named "Elasticity"
          LinearImplicitSystem& system =
            equation_systems.add_system<LinearImplicitSystem> ("Elasticity");
        
        
Add two displacement variables, u and v, to the system
          unsigned int u_var = system.add_variable("u", SECOND, LAGRANGE);
          unsigned int v_var = system.add_variable("v", SECOND, LAGRANGE);
        
        
          system.attach_assemble_function (assemble_elasticity);
        
Construct a Dirichlet boundary condition object We impose a "clamped" boundary condition on the "left" boundary, i.e. bc_id = 3
          std::set<boundary_id_type> boundary_ids;
          boundary_ids.insert(3);
        
Create a vector storing the variable numbers which the BC applies to
          std::vector<unsigned int> variables(2);
          variables[0] = u_var; variables[1] = v_var;
        
Create a ZeroFunction to initialize dirichlet_bc
          ZeroFunction<> zf;
        
          DirichletBoundary dirichlet_bc(boundary_ids,
                                         variables,
                                         &zf);
        
We must add the Dirichlet boundary condition _before_ we call equation_systems.init()
          system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
        
Initialize the data structures for the equation system.
          equation_systems.init();
        
Print information about the system to the screen.
          equation_systems.print_info();
        
Solve the system
          system.solve();
        
Plot the solution
        #ifdef LIBMESH_HAVE_EXODUS_API
          ExodusII_IO (mesh).write_equation_systems("displacement.e",equation_systems);
        #endif // #ifdef LIBMESH_HAVE_EXODUS_API
        
All done.
          return 0;
        }
        
        
        void assemble_elasticity(EquationSystems& es,
                                 const std::string& system_name)
        {
          libmesh_assert_equal_to (system_name, "Elasticity");
        
          const MeshBase& mesh = es.get_mesh();
        
          const unsigned int dim = mesh.mesh_dimension();
        
          LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Elasticity");
        
          const unsigned int u_var = system.variable_number ("u");
          const unsigned int v_var = system.variable_number ("v");
        
          const DofMap& dof_map = system.get_dof_map();
          FEType fe_type = dof_map.variable_type(0);
          AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
          QGauss qrule (dim, fe_type.default_quadrature_order());
          fe->attach_quadrature_rule (&qrule);
        
          AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
          QGauss qface(dim-1, fe_type.default_quadrature_order());
          fe_face->attach_quadrature_rule (&qface);
        
          const std::vector<Real>& JxW = fe->get_JxW();
          const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
        
          DenseMatrix<Number> Ke;
          DenseVector<Number> Fe;
        
          DenseSubMatrix<Number>
            Kuu(Ke), Kuv(Ke),
            Kvu(Ke), Kvv(Ke);
        
          DenseSubVector<Number>
            Fu(Fe),
            Fv(Fe);
        
          std::vector<dof_id_type> dof_indices;
          std::vector<dof_id_type> dof_indices_u;
          std::vector<dof_id_type> dof_indices_v;
        
          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)
            {
              const Elem* elem = *el;
        
              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();
        
              fe->reinit (elem);
        
              Ke.resize (n_dofs, n_dofs);
              Fe.resize (n_dofs);
        
              Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
              Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
        
              Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
              Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
        
              Fu.reposition (u_var*n_u_dofs, n_u_dofs);
              Fv.reposition (v_var*n_u_dofs, n_v_dofs);
        
              for (unsigned int qp=0; qp<qrule.n_points(); qp++)
              {
                  for (unsigned int i=0; i<n_u_dofs; i++)
                    for (unsigned int j=0; j<n_u_dofs; j++)
                    {
Tensor indices
                      unsigned int C_i, C_j, C_k, C_l;
                      C_i=0, C_k=0;
        
        
                      C_j=0, C_l=0;
                      Kuu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
        
                      C_j=1, C_l=0;
                      Kuu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
        
                      C_j=0, C_l=1;
                      Kuu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
        
                      C_j=1, C_l=1;
                      Kuu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                    }
        
                  for (unsigned int i=0; i<n_u_dofs; i++)
                    for (unsigned int j=0; j<n_v_dofs; j++)
                    {
Tensor indices
                      unsigned int C_i, C_j, C_k, C_l;
                      C_i=0, C_k=1;
        
        
                      C_j=0, C_l=0;
                      Kuv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
        
                      C_j=1, C_l=0;
                      Kuv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
        
                      C_j=0, C_l=1;
                      Kuv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
        
                      C_j=1, C_l=1;
                      Kuv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                    }
        
                  for (unsigned int i=0; i<n_v_dofs; i++)
                    for (unsigned int j=0; j<n_u_dofs; j++)
                    {
Tensor indices
                      unsigned int C_i, C_j, C_k, C_l;
                      C_i=1, C_k=0;
        
        
                      C_j=0, C_l=0;
                      Kvu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
        
                      C_j=1, C_l=0;
                      Kvu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
        
                      C_j=0, C_l=1;
                      Kvu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
        
                      C_j=1, C_l=1;
                      Kvu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                    }
        
                  for (unsigned int i=0; i<n_v_dofs; i++)
                    for (unsigned int j=0; j<n_v_dofs; j++)
                    {
Tensor indices
                      unsigned int C_i, C_j, C_k, C_l;
                      C_i=1, C_k=1;
        
        
                      C_j=0, C_l=0;
                      Kvv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
        
                      C_j=1, C_l=0;
                      Kvv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
        
                      C_j=0, C_l=1;
                      Kvv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
        
                      C_j=1, C_l=1;
                      Kvv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
                    }
              }
        
              {
                for (unsigned int side=0; side<elem->n_sides(); side++)
                  if (elem->neighbor(side) == NULL)
                    {
                      const std::vector<std::vector<Real> >&  phi_face = fe_face->get_phi();
                      const std::vector<Real>& JxW_face = fe_face->get_JxW();
        
                      fe_face->reinit(elem, side);
        
                      if( mesh.boundary_info->has_boundary_id (elem, side, 1) ) // Apply a traction on the right side
                      {
                        for (unsigned int qp=0; qp<qface.n_points(); qp++)
                        {
                          for (unsigned int i=0; i<n_v_dofs; i++)
                          {
                            Fv(i) += JxW_face[qp]* (-1.) * phi_face[i][qp];
                          }
                        }
                      }
                    }
              }
        
              dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
        
              system.matrix->add_matrix (Ke, dof_indices);
              system.rhs->add_vector    (Fe, dof_indices);
            }
        }
        
        Real eval_elasticity_tensor(unsigned int i,
                                    unsigned int j,
                                    unsigned int k,
                                    unsigned int l)
        {
Define the Poisson ratio
          const Real nu = 0.3;
        
Define the Lame constants (lambda_1 and lambda_2) based on Poisson ratio
          const Real lambda_1 = nu / ( (1. + nu) * (1. - 2.*nu) );
          const Real lambda_2 = 0.5 / (1 + nu);
        
Define the Kronecker delta functions that we need here
          Real delta_ij = (i == j) ? 1. : 0.;
          Real delta_il = (i == l) ? 1. : 0.;
          Real delta_ik = (i == k) ? 1. : 0.;
          Real delta_jl = (j == l) ? 1. : 0.;
          Real delta_jk = (j == k) ? 1. : 0.;
          Real delta_kl = (k == l) ? 1. : 0.;
        
          return lambda_1 * delta_ij * delta_kl + lambda_2 * (delta_ik * delta_jl + delta_il * delta_jk);
        }



The source file systems_of_equations_ex4.C without comments:

 
  
  
  #include <iostream>
  #include <algorithm>
  #include <math.h>
  
  #include "libmesh/libmesh.h"
  #include "libmesh/mesh.h"
  #include "libmesh/mesh_generation.h"
  #include "libmesh/exodusII_io.h"
  #include "libmesh/gnuplot_io.h"
  #include "libmesh/linear_implicit_system.h"
  #include "libmesh/equation_systems.h"
  #include "libmesh/fe.h"
  #include "libmesh/quadrature_gauss.h"
  #include "libmesh/dof_map.h"
  #include "libmesh/sparse_matrix.h"
  #include "libmesh/numeric_vector.h"
  #include "libmesh/dense_matrix.h"
  #include "libmesh/dense_submatrix.h"
  #include "libmesh/dense_vector.h"
  #include "libmesh/dense_subvector.h"
  #include "libmesh/perf_log.h"
  #include "libmesh/elem.h"
  #include "libmesh/boundary_info.h"
  #include "libmesh/zero_function.h"
  #include "libmesh/dirichlet_boundaries.h"
  #include "libmesh/string_to_enum.h"
  #include "libmesh/getpot.h"
  
  using namespace libMesh;
  
  void assemble_elasticity(EquationSystems& es,
                           const std::string& system_name);
  
  Real eval_elasticity_tensor(unsigned int i,
                              unsigned int j,
                              unsigned int k,
                              unsigned int l);
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
  
    const unsigned int dim = 2;
  
    libmesh_example_assert(dim <= LIBMESH_DIM, "2D support");
  
    Mesh mesh(init.comm(), dim);
    MeshTools::Generation::build_square (mesh,
                                         50, 10,
                                         0., 1.,
                                         0., 0.2,
                                         QUAD9);
  
  
    mesh.print_info();
  
  
    EquationSystems equation_systems (mesh);
  
    LinearImplicitSystem& system =
      equation_systems.add_system<LinearImplicitSystem> ("Elasticity");
  
  
    unsigned int u_var = system.add_variable("u", SECOND, LAGRANGE);
    unsigned int v_var = system.add_variable("v", SECOND, LAGRANGE);
  
  
    system.attach_assemble_function (assemble_elasticity);
  
    std::set<boundary_id_type> boundary_ids;
    boundary_ids.insert(3);
  
    std::vector<unsigned int> variables(2);
    variables[0] = u_var; variables[1] = v_var;
  
    ZeroFunction<> zf;
  
    DirichletBoundary dirichlet_bc(boundary_ids,
                                   variables,
                                   &zf);
  
    system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
  
    equation_systems.init();
  
    equation_systems.print_info();
  
    system.solve();
  
  #ifdef LIBMESH_HAVE_EXODUS_API
    ExodusII_IO (mesh).write_equation_systems("displacement.e",equation_systems);
  #endif // #ifdef LIBMESH_HAVE_EXODUS_API
  
    return 0;
  }
  
  
  void assemble_elasticity(EquationSystems& es,
                           const std::string& system_name)
  {
    libmesh_assert_equal_to (system_name, "Elasticity");
  
    const MeshBase& mesh = es.get_mesh();
  
    const unsigned int dim = mesh.mesh_dimension();
  
    LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Elasticity");
  
    const unsigned int u_var = system.variable_number ("u");
    const unsigned int v_var = system.variable_number ("v");
  
    const DofMap& dof_map = system.get_dof_map();
    FEType fe_type = dof_map.variable_type(0);
    AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
    QGauss qrule (dim, fe_type.default_quadrature_order());
    fe->attach_quadrature_rule (&qrule);
  
    AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
    QGauss qface(dim-1, fe_type.default_quadrature_order());
    fe_face->attach_quadrature_rule (&qface);
  
    const std::vector<Real>& JxW = fe->get_JxW();
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
    DenseMatrix<Number> Ke;
    DenseVector<Number> Fe;
  
    DenseSubMatrix<Number>
      Kuu(Ke), Kuv(Ke),
      Kvu(Ke), Kvv(Ke);
  
    DenseSubVector<Number>
      Fu(Fe),
      Fv(Fe);
  
    std::vector<dof_id_type> dof_indices;
    std::vector<dof_id_type> dof_indices_u;
    std::vector<dof_id_type> dof_indices_v;
  
    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)
      {
        const Elem* elem = *el;
  
        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();
  
        fe->reinit (elem);
  
        Ke.resize (n_dofs, n_dofs);
        Fe.resize (n_dofs);
  
        Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
        Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
  
        Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
        Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
  
        Fu.reposition (u_var*n_u_dofs, n_u_dofs);
        Fv.reposition (v_var*n_u_dofs, n_v_dofs);
  
        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        {
            for (unsigned int i=0; i<n_u_dofs; i++)
              for (unsigned int j=0; j<n_u_dofs; j++)
              {
                unsigned int C_i, C_j, C_k, C_l;
                C_i=0, C_k=0;
  
  
                C_j=0, C_l=0;
                Kuu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
  
                C_j=1, C_l=0;
                Kuu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
  
                C_j=0, C_l=1;
                Kuu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
  
                C_j=1, C_l=1;
                Kuu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
              }
  
            for (unsigned int i=0; i<n_u_dofs; i++)
              for (unsigned int j=0; j<n_v_dofs; j++)
              {
                unsigned int C_i, C_j, C_k, C_l;
                C_i=0, C_k=1;
  
  
                C_j=0, C_l=0;
                Kuv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
  
                C_j=1, C_l=0;
                Kuv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
  
                C_j=0, C_l=1;
                Kuv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
  
                C_j=1, C_l=1;
                Kuv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
              }
  
            for (unsigned int i=0; i<n_v_dofs; i++)
              for (unsigned int j=0; j<n_u_dofs; j++)
              {
                unsigned int C_i, C_j, C_k, C_l;
                C_i=1, C_k=0;
  
  
                C_j=0, C_l=0;
                Kvu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
  
                C_j=1, C_l=0;
                Kvu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
  
                C_j=0, C_l=1;
                Kvu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
  
                C_j=1, C_l=1;
                Kvu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
              }
  
            for (unsigned int i=0; i<n_v_dofs; i++)
              for (unsigned int j=0; j<n_v_dofs; j++)
              {
                unsigned int C_i, C_j, C_k, C_l;
                C_i=1, C_k=1;
  
  
                C_j=0, C_l=0;
                Kvv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
  
                C_j=1, C_l=0;
                Kvv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
  
                C_j=0, C_l=1;
                Kvv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
  
                C_j=1, C_l=1;
                Kvv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
              }
        }
  
        {
          for (unsigned int side=0; side<elem->n_sides(); side++)
            if (elem->neighbor(side) == NULL)
              {
                const std::vector<std::vector<Real> >&  phi_face = fe_face->get_phi();
                const std::vector<Real>& JxW_face = fe_face->get_JxW();
  
                fe_face->reinit(elem, side);
  
                if( mesh.boundary_info->has_boundary_id (elem, side, 1) ) // Apply a traction on the right side
                {
                  for (unsigned int qp=0; qp<qface.n_points(); qp++)
                  {
                    for (unsigned int i=0; i<n_v_dofs; i++)
                    {
                      Fv(i) += JxW_face[qp]* (-1.) * phi_face[i][qp];
                    }
                  }
                }
              }
        }
  
        dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
  
        system.matrix->add_matrix (Ke, dof_indices);
        system.rhs->add_vector    (Fe, dof_indices);
      }
  }
  
  Real eval_elasticity_tensor(unsigned int i,
                              unsigned int j,
                              unsigned int k,
                              unsigned int l)
  {
    const Real nu = 0.3;
  
    const Real lambda_1 = nu / ( (1. + nu) * (1. - 2.*nu) );
    const Real lambda_2 = 0.5 / (1 + nu);
  
    Real delta_ij = (i == j) ? 1. : 0.;
    Real delta_il = (i == l) ? 1. : 0.;
    Real delta_ik = (i == k) ? 1. : 0.;
    Real delta_jl = (j == l) ? 1. : 0.;
    Real delta_jk = (j == k) ? 1. : 0.;
    Real delta_kl = (k == l) ? 1. : 0.;
  
    return lambda_1 * delta_ij * delta_kl + lambda_2 * (delta_ik * delta_jl + delta_il * delta_jk);
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/systems_of_equations/systems_of_equations_ex4'
***************************************************************
* Running Example systems_of_equations_ex4:
*  mpirun -np 4 example-devel -ksp_type cg -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
 
 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=2121
    n_local_nodes()=547
  n_elem()=500
    n_local_elem()=125
    n_active_elem()=500
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "Elasticity"
    Type "LinearImplicit"
    Variables={ "u" "v" } 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="SECOND", "THIRD" 
    n_dofs()=4242
    n_local_dofs()=1094
    n_constrained_dofs()=42
    n_local_constrained_dofs()=0
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 29.9255
      Average Off-Processor Bandwidth <= 0.946723
      Maximum  On-Processor Bandwidth <= 52
      Maximum Off-Processor Bandwidth <= 28
    DofMap Constraints
      Number of DoF Constraints = 42
      Average DoF Constraint Length= 0
      Number of Node Constraints = 0


 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:54:06 2013                                                                          |
| OS:             Linux                                                                                             |
| HostName:       spark.ices.utexas.edu                                                                             |
| OS Release:     2.6.32-279.22.1.el6.x86_64                                                                        |
| OS Version:     #1 SMP Tue Feb 5 14:33:39 CST 2013                                                                |
| Machine:        x86_64                                                                                            |
| Username:       roystgnr                                                                                          |
| Configuration:  ../configure  '--enable-everything'                                                               |
|  'METHODS=devel'                                                                                                  |
|  '--prefix=/h2/roystgnr/libmesh-test'                                                                             |
|  'CXX=distcc /usr/bin/g++'                                                                                        |
|  'CC=distcc /usr/bin/gcc'                                                                                         |
|  'FC=distcc /usr/bin/gfortran'                                                                                    |
|  'F77=distcc /usr/bin/gfortran'                                                                                   |
|  'PETSC_DIR=/opt/apps/ossw/libraries/petsc/petsc-3.3-p2'                                                          |
|  'PETSC_ARCH=gcc-system-mkl-gf-10.3.12.361-mpich2-1.4.1p1-cxx-opt'                                                |
|  'SLEPC_DIR=/opt/apps/ossw/libraries/slepc/slepc-3.3-p2-petsc-3.3-p2-cxx-opt'                                     |
|  'TRILINOS_DIR=/opt/apps/ossw/libraries/trilinos/trilinos-10.12.2/sl6/gcc-system/mpich2-1.4.1p1/mkl-gf-10.3.12.361'|
|  'VTK_DIR=/opt/apps/ossw/libraries/vtk/vtk-5.10.0/sl6/gcc-system'                                                 |
|  'HDF5_DIR=/opt/apps/ossw/libraries/hdf5/hdf5-1.8.9/sl6/gcc-system'                                               |
 -------------------------------------------------------------------------------------------------------------------
 ----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=0.397633, Active time=0.294797                                                 |
 ----------------------------------------------------------------------------------------------------------------
| Event                              nCalls    Total Time  Avg Time    Total Time  Avg Time    % of Active Time  |
|                                              w/o Sub     w/o Sub     With Sub    With Sub    w/o S    With S   |
|----------------------------------------------------------------------------------------------------------------|
|                                                                                                                |
|                                                                                                                |
| DofMap                                                                                                         |
|   add_neighbors_to_send_list()     1         0.0013      0.001261    0.0015      0.001523    0.43     0.52     |
|   build_constraint_matrix()        125       0.0001      0.000001    0.0001      0.000001    0.05     0.05     |
|   build_sparsity()                 1         0.0012      0.001171    0.0048      0.004793    0.40     1.63     |
|   cnstrn_elem_mat_vec()            125       0.0001      0.000001    0.0001      0.000001    0.03     0.03     |
|   create_dof_constraints()         1         0.0020      0.001966    0.0090      0.009022    0.67     3.06     |
|   distribute_dofs()                1         0.0023      0.002314    0.0084      0.008429    0.78     2.86     |
|   dof_indices()                    1772      0.0138      0.000008    0.0138      0.000008    4.68     4.68     |
|   prepare_send_list()              1         0.0000      0.000011    0.0000      0.000011    0.00     0.00     |
|   reinit()                         1         0.0028      0.002783    0.0028      0.002783    0.94     0.94     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          1         0.0008      0.000797    0.0040      0.004002    0.27     1.36     |
|                                                                                                                |
| ExodusII_IO                                                                                                    |
|   write_nodal_data()               1         0.0665      0.066522    0.0665      0.066522    22.57    22.57    |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        180       0.0005      0.000003    0.0005      0.000003    0.18     0.18     |
|   init_shape_functions()           56        0.0001      0.000002    0.0001      0.000002    0.03     0.03     |
|   inverse_map()                    165       0.0006      0.000004    0.0006      0.000004    0.21     0.21     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             180       0.0006      0.000003    0.0006      0.000003    0.19     0.19     |
|   compute_face_map()               55        0.0004      0.000007    0.0010      0.000019    0.14     0.35     |
|   init_face_shape_functions()      21        0.0000      0.000002    0.0000      0.000002    0.02     0.02     |
|   init_reference_to_physical_map() 56        0.0004      0.000008    0.0004      0.000008    0.15     0.15     |
|                                                                                                                |
| Mesh                                                                                                           |
|   find_neighbors()                 1         0.0008      0.000828    0.0015      0.001540    0.28     0.52     |
|   renumber_nodes_and_elem()        2         0.0002      0.000081    0.0002      0.000081    0.05     0.05     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   compute_hilbert_indices()        2         0.0026      0.001299    0.0026      0.001299    0.88     0.88     |
|   find_global_indices()            2         0.0004      0.000190    0.0056      0.002789    0.13     1.89     |
|   parallel_sort()                  2         0.0003      0.000166    0.0022      0.001115    0.11     0.76     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         1         0.0001      0.000067    0.0707      0.070662    0.02     23.97    |
|                                                                                                                |
| MeshTools::Generation                                                                                          |
|   build_cube()                     1         0.0006      0.000628    0.0006      0.000628    0.21     0.21     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      1         0.0031      0.003088    0.0061      0.006114    1.05     2.07     |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      9         0.0024      0.000262    0.0024      0.000269    0.80     0.82     |
|   max(bool)                        1         0.0000      0.000004    0.0000      0.000004    0.00     0.00     |
|   max(scalar)                      105       0.0005      0.000005    0.0005      0.000005    0.18     0.18     |
|   max(vector)                      24        0.0001      0.000006    0.0004      0.000018    0.05     0.15     |
|   min(bool)                        121       0.0004      0.000004    0.0004      0.000004    0.15     0.15     |
|   min(scalar)                      99        0.0138      0.000140    0.0138      0.000140    4.70     4.70     |
|   min(vector)                      24        0.0002      0.000008    0.0009      0.000039    0.07     0.32     |
|   probe()                          36        0.0009      0.000025    0.0009      0.000025    0.30     0.30     |
|   receive()                        36        0.0001      0.000003    0.0010      0.000028    0.04     0.34     |
|   send()                           36        0.0001      0.000003    0.0001      0.000003    0.03     0.03     |
|   send_receive()                   40        0.0002      0.000005    0.0013      0.000034    0.07     0.46     |
|   sum()                            20        0.0021      0.000107    0.0024      0.000120    0.72     0.81     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           36        0.0000      0.000001    0.0000      0.000001    0.02     0.02     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         1         0.0005      0.000506    0.0023      0.002324    0.17     0.79     |
|   set_parent_processor_ids()       1         0.0001      0.000078    0.0001      0.000078    0.03     0.03     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          1         0.1384      0.138356    0.1384      0.138356    46.93    46.93    |
|                                                                                                                |
| System                                                                                                         |
|   assemble()                       1         0.0332      0.033209    0.0388      0.038782    11.27    13.16    |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            3346      0.2948                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example systems_of_equations_ex4:
*  mpirun -np 4 example-devel -ksp_type cg -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
make[4]: Leaving directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/systems_of_equations/systems_of_equations_ex4'

Site Created By: libMesh Developers
Last modified: April 23 2013 04:19:30 UTC

Hosted By:
SourceForge.net Logo