The source file systems_of_equations_ex6.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"
#define x_scaling 1.3
#define x_load 0.0
#define y_load 0.0
#define z_load -1.0
boundary IDs
#define BOUNDARY_ID_MIN_Z 0
#define BOUNDARY_ID_MIN_Y 1
#define BOUNDARY_ID_MAX_X 2
#define BOUNDARY_ID_MAX_Y 3
#define BOUNDARY_ID_MIN_X 4
#define BOUNDARY_ID_MAX_Z 5
Bring in everything from the libMesh namespace
using namespace libMesh;
Matrix and right-hand side assembly
void assemble_elasticity(EquationSystems& es,
const std::string& system_name);
Post-process the solution to compute stresses
void compute_stresses(EquationSystems& es);
The Kronecker delta function, used in eval_elasticity_tensor
Real kronecker_delta(unsigned int i,
unsigned int j);
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 = 3;
Make sure libMesh was compiled for 3D
libmesh_example_assert(dim == LIBMESH_DIM, "3D support");
Mesh mesh(dim);
MeshTools::Generation::build_cube (mesh,
40,
8,
4,
0., 1.*x_scaling,
0., 0.3,
0., 0.1,
HEX8);
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 three displacement variables, u and v, to the system
unsigned int u_var = system.add_variable("u", FIRST, LAGRANGE);
unsigned int v_var = system.add_variable("v", FIRST, LAGRANGE);
unsigned int w_var = system.add_variable("w", FIRST, LAGRANGE);
system.attach_assemble_function (assemble_elasticity);
std::set<boundary_id_type> boundary_ids;
boundary_ids.insert(BOUNDARY_ID_MIN_X);
Create a vector storing the variable numbers which the BC applies to
std::vector<unsigned int> variables;
variables.push_back(u_var);
variables.push_back(v_var);
variables.push_back(w_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);
Also, initialize an ExplicitSystem to store stresses
ExplicitSystem& stress_system =
equation_systems.add_system<ExplicitSystem> ("StressSystem");
stress_system.add_variable("sigma_00", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_01", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_02", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_10", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_11", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_12", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_20", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_21", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_22", CONSTANT, MONOMIAL);
stress_system.add_variable("vonMises", CONSTANT, MONOMIAL);
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();
Post-process the solution to compute the stresses
compute_stresses(equation_systems);
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 n_components = 3;
const unsigned int u_var = system.variable_number ("u");
const unsigned int v_var = system.variable_number ("v");
const unsigned int w_var = system.variable_number ("w");
const DofMap& dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(u_var);
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), Kuw(Ke),
Kvu(Ke), Kvv(Ke), Kvw(Ke),
Kwu(Ke), Kwv(Ke), Kww(Ke);
DenseSubVector<Number>
Fu(Fe),
Fv(Fe),
Fw(Fe);
std::vector<dof_id_type> dof_indices;
std::vector<dof_id_type> dof_indices_u;
std::vector<dof_id_type> dof_indices_v;
std::vector<dof_id_type> dof_indices_w;
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);
dof_map.dof_indices (elem, dof_indices_w, w_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();
const unsigned int n_w_dofs = dof_indices_w.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);
Kuw.reposition (u_var*n_u_dofs, w_var*n_u_dofs, n_u_dofs, n_w_dofs);
Kvu.reposition (v_var*n_u_dofs, u_var*n_u_dofs, n_v_dofs, n_u_dofs);
Kvv.reposition (v_var*n_u_dofs, v_var*n_u_dofs, n_v_dofs, n_v_dofs);
Kvw.reposition (v_var*n_u_dofs, w_var*n_u_dofs, n_v_dofs, n_w_dofs);
Kwu.reposition (w_var*n_u_dofs, u_var*n_u_dofs, n_w_dofs, n_u_dofs);
Kwv.reposition (w_var*n_u_dofs, v_var*n_u_dofs, n_w_dofs, n_v_dofs);
Kww.reposition (w_var*n_u_dofs, w_var*n_u_dofs, n_w_dofs, n_w_dofs);
Fu.reposition (u_var*n_u_dofs, n_u_dofs);
Fv.reposition (v_var*n_u_dofs, n_v_dofs);
Fw.reposition (w_var*n_u_dofs, n_w_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=0, C_k=0;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
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=0, C_k=1;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
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_u_dofs; i++)
for (unsigned int j=0; j<n_w_dofs; j++)
{
Tensor indices
unsigned int C_i=0, C_k=2;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
Kuw(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 = 1, C_k = 0;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
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 = 1, C_k = 1;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
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 i=0; i<n_v_dofs; i++)
for (unsigned int j=0; j<n_w_dofs; j++)
{
Tensor indices
unsigned int C_i = 1, C_k = 2;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
Kvw(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_w_dofs; i++)
for (unsigned int j=0; j<n_u_dofs; j++)
{
Tensor indices
unsigned int C_i = 2, C_k = 0;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
Kwu(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_w_dofs; i++)
for (unsigned int j=0; j<n_v_dofs; j++)
{
Tensor indices
unsigned int C_i = 2, C_k = 1;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
Kwv(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_w_dofs; i++)
for (unsigned int j=0; j<n_w_dofs; j++)
{
Tensor indices
unsigned int C_i = 2, C_k = 2;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
Kww(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);
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
Apply a traction
if( mesh.boundary_info->has_boundary_id(elem, side, BOUNDARY_ID_MAX_X) )
{
for (unsigned int i=0; i<n_v_dofs; i++)
{
Fu(i) += JxW_face[qp] * x_load * phi_face[i][qp];
}
for (unsigned int i=0; i<n_v_dofs; i++)
{
Fv(i) += JxW_face[qp] * y_load * phi_face[i][qp];
}
for (unsigned int i=0; i<n_v_dofs; i++)
{
Fw(i) += JxW_face[qp] * z_load * 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);
}
}
void compute_stresses(EquationSystems& es)
{
const MeshBase& mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Elasticity");
unsigned int displacement_vars[3];
displacement_vars[0] = system.variable_number ("u");
displacement_vars[1] = system.variable_number ("v");
displacement_vars[2] = system.variable_number ("w");
const unsigned int u_var = system.variable_number ("u");
const DofMap& dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(u_var);
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
QGauss qrule (dim, fe_type.default_quadrature_order());
fe->attach_quadrature_rule (&qrule);
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
Also, get a reference to the ExplicitSystem
ExplicitSystem& stress_system = es.get_system<ExplicitSystem>("StressSystem");
const DofMap& stress_dof_map = stress_system.get_dof_map();
unsigned int sigma_vars[3][3];
sigma_vars[0][0] = stress_system.variable_number ("sigma_00");
sigma_vars[0][1] = stress_system.variable_number ("sigma_01");
sigma_vars[0][2] = stress_system.variable_number ("sigma_02");
sigma_vars[1][0] = stress_system.variable_number ("sigma_10");
sigma_vars[1][1] = stress_system.variable_number ("sigma_11");
sigma_vars[1][2] = stress_system.variable_number ("sigma_12");
sigma_vars[2][0] = stress_system.variable_number ("sigma_20");
sigma_vars[2][1] = stress_system.variable_number ("sigma_21");
sigma_vars[2][2] = stress_system.variable_number ("sigma_22");
unsigned int vonMises_var = stress_system.variable_number ("vonMises");
Storage for the stress dof indices on each element
std::vector< std::vector<dof_id_type> > dof_indices_var(system.n_vars());
std::vector<dof_id_type> stress_dof_indices_var;
To store the stress tensor on each element
DenseMatrix<Number> elem_sigma;
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;
for(unsigned int var=0; var<3; var++)
{
dof_map.dof_indices (elem, dof_indices_var[var], displacement_vars[var]);
}
fe->reinit (elem);
elem_sigma.resize(3,3);
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
for(unsigned int C_i=0; C_i<3; C_i++)
for(unsigned int C_j=0; C_j<3; C_j++)
for(unsigned int C_k=0; C_k<3; C_k++)
{
const unsigned int n_var_dofs = dof_indices_var[C_k].size();
Get the gradient at this quadrature point
Gradient displacement_gradient;
for(unsigned int l=0; l<n_var_dofs; l++)
{
displacement_gradient.add_scaled(dphi[l][qp], system.current_solution(dof_indices_var[C_k][l]));
}
for(unsigned int C_l=0; C_l<3; C_l++)
{
elem_sigma(C_i,C_j) += JxW[qp]*( eval_elasticity_tensor(C_i,C_j,C_k,C_l) * displacement_gradient(C_l) );
}
}
}
Get the average stresses by dividing by the element volume
elem_sigma.scale(1./elem->volume());
load elem_sigma data into stress_system
for(unsigned int i=0; i<3; i++)
for(unsigned int j=0; j<3; j++)
{
stress_dof_map.dof_indices (elem, stress_dof_indices_var, sigma_vars[i][j]);
We are using CONSTANT MONOMIAL basis functions, hence we only need to get
one dof index per variable
dof_id_type dof_index = stress_dof_indices_var[0];
if( (stress_system.solution->first_local_index() <= dof_index) &&
(dof_index < stress_system.solution->last_local_index()) )
{
stress_system.solution->set(dof_index, elem_sigma(i,j));
}
}
Also, the von Mises stress
Number vonMises_value = std::sqrt( 0.5*( pow(elem_sigma(0,0) - elem_sigma(1,1),2.) +
pow(elem_sigma(1,1) - elem_sigma(2,2),2.) +
pow(elem_sigma(2,2) - elem_sigma(0,0),2.) +
6.*(pow(elem_sigma(0,1),2.) + pow(elem_sigma(1,2),2.) + pow(elem_sigma(2,0),2.))
) );
stress_dof_map.dof_indices (elem, stress_dof_indices_var, vonMises_var);
dof_id_type dof_index = stress_dof_indices_var[0];
if( (stress_system.solution->first_local_index() <= dof_index) &&
(dof_index < stress_system.solution->last_local_index()) )
{
stress_system.solution->set(dof_index, vonMises_value);
}
}
Should call close and update when we set vector entries directly
stress_system.solution->close();
stress_system.update();
}
Real kronecker_delta(unsigned int i,
unsigned int j)
{
return i == j ? 1. : 0.;
}
Real eval_elasticity_tensor(unsigned int i,
unsigned int j,
unsigned int k,
unsigned int l)
{
Define the Poisson ratio and Young's modulus
const Real nu = 0.3;
const Real E = 1.;
Define the Lame constants (lambda_1 and lambda_2) based on nu and E
const Real lambda_1 = E * nu / ( (1. + nu) * (1. - 2.*nu) );
const Real lambda_2 = 0.5 * E / (1. + nu);
return lambda_1 * kronecker_delta(i,j) * kronecker_delta(k,l)
+ lambda_2 * (kronecker_delta(i,k) * kronecker_delta(j,l) + kronecker_delta(i,l) * kronecker_delta(j,k));
}
The source file systems_of_equations_ex6.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"
#define x_scaling 1.3
#define x_load 0.0
#define y_load 0.0
#define z_load -1.0
#define BOUNDARY_ID_MIN_Z 0
#define BOUNDARY_ID_MIN_Y 1
#define BOUNDARY_ID_MAX_X 2
#define BOUNDARY_ID_MAX_Y 3
#define BOUNDARY_ID_MIN_X 4
#define BOUNDARY_ID_MAX_Z 5
using namespace libMesh;
void assemble_elasticity(EquationSystems& es,
const std::string& system_name);
void compute_stresses(EquationSystems& es);
Real kronecker_delta(unsigned int i,
unsigned int j);
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 = 3;
libmesh_example_assert(dim == LIBMESH_DIM, "3D support");
Mesh mesh(dim);
MeshTools::Generation::build_cube (mesh,
40,
8,
4,
0., 1.*x_scaling,
0., 0.3,
0., 0.1,
HEX8);
mesh.print_info();
EquationSystems equation_systems (mesh);
LinearImplicitSystem& system =
equation_systems.add_system<LinearImplicitSystem> ("Elasticity");
unsigned int u_var = system.add_variable("u", FIRST, LAGRANGE);
unsigned int v_var = system.add_variable("v", FIRST, LAGRANGE);
unsigned int w_var = system.add_variable("w", FIRST, LAGRANGE);
system.attach_assemble_function (assemble_elasticity);
std::set<boundary_id_type> boundary_ids;
boundary_ids.insert(BOUNDARY_ID_MIN_X);
std::vector<unsigned int> variables;
variables.push_back(u_var);
variables.push_back(v_var);
variables.push_back(w_var);
ZeroFunction<> zf;
DirichletBoundary dirichlet_bc(boundary_ids,
variables,
&zf);
system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
ExplicitSystem& stress_system =
equation_systems.add_system<ExplicitSystem> ("StressSystem");
stress_system.add_variable("sigma_00", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_01", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_02", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_10", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_11", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_12", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_20", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_21", CONSTANT, MONOMIAL);
stress_system.add_variable("sigma_22", CONSTANT, MONOMIAL);
stress_system.add_variable("vonMises", CONSTANT, MONOMIAL);
equation_systems.init();
equation_systems.print_info();
system.solve();
compute_stresses(equation_systems);
#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 n_components = 3;
const unsigned int u_var = system.variable_number ("u");
const unsigned int v_var = system.variable_number ("v");
const unsigned int w_var = system.variable_number ("w");
const DofMap& dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(u_var);
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), Kuw(Ke),
Kvu(Ke), Kvv(Ke), Kvw(Ke),
Kwu(Ke), Kwv(Ke), Kww(Ke);
DenseSubVector<Number>
Fu(Fe),
Fv(Fe),
Fw(Fe);
std::vector<dof_id_type> dof_indices;
std::vector<dof_id_type> dof_indices_u;
std::vector<dof_id_type> dof_indices_v;
std::vector<dof_id_type> dof_indices_w;
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);
dof_map.dof_indices (elem, dof_indices_w, w_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();
const unsigned int n_w_dofs = dof_indices_w.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);
Kuw.reposition (u_var*n_u_dofs, w_var*n_u_dofs, n_u_dofs, n_w_dofs);
Kvu.reposition (v_var*n_u_dofs, u_var*n_u_dofs, n_v_dofs, n_u_dofs);
Kvv.reposition (v_var*n_u_dofs, v_var*n_u_dofs, n_v_dofs, n_v_dofs);
Kvw.reposition (v_var*n_u_dofs, w_var*n_u_dofs, n_v_dofs, n_w_dofs);
Kwu.reposition (w_var*n_u_dofs, u_var*n_u_dofs, n_w_dofs, n_u_dofs);
Kwv.reposition (w_var*n_u_dofs, v_var*n_u_dofs, n_w_dofs, n_v_dofs);
Kww.reposition (w_var*n_u_dofs, w_var*n_u_dofs, n_w_dofs, n_w_dofs);
Fu.reposition (u_var*n_u_dofs, n_u_dofs);
Fv.reposition (v_var*n_u_dofs, n_v_dofs);
Fw.reposition (w_var*n_u_dofs, n_w_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=0, C_k=0;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
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=0, C_k=1;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
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_u_dofs; i++)
for (unsigned int j=0; j<n_w_dofs; j++)
{
unsigned int C_i=0, C_k=2;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
Kuw(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 = 1, C_k = 0;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
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 = 1, C_k = 1;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
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 i=0; i<n_v_dofs; i++)
for (unsigned int j=0; j<n_w_dofs; j++)
{
unsigned int C_i = 1, C_k = 2;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
Kvw(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_w_dofs; i++)
for (unsigned int j=0; j<n_u_dofs; j++)
{
unsigned int C_i = 2, C_k = 0;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
Kwu(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_w_dofs; i++)
for (unsigned int j=0; j<n_v_dofs; j++)
{
unsigned int C_i = 2, C_k = 1;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
Kwv(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_w_dofs; i++)
for (unsigned int j=0; j<n_w_dofs; j++)
{
unsigned int C_i = 2, C_k = 2;
for(unsigned int C_j=0; C_j<n_components; C_j++)
for(unsigned int C_l=0; C_l<n_components; C_l++)
{
Kww(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);
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
if( mesh.boundary_info->has_boundary_id(elem, side, BOUNDARY_ID_MAX_X) )
{
for (unsigned int i=0; i<n_v_dofs; i++)
{
Fu(i) += JxW_face[qp] * x_load * phi_face[i][qp];
}
for (unsigned int i=0; i<n_v_dofs; i++)
{
Fv(i) += JxW_face[qp] * y_load * phi_face[i][qp];
}
for (unsigned int i=0; i<n_v_dofs; i++)
{
Fw(i) += JxW_face[qp] * z_load * 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);
}
}
void compute_stresses(EquationSystems& es)
{
const MeshBase& mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Elasticity");
unsigned int displacement_vars[3];
displacement_vars[0] = system.variable_number ("u");
displacement_vars[1] = system.variable_number ("v");
displacement_vars[2] = system.variable_number ("w");
const unsigned int u_var = system.variable_number ("u");
const DofMap& dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(u_var);
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
QGauss qrule (dim, fe_type.default_quadrature_order());
fe->attach_quadrature_rule (&qrule);
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
ExplicitSystem& stress_system = es.get_system<ExplicitSystem>("StressSystem");
const DofMap& stress_dof_map = stress_system.get_dof_map();
unsigned int sigma_vars[3][3];
sigma_vars[0][0] = stress_system.variable_number ("sigma_00");
sigma_vars[0][1] = stress_system.variable_number ("sigma_01");
sigma_vars[0][2] = stress_system.variable_number ("sigma_02");
sigma_vars[1][0] = stress_system.variable_number ("sigma_10");
sigma_vars[1][1] = stress_system.variable_number ("sigma_11");
sigma_vars[1][2] = stress_system.variable_number ("sigma_12");
sigma_vars[2][0] = stress_system.variable_number ("sigma_20");
sigma_vars[2][1] = stress_system.variable_number ("sigma_21");
sigma_vars[2][2] = stress_system.variable_number ("sigma_22");
unsigned int vonMises_var = stress_system.variable_number ("vonMises");
std::vector< std::vector<dof_id_type> > dof_indices_var(system.n_vars());
std::vector<dof_id_type> stress_dof_indices_var;
DenseMatrix<Number> elem_sigma;
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;
for(unsigned int var=0; var<3; var++)
{
dof_map.dof_indices (elem, dof_indices_var[var], displacement_vars[var]);
}
fe->reinit (elem);
elem_sigma.resize(3,3);
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
for(unsigned int C_i=0; C_i<3; C_i++)
for(unsigned int C_j=0; C_j<3; C_j++)
for(unsigned int C_k=0; C_k<3; C_k++)
{
const unsigned int n_var_dofs = dof_indices_var[C_k].size();
Gradient displacement_gradient;
for(unsigned int l=0; l<n_var_dofs; l++)
{
displacement_gradient.add_scaled(dphi[l][qp], system.current_solution(dof_indices_var[C_k][l]));
}
for(unsigned int C_l=0; C_l<3; C_l++)
{
elem_sigma(C_i,C_j) += JxW[qp]*( eval_elasticity_tensor(C_i,C_j,C_k,C_l) * displacement_gradient(C_l) );
}
}
}
elem_sigma.scale(1./elem->volume());
for(unsigned int i=0; i<3; i++)
for(unsigned int j=0; j<3; j++)
{
stress_dof_map.dof_indices (elem, stress_dof_indices_var, sigma_vars[i][j]);
dof_id_type dof_index = stress_dof_indices_var[0];
if( (stress_system.solution->first_local_index() <= dof_index) &&
(dof_index < stress_system.solution->last_local_index()) )
{
stress_system.solution->set(dof_index, elem_sigma(i,j));
}
}
Number vonMises_value = std::sqrt( 0.5*( pow(elem_sigma(0,0) - elem_sigma(1,1),2.) +
pow(elem_sigma(1,1) - elem_sigma(2,2),2.) +
pow(elem_sigma(2,2) - elem_sigma(0,0),2.) +
6.*(pow(elem_sigma(0,1),2.) + pow(elem_sigma(1,2),2.) + pow(elem_sigma(2,0),2.))
) );
stress_dof_map.dof_indices (elem, stress_dof_indices_var, vonMises_var);
dof_id_type dof_index = stress_dof_indices_var[0];
if( (stress_system.solution->first_local_index() <= dof_index) &&
(dof_index < stress_system.solution->last_local_index()) )
{
stress_system.solution->set(dof_index, vonMises_value);
}
}
stress_system.solution->close();
stress_system.update();
}
Real kronecker_delta(unsigned int i,
unsigned int j)
{
return i == j ? 1. : 0.;
}
Real eval_elasticity_tensor(unsigned int i,
unsigned int j,
unsigned int k,
unsigned int l)
{
const Real nu = 0.3;
const Real E = 1.;
const Real lambda_1 = E * nu / ( (1. + nu) * (1. - 2.*nu) );
const Real lambda_2 = 0.5 * E / (1. + nu);
return lambda_1 * kronecker_delta(i,j) * kronecker_delta(k,l)
+ lambda_2 * (kronecker_delta(i,k) * kronecker_delta(j,l) + kronecker_delta(i,l) * kronecker_delta(j,k));
}
The console output of the program:
***************************************************************
* Running Example systems_of_equations_ex6:
* mpirun -np 2 example-devel -ksp_type cg -pc_type bjacobi -sub_pc_type icc -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()=3
spatial_dimension()=3
n_nodes()=1845
n_local_nodes()=945
n_elem()=1280
n_local_elem()=640
n_active_elem()=1280
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
EquationSystems
n_systems()=2
System #0, "Elasticity"
Type "LinearImplicit"
Variables={ "u" "v" "w" }
Finite Element Types="LAGRANGE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="FIRST", "THIRD"
n_dofs()=5535
n_local_dofs()=2835
n_constrained_dofs()=135
n_local_constrained_dofs()=135
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 63.4146
Average Off-Processor Bandwidth <= 1.05691
Maximum On-Processor Bandwidth <= 81
Maximum Off-Processor Bandwidth <= 27
DofMap Constraints
Number of DoF Constraints = 135
Average DoF Constraint Length= 0
Number of Node Constraints = 0
System #1, "StressSystem"
Type "Explicit"
Variables={ "sigma_00" "sigma_01" "sigma_02" "sigma_10" "sigma_11" "sigma_12" "sigma_20" "sigma_21" "sigma_22" "vonMises" }
Finite Element Types="MONOMIAL", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="CONSTANT", "THIRD"
n_dofs()=12800
n_local_dofs()=6400
n_constrained_dofs()=0
n_local_constrained_dofs()=0
n_vectors()=1
n_matrices()=0
DofMap Sparsity
Average On-Processor Bandwidth <= 0
Average Off-Processor Bandwidth <= 0
Maximum On-Processor Bandwidth <= 0
Maximum Off-Processor Bandwidth <= 0
DofMap Constraints
Number of DoF Constraints = 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_ex6/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:45:40 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): 4.198e+00 1.00000 4.198e+00
Objects: 4.100e+01 1.00000 4.100e+01
Flops: 2.617e+08 1.02245 2.589e+08 5.177e+08
Flops/sec: 6.235e+07 1.02245 6.166e+07 1.233e+08
MPI Messages: 1.250e+01 1.00000 1.250e+01 2.500e+01
MPI Message Lengths: 8.502e+04 1.00000 6.801e+03 1.700e+05
MPI Reductions: 6.900e+01 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: 4.1979e+00 100.0% 5.1772e+08 100.0% 2.500e+01 100.0% 6.801e+03 100.0% 6.800e+01 98.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 3 1.0 3.4094e-05 1.2 1.70e+04 1.1 0.0e+00 0.0e+00 3.0e+00 0 0 0 0 4 0 0 0 0 4 974
VecNorm 3 1.0 6.9370e-03 1.0 1.70e+04 1.1 0.0e+00 0.0e+00 3.0e+00 0 0 0 0 4 0 0 0 0 4 5
VecCopy 3 1.0 1.2159e-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
VecSet 11 1.0 1.8835e-05 1.1 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 2 1.0 3.3140e-05 1.1 1.13e+04 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 668
VecAYPX 1 1.0 5.0068e-06 1.0 2.84e+03 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 1105
VecAssemblyBegin 5 1.0 6.3801e-04 8.5 0.00e+00 0.0 2.0e+00 2.3e+03 1.5e+01 0 0 8 3 22 0 0 8 3 22 0
VecAssemblyEnd 5 1.0 2.3127e-05 1.1 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 4 1.0 6.4850e-05 1.2 0.00e+00 0.0 8.0e+00 1.6e+03 0.0e+00 0 0 32 7 0 0 0 32 7 0 0
VecScatterEnd 4 1.0 1.7468e-021356.8 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 2 1.0 1.7936e-0237.3 7.20e+05 1.1 4.0e+00 1.1e+03 0.0e+00 0 0 16 3 0 0 0 16 3 0 78
MatSolve 3 1.0 3.7189e-03 1.1 6.82e+06 1.1 0.0e+00 0.0e+00 0.0e+00 0 3 0 0 0 0 3 0 0 0 3566
MatLUFactorNum 1 1.0 9.4637e-02 1.0 2.54e+08 1.0 0.0e+00 0.0e+00 0.0e+00 2 97 0 0 0 2 97 0 0 0 5315
MatILUFactorSym 1 1.0 2.7927e-01 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 3.0e+00 6 0 0 0 4 6 0 0 0 4 0
MatAssemblyBegin 2 1.0 1.5917e-02137.4 0.00e+00 0.0 3.0e+00 4.9e+04 4.0e+00 0 0 12 87 6 0 0 12 87 6 0
MatAssemblyEnd 2 1.0 1.2932e-03 1.1 0.00e+00 0.0 4.0e+00 2.7e+02 8.0e+00 0 0 16 1 12 0 0 16 1 12 0
MatGetRowIJ 1 1.0 5.0068e-06 1.2 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.1802e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 4.0e+00 0 0 0 0 6 0 0 0 0 6 0
MatZeroEntries 3 1.0 5.0831e-04 1.3 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 9.2983e-05 1.1 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 3.8583e-01 1.0 2.62e+08 1.0 4.0e+00 1.1e+03 1.5e+01 9100 16 3 22 9100 16 3 22 1342
PCSetUp 2 1.0 3.7444e-01 1.0 2.54e+08 1.0 0.0e+00 0.0e+00 9.0e+00 9 97 0 0 13 9 97 0 0 13 1343
PCSetUpOnBlocks 1 1.0 3.7415e-01 1.0 2.54e+08 1.0 0.0e+00 0.0e+00 7.0e+00 9 97 0 0 10 9 97 0 0 10 1344
PCApply 3 1.0 3.8271e-03 1.0 6.82e+06 1.1 0.0e+00 0.0e+00 0.0e+00 0 3 0 0 0 0 3 0 0 0 3465
------------------------------------------------------------------------------------------------------------------------
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 16 16 340696 0
Vector Scatter 3 3 3108 0
Index Set 11 11 36980 0
IS L to G Mapping 2 2 1128 0
Matrix 4 4 15965512 0
Krylov Solver 2 2 2368 0
Preconditioner 2 2 1784 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 0
Average time for MPI_Barrier(): 2.38419e-06
Average time for zero size MPI_Send(): 1.9908e-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:45:40 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=4.20892, Active time=4.04758 |
----------------------------------------------------------------------------------------------------------------
| 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() 2 0.1721 0.086031 0.1884 0.094195 4.25 4.65 |
| build_constraint_matrix() 640 0.0063 0.000010 0.0063 0.000010 0.16 0.16 |
| build_sparsity() 1 0.1789 0.178855 0.2855 0.285462 4.42 7.05 |
| cnstrn_elem_mat_vec() 640 0.0086 0.000013 0.0086 0.000013 0.21 0.21 |
| create_dof_constraints() 2 0.0442 0.022083 0.2601 0.130057 1.09 6.43 |
| distribute_dofs() 2 0.0364 0.018188 0.0983 0.049144 0.90 2.43 |
| dof_indices() 23808 0.9057 0.000038 0.9057 0.000038 22.38 22.38 |
| prepare_send_list() 2 0.0004 0.000224 0.0004 0.000224 0.01 0.01 |
| reinit() 2 0.0615 0.030758 0.0615 0.030758 1.52 1.52 |
| |
| EquationSystems |
| build_solution_vector() 1 0.0248 0.024809 0.2039 0.203928 0.61 5.04 |
| |
| ExodusII_IO |
| write_nodal_data() 1 0.0097 0.009669 0.0097 0.009669 0.24 0.24 |
| |
| FE |
| compute_shape_functions() 1792 0.0242 0.000014 0.0242 0.000014 0.60 0.60 |
| init_shape_functions() 514 0.0014 0.000003 0.0014 0.000003 0.03 0.03 |
| |
| FEMap |
| compute_affine_map() 1792 0.0182 0.000010 0.0182 0.000010 0.45 0.45 |
| compute_face_map() 512 0.0054 0.000011 0.0054 0.000011 0.13 0.13 |
| init_face_shape_functions() 1 0.0000 0.000032 0.0000 0.000032 0.00 0.00 |
| init_reference_to_physical_map() 514 0.0133 0.000026 0.0133 0.000026 0.33 0.33 |
| |
| Mesh |
| find_neighbors() 1 0.0232 0.023213 0.0241 0.024098 0.57 0.60 |
| renumber_nodes_and_elem() 2 0.0011 0.000533 0.0011 0.000533 0.03 0.03 |
| |
| MeshCommunication |
| compute_hilbert_indices() 2 0.0114 0.005683 0.0114 0.005683 0.28 0.28 |
| find_global_indices() 2 0.0042 0.002120 0.0173 0.008639 0.10 0.43 |
| parallel_sort() 2 0.0013 0.000632 0.0013 0.000669 0.03 0.03 |
| |
| MeshOutput |
| write_equation_systems() 1 0.0001 0.000084 0.2137 0.213732 0.00 5.28 |
| |
| MeshTools::Generation |
| build_cube() 1 0.0054 0.005425 0.0054 0.005425 0.13 0.13 |
| |
| MetisPartitioner |
| partition() 1 0.0167 0.016747 0.0250 0.025020 0.41 0.62 |
| |
| Parallel |
| allgather() 12 0.0002 0.000018 0.0002 0.000019 0.01 0.01 |
| max(bool) 1 0.0000 0.000004 0.0000 0.000004 0.00 0.00 |
| max(scalar) 137 0.0003 0.000003 0.0003 0.000003 0.01 0.01 |
| max(vector) 31 0.0002 0.000006 0.0004 0.000013 0.00 0.01 |
| min(bool) 157 0.0004 0.000002 0.0004 0.000002 0.01 0.01 |
| min(scalar) 128 0.0016 0.000012 0.0016 0.000012 0.04 0.04 |
| min(vector) 31 0.0002 0.000008 0.0005 0.000016 0.01 0.01 |
| probe() 16 0.0004 0.000025 0.0004 0.000025 0.01 0.01 |
| receive() 16 0.0002 0.000010 0.0006 0.000035 0.00 0.01 |
| send() 16 0.0001 0.000004 0.0001 0.000004 0.00 0.00 |
| send_receive() 20 0.0003 0.000016 0.0010 0.000050 0.01 0.02 |
| sum() 28 0.0006 0.000022 0.0007 0.000025 0.02 0.02 |
| |
| Parallel::Request |
| wait() 16 0.0000 0.000003 0.0000 0.000003 0.00 0.00 |
| |
| Partitioner |
| set_node_processor_ids() 1 0.0016 0.001613 0.0017 0.001693 0.04 0.04 |
| set_parent_processor_ids() 1 0.0010 0.000958 0.0010 0.000958 0.02 0.02 |
| |
| PetscLinearSolver |
| solve() 1 0.3880 0.387974 0.3880 0.387974 9.59 9.59 |
| |
| System |
| assemble() 1 2.0780 2.078003 2.3550 2.355004 51.34 58.18 |
----------------------------------------------------------------------------------------------------------------
| Totals: 30851 4.0476 100.00 |
----------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example systems_of_equations_ex6:
* mpirun -np 2 example-devel -ksp_type cg -pc_type bjacobi -sub_pc_type icc -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
Systems Example 6 - 3D Linear Elastic Cantilever
By David KnezevicThis is a 3D version of systems_of_equations_ex4. In this case we also compute and plot stresses.
C++ include files that we need