The source file miscellaneous_ex5.C with comments:
#include <iostream>
LibMesh include files.
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/equation_systems.h"
#include "libmesh/mesh_data.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_modification.h"
#include "libmesh/elem.h"
#include "libmesh/transient_system.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/dof_map.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/dense_submatrix.h"
#include "libmesh/dense_subvector.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/fe_interface.h"
#include "libmesh/getpot.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/error_vector.h"
#include "libmesh/kelly_error_estimator.h"
#include "libmesh/discontinuity_measure.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/exact_solution.h"
#define QORDER TWENTYSIXTH
Bring in everything from the libMesh namespace
Bring in everything from the libMesh namespace
using namespace libMesh;
Number exact_solution (const Point& p, const Parameters& parameters, const std::string&, const std::string&)
{
const Real x = p(0);
const Real y = p(1);
const Real z = p(2);
if (parameters.get<bool>("singularity"))
{
Real theta = atan2(y,x);
if (theta < 0)
theta += 2. * libMesh::pi;
return pow(x*x + y*y, 1./3.)*sin(2./3.*theta) + z;
}
else
{
return cos(x) * exp(y) * (1. - z);
}
}
We now define the gradient of the exact solution, again being careful
to obtain an angle from atan2 in the correct
quadrant.
Gradient exact_derivative(const Point& p,
const Parameters& parameters, // es parameters
const std::string&, // sys_name, not needed
const std::string&) // unk_name, not needed
{
Gradient value to be returned.
Gradient gradu;
x and y coordinates in space
const Real x = p(0);
const Real y = p(1);
const Real z = p(2);
if (parameters.get<bool>("singularity"))
{
libmesh_assert_not_equal_to (x, 0.);
For convenience...
const Real tt = 2./3.;
const Real ot = 1./3.;
The value of the radius, squared
const Real r2 = x*x + y*y;
The boundary value, given by the exact solution,
u_exact = r^(2/3)*sin(2*theta/3).
Real theta = atan2(y,x);
Make sure 0 <= theta <= 2*pi
if (theta < 0)
theta += 2. * libMesh::pi;
du/dx
gradu(0) = tt*x*pow(r2,-tt)*sin(tt*theta) - pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*y/x/x;
gradu(1) = tt*y*pow(r2,-tt)*sin(tt*theta) + pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*1./x;
gradu(2) = 1.;
}
else
{
gradu(0) = -sin(x) * exp(y) * (1. - z);
gradu(1) = cos(x) * exp(y) * (1. - z);
gradu(2) = -cos(x) * exp(y);
}
return gradu;
}
We now define the matrix assembly function for the
Laplace system. We need to first compute element volume
matrices, and then take into account the boundary
conditions and the flux integrals, which will be handled
via an interior penalty method.
void assemble_ellipticdg(EquationSystems& es, const std::string& system_name)
{
std::cout<<" assembling elliptic dg system... ";
std::cout.flush();
It is a good idea to make sure we are assembling
the proper system.
libmesh_assert_equal_to (system_name, "EllipticDG");
Get a constant reference to the mesh object.
const MeshBase& mesh = es.get_mesh();
The dimension that we are running
const unsigned int dim = mesh.mesh_dimension();
Get a reference to the LinearImplicitSystem we are solving
LinearImplicitSystem & ellipticdg_system = es.get_system<LinearImplicitSystem> ("EllipticDG");
Get some parameters that we need during assembly
const Real penalty = es.parameters.get<Real> ("penalty");
std::string refinement_type = es.parameters.get<std::string> ("refinement");
A reference to the \p DofMap object for this system. The \p DofMap
object handles the index translation from node and element numbers
to degree of freedom numbers. We will talk more about the \p DofMap
const DofMap & dof_map = ellipticdg_system.get_dof_map();
Get a constant reference to the Finite Element type
for the first (and only) variable in the system.
FEType fe_type = ellipticdg_system.variable_type(0);
Build a Finite Element object of the specified type. Since the
\p FEBase::build() member dynamically creates memory we will
store the object as an \p AutoPtr. This can be thought
of as a pointer that will clean up after itself.
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
AutoPtr<FEBase> fe_elem_face(FEBase::build(dim, fe_type));
AutoPtr<FEBase> fe_neighbor_face(FEBase::build(dim, fe_type));
Quadrature rules for numerical integration.
#ifdef QORDER
QGauss qrule (dim, QORDER);
#else
QGauss qrule (dim, fe_type.default_quadrature_order());
#endif
fe->attach_quadrature_rule (&qrule);
#ifdef QORDER
QGauss qface(dim-1, QORDER);
#else
QGauss qface(dim-1, fe_type.default_quadrature_order());
#endif
Tell the finite element object to use our quadrature rule.
fe_elem_face->attach_quadrature_rule(&qface);
fe_neighbor_face->attach_quadrature_rule(&qface);
Here we define some references to cell-specific data that
will be used to assemble the linear system.
Data for interior volume integrals
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
Data for surface integrals on the element boundary
const std::vector<std::vector<Real> >& phi_face = fe_elem_face->get_phi();
const std::vector<std::vector<RealGradient> >& dphi_face = fe_elem_face->get_dphi();
const std::vector<Real>& JxW_face = fe_elem_face->get_JxW();
const std::vector<Point>& qface_normals = fe_elem_face->get_normals();
const std::vector<Point>& qface_points = fe_elem_face->get_xyz();
Data for surface integrals on the neighbor boundary
const std::vector<std::vector<Real> >& phi_neighbor_face = fe_neighbor_face->get_phi();
const std::vector<std::vector<RealGradient> >& dphi_neighbor_face = fe_neighbor_face->get_dphi();
Define data structures to contain the element interior matrix
and right-hand-side vector contribution. Following
basic finite element terminology we will denote these
"Ke" and "Fe".
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
Data structures to contain the element and neighbor boundary matrix
contribution. This matrices will do the coupling beetwen the dofs of
the element and those of his neighbors.
Ken: matrix coupling elem and neighbor dofs
DenseMatrix<Number> Kne;
DenseMatrix<Number> Ken;
DenseMatrix<Number> Kee;
DenseMatrix<Number> Knn;
This vector will hold the degree of freedom indices for
the element. These define where in the global system
the element degrees of freedom get mapped.
std::vector<dof_id_type> dof_indices;
Now we will loop over all the elements in the mesh. We will
compute first the element interior matrix and right-hand-side contribution
and then the element and neighbors boundary matrix contributions.
MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
for ( ; el != end_el; ++el)
{
Store a pointer to the element we are currently
working on. This allows for nicer syntax later.
const Elem* elem = *el;
Get the degree of freedom indices for the
current element. These define where in the global
matrix and right-hand-side this element will
contribute to.
dof_map.dof_indices (elem, dof_indices);
const unsigned int n_dofs = dof_indices.size();
Compute the element-specific data for the current
element. This involves computing the location of the
quadrature points (q_point) and the shape functions
(phi, dphi) for the current element.
fe->reinit (elem);
Zero the element matrix and right-hand side before
summing them. We use the resize member here because
the number of degrees of freedom might have changed from
the last element.
Ke.resize (n_dofs, n_dofs);
Fe.resize (n_dofs);
Now we will build the element interior matrix. This involves
a double loop to integrate the test funcions (i) against
the trial functions (j).
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
for (unsigned int i=0; i<n_dofs; i++)
{
for (unsigned int j=0; j<n_dofs; j++)
{
Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
}
}
}
Now we adress boundary conditions.
We consider Dirichlet bc imposed via the interior penalty method
The following loops over the sides of the element.
If the element has no neighbor on a side then that
side MUST live on a boundary of the domain.
for (unsigned int side=0; side<elem->n_sides(); side++)
{
if (elem->neighbor(side) == NULL)
{
Pointer to the element face
fe_elem_face->reinit(elem, side);
AutoPtr<Elem> elem_side (elem->build_side(side));
h elemet dimension to compute the interior penalty penalty parameter
const unsigned int elem_b_order = static_cast<unsigned int> (fe_elem_face->get_order());
const double h_elem = elem->volume()/elem_side->volume() * 1./pow(elem_b_order, 2.);
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
Number bc_value = exact_solution(qface_points[qp], es.parameters,"null","void");
for (unsigned int i=0; i<n_dofs; i++)
{
Matrix contribution
for (unsigned int j=0; j<n_dofs; j++)
{
Ke(i,j) += JxW_face[qp] * penalty/h_elem * phi_face[i][qp] * phi_face[j][qp]; // stability
Ke(i,j) -= JxW_face[qp] * (phi_face[i][qp] * (dphi_face[j][qp]*qface_normals[qp]) + phi_face[j][qp] * (dphi_face[i][qp]*qface_normals[qp])); // consistency
}
RHS contribution
Fe(i) += JxW_face[qp] * bc_value * penalty/h_elem * phi_face[i][qp]; // stability
Fe(i) -= JxW_face[qp] * dphi_face[i][qp] * (bc_value*qface_normals[qp]); // consistency
}
}
}
If the element is not on a boundary of the domain
we loop over his neighbors to compute the element
and neighbor boundary matrix contributions
else
{
Store a pointer to the neighbor we are currently
working on.
const Elem* neighbor = elem->neighbor(side);
Get the global id of the element and the neighbor
const unsigned int elem_id = elem->id();
const unsigned int neighbor_id = neighbor->id();
If the neighbor has the same h level and is active
perform integration only if our global id is bigger than our neighbor id.
We don't want to compute twice the same contributions.
If the neighbor has a different h level perform integration
only if the neighbor is at a lower level.
if ((neighbor->active() && (neighbor->level() == elem->level()) && (elem_id < neighbor_id)) || (neighbor->level() < elem->level()))
{
Pointer to the element side
AutoPtr<Elem> elem_side (elem->build_side(side));
h dimension to compute the interior penalty penalty parameter
const unsigned int elem_b_order = static_cast<unsigned int>(fe_elem_face->get_order());
const unsigned int neighbor_b_order = static_cast<unsigned int>(fe_neighbor_face->get_order());
const double side_order = (elem_b_order + neighbor_b_order)/2.;
const double h_elem = (elem->volume()/elem_side->volume()) * 1./pow(side_order,2.);
The quadrature point locations on the neighbor side
std::vector<Point> qface_neighbor_point;
The quadrature point locations on the element side
std::vector<Point > qface_point;
Reinitialize shape functions on the element side
fe_elem_face->reinit(elem, side);
Get the physical locations of the element quadrature points
qface_point = fe_elem_face->get_xyz();
Find their locations on the neighbor
unsigned int side_neighbor = neighbor->which_neighbor_am_i(elem);
if (refinement_type == "p")
fe_neighbor_face->side_map (neighbor, elem_side.get(), side_neighbor, qface.get_points(), qface_neighbor_point);
else
FEInterface::inverse_map (elem->dim(), fe->get_fe_type(), neighbor, qface_point, qface_neighbor_point);
Calculate the neighbor element shape functions at those locations
fe_neighbor_face->reinit(neighbor, &qface_neighbor_point);
Get the degree of freedom indices for the
neighbor. These define where in the global
matrix this neighbor will contribute to.
std::vector<dof_id_type> neighbor_dof_indices;
dof_map.dof_indices (neighbor, neighbor_dof_indices);
const unsigned int n_neighbor_dofs = neighbor_dof_indices.size();
Zero the element and neighbor side matrix before
summing them. We use the resize member here because
the number of degrees of freedom might have changed from
the last element or neighbor.
Note that Kne and Ken are not square matrices if neighbor
and element have a different p level
Kne.resize (n_neighbor_dofs, n_dofs);
Ken.resize (n_dofs, n_neighbor_dofs);
Kee.resize (n_dofs, n_dofs);
Knn.resize (n_neighbor_dofs, n_neighbor_dofs);
Now we will build the element and neighbor
boundary matrices. This involves
a double loop to integrate the test funcions
(i) against the trial functions (j).
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
Kee Matrix. Integrate the element test function i
against the element test function j
for (unsigned int i=0; i<n_dofs; i++)
{
for (unsigned int j=0; j<n_dofs; j++)
{
Kee(i,j) -= 0.5 * JxW_face[qp] * (phi_face[j][qp]*(qface_normals[qp]*dphi_face[i][qp]) + phi_face[i][qp]*(qface_normals[qp]*dphi_face[j][qp])); // consistency
Kee(i,j) += JxW_face[qp] * penalty/h_elem * phi_face[j][qp]*phi_face[i][qp]; // stability
}
}
Knn Matrix. Integrate the neighbor test function i
against the neighbor test function j
for (unsigned int i=0; i<n_neighbor_dofs; i++)
{
for (unsigned int j=0; j<n_neighbor_dofs; j++)
{
Knn(i,j) += 0.5 * JxW_face[qp] * (phi_neighbor_face[j][qp]*(qface_normals[qp]*dphi_neighbor_face[i][qp]) + phi_neighbor_face[i][qp]*(qface_normals[qp]*dphi_neighbor_face[j][qp])); // consistency
Knn(i,j) += JxW_face[qp] * penalty/h_elem * phi_neighbor_face[j][qp]*phi_neighbor_face[i][qp]; // stability
}
}
Kne Matrix. Integrate the neighbor test function i
against the element test function j
for (unsigned int i=0; i<n_neighbor_dofs; i++)
{
for (unsigned int j=0; j<n_dofs; j++)
{
Kne(i,j) += 0.5 * JxW_face[qp] * (phi_neighbor_face[i][qp]*(qface_normals[qp]*dphi_face[j][qp]) - phi_face[j][qp]*(qface_normals[qp]*dphi_neighbor_face[i][qp])); // consistency
Kne(i,j) -= JxW_face[qp] * penalty/h_elem * phi_face[j][qp]*phi_neighbor_face[i][qp]; // stability
}
}
Ken Matrix. Integrate the element test function i
against the neighbor test function j
for (unsigned int i=0; i<n_dofs; i++)
{
for (unsigned int j=0; j<n_neighbor_dofs; j++)
{
Ken(i,j) += 0.5 * JxW_face[qp] * (phi_neighbor_face[j][qp]*(qface_normals[qp]*dphi_face[i][qp]) - phi_face[i][qp]*(qface_normals[qp]*dphi_neighbor_face[j][qp])); // consistency
Ken(i,j) -= JxW_face[qp] * penalty/h_elem * phi_face[i][qp]*phi_neighbor_face[j][qp]; // stability
}
}
}
The element and neighbor boundary matrix are now built
for this side. Add them to the global matrix
The \p SparseMatrix::add_matrix() members do this for us.
ellipticdg_system.matrix->add_matrix(Kne,neighbor_dof_indices,dof_indices);
ellipticdg_system.matrix->add_matrix(Ken,dof_indices,neighbor_dof_indices);
ellipticdg_system.matrix->add_matrix(Kee,dof_indices);
ellipticdg_system.matrix->add_matrix(Knn,neighbor_dof_indices);
}
}
}
The element interior matrix and right-hand-side are now built
for this element. Add them to the global matrix and
right-hand-side vector. The \p SparseMatrix::add_matrix()
and \p NumericVector::add_vector() members do this for us.
ellipticdg_system.matrix->add_matrix(Ke, dof_indices);
ellipticdg_system.rhs->add_vector(Fe, dof_indices);
}
std::cout << "done" << std::endl;
}
int main (int argc, char** argv)
{
LibMeshInit init(argc, argv);
Parse the input file
GetPot input_file("miscellaneous_ex5.in");
Read in parameters from the input file
const unsigned int adaptive_refinement_steps = input_file("max_adaptive_r_steps",3);
const unsigned int uniform_refinement_steps = input_file("uniform_h_r_steps",3);
const Real refine_fraction = input_file("refine_fraction",0.5);
const Real coarsen_fraction = input_file("coarsen_fraction",0.);
const unsigned int max_h_level = input_file("max_h_level", 10);
const std::string refinement_type = input_file("refinement_type","p");
Order p_order = static_cast<Order>(input_file("p_order",1));
const std::string element_type = input_file("element_type", "tensor");
const Real penalty = input_file("ip_penalty", 10.);
const bool singularity = input_file("singularity", true);
const unsigned int dim = input_file("dimension", 3);
Skip higher-dimensional examples on a lower-dimensional libMesh build
libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");
Skip adaptive examples on a non-adaptive libMesh build
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_assert(false, "--enable-amr");
#else
Create or read the mesh
Mesh mesh;
if (dim == 1)
MeshTools::Generation::build_line(mesh,1,-1.,0.);
else if (dim == 2)
mesh.read("lshaped.xda");
else
mesh.read ("lshaped3D.xda");
Use triangles if the config file says so
if (element_type == "simplex")
MeshTools::Modification::all_tri(mesh);
Mesh Refinement object
MeshRefinement mesh_refinement(mesh);
mesh_refinement.refine_fraction() = refine_fraction;
mesh_refinement.coarsen_fraction() = coarsen_fraction;
mesh_refinement.max_h_level() = max_h_level;
Do uniform refinement
for (unsigned int rstep=0; rstep<uniform_refinement_steps; rstep++)
{
mesh_refinement.uniformly_refine(1);
}
Crate an equation system object
EquationSystems equation_system (mesh);
Set parameters for the equation system and the solver
equation_system.parameters.set<Real>("linear solver tolerance") = TOLERANCE * TOLERANCE;
equation_system.parameters.set<unsigned int>("linear solver maximum iterations") = 1000;
equation_system.parameters.set<Real>("penalty") = penalty;
equation_system.parameters.set<bool>("singularity") = singularity;
equation_system.parameters.set<std::string>("refinement") = refinement_type;
Create a system named ellipticdg
LinearImplicitSystem& ellipticdg_system = equation_system.add_system<LinearImplicitSystem> ("EllipticDG");
Add a variable "u" to "ellipticdg" using the p_order specified in the config file
if( on_command_line( "element_type" ) )
{
std::string fe_str = command_line_value( std::string("element_type"),
std::string("MONOMIAL") );
if( fe_str != "MONOMIAL" || fe_str != "XYZ" )
{
std::cerr << "Error: This example must be run with MONOMIAL or XYZ element types." << std::endl;
libmesh_error();
}
ellipticdg_system.add_variable ("u", p_order, Utility::string_to_enum<FEFamily>(fe_str) );
}
else
ellipticdg_system.add_variable ("u", p_order, MONOMIAL);
Give the system a pointer to the matrix assembly function
ellipticdg_system.attach_assemble_function (assemble_ellipticdg);
Initialize the data structures for the equation system
equation_system.init();
Construct ExactSolution object and attach solution functions
ExactSolution exact_sol(equation_system);
exact_sol.attach_exact_value(exact_solution);
exact_sol.attach_exact_deriv(exact_derivative);
A refinement loop.
for (unsigned int rstep=0; rstep<adaptive_refinement_steps; ++rstep)
{
std::cout << " Beginning Solve " << rstep << std::endl;
std::cout << "Number of elements: " << mesh.n_elem() << std::endl;
Solve the system
ellipticdg_system.solve();
std::cout << "System has: " << equation_system.n_active_dofs()
<< " degrees of freedom."
<< std::endl;
std::cout << "Linear solver converged at step: "
<< ellipticdg_system.n_linear_iterations()
<< ", final residual: "
<< ellipticdg_system.final_linear_residual()
<< std::endl;
Compute the error
exact_sol.compute_error("EllipticDG", "u");
Print out the error values
std::cout << "L2-Error is: "
<< exact_sol.l2_error("EllipticDG", "u")
<< std::endl;
Possibly refine the mesh
if (rstep+1 < adaptive_refinement_steps)
{
The ErrorVector is a particular StatisticsVector
for computing error information on a finite element mesh.
ErrorVector error;
The discontinuity error estimator
evaluate the jump of the solution
on elements faces
DiscontinuityMeasure error_estimator;
error_estimator.estimate_error(ellipticdg_system,error);
Take the error in error and decide which elements will be coarsened or refined
mesh_refinement.flag_elements_by_error_fraction(error);
if (refinement_type == "p")
mesh_refinement.switch_h_to_p_refinement();
if (refinement_type == "hp")
mesh_refinement.add_p_to_h_refinement();
Refine and coarsen the flagged elements
mesh_refinement.refine_and_coarsen_elements();
equation_system.reinit();
}
}
Write out the solution
After solving the system write the solution
to a ExodusII-formatted plot file.
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO (mesh).write_discontinuous_exodusII("lshaped_dg.e",equation_system);
#endif
#endif // #ifndef LIBMESH_ENABLE_AMR
All done.
return 0;
}
The source file miscellaneous_ex5.C without comments:
#include <iostream>
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/equation_systems.h"
#include "libmesh/mesh_data.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_modification.h"
#include "libmesh/elem.h"
#include "libmesh/transient_system.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/dof_map.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/dense_submatrix.h"
#include "libmesh/dense_subvector.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/fe_interface.h"
#include "libmesh/getpot.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/error_vector.h"
#include "libmesh/kelly_error_estimator.h"
#include "libmesh/discontinuity_measure.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/exact_solution.h"
using namespace libMesh;
Number exact_solution (const Point& p, const Parameters& parameters, const std::string&, const std::string&)
{
const Real x = p(0);
const Real y = p(1);
const Real z = p(2);
if (parameters.get<bool>("singularity"))
{
Real theta = atan2(y,x);
if (theta < 0)
theta += 2. * libMesh::pi;
return pow(x*x + y*y, 1./3.)*sin(2./3.*theta) + z;
}
else
{
return cos(x) * exp(y) * (1. - z);
}
}
Gradient exact_derivative(const Point& p,
const Parameters& parameters, // es parameters
const std::string&, // sys_name, not needed
const std::string&) // unk_name, not needed
{
Gradient gradu;
const Real x = p(0);
const Real y = p(1);
const Real z = p(2);
if (parameters.get<bool>("singularity"))
{
libmesh_assert_not_equal_to (x, 0.);
const Real tt = 2./3.;
const Real ot = 1./3.;
const Real r2 = x*x + y*y;
Real theta = atan2(y,x);
if (theta < 0)
theta += 2. * libMesh::pi;
gradu(0) = tt*x*pow(r2,-tt)*sin(tt*theta) - pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*y/x/x;
gradu(1) = tt*y*pow(r2,-tt)*sin(tt*theta) + pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*1./x;
gradu(2) = 1.;
}
else
{
gradu(0) = -sin(x) * exp(y) * (1. - z);
gradu(1) = cos(x) * exp(y) * (1. - z);
gradu(2) = -cos(x) * exp(y);
}
return gradu;
}
void assemble_ellipticdg(EquationSystems& es, const std::string& system_name)
{
std::cout<<" assembling elliptic dg system... ";
std::cout.flush();
libmesh_assert_equal_to (system_name, "EllipticDG");
const MeshBase& mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
LinearImplicitSystem & ellipticdg_system = es.get_system<LinearImplicitSystem> ("EllipticDG");
const Real penalty = es.parameters.get<Real> ("penalty");
std::string refinement_type = es.parameters.get<std::string> ("refinement");
const DofMap & dof_map = ellipticdg_system.get_dof_map();
FEType fe_type = ellipticdg_system.variable_type(0);
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
AutoPtr<FEBase> fe_elem_face(FEBase::build(dim, fe_type));
AutoPtr<FEBase> fe_neighbor_face(FEBase::build(dim, fe_type));
#ifdef QORDER
QGauss qrule (dim, QORDER);
#else
QGauss qrule (dim, fe_type.default_quadrature_order());
#endif
fe->attach_quadrature_rule (&qrule);
#ifdef QORDER
QGauss qface(dim-1, QORDER);
#else
QGauss qface(dim-1, fe_type.default_quadrature_order());
#endif
fe_elem_face->attach_quadrature_rule(&qface);
fe_neighbor_face->attach_quadrature_rule(&qface);
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
const std::vector<std::vector<Real> >& phi_face = fe_elem_face->get_phi();
const std::vector<std::vector<RealGradient> >& dphi_face = fe_elem_face->get_dphi();
const std::vector<Real>& JxW_face = fe_elem_face->get_JxW();
const std::vector<Point>& qface_normals = fe_elem_face->get_normals();
const std::vector<Point>& qface_points = fe_elem_face->get_xyz();
const std::vector<std::vector<Real> >& phi_neighbor_face = fe_neighbor_face->get_phi();
const std::vector<std::vector<RealGradient> >& dphi_neighbor_face = fe_neighbor_face->get_dphi();
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
DenseMatrix<Number> Kne;
DenseMatrix<Number> Ken;
DenseMatrix<Number> Kee;
DenseMatrix<Number> Knn;
std::vector<dof_id_type> dof_indices;
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);
const unsigned int n_dofs = dof_indices.size();
fe->reinit (elem);
Ke.resize (n_dofs, n_dofs);
Fe.resize (n_dofs);
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
for (unsigned int i=0; i<n_dofs; i++)
{
for (unsigned int j=0; j<n_dofs; j++)
{
Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
}
}
}
for (unsigned int side=0; side<elem->n_sides(); side++)
{
if (elem->neighbor(side) == NULL)
{
fe_elem_face->reinit(elem, side);
AutoPtr<Elem> elem_side (elem->build_side(side));
const unsigned int elem_b_order = static_cast<unsigned int> (fe_elem_face->get_order());
const double h_elem = elem->volume()/elem_side->volume() * 1./pow(elem_b_order, 2.);
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
Number bc_value = exact_solution(qface_points[qp], es.parameters,"null","void");
for (unsigned int i=0; i<n_dofs; i++)
{
for (unsigned int j=0; j<n_dofs; j++)
{
Ke(i,j) += JxW_face[qp] * penalty/h_elem * phi_face[i][qp] * phi_face[j][qp]; // stability
Ke(i,j) -= JxW_face[qp] * (phi_face[i][qp] * (dphi_face[j][qp]*qface_normals[qp]) + phi_face[j][qp] * (dphi_face[i][qp]*qface_normals[qp])); // consistency
}
Fe(i) += JxW_face[qp] * bc_value * penalty/h_elem * phi_face[i][qp]; // stability
Fe(i) -= JxW_face[qp] * dphi_face[i][qp] * (bc_value*qface_normals[qp]); // consistency
}
}
}
else
{
const Elem* neighbor = elem->neighbor(side);
const unsigned int elem_id = elem->id();
const unsigned int neighbor_id = neighbor->id();
if ((neighbor->active() && (neighbor->level() == elem->level()) && (elem_id < neighbor_id)) || (neighbor->level() < elem->level()))
{
AutoPtr<Elem> elem_side (elem->build_side(side));
const unsigned int elem_b_order = static_cast<unsigned int>(fe_elem_face->get_order());
const unsigned int neighbor_b_order = static_cast<unsigned int>(fe_neighbor_face->get_order());
const double side_order = (elem_b_order + neighbor_b_order)/2.;
const double h_elem = (elem->volume()/elem_side->volume()) * 1./pow(side_order,2.);
std::vector<Point> qface_neighbor_point;
std::vector<Point > qface_point;
fe_elem_face->reinit(elem, side);
qface_point = fe_elem_face->get_xyz();
unsigned int side_neighbor = neighbor->which_neighbor_am_i(elem);
if (refinement_type == "p")
fe_neighbor_face->side_map (neighbor, elem_side.get(), side_neighbor, qface.get_points(), qface_neighbor_point);
else
FEInterface::inverse_map (elem->dim(), fe->get_fe_type(), neighbor, qface_point, qface_neighbor_point);
fe_neighbor_face->reinit(neighbor, &qface_neighbor_point);
std::vector<dof_id_type> neighbor_dof_indices;
dof_map.dof_indices (neighbor, neighbor_dof_indices);
const unsigned int n_neighbor_dofs = neighbor_dof_indices.size();
Kne.resize (n_neighbor_dofs, n_dofs);
Ken.resize (n_dofs, n_neighbor_dofs);
Kee.resize (n_dofs, n_dofs);
Knn.resize (n_neighbor_dofs, n_neighbor_dofs);
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
for (unsigned int i=0; i<n_dofs; i++)
{
for (unsigned int j=0; j<n_dofs; j++)
{
Kee(i,j) -= 0.5 * JxW_face[qp] * (phi_face[j][qp]*(qface_normals[qp]*dphi_face[i][qp]) + phi_face[i][qp]*(qface_normals[qp]*dphi_face[j][qp])); // consistency
Kee(i,j) += JxW_face[qp] * penalty/h_elem * phi_face[j][qp]*phi_face[i][qp]; // stability
}
}
for (unsigned int i=0; i<n_neighbor_dofs; i++)
{
for (unsigned int j=0; j<n_neighbor_dofs; j++)
{
Knn(i,j) += 0.5 * JxW_face[qp] * (phi_neighbor_face[j][qp]*(qface_normals[qp]*dphi_neighbor_face[i][qp]) + phi_neighbor_face[i][qp]*(qface_normals[qp]*dphi_neighbor_face[j][qp])); // consistency
Knn(i,j) += JxW_face[qp] * penalty/h_elem * phi_neighbor_face[j][qp]*phi_neighbor_face[i][qp]; // stability
}
}
for (unsigned int i=0; i<n_neighbor_dofs; i++)
{
for (unsigned int j=0; j<n_dofs; j++)
{
Kne(i,j) += 0.5 * JxW_face[qp] * (phi_neighbor_face[i][qp]*(qface_normals[qp]*dphi_face[j][qp]) - phi_face[j][qp]*(qface_normals[qp]*dphi_neighbor_face[i][qp])); // consistency
Kne(i,j) -= JxW_face[qp] * penalty/h_elem * phi_face[j][qp]*phi_neighbor_face[i][qp]; // stability
}
}
for (unsigned int i=0; i<n_dofs; i++)
{
for (unsigned int j=0; j<n_neighbor_dofs; j++)
{
Ken(i,j) += 0.5 * JxW_face[qp] * (phi_neighbor_face[j][qp]*(qface_normals[qp]*dphi_face[i][qp]) - phi_face[i][qp]*(qface_normals[qp]*dphi_neighbor_face[j][qp])); // consistency
Ken(i,j) -= JxW_face[qp] * penalty/h_elem * phi_face[i][qp]*phi_neighbor_face[j][qp]; // stability
}
}
}
ellipticdg_system.matrix->add_matrix(Kne,neighbor_dof_indices,dof_indices);
ellipticdg_system.matrix->add_matrix(Ken,dof_indices,neighbor_dof_indices);
ellipticdg_system.matrix->add_matrix(Kee,dof_indices);
ellipticdg_system.matrix->add_matrix(Knn,neighbor_dof_indices);
}
}
}
ellipticdg_system.matrix->add_matrix(Ke, dof_indices);
ellipticdg_system.rhs->add_vector(Fe, dof_indices);
}
std::cout << "done" << std::endl;
}
int main (int argc, char** argv)
{
LibMeshInit init(argc, argv);
GetPot input_file("miscellaneous_ex5.in");
const unsigned int adaptive_refinement_steps = input_file("max_adaptive_r_steps",3);
const unsigned int uniform_refinement_steps = input_file("uniform_h_r_steps",3);
const Real refine_fraction = input_file("refine_fraction",0.5);
const Real coarsen_fraction = input_file("coarsen_fraction",0.);
const unsigned int max_h_level = input_file("max_h_level", 10);
const std::string refinement_type = input_file("refinement_type","p");
Order p_order = static_cast<Order>(input_file("p_order",1));
const std::string element_type = input_file("element_type", "tensor");
const Real penalty = input_file("ip_penalty", 10.);
const bool singularity = input_file("singularity", true);
const unsigned int dim = input_file("dimension", 3);
libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_assert(false, "--enable-amr");
#else
Mesh mesh;
if (dim == 1)
MeshTools::Generation::build_line(mesh,1,-1.,0.);
else if (dim == 2)
mesh.read("lshaped.xda");
else
mesh.read ("lshaped3D.xda");
if (element_type == "simplex")
MeshTools::Modification::all_tri(mesh);
MeshRefinement mesh_refinement(mesh);
mesh_refinement.refine_fraction() = refine_fraction;
mesh_refinement.coarsen_fraction() = coarsen_fraction;
mesh_refinement.max_h_level() = max_h_level;
for (unsigned int rstep=0; rstep<uniform_refinement_steps; rstep++)
{
mesh_refinement.uniformly_refine(1);
}
EquationSystems equation_system (mesh);
equation_system.parameters.set<Real>("linear solver tolerance") = TOLERANCE * TOLERANCE;
equation_system.parameters.set<unsigned int>("linear solver maximum iterations") = 1000;
equation_system.parameters.set<Real>("penalty") = penalty;
equation_system.parameters.set<bool>("singularity") = singularity;
equation_system.parameters.set<std::string>("refinement") = refinement_type;
LinearImplicitSystem& ellipticdg_system = equation_system.add_system<LinearImplicitSystem> ("EllipticDG");
if( on_command_line( "element_type" ) )
{
std::string fe_str = command_line_value( std::string("element_type"),
std::string("MONOMIAL") );
if( fe_str != "MONOMIAL" || fe_str != "XYZ" )
{
std::cerr << "Error: This example must be run with MONOMIAL or XYZ element types." << std::endl;
libmesh_error();
}
ellipticdg_system.add_variable ("u", p_order, Utility::string_to_enum<FEFamily>(fe_str) );
}
else
ellipticdg_system.add_variable ("u", p_order, MONOMIAL);
ellipticdg_system.attach_assemble_function (assemble_ellipticdg);
equation_system.init();
ExactSolution exact_sol(equation_system);
exact_sol.attach_exact_value(exact_solution);
exact_sol.attach_exact_deriv(exact_derivative);
for (unsigned int rstep=0; rstep<adaptive_refinement_steps; ++rstep)
{
std::cout << " Beginning Solve " << rstep << std::endl;
std::cout << "Number of elements: " << mesh.n_elem() << std::endl;
ellipticdg_system.solve();
std::cout << "System has: " << equation_system.n_active_dofs()
<< " degrees of freedom."
<< std::endl;
std::cout << "Linear solver converged at step: "
<< ellipticdg_system.n_linear_iterations()
<< ", final residual: "
<< ellipticdg_system.final_linear_residual()
<< std::endl;
exact_sol.compute_error("EllipticDG", "u");
std::cout << "L2-Error is: "
<< exact_sol.l2_error("EllipticDG", "u")
<< std::endl;
if (rstep+1 < adaptive_refinement_steps)
{
ErrorVector error;
DiscontinuityMeasure error_estimator;
error_estimator.estimate_error(ellipticdg_system,error);
mesh_refinement.flag_elements_by_error_fraction(error);
if (refinement_type == "p")
mesh_refinement.switch_h_to_p_refinement();
if (refinement_type == "hp")
mesh_refinement.add_p_to_h_refinement();
mesh_refinement.refine_and_coarsen_elements();
equation_system.reinit();
}
}
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO (mesh).write_discontinuous_exodusII("lshaped_dg.e",equation_system);
#endif
#endif // #ifndef LIBMESH_ENABLE_AMR
return 0;
}
The console output of the program:
***************************************************************
* Running Example miscellaneous_ex5:
* mpirun -np 2 example-devel -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
Beginning Solve 0
Number of elements: 219
assembling elliptic dg system... done
System has: 768 degrees of freedom.
Linear solver converged at step: 25, final residual: 9.49449e-12
L2-Error is: 0.00666744
Beginning Solve 1
Number of elements: 827
assembling elliptic dg system... done
System has: 2896 degrees of freedom.
Linear solver converged at step: 32, final residual: 1.90449e-11
L2-Error is: 0.00264921
Beginning Solve 2
Number of elements: 3003
assembling elliptic dg system... done
System has: 10512 degrees of freedom.
Linear solver converged at step: 46, final residual: 3.57665e-11
L2-Error is: 0.0016323
************************************************************************************************************************
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***
************************************************************************************************************************
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------
/workspace/libmesh/examples/miscellaneous/miscellaneous_ex5/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:41:23 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): 5.963e+00 1.00001 5.963e+00
Objects: 2.020e+02 1.00000 2.020e+02
Flops: 8.410e+08 1.08167 8.093e+08 1.619e+09
Flops/sec: 1.410e+08 1.08169 1.357e+08 2.714e+08
MPI Messages: 1.405e+02 1.00000 1.405e+02 2.810e+02
MPI Message Lengths: 4.403e+05 1.00000 3.134e+03 8.807e+05
MPI Reductions: 3.990e+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: 5.9633e+00 100.0% 1.6186e+09 100.0% 2.810e+02 100.0% 3.134e+03 100.0% 3.980e+02 99.7%
------------------------------------------------------------------------------------------------------------------------
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
VecMDot 103 1.0 7.6683e-03 1.5 7.92e+06 1.0 0.0e+00 0.0e+00 1.0e+02 0 1 0 0 26 0 1 0 0 26 2066
VecNorm 111 1.0 9.8920e-04 1.4 6.37e+05 1.0 0.0e+00 0.0e+00 1.1e+02 0 0 0 0 28 0 0 0 0 28 1288
VecScale 108 1.0 1.6475e-04 1.0 3.12e+05 1.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 3782
VecCopy 14 1.0 3.2425e-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 139 1.0 2.4986e-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
VecAXPY 10 1.0 6.2091e-03 1.0 5.52e+04 1.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 18
VecMAXPY 108 1.0 2.3346e-03 1.0 8.52e+06 1.0 0.0e+00 0.0e+00 0.0e+00 0 1 0 0 0 0 1 0 0 0 7297
VecAssemblyBegin 23 1.0 1.7472e-0289.3 0.00e+00 0.0 0.0e+00 0.0e+00 5.7e+01 0 0 0 0 14 0 0 0 0 14 0
VecAssemblyEnd 23 1.0 2.8610e-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 117 1.0 4.2748e-04 1.0 0.00e+00 0.0 2.3e+02 2.9e+03 0.0e+00 0 0 83 78 0 0 0 83 78 0 0
VecScatterEnd 117 1.0 6.2580e-02 2.7 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
VecNormalize 108 1.0 1.1902e-03 1.3 9.35e+05 1.0 0.0e+00 0.0e+00 1.1e+02 0 0 0 0 27 0 0 0 0 27 1570
MatMult 108 1.0 7.1138e-02 2.2 1.68e+07 1.0 2.2e+02 2.9e+03 0.0e+00 1 2 77 71 0 1 2 77 71 0 471
MatSolve 111 1.0 1.1712e-01 1.0 2.15e+08 1.0 0.0e+00 0.0e+00 0.0e+00 2 26 0 0 0 2 26 0 0 0 3635
MatLUFactorNum 3 1.0 1.9490e-01 1.1 5.92e+08 1.1 0.0e+00 0.0e+00 0.0e+00 3 69 0 0 0 3 69 0 0 0 5770
MatILUFactorSym 3 1.0 3.4527e-01 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 9.0e+00 6 0 0 0 2 6 0 0 0 2 0
MatAssemblyBegin 6 1.0 3.4488e-02110.5 0.00e+00 0.0 1.5e+01 1.1e+04 1.2e+01 0 0 5 19 3 0 0 5 19 3 0
MatAssemblyEnd 6 1.0 1.2739e-03 1.1 0.00e+00 0.0 1.2e+01 6.3e+02 2.4e+01 0 0 4 1 6 0 0 4 1 6 0
MatGetRowIJ 3 1.0 3.3379e-06 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
MatGetOrdering 3 1.0 2.0790e-04 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 1.2e+01 0 0 0 0 3 0 0 0 0 3 0
MatZeroEntries 9 1.0 4.9257e-04 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
KSPGMRESOrthog 103 1.0 9.9757e-03 1.4 1.58e+07 1.0 0.0e+00 0.0e+00 1.0e+02 0 2 0 0 26 0 2 0 0 26 3177
KSPSetUp 6 1.0 1.2803e-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
KSPSolve 3 1.0 7.0725e-01 1.0 8.41e+08 1.1 2.2e+02 2.9e+03 2.4e+02 12100 77 71 60 12100 77 71 61 2289
PCSetUp 6 1.0 5.4099e-01 1.1 5.92e+08 1.1 0.0e+00 0.0e+00 2.7e+01 9 69 0 0 7 9 69 0 0 7 2079
PCSetUpOnBlocks 3 1.0 5.4058e-01 1.1 5.92e+08 1.1 0.0e+00 0.0e+00 2.1e+01 9 69 0 0 5 9 69 0 0 5 2080
PCApply 111 1.0 1.1837e-01 1.0 2.15e+08 1.0 0.0e+00 0.0e+00 0.0e+00 2 26 0 0 0 2 26 0 0 0 3596
------------------------------------------------------------------------------------------------------------------------
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 137 137 2580728 0
Vector Scatter 8 8 8288 0
Index Set 29 29 93912 0
IS L to G Mapping 3 3 1692 0
Matrix 12 12 30013852 0
Krylov Solver 6 6 58080 0
Preconditioner 6 6 5352 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 0
Average time for MPI_Barrier(): 1.38283e-06
Average time for zero size MPI_Send(): 1.84774e-05
#PETSc Option Table entries:
-ksp_right_pc
-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:41:23 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=5.9737, Active time=5.7766 |
---------------------------------------------------------------------------------------------------------------------
| 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() 3 0.0355 0.011821 0.0539 0.017951 0.61 0.93 |
| build_sparsity() 3 0.2930 0.097666 1.5612 0.520393 5.07 27.03 |
| create_dof_constraints() 3 0.0485 0.016166 0.0795 0.026511 0.84 1.38 |
| distribute_dofs() 3 0.0408 0.013585 0.1042 0.034731 0.71 1.80 |
| dof_indices() 58964 1.7292 0.000029 1.7292 0.000029 29.93 29.93 |
| old_dof_indices() 3512 0.1004 0.000029 0.1004 0.000029 1.74 1.74 |
| prepare_send_list() 3 0.0006 0.000192 0.0006 0.000192 0.01 0.01 |
| reinit() 3 0.0625 0.020839 0.0625 0.020839 1.08 1.08 |
| |
| EquationSystems |
| build_discontinuous_solution_vector() 1 0.0188 0.018767 0.0957 0.095707 0.32 1.66 |
| |
| ExodusII_IO |
| write_nodal_data_discontinuous() 1 0.0102 0.010223 0.0102 0.010223 0.18 0.18 |
| |
| FE |
| compute_shape_functions() 20120 0.2358 0.000012 0.2358 0.000012 4.08 4.08 |
| init_shape_functions() 16582 0.1045 0.000006 0.1045 0.000006 1.81 1.81 |
| inverse_map() 43428 0.4920 0.000011 0.4920 0.000011 8.52 8.52 |
| |
| FEMap |
| compute_affine_map() 20120 0.1891 0.000009 0.1891 0.000009 3.27 3.27 |
| compute_face_map() 7191 0.0757 0.000011 0.0757 0.000011 1.31 1.31 |
| init_face_shape_functions() 5 0.0001 0.000025 0.0001 0.000025 0.00 0.00 |
| init_reference_to_physical_map() 16582 0.5164 0.000031 0.5164 0.000031 8.94 8.94 |
| |
| JumpErrorEstimator |
| estimate_error() 2 0.0662 0.033098 0.3373 0.168633 1.15 5.84 |
| |
| LocationMap |
| find() 22344 0.0486 0.000002 0.0486 0.000002 0.84 0.84 |
| init() 6 0.0053 0.000885 0.0053 0.000885 0.09 0.09 |
| |
| Mesh |
| contract() 2 0.0010 0.000485 0.0028 0.001389 0.02 0.05 |
| find_neighbors() 5 0.0815 0.016298 0.0825 0.016493 1.41 1.43 |
| renumber_nodes_and_elem() 12 0.0050 0.000420 0.0050 0.000420 0.09 0.09 |
| |
| MeshCommunication |
| compute_hilbert_indices() 5 0.0155 0.003103 0.0155 0.003103 0.27 0.27 |
| find_global_indices() 5 0.0060 0.001202 0.0243 0.004860 0.10 0.42 |
| parallel_sort() 5 0.0016 0.000322 0.0020 0.000401 0.03 0.03 |
| |
| MeshRefinement |
| _coarsen_elements() 4 0.0010 0.000244 0.0010 0.000251 0.02 0.02 |
| _refine_elements() 6 0.0714 0.011908 0.1948 0.032462 1.24 3.37 |
| add_point() 22344 0.0566 0.000003 0.1093 0.000005 0.98 1.89 |
| make_coarsening_compatible() 8 0.0186 0.002324 0.0186 0.002324 0.32 0.32 |
| make_refinement_compatible() 8 0.0010 0.000127 0.0011 0.000132 0.02 0.02 |
| |
| MetisPartitioner |
| partition() 5 0.0761 0.015213 0.1009 0.020177 1.32 1.75 |
| |
| Parallel |
| allgather() 19 0.0005 0.000025 0.0005 0.000028 0.01 0.01 |
| broadcast() 20 0.0003 0.000015 0.0002 0.000012 0.01 0.00 |
| max(bool) 23 0.0105 0.000459 0.0105 0.000459 0.18 0.18 |
| max(scalar) 494 0.0106 0.000021 0.0106 0.000021 0.18 0.18 |
| max(vector) 119 0.0008 0.000007 0.0019 0.000016 0.01 0.03 |
| min(bool) 611 0.0017 0.000003 0.0017 0.000003 0.03 0.03 |
| min(scalar) 483 0.0435 0.000090 0.0435 0.000090 0.75 0.75 |
| min(vector) 119 0.0009 0.000007 0.0021 0.000017 0.02 0.04 |
| probe() 40 0.0002 0.000005 0.0002 0.000005 0.00 0.00 |
| receive() 40 0.0003 0.000007 0.0005 0.000012 0.00 0.01 |
| send() 40 0.0001 0.000003 0.0001 0.000003 0.00 0.00 |
| send_receive() 50 0.0004 0.000007 0.0011 0.000021 0.01 0.02 |
| sum() 25 0.0004 0.000017 0.0005 0.000021 0.01 0.01 |
| |
| Parallel::Request |
| wait() 40 0.0001 0.000002 0.0001 0.000002 0.00 0.00 |
| |
| Partitioner |
| set_node_processor_ids() 6 0.0057 0.000953 0.0077 0.001276 0.10 0.13 |
| set_parent_processor_ids() 5 0.0058 0.001169 0.0058 0.001169 0.10 0.10 |
| |
| PetscLinearSolver |
| solve() 3 0.7440 0.247995 0.7440 0.247995 12.88 12.88 |
| |
| ProjectVector |
| operator() 2 0.0630 0.031521 0.5182 0.259119 1.09 8.97 |
| |
| System |
| assemble() 3 0.4545 0.151510 1.6104 0.536811 7.87 27.88 |
| project_vector() 2 0.0239 0.011971 0.5922 0.296083 0.41 10.25 |
| |
| XdrIO |
| read() 1 0.0007 0.000654 0.0007 0.000721 0.01 0.01 |
---------------------------------------------------------------------------------------------------------------------
| Totals: 233430 5.7766 100.00 |
---------------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example miscellaneous_ex5:
* mpirun -np 2 example-devel -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
Miscellaneous Example 5 - Interior Penalty Discontinuous Galerkin
By Lorenzo Botti
This example is based on Adaptivity Example 3, but uses an Interior Penalty Discontinuous Galerkin formulation.