The source file systems_of_equations_ex4.C with comments:
#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");
Mesh mesh(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(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:
***************************************************************
* Running Example systems_of_equations_ex4:
* mpirun -np 2 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 -log_summary
***************************************************************
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=2121
n_local_nodes()=1071
n_elem()=500
n_local_elem()=250
n_active_elem()=500
n_subdomains()=1
n_partitions()=2
n_processors()=2
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()=2142
n_constrained_dofs()=42
n_local_constrained_dofs()=0
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 30.3989
Average Off-Processor Bandwidth <= 0.305516
Maximum On-Processor Bandwidth <= 50
Maximum Off-Processor Bandwidth <= 20
DofMap Constraints
Number of DoF Constraints = 42
Average DoF Constraint Length= 0
Number of Node Constraints = 0
************************************************************************************************************************
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***
************************************************************************************************************************
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------
/workspace/libmesh/examples/systems_of_equations/systems_of_equations_ex4/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:44:43 2013
Using Petsc Release Version 3.3.0, Patch 2, Fri Jul 13 15:42:00 CDT 2012
Max Max/Min Avg Total
Time (sec): 8.923e-01 1.00000 8.923e-01
Objects: 3.300e+01 1.00000 3.300e+01
Flops: 7.911e+07 1.02910 7.799e+07 1.560e+08
Flops/sec: 8.865e+07 1.02910 8.740e+07 1.748e+08
MPI Messages: 7.150e+01 1.00000 7.150e+01 1.430e+02
MPI Message Lengths: 4.278e+04 1.00000 5.983e+02 8.556e+04
MPI Reductions: 2.380e+02 1.00000
Flop counting convention: 1 flop = 1 real number operation of type (multiply/divide/add/subtract)
e.g., VecAXPY() for real vectors of length N --> 2N flops
and VecAXPY() for complex vectors of length N --> 8N flops
Summary of Stages: ----- Time ------ ----- Flops ----- --- Messages --- -- Message Lengths -- -- Reductions --
Avg %Total Avg %Total counts %Total Avg %Total counts %Total
0: Main Stage: 8.9230e-01 100.0% 1.5598e+08 100.0% 1.430e+02 100.0% 5.983e+02 100.0% 2.370e+02 99.6%
------------------------------------------------------------------------------------------------------------------------
See the 'Profiling' chapter of the users' manual for details on interpreting output.
Phase summary info:
Count: number of times phase was executed
Time and Flops: Max - maximum over all processors
Ratio - ratio of maximum to minimum over all processors
Mess: number of messages sent
Avg. len: average message length
Reduct: number of global reductions
Global: entire computation
Stage: stages of a computation. Set stages with PetscLogStagePush() and PetscLogStagePop().
%T - percent time in this phase %f - percent flops in this phase
%M - percent messages in this phase %L - percent message lengths in this phase
%R - percent reductions in this phase
Total Mflop/s: 10e-6 * (sum of flops over all processors)/(max time over all processors)
------------------------------------------------------------------------------------------------------------------------
Event Count Time (sec) Flops --- Global --- --- Stage --- Total
Max Ratio Max Ratio Max Ratio Mess Avg len Reduct %T %f %M %L %R %T %f %M %L %R Mflop/s
------------------------------------------------------------------------------------------------------------------------
--- Event Stage 0: Main Stage
VecTDot 126 1.0 5.8198e-04 1.0 5.40e+05 1.0 0.0e+00 0.0e+00 1.3e+02 0 1 0 0 53 0 1 0 0 53 1836
VecNorm 65 1.0 9.0971e-03 1.0 2.78e+05 1.0 0.0e+00 0.0e+00 6.5e+01 1 0 0 0 27 1 0 0 0 27 61
VecCopy 2 1.0 8.1062e-06 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecSet 70 1.0 8.7738e-05 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecAXPY 126 1.0 3.9697e-04 1.1 5.40e+05 1.0 0.0e+00 0.0e+00 0.0e+00 0 1 0 0 0 0 1 0 0 0 2693
VecAYPX 63 1.0 2.5129e-04 1.0 2.68e+05 1.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 2110
VecAssemblyBegin 3 1.0 7.8201e-05 1.0 0.00e+00 0.0 2.0e+00 3.6e+02 9.0e+00 0 0 1 1 4 0 0 1 1 4 0
VecAssemblyEnd 3 1.0 2.4796e-05 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecScatterBegin 65 1.0 2.0361e-04 1.2 0.00e+00 0.0 1.3e+02 5.1e+02 0.0e+00 0 0 91 77 0 0 0 91 77 0 0
VecScatterEnd 65 1.0 1.3089e-0314.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
MatMult 64 1.0 7.5910e-03 1.2 8.28e+06 1.0 1.3e+02 5.0e+02 0.0e+00 1 10 90 75 0 1 10 90 75 0 2155
MatSolve 65 1.0 2.6434e-02 1.0 4.25e+07 1.0 0.0e+00 0.0e+00 0.0e+00 3 54 0 0 0 3 54 0 0 0 3181
MatLUFactorNum 1 1.0 2.1129e-02 1.0 2.67e+07 1.0 0.0e+00 0.0e+00 0.0e+00 2 34 0 0 0 2 34 0 0 0 2476
MatILUFactorSym 1 1.0 5.8221e-02 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 3.0e+00 6 0 0 0 1 6 0 0 0 1 0
MatAssemblyBegin 2 1.0 1.2331e-0264.2 0.00e+00 0.0 3.0e+00 5.8e+03 4.0e+00 1 0 2 20 2 1 0 2 20 2 0
MatAssemblyEnd 2 1.0 6.5923e-04 1.0 0.00e+00 0.0 4.0e+00 1.3e+02 8.0e+00 0 0 3 1 3 0 0 3 1 3 0
MatGetRowIJ 1 1.0 5.9605e-06 1.5 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
MatGetOrdering 1 1.0 1.8287e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 4.0e+00 0 0 0 0 2 0 0 0 0 2 0
MatZeroEntries 3 1.0 1.9717e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
KSPSetUp 2 1.0 1.2994e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
KSPSolve 1 1.0 1.2408e-01 1.0 7.91e+07 1.0 1.3e+02 5.0e+02 2.0e+02 14100 90 75 84 14100 90 75 84 1257
PCSetUp 2 1.0 8.0088e-02 1.0 2.67e+07 1.0 0.0e+00 0.0e+00 9.0e+00 9 34 0 0 4 9 34 0 0 4 653
PCSetUpOnBlocks 1 1.0 7.9685e-02 1.0 2.67e+07 1.0 0.0e+00 0.0e+00 7.0e+00 9 34 0 0 3 9 34 0 0 3 656
PCApply 65 1.0 2.7169e-02 1.0 4.25e+07 1.0 0.0e+00 0.0e+00 0.0e+00 3 54 0 0 0 3 54 0 0 0 3095
------------------------------------------------------------------------------------------------------------------------
Memory usage is given in bytes:
Object Type Creations Destructions Memory Descendants' Mem.
Reports information only for process 0.
--- Event Stage 0: Main Stage
Vector 12 12 138992 0
Vector Scatter 2 2 2072 0
Index Set 9 9 28536 0
IS L to G Mapping 1 1 564 0
Matrix 4 4 4811284 0
Krylov Solver 2 2 2368 0
Preconditioner 2 2 1784 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 1.90735e-07
Average time for MPI_Barrier(): 3.38554e-06
Average time for zero size MPI_Send(): 1.69277e-05
#PETSc Option Table entries:
-ksp_right_pc
-ksp_type cg
-log_summary
-pc_type bjacobi
-sub_pc_factor_levels 4
-sub_pc_factor_zeropivot 0
-sub_pc_type ilu
#End of PETSc Option Table entries
Compiled without FORTRAN kernels
Compiled with full precision matrices (default)
sizeof(short) 2 sizeof(int) 4 sizeof(long) 8 sizeof(void*) 8 sizeof(PetscScalar) 8 sizeof(PetscInt) 4
Configure run at: Thu Nov 8 11:21:02 2012
Configure options: --with-debugging=false --COPTFLAGS=-O3 --CXXOPTFLAGS=-O3 --FOPTFLAGS=-O3 --with-clanguage=C++ --with-shared-libraries=1 --with-mpi-dir=/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1 --with-mumps=true --download-mumps=1 --with-metis=true --download-metis=1 --with-parmetis=true --download-parmetis=1 --with-superlu=true --download-superlu=1 --with-superludir=true --download-superlu_dist=1 --with-blacs=true --download-blacs=1 --with-scalapack=true --download-scalapack=1 --with-hypre=true --download-hypre=1 --with-blas-lib="[/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_intel_lp64.so,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_sequential.so,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_core.so]" --with-lapack-lib="[/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_lapack95_lp64.a]"
-----------------------------------------
Libraries compiled on Thu Nov 8 11:21:02 2012 on daedalus.ices.utexas.edu
Machine characteristics: Linux-2.6.32-279.1.1.el6.x86_64-x86_64-with-redhat-6.3-Carbon
Using PETSc directory: /opt/apps/ossw/libraries/petsc/petsc-3.3-p2
Using PETSc arch: intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt
-----------------------------------------
Using C compiler: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpicxx -wd1572 -O3 -fPIC ${COPTFLAGS} ${CFLAGS}
Using Fortran compiler: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpif90 -fPIC -O3 ${FOPTFLAGS} ${FFLAGS}
-----------------------------------------
Using include paths: -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/include -I/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/include
-----------------------------------------
Using C linker: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpicxx
Using Fortran linker: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpif90
Using libraries: -Wl,-rpath,/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -L/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -lpetsc -lX11 -Wl,-rpath,/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -L/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -lcmumps -ldmumps -lsmumps -lzmumps -lmumps_common -lpord -lHYPRE -lpthread -lsuperlu_dist_3.0 -lparmetis -lmetis -lscalapack -lblacs -lsuperlu_4.3 -Wl,-rpath,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64 -L/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64 -lmkl_lapack95_lp64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -Wl,-rpath,/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/lib -L/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/lib -Wl,-rpath,/opt/apps/sysnet/intel/12.1/composer_xe_2011_sp1.7.256/compiler/lib/intel64 -L/opt/apps/sysnet/intel/12.1/composer_xe_2011_sp1.7.256/compiler/lib/intel64 -Wl,-rpath,/usr/lib/gcc/x86_64-redhat-linux/4.4.6 -L/usr/lib/gcc/x86_64-redhat-linux/4.4.6 -lmpichf90 -lifport -lifcore -lm -lm -lmpichcxx -ldl -lmpich -lopa -lmpl -lrt -lpthread -limf -lsvml -lipgo -ldecimal -lcilkrts -lstdc++ -lgcc_s -lirc -lirc_s -ldl
-----------------------------------------
----------------------------------------------------------------------------------------------------------------------
| Processor id: 0 |
| Num Processors: 2 |
| Time: Tue Feb 5 13:44:43 2013 |
| OS: Linux |
| HostName: hbar.ices.utexas.edu |
| OS Release: 2.6.32-279.1.1.el6.x86_64 |
| OS Version: #1 SMP Tue Jul 10 11:24:23 CDT 2012 |
| Machine: x86_64 |
| Username: benkirk |
| Configuration: ./configure '--enable-everything' |
| '--prefix=/workspace/libmesh/install' |
| 'CXX=mpicxx' |
| 'CC=mpicc' |
| 'F77=mpif77' |
| 'FC=mpif90' |
| 'PETSC_DIR=/opt/apps/ossw/libraries/petsc/petsc-3.3-p2' |
| 'PETSC_ARCH=intel-12.1-mkl-intel-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/intel-12.1/mpich2-1.4.1p1/mkl-intel-10.3.12.361'|
| 'VTK_DIR=/opt/apps/ossw/libraries/vtk/vtk-5.10.0/sl6/intel-12.1' |
----------------------------------------------------------------------------------------------------------------------
----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=0.905247, Active time=0.877297 |
----------------------------------------------------------------------------------------------------------------
| 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.0220 0.022022 0.0243 0.024255 2.51 2.76 |
| build_constraint_matrix() 250 0.0020 0.000008 0.0020 0.000008 0.23 0.23 |
| build_sparsity() 1 0.0253 0.025256 0.0543 0.054286 2.88 6.19 |
| cnstrn_elem_mat_vec() 250 0.0004 0.000002 0.0004 0.000002 0.05 0.05 |
| create_dof_constraints() 1 0.0095 0.009541 0.0685 0.068484 1.09 7.81 |
| distribute_dofs() 1 0.0144 0.014438 0.0346 0.034581 1.65 3.94 |
| dof_indices() 2520 0.2372 0.000094 0.2372 0.000094 27.04 27.04 |
| prepare_send_list() 1 0.0001 0.000064 0.0001 0.000064 0.01 0.01 |
| reinit() 1 0.0198 0.019768 0.0198 0.019768 2.25 2.25 |
| |
| EquationSystems |
| build_solution_vector() 1 0.0037 0.003691 0.0548 0.054788 0.42 6.25 |
| |
| ExodusII_IO |
| write_nodal_data() 1 0.0066 0.006577 0.0066 0.006577 0.75 0.75 |
| |
| FE |
| compute_shape_functions() 330 0.0059 0.000018 0.0059 0.000018 0.68 0.68 |
| init_shape_functions() 81 0.0004 0.000005 0.0004 0.000005 0.04 0.04 |
| inverse_map() 240 0.0034 0.000014 0.0034 0.000014 0.39 0.39 |
| |
| FEMap |
| compute_affine_map() 330 0.0049 0.000015 0.0049 0.000015 0.56 0.56 |
| compute_face_map() 80 0.0022 0.000028 0.0057 0.000071 0.26 0.65 |
| init_face_shape_functions() 21 0.0001 0.000006 0.0001 0.000006 0.02 0.02 |
| init_reference_to_physical_map() 81 0.0024 0.000029 0.0024 0.000029 0.27 0.27 |
| |
| Mesh |
| find_neighbors() 1 0.0056 0.005628 0.0057 0.005743 0.64 0.65 |
| renumber_nodes_and_elem() 2 0.0006 0.000300 0.0006 0.000300 0.07 0.07 |
| |
| MeshCommunication |
| compute_hilbert_indices() 2 0.0044 0.002205 0.0044 0.002205 0.50 0.50 |
| find_global_indices() 2 0.0018 0.000899 0.0073 0.003670 0.20 0.84 |
| parallel_sort() 2 0.0008 0.000392 0.0009 0.000436 0.09 0.10 |
| |
| MeshOutput |
| write_equation_systems() 1 0.0001 0.000107 0.0615 0.061541 0.01 7.01 |
| |
| MeshTools::Generation |
| build_cube() 1 0.0029 0.002870 0.0029 0.002870 0.33 0.33 |
| |
| MetisPartitioner |
| partition() 1 0.0063 0.006261 0.0096 0.009602 0.71 1.09 |
| |
| Parallel |
| allgather() 9 0.0003 0.000030 0.0003 0.000032 0.03 0.03 |
| max(bool) 1 0.0000 0.000003 0.0000 0.000003 0.00 0.00 |
| max(scalar) 105 0.0003 0.000003 0.0003 0.000003 0.03 0.03 |
| max(vector) 24 0.0001 0.000006 0.0003 0.000014 0.02 0.04 |
| min(bool) 121 0.0003 0.000002 0.0003 0.000002 0.03 0.03 |
| min(scalar) 99 0.0025 0.000025 0.0025 0.000025 0.28 0.28 |
| min(vector) 24 0.0002 0.000008 0.0004 0.000017 0.02 0.05 |
| probe() 12 0.0002 0.000016 0.0002 0.000016 0.02 0.02 |
| receive() 12 0.0001 0.000007 0.0003 0.000024 0.01 0.03 |
| send() 12 0.0001 0.000005 0.0001 0.000005 0.01 0.01 |
| send_receive() 16 0.0002 0.000013 0.0006 0.000036 0.02 0.07 |
| sum() 20 0.0003 0.000013 0.0009 0.000045 0.03 0.10 |
| |
| Parallel::Request |
| wait() 12 0.0000 0.000003 0.0000 0.000003 0.00 0.00 |
| |
| Partitioner |
| set_node_processor_ids() 1 0.0011 0.001107 0.0012 0.001205 0.13 0.14 |
| set_parent_processor_ids() 1 0.0004 0.000385 0.0004 0.000385 0.04 0.04 |
| |
| PetscLinearSolver |
| solve() 1 0.1382 0.138174 0.1382 0.138174 15.75 15.75 |
| |
| System |
| assemble() 1 0.3503 0.350335 0.4706 0.470622 39.93 53.64 |
----------------------------------------------------------------------------------------------------------------
| Totals: 4674 0.8773 100.00 |
----------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example systems_of_equations_ex4:
* mpirun -np 2 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 -log_summary
***************************************************************
Systems Example 4 - Linear Elastic Cantilever
By David KnezevicIn 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