The source file systems_of_equations_ex5.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/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);
Add a SCALAR variable for the Lagrange multiplier to enforce our constraint
unsigned int lambda_var = system.add_variable("lambda", FIRST, SCALAR);
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 unsigned int lambda_var = system.variable_number ("lambda");
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);
DenseSubMatrix<Number> Klambda_v(Ke), Kv_lambda(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;
std::vector<dof_id_type> dof_indices_lambda;
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_lambda, lambda_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_lambda_dofs = dof_indices_lambda.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);
Also, add a row and a column to enforce the constraint
Kv_lambda.reposition (v_var*n_u_dofs, v_var*n_u_dofs+n_v_dofs, n_v_dofs, 1);
Klambda_v.reposition (v_var*n_v_dofs+n_v_dofs, v_var*n_v_dofs, 1, 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<boundary_id_type> bc_ids =
mesh.boundary_info->boundary_ids (elem,side);
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 (std::vector<boundary_id_type>::const_iterator b =
bc_ids.begin(); b != bc_ids.end(); ++b)
{
const boundary_id_type bc_id = *b;
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
Add the loading
if( bc_id == 2 )
{
for (unsigned int i=0; i<n_v_dofs; i++)
{
Fv(i) += JxW_face[qp]* (-1.) * phi_face[i][qp];
}
}
Add the constraint contributions
if( bc_id == 1 )
{
for (unsigned int i=0; i<n_v_dofs; i++)
for (unsigned int j=0; j<n_lambda_dofs; j++)
{
Kv_lambda(i,j) += JxW_face[qp]* (-1.) * phi_face[i][qp];
}
for (unsigned int i=0; i<n_lambda_dofs; i++)
for (unsigned int j=0; j<n_v_dofs; j++)
{
Klambda_v(i,j) += JxW_face[qp]* (-1.) * phi_face[j][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_ex5.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/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);
unsigned int lambda_var = system.add_variable("lambda", FIRST, SCALAR);
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 unsigned int lambda_var = system.variable_number ("lambda");
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);
DenseSubMatrix<Number> Klambda_v(Ke), Kv_lambda(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;
std::vector<dof_id_type> dof_indices_lambda;
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_lambda, lambda_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_lambda_dofs = dof_indices_lambda.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);
Kv_lambda.reposition (v_var*n_u_dofs, v_var*n_u_dofs+n_v_dofs, n_v_dofs, 1);
Klambda_v.reposition (v_var*n_v_dofs+n_v_dofs, v_var*n_v_dofs, 1, 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<boundary_id_type> bc_ids =
mesh.boundary_info->boundary_ids (elem,side);
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 (std::vector<boundary_id_type>::const_iterator b =
bc_ids.begin(); b != bc_ids.end(); ++b)
{
const boundary_id_type bc_id = *b;
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
if( bc_id == 2 )
{
for (unsigned int i=0; i<n_v_dofs; i++)
{
Fv(i) += JxW_face[qp]* (-1.) * phi_face[i][qp];
}
}
if( bc_id == 1 )
{
for (unsigned int i=0; i<n_v_dofs; i++)
for (unsigned int j=0; j<n_lambda_dofs; j++)
{
Kv_lambda(i,j) += JxW_face[qp]* (-1.) * phi_face[i][qp];
}
for (unsigned int i=0; i<n_lambda_dofs; i++)
for (unsigned int j=0; j<n_v_dofs; j++)
{
Klambda_v(i,j) += JxW_face[qp]* (-1.) * phi_face[j][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_ex5:
* 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" } "lambda"
Finite Element Types="LAGRANGE", "JACOBI_20_00" "SCALAR", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN" "CARTESIAN"
Approximation Orders="SECOND", "THIRD" "FIRST", "THIRD"
n_dofs()=4243
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 <= 31.3818
Average Off-Processor Bandwidth <= 1.32501
Maximum On-Processor Bandwidth <= 2101
Maximum Off-Processor Bandwidth <= 2142
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_ex5/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:45:09 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): 6.836e-01 1.00000 6.836e-01
Objects: 3.300e+01 1.00000 3.300e+01
Flops: 5.416e+07 1.87204 4.154e+07 8.309e+07
Flops/sec: 7.922e+07 1.87204 6.077e+07 1.215e+08
MPI Messages: 1.200e+01 1.00000 1.200e+01 2.400e+01
MPI Message Lengths: 7.253e+04 1.00000 6.044e+03 1.451e+05
MPI Reductions: 5.200e+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: 6.8356e-01 100.0% 8.3087e+07 100.0% 2.400e+01 100.0% 6.044e+03 100.0% 5.100e+01 98.1%
------------------------------------------------------------------------------------------------------------------------
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 2 1.0 1.9073e-05 1.0 8.57e+03 1.0 0.0e+00 0.0e+00 2.0e+00 0 0 0 0 4 0 0 0 0 4 890
VecNorm 3 1.0 6.7251e-03 1.0 1.29e+04 1.0 0.0e+00 0.0e+00 3.0e+00 1 0 0 0 6 1 0 0 0 6 4
VecCopy 2 1.0 4.0531e-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 8 1.0 1.1444e-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 2.1935e-05 1.0 8.57e+03 1.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 774
VecAYPX 1 1.0 5.0068e-06 1.2 2.14e+03 1.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 847
VecAssemblyBegin 3 1.0 6.3896e-05 1.0 0.00e+00 0.0 4.0e+00 9.3e+02 9.0e+00 0 0 17 3 17 0 0 17 3 18 0
VecAssemblyEnd 3 1.0 2.0266e-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 3 1.0 6.7949e-05 1.7 0.00e+00 0.0 6.0e+00 6.2e+03 0.0e+00 0 0 25 26 0 0 0 25 26 0 0
VecScatterEnd 3 1.0 1.4280e-02820.5 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 1 0 0 0 0 1 0 0 0 0 0
MatMult 2 1.0 1.4525e-0266.9 2.78e+05 1.0 4.0e+00 8.9e+03 0.0e+00 1 1 17 25 0 1 1 17 25 0 38
MatSolve 3 1.0 8.5425e-04 1.1 1.96e+06 1.0 0.0e+00 0.0e+00 0.0e+00 0 5 0 0 0 0 5 0 0 0 4573
MatLUFactorNum 1 1.0 2.2945e-02 1.8 5.19e+07 1.9 0.0e+00 0.0e+00 0.0e+00 3 95 0 0 0 3 95 0 0 0 3424
MatILUFactorSym 1 1.0 4.1117e-02 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 3.0e+00 6 0 0 0 6 6 0 0 0 6 0
MatAssemblyBegin 2 1.0 7.3969e-0344.3 0.00e+00 0.0 6.0e+00 1.6e+04 4.0e+00 1 0 25 65 8 1 0 25 65 8 0
MatAssemblyEnd 2 1.0 8.4805e-04 1.1 0.00e+00 0.0 4.0e+00 2.2e+03 8.0e+00 0 0 17 6 15 0 0 17 6 16 0
MatGetRowIJ 1 1.0 3.0994e-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
MatGetOrdering 1 1.0 1.1611e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 4.0e+00 0 0 0 0 8 0 0 0 0 8 0
MatZeroEntries 3 1.0 2.2197e-04 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
KSPSetUp 2 1.0 8.6069e-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 7.2719e-02 1.0 5.42e+07 1.9 4.0e+00 8.9e+03 1.4e+01 11100 17 25 27 11100 17 25 27 1143
PCSetUp 2 1.0 6.4539e-02 1.3 5.19e+07 1.9 0.0e+00 0.0e+00 9.0e+00 8 95 0 0 17 8 95 0 0 18 1217
PCSetUpOnBlocks 1 1.0 6.4283e-02 1.3 5.19e+07 1.9 0.0e+00 0.0e+00 7.0e+00 8 95 0 0 13 8 95 0 0 14 1222
PCApply 3 1.0 9.5892e-04 1.1 1.96e+06 1.0 0.0e+00 0.0e+00 0.0e+00 0 5 0 0 0 0 5 0 0 0 4074
------------------------------------------------------------------------------------------------------------------------
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 139008 0
Vector Scatter 2 2 2072 0
Index Set 9 9 28544 0
IS L to G Mapping 1 1 564 0
Matrix 4 4 4837496 0
Krylov Solver 2 2 2368 0
Preconditioner 2 2 1784 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 9.53674e-08
Average time for MPI_Barrier(): 2.6226e-06
Average time for zero size MPI_Send(): 1.60933e-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:09 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.69837, Active time=0.668448 |
----------------------------------------------------------------------------------------------------------------
| 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 |
| SCALAR_dof_indices() 1020 0.0035 0.000003 0.0035 0.000003 0.52 0.52 |
| add_neighbors_to_send_list() 1 0.0317 0.031698 0.0346 0.034579 4.74 5.17 |
| build_constraint_matrix() 250 0.0013 0.000005 0.0013 0.000005 0.20 0.20 |
| build_sparsity() 1 0.0439 0.043869 0.0812 0.081166 6.56 12.14 |
| cnstrn_elem_mat_vec() 250 0.0003 0.000001 0.0003 0.000001 0.04 0.04 |
| create_dof_constraints() 1 0.0094 0.009353 0.0789 0.078917 1.40 11.81 |
| distribute_dofs() 1 0.0171 0.017116 0.0404 0.040362 2.56 6.04 |
| dof_indices() 3020 0.2102 0.000070 0.2138 0.000071 31.44 31.98 |
| prepare_send_list() 1 0.0001 0.000069 0.0001 0.000069 0.01 0.01 |
| reinit() 1 0.0211 0.021106 0.0211 0.021106 3.16 3.16 |
| |
| EquationSystems |
| build_solution_vector() 1 0.0032 0.003154 0.0400 0.040000 0.47 5.98 |
| |
| ExodusII_IO |
| write_nodal_data() 1 0.0043 0.004331 0.0043 0.004331 0.65 0.65 |
| |
| FE |
| compute_shape_functions() 330 0.0038 0.000012 0.0038 0.000012 0.58 0.58 |
| init_shape_functions() 81 0.0003 0.000004 0.0003 0.000004 0.05 0.05 |
| inverse_map() 240 0.0020 0.000008 0.0020 0.000008 0.29 0.29 |
| |
| FEMap |
| compute_affine_map() 330 0.0031 0.000009 0.0031 0.000009 0.46 0.46 |
| compute_face_map() 80 0.0014 0.000018 0.0034 0.000043 0.21 0.51 |
| init_face_shape_functions() 21 0.0001 0.000007 0.0001 0.000007 0.02 0.02 |
| init_reference_to_physical_map() 81 0.0017 0.000021 0.0017 0.000021 0.25 0.25 |
| |
| Mesh |
| find_neighbors() 1 0.0057 0.005710 0.0060 0.005966 0.85 0.89 |
| renumber_nodes_and_elem() 2 0.0006 0.000314 0.0006 0.000314 0.09 0.09 |
| |
| MeshCommunication |
| compute_hilbert_indices() 2 0.0046 0.002307 0.0046 0.002307 0.69 0.69 |
| find_global_indices() 2 0.0019 0.000932 0.0076 0.003815 0.28 1.14 |
| parallel_sort() 2 0.0008 0.000401 0.0009 0.000442 0.12 0.13 |
| |
| MeshOutput |
| write_equation_systems() 1 0.0001 0.000083 0.0445 0.044471 0.01 6.65 |
| |
| MeshTools::Generation |
| build_cube() 1 0.0032 0.003226 0.0032 0.003226 0.48 0.48 |
| |
| MetisPartitioner |
| partition() 1 0.0062 0.006247 0.0097 0.009736 0.93 1.46 |
| |
| Parallel |
| allgather() 9 0.0021 0.000233 0.0021 0.000235 0.31 0.32 |
| 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.04 0.04 |
| max(vector) 24 0.0002 0.000006 0.0003 0.000014 0.02 0.05 |
| min(bool) 121 0.0003 0.000002 0.0003 0.000002 0.04 0.04 |
| min(scalar) 99 0.0035 0.000035 0.0035 0.000035 0.52 0.52 |
| min(vector) 24 0.0002 0.000009 0.0005 0.000019 0.03 0.07 |
| probe() 12 0.0001 0.000009 0.0001 0.000009 0.02 0.02 |
| receive() 12 0.0001 0.000008 0.0002 0.000018 0.02 0.03 |
| send() 12 0.0001 0.000005 0.0001 0.000005 0.01 0.01 |
| send_receive() 16 0.0002 0.000013 0.0005 0.000033 0.03 0.08 |
| sum() 20 0.0002 0.000012 0.0004 0.000022 0.04 0.06 |
| |
| 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.001110 0.0012 0.001185 0.17 0.18 |
| set_parent_processor_ids() 1 0.0004 0.000422 0.0004 0.000422 0.06 0.06 |
| |
| PetscLinearSolver |
| solve() 1 0.0817 0.081707 0.0817 0.081707 12.22 12.22 |
| |
| System |
| assemble() 1 0.1963 0.196319 0.2808 0.280796 29.37 42.01 |
----------------------------------------------------------------------------------------------------------------
| Totals: 6194 0.6684 100.00 |
----------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example systems_of_equations_ex5:
* 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 5 - Linear Elastic Cantilever with Constraint
By David KnezevicIn this example we extend systems_of_equations_ex4 to enforce a constraint. We apply a uniform load on the top surface of the cantilever, and we determine the traction on the right boundary in order to obtain zero average vertical displacement on the right boundary of the domain.
This constraint is enforced via a Lagrange multiplier (SCALAR variable). The system we solve, therefore, is of the form: a(u,v) + \lambda g(v) = f(v) g(u) = 0 Here \lambda tells us the traction required to satisfy the constraint.
C++ include files that we need