The source file adaptivity_ex3.C with comments:
#include "libmesh/mesh.h"
#include "libmesh/equation_systems.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/tecplot_io.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/error_vector.h"
#include "libmesh/exact_error_estimator.h"
#include "libmesh/kelly_error_estimator.h"
#include "libmesh/patch_recovery_error_estimator.h"
#include "libmesh/uniform_refinement_estimator.h"
#include "libmesh/hp_coarsentest.h"
#include "libmesh/hp_singular.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_modification.h"
#include "libmesh/perf_log.h"
#include "libmesh/getpot.h"
#include "libmesh/exact_solution.h"
#include "libmesh/dof_map.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/elem.h"
#include "libmesh/string_to_enum.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
Function prototype. This is the function that will assemble
the linear system for our Laplace problem. Note that the
function will take the \p EquationSystems object and the
name of the system we are assembling as input. From the
\p EquationSystems object we have acess to the \p Mesh and
other objects we might need.
void assemble_laplace(EquationSystems& es,
const std::string& system_name);
Prototype for calculation of the exact solution. Useful
for setting boundary conditions.
Number exact_solution(const Point& p,
const Parameters&, // EquationSystem parameters, not needed
const std::string&, // sys_name, not needed
const std::string&); // unk_name, not needed);
Prototype for calculation of the gradient of the exact solution.
Gradient exact_derivative(const Point& p,
const Parameters&, // EquationSystems parameters, not needed
const std::string&, // sys_name, not needed
const std::string&); // unk_name, not needed);
These are non-const because the input file may change it,
It is global because our exact_* functions use it.
Set the dimensionality of the mesh
Set the dimensionality of the mesh
unsigned int dim = 2;
Choose whether or not to use the singular solution
bool singularity = true;
int main(int argc, char** argv)
{
Initialize libMesh.
LibMeshInit init (argc, argv);
Skip adaptive examples on a non-adaptive libMesh build
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_assert(false, "--enable-amr");
#else
Parse the input file
GetPot input_file("adaptivity_ex3.in");
Read in parameters from the input file
const unsigned int max_r_steps = input_file("max_r_steps", 3);
const unsigned int max_r_level = input_file("max_r_level", 3);
const Real refine_percentage = input_file("refine_percentage", 0.5);
const Real coarsen_percentage = input_file("coarsen_percentage", 0.5);
const unsigned int uniform_refine = input_file("uniform_refine",0);
const std::string refine_type = input_file("refinement_type", "h");
const std::string approx_type = input_file("approx_type", "LAGRANGE");
const unsigned int approx_order = input_file("approx_order", 1);
const std::string element_type = input_file("element_type", "tensor");
const int extra_error_quadrature = input_file("extra_error_quadrature", 0);
const int max_linear_iterations = input_file("max_linear_iterations", 5000);
const bool output_intermediate = input_file("output_intermediate", false);
dim = input_file("dimension", 2);
const std::string indicator_type = input_file("indicator_type", "kelly");
singularity = input_file("singularity", true);
Skip higher-dimensional examples on a lower-dimensional libMesh build
libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");
Output file for plotting the error as a function of
the number of degrees of freedom.
std::string approx_name = "";
if (element_type == "tensor")
approx_name += "bi";
if (approx_order == 1)
approx_name += "linear";
else if (approx_order == 2)
approx_name += "quadratic";
else if (approx_order == 3)
approx_name += "cubic";
else if (approx_order == 4)
approx_name += "quartic";
std::string output_file = approx_name;
output_file += "_";
output_file += refine_type;
if (uniform_refine == 0)
output_file += "_adaptive.m";
else
output_file += "_uniform.m";
std::ofstream out (output_file.c_str());
out << "% dofs L2-error H1-error" << std::endl;
out << "e = [" << std::endl;
Create a mesh.
Mesh mesh;
Read in the 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);
We used first order elements to describe the geometry,
but we may need second order elements to hold the degrees
of freedom
if (approx_order > 1 || refine_type != "h")
mesh.all_second_order();
Mesh Refinement object
MeshRefinement mesh_refinement(mesh);
mesh_refinement.refine_fraction() = refine_percentage;
mesh_refinement.coarsen_fraction() = coarsen_percentage;
mesh_refinement.max_h_level() = max_r_level;
Create an equation systems object.
EquationSystems equation_systems (mesh);
Declare the system and its variables.
Creates a system named "Laplace"
LinearImplicitSystem& system =
equation_systems.add_system<LinearImplicitSystem> ("Laplace");
Adds the variable "u" to "Laplace", using
the finite element type and order specified
in the config file
system.add_variable("u", static_cast<Order>(approx_order),
Utility::string_to_enum<FEFamily>(approx_type));
Give the system a pointer to the matrix assembly
function.
system.attach_assemble_function (assemble_laplace);
Initialize the data structures for the equation system.
equation_systems.init();
Set linear solver max iterations
equation_systems.parameters.set<unsigned int>("linear solver maximum iterations")
= max_linear_iterations;
Linear solver tolerance.
equation_systems.parameters.set<Real>("linear solver tolerance") =
std::pow(TOLERANCE, 2.5);
Prints information about the system to the screen.
equation_systems.print_info();
Construct ExactSolution object and attach solution functions
ExactSolution exact_sol(equation_systems);
exact_sol.attach_exact_value(exact_solution);
exact_sol.attach_exact_deriv(exact_derivative);
Use higher quadrature order for more accurate error results
exact_sol.extra_quadrature_order(extra_error_quadrature);
A refinement loop.
for (unsigned int r_step=0; r_step<max_r_steps; r_step++)
{
std::cout << "Beginning Solve " << r_step << std::endl;
Solve the system "Laplace", just like example 2.
system.solve();
std::cout << "System has: " << equation_systems.n_active_dofs()
<< " degrees of freedom."
<< std::endl;
std::cout << "Linear solver converged at step: "
<< system.n_linear_iterations()
<< ", final residual: "
<< system.final_linear_residual()
<< std::endl;
#ifdef LIBMESH_HAVE_EXODUS_API
After solving the system write the solution
to a ExodusII-formatted plot file.
if (output_intermediate)
{
std::ostringstream outfile;
outfile << "lshaped_" << r_step << ".e";
ExodusII_IO (mesh).write_equation_systems (outfile.str(),
equation_systems);
}
#endif // #ifdef LIBMESH_HAVE_EXODUS_API
Compute the error.
exact_sol.compute_error("Laplace", "u");
Print out the error values
std::cout << "L2-Error is: "
<< exact_sol.l2_error("Laplace", "u")
<< std::endl;
std::cout << "H1-Error is: "
<< exact_sol.h1_error("Laplace", "u")
<< std::endl;
Print to output file
out << equation_systems.n_active_dofs() << " "
<< exact_sol.l2_error("Laplace", "u") << " "
<< exact_sol.h1_error("Laplace", "u") << std::endl;
Possibly refine the mesh
if (r_step+1 != max_r_steps)
{
std::cout << " Refining the mesh..." << std::endl;
if (uniform_refine == 0)
{
The \p ErrorVector is a particular \p StatisticsVector
for computing error information on a finite element mesh.
ErrorVector error;
if (indicator_type == "exact")
{
The \p ErrorEstimator class interrogates a
finite element solution and assigns to each
element a positive error value.
This value is used for deciding which elements to
refine and which to coarsen.
For these simple test problems, we can use
numerical quadrature of the exact error between
the approximate and analytic solutions.
However, for real problems, we would need an error
indicator which only relies on the approximate
solution.
ExactErrorEstimator error_estimator;
error_estimator.attach_exact_value(exact_solution);
error_estimator.attach_exact_deriv(exact_derivative);
We optimize in H1 norm, the default
error_estimator.error_norm = H1;
Compute the error for each active element using the provided indicator. Note in general you will need to provide an error estimator specifically designed for your application.
Compute the error for each active element using the provided indicator. Note in general you will need to provide an error estimator specifically designed for your application.
error_estimator.estimate_error (system, error);
}
else if (indicator_type == "patch")
{
The patch recovery estimator should give a
good estimate of the solution interpolation
error.
PatchRecoveryErrorEstimator error_estimator;
error_estimator.estimate_error (system, error);
}
else if (indicator_type == "uniform")
{
Error indication based on uniform refinement
is reliable, but very expensive.
UniformRefinementEstimator error_estimator;
error_estimator.estimate_error (system, error);
}
else
{
libmesh_assert_equal_to (indicator_type, "kelly");
The Kelly error estimator is based on
an error bound for the Poisson problem
on linear elements, but is useful for
driving adaptive refinement in many problems
KellyErrorEstimator error_estimator;
error_estimator.estimate_error (system, error);
}
Write out the error distribution
std::ostringstream ss;
ss << r_step;
#ifdef LIBMESH_HAVE_EXODUS_API
std::string error_output = "error_"+ss.str()+".e";
#else
std::string error_output = "error_"+ss.str()+".gmv";
#endif
error.plot_error( error_output, mesh );
This takes the error in \p error and decides which elements
will be coarsened or refined. Any element within 20% of the
maximum error on any element will be refined, and any
element within 10% of the minimum error on any element might
be coarsened. Note that the elements flagged for refinement
will be refined, but those flagged for coarsening _might_ be
coarsened.
mesh_refinement.flag_elements_by_error_fraction (error);
If we are doing adaptive p refinement, we want
elements flagged for that instead.
if (refine_type == "p")
mesh_refinement.switch_h_to_p_refinement();
If we are doing "matched hp" refinement, we
flag elements for both h and p
if (refine_type == "matchedhp")
mesh_refinement.add_p_to_h_refinement();
If we are doing hp refinement, we
try switching some elements from h to p
if (refine_type == "hp")
{
HPCoarsenTest hpselector;
hpselector.select_refinement(system);
}
If we are doing "singular hp" refinement, we
try switching most elements from h to p
if (refine_type == "singularhp")
{
This only differs from p refinement for
the singular problem
libmesh_assert (singularity);
HPSingularity hpselector;
Our only singular point is at the origin
hpselector.singular_points.push_back(Point());
hpselector.select_refinement(system);
}
This call actually refines and coarsens the flagged
elements.
mesh_refinement.refine_and_coarsen_elements();
}
else if (uniform_refine == 1)
{
if (refine_type == "h" || refine_type == "hp" ||
refine_type == "matchedhp")
mesh_refinement.uniformly_refine(1);
if (refine_type == "p" || refine_type == "hp" ||
refine_type == "matchedhp")
mesh_refinement.uniformly_p_refine(1);
}
This call reinitializes the \p EquationSystems object for
the newly refined mesh. One of the steps in the
reinitialization is projecting the \p solution,
\p old_solution, etc... vectors from the old mesh to
the current one.
equation_systems.reinit ();
}
}
#ifdef LIBMESH_HAVE_EXODUS_API
Write out the solution
After solving the system write the solution
to a ExodusII-formatted plot file.
ExodusII_IO (mesh).write_equation_systems ("lshaped.e",
equation_systems);
#endif // #ifdef LIBMESH_HAVE_EXODUS_API
Close up the output file.
out << "];" << std::endl;
out << "hold on" << std::endl;
out << "plot(e(:,1), e(:,2), 'bo-');" << std::endl;
out << "plot(e(:,1), e(:,3), 'ro-');" << std::endl;
out << "set(gca,'XScale', 'Log');" << std::endl;
out << "set(gca,'YScale', 'Log');" << std::endl;
out << "xlabel('dofs');" << std::endl;
out << "title('" << approx_name << " elements');" << std::endl;
out << "legend('L2-error', 'H1-error');" << std::endl;
out << "disp('L2-error linear fit');" << std::endl;
out << "polyfit(log10(e(:,1)), log10(e(:,2)), 1)" << std::endl;
out << "disp('H1-error linear fit');" << std::endl;
out << "polyfit(log10(e(:,1)), log10(e(:,3)), 1)" << std::endl;
#endif // #ifndef LIBMESH_ENABLE_AMR
All done.
return 0;
}
We now define the exact solution, being careful
to obtain an angle from atan2 in the correct
quadrant.
Number exact_solution(const Point& p,
const Parameters&, // parameters, not needed
const std::string&, // sys_name, not needed
const std::string&) // unk_name, not needed
{
const Real x = p(0);
const Real y = (dim > 1) ? p(1) : 0.;
if (singularity)
{
The exact solution to the singular problem,
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;
Make the 3D solution similar
const Real z = (dim > 2) ? p(2) : 0;
return pow(x*x + y*y, 1./3.)*sin(2./3.*theta) + z;
}
else
{
The exact solution to a nonsingular problem,
good for testing ideal convergence rates
const Real z = (dim > 2) ? p(2) : 0;
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, not needed
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 = dim > 1 ? p(1) : 0.;
if (singularity)
{
We can't compute the gradient at x=0, it is not defined.
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;
du/dy
if (dim > 1)
gradu(1) = tt*y*pow(r2,-tt)*sin(tt*theta) + pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*1./x;
if (dim > 2)
gradu(2) = 1.;
}
else
{
const Real z = (dim > 2) ? p(2) : 0;
gradu(0) = -sin(x) * exp(y) * (1. - z);
if (dim > 1)
gradu(1) = cos(x) * exp(y) * (1. - z);
if (dim > 2)
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
matrices and right-hand sides, and then take into
account the boundary conditions, which will be handled
via a penalty method.
void assemble_laplace(EquationSystems& es,
const std::string& system_name)
{
#ifdef LIBMESH_ENABLE_AMR
It is a good idea to make sure we are assembling
the proper system.
libmesh_assert_equal_to (system_name, "Laplace");
Declare a performance log. Give it a descriptive
string to identify what part of the code we are
logging, since there may be many PerfLogs in an
application.
PerfLog perf_log ("Matrix Assembly",false);
Get a constant reference to the mesh object.
const MeshBase& mesh = es.get_mesh();
The dimension that we are running
const unsigned int mesh_dim = mesh.mesh_dimension();
Get a reference to the LinearImplicitSystem we are solving
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Laplace");
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
in future examples.
const DofMap& dof_map = 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 = dof_map.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(mesh_dim, fe_type));
AutoPtr<FEBase> fe_face (FEBase::build(mesh_dim, fe_type));
Quadrature rules for numerical integration.
AutoPtr<QBase> qrule(fe_type.default_quadrature_rule(mesh_dim));
AutoPtr<QBase> qface(fe_type.default_quadrature_rule(mesh_dim-1));
Tell the finite element object to use our quadrature rule.
fe->attach_quadrature_rule (qrule.get());
fe_face->attach_quadrature_rule (qface.get());
Here we define some references to cell-specific data that
will be used to assemble the linear system.
We begin with the element Jacobian * quadrature weight at each
integration point.
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<Real>& JxW_face = fe_face->get_JxW();
The physical XY locations of the quadrature points on the element.
These might be useful for evaluating spatially varying material
properties or forcing functions at the quadrature points.
const std::vector<Point>& q_point = fe->get_xyz();
The element shape functions evaluated at the quadrature points.
For this simple problem we usually only need them on element
boundaries.
const std::vector<std::vector<Real> >& phi = fe->get_phi();
const std::vector<std::vector<Real> >& psi = fe_face->get_phi();
The element shape function gradients evaluated at the quadrature
points.
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
The XY locations of the quadrature points used for face integration
const std::vector<Point>& qface_points = fe_face->get_xyz();
Define data structures to contain the element matrix
and right-hand-side vector contribution. Following
basic finite element terminology we will denote these
"Ke" and "Fe". More detail is in example 3.
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
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 the element matrix and right-hand-side contribution. See
example 3 for a discussion of the element iterators. Here we use
the \p const_active_local_elem_iterator to indicate we only want
to loop over elements that are assigned to the local processor
which are "active" in the sense of AMR. This allows each
processor to compute its components of the global matrix for
active elements while ignoring parent elements which have been
refined.
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)
{
Start logging the shape function initialization.
This is done through a simple function call with
the name of the event to log.
perf_log.push("elem init");
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);
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. Note that this will be the case if the
element type is different (i.e. the last element was a
triangle, now we are on a quadrilateral).
Ke.resize (dof_indices.size(),
dof_indices.size());
Fe.resize (dof_indices.size());
Stop logging the shape function initialization.
If you forget to stop logging an event the PerfLog
object will probably catch the error and abort.
perf_log.pop("elem init");
Now we will build the element matrix. This involves
a double loop to integrate the test funcions (i) against
the trial functions (j).
Now start logging the element matrix computation
Now start logging the element matrix computation
perf_log.push ("Ke");
for (unsigned int qp=0; qp<qrule->n_points(); qp++)
for (unsigned int i=0; i<dphi.size(); i++)
for (unsigned int j=0; j<dphi.size(); j++)
Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
We need a forcing function to make the 1D case interesting
if (mesh_dim == 1)
for (unsigned int qp=0; qp<qrule->n_points(); qp++)
{
Real x = q_point[qp](0);
Real f = singularity ? sqrt(3.)/9.*pow(-x, -4./3.) :
cos(x);
for (unsigned int i=0; i<dphi.size(); ++i)
Fe(i) += JxW[qp]*phi[i][qp]*f;
}
Stop logging the matrix computation
perf_log.pop ("Ke");
At this point the interior element integration has
been completed. However, we have not yet addressed
boundary conditions. For this example we will only
consider simple Dirichlet boundary conditions imposed
via the penalty method.
This approach adds the L2 projection of the boundary data in penalty form to the weak statement. This is a more generic approach for applying Dirichlet BCs which is applicable to non-Lagrange finite element discretizations.
This approach adds the L2 projection of the boundary data in penalty form to the weak statement. This is a more generic approach for applying Dirichlet BCs which is applicable to non-Lagrange finite element discretizations.
{
Start logging the boundary condition computation
perf_log.push ("BCs");
The penalty value.
const Real penalty = 1.e10;
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 s=0; s<elem->n_sides(); s++)
if (elem->neighbor(s) == NULL)
{
fe_face->reinit(elem,s);
for (unsigned int qp=0; qp<qface->n_points(); qp++)
{
const Number value = exact_solution (qface_points[qp],
es.parameters,
"null",
"void");
RHS contribution
for (unsigned int i=0; i<psi.size(); i++)
Fe(i) += penalty*JxW_face[qp]*value*psi[i][qp];
Matrix contribution
for (unsigned int i=0; i<psi.size(); i++)
for (unsigned int j=0; j<psi.size(); j++)
Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
}
}
Stop logging the boundary condition computation
perf_log.pop ("BCs");
}
The element 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.
Start logging the insertion of the local (element)
matrix and vector into the global matrix and vector
perf_log.push ("matrix insertion");
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);
Start logging the insertion of the local (element)
matrix and vector into the global matrix and vector
perf_log.pop ("matrix insertion");
}
That's it. We don't need to do anything else to the
PerfLog. When it goes out of scope (at this function return)
it will print its log to the screen. Pretty easy, huh?
#endif // #ifdef LIBMESH_ENABLE_AMR
}
The source file adaptivity_ex3.C without comments:
#include "libmesh/mesh.h"
#include "libmesh/equation_systems.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/tecplot_io.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/error_vector.h"
#include "libmesh/exact_error_estimator.h"
#include "libmesh/kelly_error_estimator.h"
#include "libmesh/patch_recovery_error_estimator.h"
#include "libmesh/uniform_refinement_estimator.h"
#include "libmesh/hp_coarsentest.h"
#include "libmesh/hp_singular.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_modification.h"
#include "libmesh/perf_log.h"
#include "libmesh/getpot.h"
#include "libmesh/exact_solution.h"
#include "libmesh/dof_map.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/elem.h"
#include "libmesh/string_to_enum.h"
using namespace libMesh;
void assemble_laplace(EquationSystems& es,
const std::string& system_name);
Number exact_solution(const Point& p,
const Parameters&, // EquationSystem parameters, not needed
const std::string&, // sys_name, not needed
const std::string&); // unk_name, not needed);
Gradient exact_derivative(const Point& p,
const Parameters&, // EquationSystems parameters, not needed
const std::string&, // sys_name, not needed
const std::string&); // unk_name, not needed);
unsigned int dim = 2;
bool singularity = true;
int main(int argc, char** argv)
{
LibMeshInit init (argc, argv);
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_assert(false, "--enable-amr");
#else
GetPot input_file("adaptivity_ex3.in");
const unsigned int max_r_steps = input_file("max_r_steps", 3);
const unsigned int max_r_level = input_file("max_r_level", 3);
const Real refine_percentage = input_file("refine_percentage", 0.5);
const Real coarsen_percentage = input_file("coarsen_percentage", 0.5);
const unsigned int uniform_refine = input_file("uniform_refine",0);
const std::string refine_type = input_file("refinement_type", "h");
const std::string approx_type = input_file("approx_type", "LAGRANGE");
const unsigned int approx_order = input_file("approx_order", 1);
const std::string element_type = input_file("element_type", "tensor");
const int extra_error_quadrature = input_file("extra_error_quadrature", 0);
const int max_linear_iterations = input_file("max_linear_iterations", 5000);
const bool output_intermediate = input_file("output_intermediate", false);
dim = input_file("dimension", 2);
const std::string indicator_type = input_file("indicator_type", "kelly");
singularity = input_file("singularity", true);
libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");
std::string approx_name = "";
if (element_type == "tensor")
approx_name += "bi";
if (approx_order == 1)
approx_name += "linear";
else if (approx_order == 2)
approx_name += "quadratic";
else if (approx_order == 3)
approx_name += "cubic";
else if (approx_order == 4)
approx_name += "quartic";
std::string output_file = approx_name;
output_file += "_";
output_file += refine_type;
if (uniform_refine == 0)
output_file += "_adaptive.m";
else
output_file += "_uniform.m";
std::ofstream out (output_file.c_str());
out << "% dofs L2-error H1-error" << std::endl;
out << "e = [" << std::endl;
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);
if (approx_order > 1 || refine_type != "h")
mesh.all_second_order();
MeshRefinement mesh_refinement(mesh);
mesh_refinement.refine_fraction() = refine_percentage;
mesh_refinement.coarsen_fraction() = coarsen_percentage;
mesh_refinement.max_h_level() = max_r_level;
EquationSystems equation_systems (mesh);
LinearImplicitSystem& system =
equation_systems.add_system<LinearImplicitSystem> ("Laplace");
system.add_variable("u", static_cast<Order>(approx_order),
Utility::string_to_enum<FEFamily>(approx_type));
system.attach_assemble_function (assemble_laplace);
equation_systems.init();
equation_systems.parameters.set<unsigned int>("linear solver maximum iterations")
= max_linear_iterations;
equation_systems.parameters.set<Real>("linear solver tolerance") =
std::pow(TOLERANCE, 2.5);
equation_systems.print_info();
ExactSolution exact_sol(equation_systems);
exact_sol.attach_exact_value(exact_solution);
exact_sol.attach_exact_deriv(exact_derivative);
exact_sol.extra_quadrature_order(extra_error_quadrature);
for (unsigned int r_step=0; r_step<max_r_steps; r_step++)
{
std::cout << "Beginning Solve " << r_step << std::endl;
system.solve();
std::cout << "System has: " << equation_systems.n_active_dofs()
<< " degrees of freedom."
<< std::endl;
std::cout << "Linear solver converged at step: "
<< system.n_linear_iterations()
<< ", final residual: "
<< system.final_linear_residual()
<< std::endl;
#ifdef LIBMESH_HAVE_EXODUS_API
if (output_intermediate)
{
std::ostringstream outfile;
outfile << "lshaped_" << r_step << ".e";
ExodusII_IO (mesh).write_equation_systems (outfile.str(),
equation_systems);
}
#endif // #ifdef LIBMESH_HAVE_EXODUS_API
exact_sol.compute_error("Laplace", "u");
std::cout << "L2-Error is: "
<< exact_sol.l2_error("Laplace", "u")
<< std::endl;
std::cout << "H1-Error is: "
<< exact_sol.h1_error("Laplace", "u")
<< std::endl;
out << equation_systems.n_active_dofs() << " "
<< exact_sol.l2_error("Laplace", "u") << " "
<< exact_sol.h1_error("Laplace", "u") << std::endl;
if (r_step+1 != max_r_steps)
{
std::cout << " Refining the mesh..." << std::endl;
if (uniform_refine == 0)
{
ErrorVector error;
if (indicator_type == "exact")
{
ExactErrorEstimator error_estimator;
error_estimator.attach_exact_value(exact_solution);
error_estimator.attach_exact_deriv(exact_derivative);
error_estimator.estimate_error (system, error);
}
else if (indicator_type == "patch")
{
PatchRecoveryErrorEstimator error_estimator;
error_estimator.estimate_error (system, error);
}
else if (indicator_type == "uniform")
{
UniformRefinementEstimator error_estimator;
error_estimator.estimate_error (system, error);
}
else
{
libmesh_assert_equal_to (indicator_type, "kelly");
KellyErrorEstimator error_estimator;
error_estimator.estimate_error (system, error);
}
std::ostringstream ss;
ss << r_step;
#ifdef LIBMESH_HAVE_EXODUS_API
std::string error_output = "error_"+ss.str()+".e";
#else
std::string error_output = "error_"+ss.str()+".gmv";
#endif
error.plot_error( error_output, mesh );
mesh_refinement.flag_elements_by_error_fraction (error);
if (refine_type == "p")
mesh_refinement.switch_h_to_p_refinement();
if (refine_type == "matchedhp")
mesh_refinement.add_p_to_h_refinement();
if (refine_type == "hp")
{
HPCoarsenTest hpselector;
hpselector.select_refinement(system);
}
if (refine_type == "singularhp")
{
libmesh_assert (singularity);
HPSingularity hpselector;
hpselector.singular_points.push_back(Point());
hpselector.select_refinement(system);
}
mesh_refinement.refine_and_coarsen_elements();
}
else if (uniform_refine == 1)
{
if (refine_type == "h" || refine_type == "hp" ||
refine_type == "matchedhp")
mesh_refinement.uniformly_refine(1);
if (refine_type == "p" || refine_type == "hp" ||
refine_type == "matchedhp")
mesh_refinement.uniformly_p_refine(1);
}
equation_systems.reinit ();
}
}
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO (mesh).write_equation_systems ("lshaped.e",
equation_systems);
#endif // #ifdef LIBMESH_HAVE_EXODUS_API
out << "];" << std::endl;
out << "hold on" << std::endl;
out << "plot(e(:,1), e(:,2), 'bo-');" << std::endl;
out << "plot(e(:,1), e(:,3), 'ro-');" << std::endl;
out << "xlabel('dofs');" << std::endl;
out << "title('" << approx_name << " elements');" << std::endl;
out << "legend('L2-error', 'H1-error');" << std::endl;
#endif // #ifndef LIBMESH_ENABLE_AMR
return 0;
}
Number exact_solution(const Point& p,
const Parameters&, // parameters, not needed
const std::string&, // sys_name, not needed
const std::string&) // unk_name, not needed
{
const Real x = p(0);
const Real y = (dim > 1) ? p(1) : 0.;
if (singularity)
{
Real theta = atan2(y,x);
if (theta < 0)
theta += 2. * libMesh::pi;
const Real z = (dim > 2) ? p(2) : 0;
return pow(x*x + y*y, 1./3.)*sin(2./3.*theta) + z;
}
else
{
const Real z = (dim > 2) ? p(2) : 0;
return cos(x) * exp(y) * (1. - z);
}
}
Gradient exact_derivative(const Point& p,
const Parameters&, // parameters, not needed
const std::string&, // sys_name, not needed
const std::string&) // unk_name, not needed
{
Gradient gradu;
const Real x = p(0);
const Real y = dim > 1 ? p(1) : 0.;
if (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;
if (dim > 1)
gradu(1) = tt*y*pow(r2,-tt)*sin(tt*theta) + pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*1./x;
if (dim > 2)
gradu(2) = 1.;
}
else
{
const Real z = (dim > 2) ? p(2) : 0;
gradu(0) = -sin(x) * exp(y) * (1. - z);
if (dim > 1)
gradu(1) = cos(x) * exp(y) * (1. - z);
if (dim > 2)
gradu(2) = -cos(x) * exp(y);
}
return gradu;
}
void assemble_laplace(EquationSystems& es,
const std::string& system_name)
{
#ifdef LIBMESH_ENABLE_AMR
libmesh_assert_equal_to (system_name, "Laplace");
PerfLog perf_log ("Matrix Assembly",false);
const MeshBase& mesh = es.get_mesh();
const unsigned int mesh_dim = mesh.mesh_dimension();
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Laplace");
const DofMap& dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(0);
AutoPtr<FEBase> fe (FEBase::build(mesh_dim, fe_type));
AutoPtr<FEBase> fe_face (FEBase::build(mesh_dim, fe_type));
AutoPtr<QBase> qrule(fe_type.default_quadrature_rule(mesh_dim));
AutoPtr<QBase> qface(fe_type.default_quadrature_rule(mesh_dim-1));
fe->attach_quadrature_rule (qrule.get());
fe_face->attach_quadrature_rule (qface.get());
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<Real>& JxW_face = fe_face->get_JxW();
const std::vector<Point>& q_point = fe->get_xyz();
const std::vector<std::vector<Real> >& phi = fe->get_phi();
const std::vector<std::vector<Real> >& psi = fe_face->get_phi();
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
const std::vector<Point>& qface_points = fe_face->get_xyz();
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
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)
{
perf_log.push("elem init");
const Elem* elem = *el;
dof_map.dof_indices (elem, dof_indices);
fe->reinit (elem);
Ke.resize (dof_indices.size(),
dof_indices.size());
Fe.resize (dof_indices.size());
perf_log.pop("elem init");
perf_log.push ("Ke");
for (unsigned int qp=0; qp<qrule->n_points(); qp++)
for (unsigned int i=0; i<dphi.size(); i++)
for (unsigned int j=0; j<dphi.size(); j++)
Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
if (mesh_dim == 1)
for (unsigned int qp=0; qp<qrule->n_points(); qp++)
{
Real x = q_point[qp](0);
Real f = singularity ? sqrt(3.)/9.*pow(-x, -4./3.) :
cos(x);
for (unsigned int i=0; i<dphi.size(); ++i)
Fe(i) += JxW[qp]*phi[i][qp]*f;
}
perf_log.pop ("Ke");
{
perf_log.push ("BCs");
const Real penalty = 1.e10;
for (unsigned int s=0; s<elem->n_sides(); s++)
if (elem->neighbor(s) == NULL)
{
fe_face->reinit(elem,s);
for (unsigned int qp=0; qp<qface->n_points(); qp++)
{
const Number value = exact_solution (qface_points[qp],
es.parameters,
"null",
"void");
for (unsigned int i=0; i<psi.size(); i++)
Fe(i) += penalty*JxW_face[qp]*value*psi[i][qp];
for (unsigned int i=0; i<psi.size(); i++)
for (unsigned int j=0; j<psi.size(); j++)
Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
}
}
perf_log.pop ("BCs");
}
perf_log.push ("matrix insertion");
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);
perf_log.pop ("matrix insertion");
}
#endif // #ifdef LIBMESH_ENABLE_AMR
}
The console output of the program:
***************************************************************
* Running Example adaptivity_ex3:
* 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
***************************************************************
EquationSystems
n_systems()=1
System #0, "Laplace"
Type "LinearImplicit"
Variables="u"
Finite Element Types="LAGRANGE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="SECOND", "THIRD"
n_dofs()=21
n_local_dofs()=9
n_constrained_dofs()=0
n_local_constrained_dofs()=0
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 8.42857
Average Off-Processor Bandwidth <= 2.28571
Maximum On-Processor Bandwidth <= 12
Maximum Off-Processor Bandwidth <= 12
DofMap Constraints
Number of DoF Constraints = 0
Number of Node Constraints = 0
Beginning Solve 0
System has: 21 degrees of freedom.
Linear solver converged at step: 7, final residual: 1.24062e-16
L2-Error is: 0.0150899
H1-Error is: 0.125332
Refining the mesh...
Beginning Solve 1
System has: 65 degrees of freedom.
Linear solver converged at step: 14, final residual: 1.43522e-15
L2-Error is: 0.00515425
H1-Error is: 0.0777803
Refining the mesh...
Beginning Solve 2
System has: 97 degrees of freedom.
Linear solver converged at step: 17, final residual: 1.17589e-15
L2-Error is: 0.00199025
H1-Error is: 0.0494318
Refining the mesh...
Beginning Solve 3
System has: 129 degrees of freedom.
Linear solver converged at step: 20, final residual: 3.95342e-16
L2-Error is: 0.000890215
H1-Error is: 0.0320056
Refining the mesh...
Beginning Solve 4
System has: 161 degrees of freedom.
Linear solver converged at step: 21, final residual: 3.23728e-15
L2-Error is: 0.000559523
H1-Error is: 0.0215069
Refining the mesh...
Beginning Solve 5
System has: 193 degrees of freedom.
Linear solver converged at step: 21, final residual: 2.55831e-15
L2-Error is: 0.000483386
H1-Error is: 0.0154837
Refining the mesh...
Beginning Solve 6
System has: 225 degrees of freedom.
Linear solver converged at step: 25, final residual: 3.92685e-15
L2-Error is: 0.000467457
H1-Error is: 0.0123018
Refining the mesh...
Beginning Solve 7
System has: 257 degrees of freedom.
Linear solver converged at step: 26, final residual: 1.73009e-15
L2-Error is: 0.000463656
H1-Error is: 0.0107815
Refining the mesh...
Beginning Solve 8
System has: 341 degrees of freedom.
Linear solver converged at step: 20, final residual: 1.52527e-15
L2-Error is: 0.000301328
H1-Error is: 0.00820153
Refining the mesh...
Beginning Solve 9
System has: 445 degrees of freedom.
Linear solver converged at step: 23, final residual: 4.3541e-15
L2-Error is: 0.000182515
H1-Error is: 0.00601168
************************************************************************************************************************
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***
************************************************************************************************************************
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------
/workspace/libmesh/examples/adaptivity/adaptivity_ex3/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:38:43 2013
Using Petsc Release Version 3.3.0, Patch 2, Fri Jul 13 15:42:00 CDT 2012
Max Max/Min Avg Total
Time (sec): 6.759e-01 1.00000 6.759e-01
Objects: 7.080e+02 1.00283 7.070e+02
Flops: 1.158e+07 1.27820 1.032e+07 2.063e+07
Flops/sec: 1.713e+07 1.27820 1.526e+07 3.053e+07
MPI Messages: 3.875e+02 1.00000 3.875e+02 7.750e+02
MPI Message Lengths: 1.128e+05 1.00000 2.910e+02 2.255e+05
MPI Reductions: 1.238e+03 1.00162
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.7583e-01 100.0% 2.0634e+07 100.0% 7.750e+02 100.0% 2.910e+02 100.0% 1.236e+03 100.0%
------------------------------------------------------------------------------------------------------------------------
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 194 1.0 1.4639e-03 2.2 6.31e+05 1.1 0.0e+00 0.0e+00 1.9e+02 0 6 0 0 16 0 6 0 0 16 825
VecNorm 214 1.0 5.1367e-03 6.8 5.92e+04 1.1 0.0e+00 0.0e+00 2.1e+02 0 1 0 0 17 0 1 0 0 17 22
VecScale 204 1.0 7.2956e-05 1.0 2.83e+04 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 745
VecCopy 47 1.0 1.9789e-05 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
VecSet 345 1.0 8.3208e-05 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecAXPY 20 1.0 1.1607e-02 1.6 4.95e+03 1.1 0.0e+00 0.0e+00 0.0e+00 1 0 0 0 0 1 0 0 0 0 1
VecMAXPY 204 1.0 2.2721e-04 1.1 6.87e+05 1.1 0.0e+00 0.0e+00 0.0e+00 0 6 0 0 0 0 6 0 0 0 5794
VecAssemblyBegin 109 1.0 1.7886e-03 1.3 0.00e+00 0.0 4.8e+01 9.4e+01 2.5e+02 0 0 6 2 20 0 0 6 2 20 0
VecAssemblyEnd 109 1.0 6.2943e-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 249 1.0 2.9969e-04 1.0 0.00e+00 0.0 5.0e+02 3.2e+02 0.0e+00 0 0 64 71 0 0 0 64 71 0 0
VecScatterEnd 249 1.0 3.3672e-03 3.9 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
VecNormalize 204 1.0 5.2562e-03 6.1 8.50e+04 1.1 0.0e+00 0.0e+00 2.0e+02 0 1 0 0 17 0 1 0 0 17 31
MatMult 204 1.0 4.1926e-03 2.5 8.95e+05 1.2 4.1e+02 2.8e+02 0.0e+00 0 8 53 50 0 0 8 53 50 0 398
MatSolve 214 1.0 1.8866e-03 1.3 4.67e+06 1.2 0.0e+00 0.0e+00 0.0e+00 0 41 0 0 0 0 41 0 0 0 4538
MatLUFactorNum 10 1.0 3.2549e-03 1.6 4.60e+06 1.5 0.0e+00 0.0e+00 0.0e+00 0 37 0 0 0 0 37 0 0 0 2366
MatILUFactorSym 10 1.0 8.6668e-03 1.3 0.00e+00 0.0 0.0e+00 0.0e+00 3.0e+01 1 0 0 0 2 1 0 0 0 2 0
MatAssemblyBegin 20 1.0 2.7285e-03 2.5 0.00e+00 0.0 4.5e+01 8.8e+02 4.0e+01 0 0 6 18 3 0 0 6 18 3 0
MatAssemblyEnd 20 1.0 1.0047e-03 1.0 0.00e+00 0.0 4.0e+01 6.4e+01 8.0e+01 0 0 5 1 6 0 0 5 1 6 0
MatGetRowIJ 10 1.0 3.8147e-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 10 1.0 3.5095e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 3.9e+01 0 0 0 0 3 0 0 0 0 3 0
MatZeroEntries 30 1.0 3.5763e-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
KSPGMRESOrthog 194 1.0 1.7748e-03 1.8 1.26e+06 1.1 0.0e+00 0.0e+00 1.9e+02 0 12 0 0 16 0 12 0 0 16 1364
KSPSetUp 20 1.0 1.6212e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
KSPSolve 10 1.0 3.2936e-02 1.0 1.16e+07 1.3 4.1e+02 2.8e+02 5.0e+02 5100 53 50 40 5100 53 50 40 626
PCSetUp 20 1.0 1.3600e-02 1.3 4.60e+06 1.5 0.0e+00 0.0e+00 8.9e+01 2 37 0 0 7 2 37 0 0 7 566
PCSetUpOnBlocks 10 1.0 1.2790e-02 1.3 4.60e+06 1.5 0.0e+00 0.0e+00 6.9e+01 2 37 0 0 6 2 37 0 0 6 602
PCApply 214 1.0 2.9938e-03 1.2 4.67e+06 1.2 0.0e+00 0.0e+00 0.0e+00 0 41 0 0 0 0 41 0 0 0 2860
------------------------------------------------------------------------------------------------------------------------
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 437 437 1059344 0
Vector Scatter 46 46 47656 0
Index Set 123 123 106420 0
IS L to G Mapping 19 19 10716 0
Matrix 40 40 1557632 0
Krylov Solver 20 20 193600 0
Preconditioner 20 20 17840 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 9.53674e-08
Average time for MPI_Barrier(): 1.23978e-06
Average time for zero size MPI_Send(): 2.15769e-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:38:43 2013 |
| OS: Linux |
| HostName: hbar.ices.utexas.edu |
| OS Release: 2.6.32-279.1.1.el6.x86_64 |
| OS Version: #1 SMP Tue Jul 10 11:24:23 CDT 2012 |
| Machine: x86_64 |
| Username: benkirk |
| Configuration: ./configure '--enable-everything' |
| '--prefix=/workspace/libmesh/install' |
| 'CXX=mpicxx' |
| 'CC=mpicc' |
| 'F77=mpif77' |
| 'FC=mpif90' |
| 'PETSC_DIR=/opt/apps/ossw/libraries/petsc/petsc-3.3-p2' |
| 'PETSC_ARCH=intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt' |
| 'SLEPC_DIR=/opt/apps/ossw/libraries/slepc/slepc-3.3-p2-petsc-3.3-p2-cxx-opt' |
| 'TRILINOS_DIR=/opt/apps/ossw/libraries/trilinos/trilinos-10.12.2/sl6/intel-12.1/mpich2-1.4.1p1/mkl-intel-10.3.12.361'|
| 'VTK_DIR=/opt/apps/ossw/libraries/vtk/vtk-5.10.0/sl6/intel-12.1' |
----------------------------------------------------------------------------------------------------------------------
----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=0.685827, Active time=0.617485 |
----------------------------------------------------------------------------------------------------------------
| 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() 19 0.0151 0.000797 0.0256 0.001345 2.45 4.14 |
| build_constraint_matrix() 229 0.0039 0.000017 0.0039 0.000017 0.63 0.63 |
| build_sparsity() 10 0.0121 0.001213 0.0290 0.002902 1.96 4.70 |
| cnstrn_elem_mat_vec() 229 0.0029 0.000013 0.0029 0.000013 0.46 0.46 |
| create_dof_constraints() 19 0.0216 0.001139 0.0472 0.002485 3.51 7.65 |
| distribute_dofs() 19 0.0231 0.001216 0.0586 0.003084 3.74 9.49 |
| dof_indices() 3236 0.1345 0.000042 0.1345 0.000042 21.79 21.79 |
| enforce_constraints_exactly() 8 0.0009 0.000117 0.0009 0.000117 0.15 0.15 |
| old_dof_indices() 470 0.0273 0.000058 0.0273 0.000058 4.43 4.43 |
| prepare_send_list() 19 0.0003 0.000016 0.0003 0.000016 0.05 0.05 |
| reinit() 19 0.0299 0.001575 0.0299 0.001575 4.85 4.85 |
| |
| EquationSystems |
| build_solution_vector() 1 0.0004 0.000410 0.0039 0.003927 0.07 0.64 |
| |
| ExodusII_IO |
| write_nodal_data() 1 0.0009 0.000856 0.0009 0.000856 0.14 0.14 |
| |
| FE |
| compute_shape_functions() 1298 0.0326 0.000025 0.0326 0.000025 5.28 5.28 |
| init_shape_functions() 846 0.0105 0.000012 0.0105 0.000012 1.70 1.70 |
| inverse_map() 5262 0.0300 0.000006 0.0300 0.000006 4.86 4.86 |
| |
| FEMap |
| compute_affine_map() 1298 0.0103 0.000008 0.0103 0.000008 1.68 1.68 |
| compute_face_map() 471 0.0080 0.000017 0.0201 0.000043 1.29 3.25 |
| init_face_shape_functions() 18 0.0002 0.000009 0.0002 0.000009 0.03 0.03 |
| init_reference_to_physical_map() 846 0.0193 0.000023 0.0193 0.000023 3.12 3.12 |
| |
| JumpErrorEstimator |
| estimate_error() 9 0.0211 0.002350 0.1302 0.014462 3.43 21.08 |
| |
| LocationMap |
| find() 740 0.0031 0.000004 0.0031 0.000004 0.50 0.50 |
| init() 18 0.0038 0.000210 0.0038 0.000210 0.61 0.61 |
| |
| Mesh |
| all_first_order() 9 0.0038 0.000424 0.0038 0.000424 0.62 0.62 |
| all_second_order() 1 0.0003 0.000265 0.0003 0.000265 0.04 0.04 |
| contract() 9 0.0003 0.000036 0.0019 0.000214 0.05 0.31 |
| find_neighbors() 29 0.0185 0.000638 0.0222 0.000766 3.00 3.60 |
| renumber_nodes_and_elem() 9 0.0016 0.000178 0.0016 0.000178 0.26 0.26 |
| |
| MeshCommunication |
| broadcast() 1 0.0007 0.000687 0.0010 0.000998 0.11 0.16 |
| compute_hilbert_indices() 21 0.0039 0.000184 0.0039 0.000184 0.63 0.63 |
| find_global_indices() 21 0.0025 0.000121 0.0118 0.000563 0.41 1.92 |
| parallel_sort() 21 0.0025 0.000118 0.0037 0.000176 0.40 0.60 |
| |
| MeshOutput |
| write_equation_systems() 1 0.0001 0.000058 0.0049 0.004890 0.01 0.79 |
| |
| MeshRefinement |
| _coarsen_elements() 18 0.0003 0.000015 0.0004 0.000021 0.04 0.06 |
| _refine_elements() 18 0.0040 0.000222 0.0112 0.000623 0.65 1.82 |
| add_point() 740 0.0026 0.000004 0.0059 0.000008 0.43 0.95 |
| make_coarsening_compatible() 36 0.0070 0.000195 0.0070 0.000195 1.13 1.13 |
| make_refinement_compatible() 36 0.0004 0.000011 0.0006 0.000016 0.07 0.09 |
| |
| MetisPartitioner |
| partition() 20 0.0152 0.000759 0.0291 0.001457 2.46 4.72 |
| |
| Parallel |
| allgather() 99 0.0017 0.000017 0.0019 0.000020 0.28 0.31 |
| broadcast() 9 0.0002 0.000025 0.0002 0.000020 0.04 0.03 |
| max(bool) 92 0.0017 0.000018 0.0017 0.000018 0.27 0.27 |
| max(scalar) 2495 0.0150 0.000006 0.0150 0.000006 2.42 2.42 |
| max(vector) 593 0.0047 0.000008 0.0122 0.000021 0.76 1.98 |
| min(bool) 3058 0.0125 0.000004 0.0125 0.000004 2.03 2.03 |
| min(scalar) 2423 0.0260 0.000011 0.0260 0.000011 4.22 4.22 |
| min(vector) 593 0.0050 0.000008 0.0131 0.000022 0.81 2.12 |
| probe() 196 0.0015 0.000007 0.0015 0.000007 0.24 0.24 |
| receive() 196 0.0009 0.000005 0.0024 0.000012 0.15 0.39 |
| send() 196 0.0005 0.000002 0.0005 0.000002 0.07 0.07 |
| send_receive() 238 0.0012 0.000005 0.0044 0.000019 0.20 0.72 |
| sum() 143 0.0014 0.000010 0.0027 0.000019 0.22 0.44 |
| |
| Parallel::Request |
| wait() 196 0.0003 0.000001 0.0003 0.000001 0.04 0.04 |
| |
| Partitioner |
| set_node_processor_ids() 29 0.0049 0.000168 0.0086 0.000297 0.79 1.40 |
| set_parent_processor_ids() 20 0.0013 0.000067 0.0013 0.000067 0.22 0.22 |
| |
| PetscLinearSolver |
| solve() 10 0.0382 0.003817 0.0382 0.003817 6.18 6.18 |
| |
| ProjectVector |
| operator() 9 0.0029 0.000318 0.0331 0.003677 0.46 5.36 |
| |
| System |
| assemble() 10 0.0180 0.001796 0.0595 0.005951 2.91 9.64 |
| project_vector() 9 0.0040 0.000448 0.0526 0.005843 0.65 8.52 |
----------------------------------------------------------------------------------------------------------------
| Totals: 26689 0.6175 100.00 |
----------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example adaptivity_ex3:
* 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
***************************************************************
Adaptivity Example 3 - Laplace Equation in the L-Shaped Domain
This example solves the Laplace equation on the classic "L-shaped" domain with adaptive mesh refinement. In this case, the exact solution is u(r,\theta) = r^{2/3} * \sin ( (2/3) * \theta), but the standard Kelly error indicator is used to estimate the error. The initial mesh contains three QUAD9 elements which represent the standard quadrants I, II, and III of the domain [-1,1]x[-1,1], i.e. Element 0: [-1,0]x[ 0,1] Element 1: [ 0,1]x[ 0,1] Element 2: [-1,0]x[-1,0] The mesh is provided in the standard libMesh ASCII format file named "lshaped.xda". In addition, an input file named "adaptivity_ex3.in" is provided which allows the user to set several parameters for the solution so that the problem can be re-run without a re-compile. The solution technique employed is to have a refinement loop with a linear solve inside followed by a refinement of the grid and projection of the solution to the new grid In the final loop iteration, there is no additional refinement after the solve. In the input file "adaptivity_ex3.in", the variable "max_r_steps" controls the number of refinement steps, "max_r_level" controls the maximum element refinement level, and "refine_percentage" / "coarsen_percentage" determine the number of elements which will be refined / coarsened at each step.
LibMesh include files.