The source file adaptivity_ex4.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/fe.h"
#include "libmesh/quadrature.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_modification.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/error_vector.h"
#include "libmesh/fourth_error_estimators.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/tensor_value.h"
#include "libmesh/perf_log.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 Biharmonic 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_biharmonic(EquationSystems& es,
const std::string& system_name);
Prototypes for calculation of the exact solution. Necessary
for setting boundary conditions.
Number exact_2D_solution(const Point& p,
const Parameters&, // parameters, not needed
const std::string&, // sys_name, not needed
const std::string&); // unk_name, not needed);
Number exact_3D_solution(const Point& p,
const Parameters&, const std::string&, const std::string&);
Prototypes for calculation of the gradient of the exact solution.
Necessary for setting boundary conditions in H^2_0 and testing
H^1 convergence of the solution
Gradient exact_2D_derivative(const Point& p,
const Parameters&, const std::string&, const std::string&);
Gradient exact_3D_derivative(const Point& p,
const Parameters&, const std::string&, const std::string&);
Tensor exact_2D_hessian(const Point& p,
const Parameters&, const std::string&, const std::string&);
Tensor exact_3D_hessian(const Point& p,
const Parameters&, const std::string&, const std::string&);
Number forcing_function_2D(const Point& p);
Number forcing_function_3D(const Point& p);
Pointers to dimension-independent functions
Number (*exact_solution)(const Point& p,
const Parameters&, const std::string&, const std::string&);
Gradient (*exact_derivative)(const Point& p,
const Parameters&, const std::string&, const std::string&);
Tensor (*exact_hessian)(const Point& p,
const Parameters&, const std::string&, const std::string&);
Number (*forcing_function)(const Point& p);
int main(int argc, char** argv)
{
Initialize libMesh.
LibMeshInit init (argc, argv);
Adaptive constraint calculations for fine Hermite elements seems
to require half-decent precision
#ifdef LIBMESH_DEFAULT_SINGLE_PRECISION
libmesh_example_assert(false, "double precision");
#endif
This example requires Adaptive Mesh Refinement support
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_assert(false, "--enable-amr");
#else
This example requires second derivative calculation support
#ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
libmesh_example_assert(false, "--enable-second");
#else
Parse the input file
GetPot input_file("adaptivity_ex4.in");
Read in parameters from the input file
const unsigned int max_r_level = input_file("max_r_level", 10);
const unsigned int max_r_steps = input_file("max_r_steps", 4);
const std::string approx_type = input_file("approx_type",
"HERMITE");
const unsigned int uniform_refine =
input_file("uniform_refine", 0);
const Real refine_percentage =
input_file("refine_percentage", 0.5);
const Real coarsen_percentage =
input_file("coarsen_percentage", 0.5);
const unsigned int dim =
input_file("dimension", 2);
const unsigned int max_linear_iterations =
input_file("max_linear_iterations", 10000);
Skip higher-dimensional examples on a lower-dimensional libMesh build
libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");
We have only defined 2 and 3 dimensional problems
libmesh_assert (dim == 2 || dim == 3);
Currently only the Hermite cubics give a 3D C^1 basis
libmesh_assert (dim == 2 || approx_type == "HERMITE");
Create a mesh.
Mesh mesh;
Output file for plotting the error
std::string output_file = "";
if (dim == 2)
output_file += "2D_";
else if (dim == 3)
output_file += "3D_";
if (approx_type == "HERMITE")
output_file += "hermite_";
else if (approx_type == "SECOND")
output_file += "reducedclough_";
else
output_file += "clough_";
if (uniform_refine == 0)
output_file += "adaptive";
else
output_file += "uniform";
std::string exd_file = output_file;
exd_file += ".e";
output_file += ".m";
std::ofstream out (output_file.c_str());
out << "% dofs L2-error H1-error H2-error\n"
<< "e = [\n";
Set up the dimension-dependent coarse mesh and solution
We build more than one cell so as to avoid bugs on fewer than
4 processors in 2D or 8 in 3D.
if (dim == 2)
{
MeshTools::Generation::build_square(mesh, 2, 2);
exact_solution = &exact_2D_solution;
exact_derivative = &exact_2D_derivative;
exact_hessian = &exact_2D_hessian;
forcing_function = &forcing_function_2D;
}
else if (dim == 3)
{
MeshTools::Generation::build_cube(mesh, 2, 2, 2);
exact_solution = &exact_3D_solution;
exact_derivative = &exact_3D_derivative;
exact_hessian = &exact_3D_hessian;
forcing_function = &forcing_function_3D;
}
Convert the mesh to second order: necessary for computing with
Clough-Tocher elements, useful for getting slightly less
broken visualization output with Hermite elements
mesh.all_second_order();
Convert it to triangles if necessary
if (approx_type != "HERMITE")
MeshTools::Modification::all_tri(mesh);
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.
Create a system named "Biharmonic"
LinearImplicitSystem& system =
equation_systems.add_system<LinearImplicitSystem> ("Biharmonic");
Adds the variable "u" to "Biharmonic". "u"
will be approximated using Hermite tensor product squares
or (possibly reduced) cubic Clough-Tocher triangles
if (approx_type == "HERMITE")
system.add_variable("u", THIRD, HERMITE);
else if (approx_type == "SECOND")
system.add_variable("u", SECOND, CLOUGH);
else if (approx_type == "CLOUGH")
system.add_variable("u", THIRD, CLOUGH);
else
libmesh_error();
Give the system a pointer to the matrix assembly
function.
system.attach_assemble_function (assemble_biharmonic);
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") = TOLERANCE * TOLERANCE;
Prints information about the system to the screen.
equation_systems.print_info();
Construct ExactSolution object and attach function to compute exact solution
ExactSolution exact_sol(equation_systems);
exact_sol.attach_exact_value(exact_solution);
exact_sol.attach_exact_deriv(exact_derivative);
exact_sol.attach_exact_hessian(exact_hessian);
Construct zero solution object, useful for computing solution norms
Attaching "zero_solution" functions is unnecessary
ExactSolution zero_sol(equation_systems);
A refinement loop.
for (unsigned int r_step=0; r_step<max_r_steps; r_step++)
{
mesh.print_info();
equation_systems.print_info();
std::cout << "Beginning Solve " << r_step << std::endl;
Solve the system "Biharmonic", just like example 2.
system.solve();
std::cout << "Linear solver converged at step: "
<< system.n_linear_iterations()
<< ", final residual: "
<< system.final_linear_residual()
<< std::endl;
Compute the error.
exact_sol.compute_error("Biharmonic", "u");
Compute the norm.
zero_sol.compute_error("Biharmonic", "u");
Print out the error values
std::cout << "L2-Norm is: "
<< zero_sol.l2_error("Biharmonic", "u")
<< std::endl;
std::cout << "H1-Norm is: "
<< zero_sol.h1_error("Biharmonic", "u")
<< std::endl;
std::cout << "H2-Norm is: "
<< zero_sol.h2_error("Biharmonic", "u")
<< std::endl
<< std::endl;
std::cout << "L2-Error is: "
<< exact_sol.l2_error("Biharmonic", "u")
<< std::endl;
std::cout << "H1-Error is: "
<< exact_sol.h1_error("Biharmonic", "u")
<< std::endl;
std::cout << "H2-Error is: "
<< exact_sol.h2_error("Biharmonic", "u")
<< std::endl
<< std::endl;
Print to output file
out << equation_systems.n_active_dofs() << " "
<< exact_sol.l2_error("Biharmonic", "u") << " "
<< exact_sol.h1_error("Biharmonic", "u") << " "
<< exact_sol.h2_error("Biharmonic", "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)
{
ErrorVector error;
LaplacianErrorEstimator error_estimator;
error_estimator.estimate_error(system, error);
mesh_refinement.flag_elements_by_elem_fraction (error);
std::cout << "Mean Error: " << error.mean() <<
std::endl;
std::cout << "Error Variance: " << error.variance() <<
std::endl;
mesh_refinement.refine_and_coarsen_elements();
}
else
{
mesh_refinement.uniformly_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 (exd_file,
equation_systems);
#endif // #ifdef LIBMESH_HAVE_EXODUS_API
Close up the output file.
out << "];\n"
<< "hold on\n"
<< "loglog(e(:,1), e(:,2), 'bo-');\n"
<< "loglog(e(:,1), e(:,3), 'ro-');\n"
<< "loglog(e(:,1), e(:,4), 'go-');\n"
<< "xlabel('log(dofs)');\n"
<< "ylabel('log(error)');\n"
<< "title('C1 " << approx_type << " elements');\n"
<< "legend('L2-error', 'H1-error', 'H2-error');\n";
All done.
return 0;
#endif // #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
#endif // #ifndef LIBMESH_ENABLE_AMR
}
Number exact_2D_solution(const Point& p,
const Parameters&, // parameters, not needed
const std::string&, // sys_name, not needed
const std::string&) // unk_name, not needed
{
x and y coordinates in space
const Real x = p(0);
const Real y = p(1);
analytic solution value
return 256.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y);
}
We now define the gradient of the exact solution
Gradient exact_2D_derivative(const Point& p,
const Parameters&, // parameters, not needed
const std::string&, // sys_name, not needed
const std::string&) // unk_name, not needed
{
x and y coordinates in space
const Real x = p(0);
const Real y = p(1);
First derivatives to be returned.
Gradient gradu;
gradu(0) = 256.*2.*(x-x*x)*(1-2*x)*(y-y*y)*(y-y*y);
gradu(1) = 256.*2.*(x-x*x)*(x-x*x)*(y-y*y)*(1-2*y);
return gradu;
}
We now define the hessian of the exact solution
Tensor exact_2D_hessian(const Point& p,
const Parameters&, // parameters, not needed
const std::string&, // sys_name, not needed
const std::string&) // unk_name, not needed
{
Second derivatives to be returned.
Tensor hessu;
x and y coordinates in space
const Real x = p(0);
const Real y = p(1);
hessu(0,0) = 256.*2.*(1-6.*x+6.*x*x)*(y-y*y)*(y-y*y);
hessu(0,1) = 256.*4.*(x-x*x)*(1.-2.*x)*(y-y*y)*(1.-2.*y);
hessu(1,1) = 256.*2.*(x-x*x)*(x-x*x)*(1.-6.*y+6.*y*y);
Hessians are always symmetric
hessu(1,0) = hessu(0,1);
return hessu;
}
Number forcing_function_2D(const Point& p)
{
x and y coordinates in space
const Real x = p(0);
const Real y = p(1);
Equals laplacian(laplacian(u))
return 256. * 8. * (3.*((y-y*y)*(y-y*y)+(x-x*x)*(x-x*x))
+ (1.-6.*x+6.*x*x)*(1.-6.*y+6.*y*y));
}
Number exact_3D_solution(const Point& p,
const Parameters&, // parameters, not needed
const std::string&, // sys_name, not needed
const std::string&) // unk_name, not needed
{
xyz coordinates in space
const Real x = p(0);
const Real y = p(1);
const Real z = p(2);
analytic solution value
return 4096.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y)*(z-z*z)*(z-z*z);
}
Gradient exact_3D_derivative(const Point& p,
const Parameters&, // parameters, not needed
const std::string&, // sys_name, not needed
const std::string&) // unk_name, not needed
{
First derivatives to be returned.
Gradient gradu;
xyz coordinates in space
const Real x = p(0);
const Real y = p(1);
const Real z = p(2);
gradu(0) = 4096.*2.*(x-x*x)*(1.-2.*x)*(y-y*y)*(y-y*y)*(z-z*z)*(z-z*z);
gradu(1) = 4096.*2.*(x-x*x)*(x-x*x)*(y-y*y)*(1.-2.*y)*(z-z*z)*(z-z*z);
gradu(2) = 4096.*2.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y)*(z-z*z)*(1.-2.*z);
return gradu;
}
We now define the hessian of the exact solution
Tensor exact_3D_hessian(const Point& p,
const Parameters&, // parameters, not needed
const std::string&, // sys_name, not needed
const std::string&) // unk_name, not needed
{
Second derivatives to be returned.
Tensor hessu;
xyz coordinates in space
const Real x = p(0);
const Real y = p(1);
const Real z = p(2);
hessu(0,0) = 4096.*(2.-12.*x+12.*x*x)*(y-y*y)*(y-y*y)*(z-z*z)*(z-z*z);
hessu(0,1) = 4096.*4.*(x-x*x)*(1.-2.*x)*(y-y*y)*(1.-2.*y)*(z-z*z)*(z-z*z);
hessu(0,2) = 4096.*4.*(x-x*x)*(1.-2.*x)*(y-y*y)*(y-y*y)*(z-z*z)*(1.-2.*z);
hessu(1,1) = 4096.*(x-x*x)*(x-x*x)*(2.-12.*y+12.*y*y)*(z-z*z)*(z-z*z);
hessu(1,2) = 4096.*4.*(x-x*x)*(x-x*x)*(y-y*y)*(1.-2.*y)*(z-z*z)*(1.-2.*z);
hessu(2,2) = 4096.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y)*(2.-12.*z+12.*z*z);
Hessians are always symmetric
hessu(1,0) = hessu(0,1);
hessu(2,0) = hessu(0,2);
hessu(2,1) = hessu(1,2);
return hessu;
}
Number forcing_function_3D(const Point& p)
{
xyz coordinates in space
const Real x = p(0);
const Real y = p(1);
const Real z = p(2);
Equals laplacian(laplacian(u))
return 4096. * 8. * (3.*((y-y*y)*(y-y*y)*(x-x*x)*(x-x*x) +
(z-z*z)*(z-z*z)*(x-x*x)*(x-x*x) +
(z-z*z)*(z-z*z)*(y-y*y)*(y-y*y)) +
(1.-6.*x+6.*x*x)*(1.-6.*y+6.*y*y)*(z-z*z)*(z-z*z) +
(1.-6.*x+6.*x*x)*(1.-6.*z+6.*z*z)*(y-y*y)*(y-y*y) +
(1.-6.*y+6.*y*y)*(1.-6.*z+6.*z*z)*(x-x*x)*(x-x*x));
}
We now define the matrix assembly function for the
Biharmonic 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_biharmonic(EquationSystems& es,
const std::string& system_name)
{
#ifdef LIBMESH_ENABLE_AMR
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
It is a good idea to make sure we are assembling
the proper system.
libmesh_assert_equal_to (system_name, "Biharmonic");
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 dim = mesh.mesh_dimension();
Get a reference to the LinearImplicitSystem we are solving
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Biharmonic");
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(dim, fe_type));
Quadrature rule for numerical integration.
With 2D triangles, the Clough quadrature rule puts a Gaussian
quadrature rule on each of the 3 subelements
AutoPtr<QBase> qrule(fe_type.default_quadrature_rule(dim));
Tell the finite element object to use our quadrature rule.
fe->attach_quadrature_rule (qrule.get());
Declare a special finite element object for
boundary integration.
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
Boundary integration requires another quadraure rule,
with dimensionality one less than the dimensionality
of the element.
In 1D, the Clough and Gauss quadrature rules are identical.
AutoPtr<QBase> qface(fe_type.default_quadrature_rule(dim-1));
Tell the finte element object to use our
quadrature rule.
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();
The physical XY locations of the quadrature points on the element.
These might be useful for evaluating spatially varying material
properties at the quadrature points.
const std::vector<Point>& q_point = fe->get_xyz();
The element shape functions evaluated at the quadrature points.
const std::vector<std::vector<Real> >& phi = fe->get_phi();
The element shape function second derivatives evaluated at the
quadrature points. Note that for the simple biharmonic, shape
function first derivatives are unnecessary.
const std::vector<std::vector<RealTensor> >& d2phi = fe->get_d2phi();
For efficiency we will compute shape function laplacians n times,
not n^2
std::vector<Real> shape_laplacian;
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.
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.
Ke.resize (dof_indices.size(),
dof_indices.size());
Fe.resize (dof_indices.size());
Make sure there is enough room in this cache
shape_laplacian.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 laplacians of the test funcions
(i) against laplacians of the trial functions (j).
This step is why we need the Clough-Tocher elements - these C1 differentiable elements have square-integrable second derivatives.
Now start logging the element matrix computation
This step is why we need the Clough-Tocher elements - these C1 differentiable elements have square-integrable second derivatives.
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<phi.size(); i++)
{
shape_laplacian[i] = d2phi[i][qp](0,0)+d2phi[i][qp](1,1);
if (dim == 3)
shape_laplacian[i] += d2phi[i][qp](2,2);
}
for (unsigned int i=0; i<phi.size(); i++)
for (unsigned int j=0; j<phi.size(); j++)
Ke(i,j) += JxW[qp]*
shape_laplacian[i]*shape_laplacian[j];
}
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. Note that this is a fourth-order
problem: Dirichlet boundary conditions include *both*
boundary values and boundary normal fluxes.
{
Start logging the boundary condition computation
perf_log.push ("BCs");
The penalty values, for solution boundary trace and flux.
const Real penalty = 1e10;
const Real penalty2 = 1e10;
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)
{
The value of the shape functions at the quadrature
points.
const std::vector<std::vector<Real> >& phi_face =
fe_face->get_phi();
The value of the shape function derivatives at the
quadrature points.
const std::vector<std::vector<RealGradient> >& dphi_face =
fe_face->get_dphi();
The Jacobian * Quadrature Weight at the quadrature
points on the face.
const std::vector<Real>& JxW_face = fe_face->get_JxW();
The XYZ locations (in physical space) of the
quadrature points on the face. This is where
we will interpolate the boundary value function.
const std::vector<Point>& qface_point = fe_face->get_xyz();
const std::vector<Point>& face_normals =
fe_face->get_normals();
Compute the shape function values on the element
face.
fe_face->reinit(elem, s);
Loop over the face quagrature points for integration.
for (unsigned int qp=0; qp<qface->n_points(); qp++)
{
The boundary value.
Number value = exact_solution(qface_point[qp],
es.parameters, "null",
"void");
Gradient flux = exact_2D_derivative(qface_point[qp],
es.parameters,
"null", "void");
Matrix contribution of the L2 projection.
Note that the basis function values are
integrated against test function values while
basis fluxes are integrated against test function
fluxes.
for (unsigned int i=0; i<phi_face.size(); i++)
for (unsigned int j=0; j<phi_face.size(); j++)
Ke(i,j) += JxW_face[qp] *
(penalty * phi_face[i][qp] *
phi_face[j][qp] + penalty2
* (dphi_face[i][qp] *
face_normals[qp]) *
(dphi_face[j][qp] *
face_normals[qp]));
Right-hand-side contribution of the L2
projection.
for (unsigned int i=0; i<phi_face.size(); i++)
Fe(i) += JxW_face[qp] *
(penalty * value * phi_face[i][qp]
+ penalty2 *
(flux * face_normals[qp])
* (dphi_face[i][qp]
* face_normals[qp]));
}
}
Stop logging the boundary condition computation
perf_log.pop ("BCs");
}
for (unsigned int qp=0; qp<qrule->n_points(); qp++)
for (unsigned int i=0; i<phi.size(); i++)
Fe(i) += JxW[qp]*phi[i][qp]*forcing_function(q_point[qp]);
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);
Stop 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?
#else
#endif // #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
#endif // #ifdef LIBMESH_ENABLE_AMR
}
The source file adaptivity_ex4.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/fe.h"
#include "libmesh/quadrature.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_modification.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/error_vector.h"
#include "libmesh/fourth_error_estimators.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/tensor_value.h"
#include "libmesh/perf_log.h"
using namespace libMesh;
void assemble_biharmonic(EquationSystems& es,
const std::string& system_name);
Number exact_2D_solution(const Point& p,
const Parameters&, // parameters, not needed
const std::string&, // sys_name, not needed
const std::string&); // unk_name, not needed);
Number exact_3D_solution(const Point& p,
const Parameters&, const std::string&, const std::string&);
Gradient exact_2D_derivative(const Point& p,
const Parameters&, const std::string&, const std::string&);
Gradient exact_3D_derivative(const Point& p,
const Parameters&, const std::string&, const std::string&);
Tensor exact_2D_hessian(const Point& p,
const Parameters&, const std::string&, const std::string&);
Tensor exact_3D_hessian(const Point& p,
const Parameters&, const std::string&, const std::string&);
Number forcing_function_2D(const Point& p);
Number forcing_function_3D(const Point& p);
Number (*exact_solution)(const Point& p,
const Parameters&, const std::string&, const std::string&);
Gradient (*exact_derivative)(const Point& p,
const Parameters&, const std::string&, const std::string&);
Tensor (*exact_hessian)(const Point& p,
const Parameters&, const std::string&, const std::string&);
Number (*forcing_function)(const Point& p);
int main(int argc, char** argv)
{
LibMeshInit init (argc, argv);
#ifdef LIBMESH_DEFAULT_SINGLE_PRECISION
libmesh_example_assert(false, "double precision");
#endif
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_assert(false, "--enable-amr");
#else
#ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
libmesh_example_assert(false, "--enable-second");
#else
GetPot input_file("adaptivity_ex4.in");
const unsigned int max_r_level = input_file("max_r_level", 10);
const unsigned int max_r_steps = input_file("max_r_steps", 4);
const std::string approx_type = input_file("approx_type",
"HERMITE");
const unsigned int uniform_refine =
input_file("uniform_refine", 0);
const Real refine_percentage =
input_file("refine_percentage", 0.5);
const Real coarsen_percentage =
input_file("coarsen_percentage", 0.5);
const unsigned int dim =
input_file("dimension", 2);
const unsigned int max_linear_iterations =
input_file("max_linear_iterations", 10000);
libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");
libmesh_assert (dim == 2 || dim == 3);
libmesh_assert (dim == 2 || approx_type == "HERMITE");
Mesh mesh;
std::string output_file = "";
if (dim == 2)
output_file += "2D_";
else if (dim == 3)
output_file += "3D_";
if (approx_type == "HERMITE")
output_file += "hermite_";
else if (approx_type == "SECOND")
output_file += "reducedclough_";
else
output_file += "clough_";
if (uniform_refine == 0)
output_file += "adaptive";
else
output_file += "uniform";
std::string exd_file = output_file;
exd_file += ".e";
output_file += ".m";
std::ofstream out (output_file.c_str());
out << "% dofs L2-error H1-error H2-error\n"
<< "e = [\n";
if (dim == 2)
{
MeshTools::Generation::build_square(mesh, 2, 2);
exact_solution = &exact_2D_solution;
exact_derivative = &exact_2D_derivative;
exact_hessian = &exact_2D_hessian;
forcing_function = &forcing_function_2D;
}
else if (dim == 3)
{
MeshTools::Generation::build_cube(mesh, 2, 2, 2);
exact_solution = &exact_3D_solution;
exact_derivative = &exact_3D_derivative;
exact_hessian = &exact_3D_hessian;
forcing_function = &forcing_function_3D;
}
mesh.all_second_order();
if (approx_type != "HERMITE")
MeshTools::Modification::all_tri(mesh);
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> ("Biharmonic");
if (approx_type == "HERMITE")
system.add_variable("u", THIRD, HERMITE);
else if (approx_type == "SECOND")
system.add_variable("u", SECOND, CLOUGH);
else if (approx_type == "CLOUGH")
system.add_variable("u", THIRD, CLOUGH);
else
libmesh_error();
system.attach_assemble_function (assemble_biharmonic);
equation_systems.init();
equation_systems.parameters.set<unsigned int>
("linear solver maximum iterations") = max_linear_iterations;
equation_systems.parameters.set<Real>
("linear solver tolerance") = TOLERANCE * TOLERANCE;
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.attach_exact_hessian(exact_hessian);
ExactSolution zero_sol(equation_systems);
for (unsigned int r_step=0; r_step<max_r_steps; r_step++)
{
mesh.print_info();
equation_systems.print_info();
std::cout << "Beginning Solve " << r_step << std::endl;
system.solve();
std::cout << "Linear solver converged at step: "
<< system.n_linear_iterations()
<< ", final residual: "
<< system.final_linear_residual()
<< std::endl;
exact_sol.compute_error("Biharmonic", "u");
zero_sol.compute_error("Biharmonic", "u");
std::cout << "L2-Norm is: "
<< zero_sol.l2_error("Biharmonic", "u")
<< std::endl;
std::cout << "H1-Norm is: "
<< zero_sol.h1_error("Biharmonic", "u")
<< std::endl;
std::cout << "H2-Norm is: "
<< zero_sol.h2_error("Biharmonic", "u")
<< std::endl
<< std::endl;
std::cout << "L2-Error is: "
<< exact_sol.l2_error("Biharmonic", "u")
<< std::endl;
std::cout << "H1-Error is: "
<< exact_sol.h1_error("Biharmonic", "u")
<< std::endl;
std::cout << "H2-Error is: "
<< exact_sol.h2_error("Biharmonic", "u")
<< std::endl
<< std::endl;
out << equation_systems.n_active_dofs() << " "
<< exact_sol.l2_error("Biharmonic", "u") << " "
<< exact_sol.h1_error("Biharmonic", "u") << " "
<< exact_sol.h2_error("Biharmonic", "u") << std::endl;
if (r_step+1 != max_r_steps)
{
std::cout << " Refining the mesh..." << std::endl;
if (uniform_refine == 0)
{
ErrorVector error;
LaplacianErrorEstimator error_estimator;
error_estimator.estimate_error(system, error);
mesh_refinement.flag_elements_by_elem_fraction (error);
std::cout << "Mean Error: " << error.mean() <<
std::endl;
std::cout << "Error Variance: " << error.variance() <<
std::endl;
mesh_refinement.refine_and_coarsen_elements();
}
else
{
mesh_refinement.uniformly_refine(1);
}
equation_systems.reinit ();
}
}
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO (mesh).write_equation_systems (exd_file,
equation_systems);
#endif // #ifdef LIBMESH_HAVE_EXODUS_API
out << "];\n"
<< "hold on\n"
<< "loglog(e(:,1), e(:,2), 'bo-');\n"
<< "loglog(e(:,1), e(:,3), 'ro-');\n"
<< "loglog(e(:,1), e(:,4), 'go-');\n"
<< "xlabel('log(dofs)');\n"
<< "ylabel('log(error)');\n"
<< "title('C1 " << approx_type << " elements');\n"
<< "legend('L2-error', 'H1-error', 'H2-error');\n";
return 0;
#endif // #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
#endif // #ifndef LIBMESH_ENABLE_AMR
}
Number exact_2D_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 = p(1);
return 256.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y);
}
Gradient exact_2D_derivative(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 = p(1);
Gradient gradu;
gradu(0) = 256.*2.*(x-x*x)*(1-2*x)*(y-y*y)*(y-y*y);
gradu(1) = 256.*2.*(x-x*x)*(x-x*x)*(y-y*y)*(1-2*y);
return gradu;
}
Tensor exact_2D_hessian(const Point& p,
const Parameters&, // parameters, not needed
const std::string&, // sys_name, not needed
const std::string&) // unk_name, not needed
{
Tensor hessu;
const Real x = p(0);
const Real y = p(1);
hessu(0,0) = 256.*2.*(1-6.*x+6.*x*x)*(y-y*y)*(y-y*y);
hessu(0,1) = 256.*4.*(x-x*x)*(1.-2.*x)*(y-y*y)*(1.-2.*y);
hessu(1,1) = 256.*2.*(x-x*x)*(x-x*x)*(1.-6.*y+6.*y*y);
hessu(1,0) = hessu(0,1);
return hessu;
}
Number forcing_function_2D(const Point& p)
{
const Real x = p(0);
const Real y = p(1);
return 256. * 8. * (3.*((y-y*y)*(y-y*y)+(x-x*x)*(x-x*x))
+ (1.-6.*x+6.*x*x)*(1.-6.*y+6.*y*y));
}
Number exact_3D_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 = p(1);
const Real z = p(2);
return 4096.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y)*(z-z*z)*(z-z*z);
}
Gradient exact_3D_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 = p(1);
const Real z = p(2);
gradu(0) = 4096.*2.*(x-x*x)*(1.-2.*x)*(y-y*y)*(y-y*y)*(z-z*z)*(z-z*z);
gradu(1) = 4096.*2.*(x-x*x)*(x-x*x)*(y-y*y)*(1.-2.*y)*(z-z*z)*(z-z*z);
gradu(2) = 4096.*2.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y)*(z-z*z)*(1.-2.*z);
return gradu;
}
Tensor exact_3D_hessian(const Point& p,
const Parameters&, // parameters, not needed
const std::string&, // sys_name, not needed
const std::string&) // unk_name, not needed
{
Tensor hessu;
const Real x = p(0);
const Real y = p(1);
const Real z = p(2);
hessu(0,0) = 4096.*(2.-12.*x+12.*x*x)*(y-y*y)*(y-y*y)*(z-z*z)*(z-z*z);
hessu(0,1) = 4096.*4.*(x-x*x)*(1.-2.*x)*(y-y*y)*(1.-2.*y)*(z-z*z)*(z-z*z);
hessu(0,2) = 4096.*4.*(x-x*x)*(1.-2.*x)*(y-y*y)*(y-y*y)*(z-z*z)*(1.-2.*z);
hessu(1,1) = 4096.*(x-x*x)*(x-x*x)*(2.-12.*y+12.*y*y)*(z-z*z)*(z-z*z);
hessu(1,2) = 4096.*4.*(x-x*x)*(x-x*x)*(y-y*y)*(1.-2.*y)*(z-z*z)*(1.-2.*z);
hessu(2,2) = 4096.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y)*(2.-12.*z+12.*z*z);
hessu(1,0) = hessu(0,1);
hessu(2,0) = hessu(0,2);
hessu(2,1) = hessu(1,2);
return hessu;
}
Number forcing_function_3D(const Point& p)
{
const Real x = p(0);
const Real y = p(1);
const Real z = p(2);
return 4096. * 8. * (3.*((y-y*y)*(y-y*y)*(x-x*x)*(x-x*x) +
(z-z*z)*(z-z*z)*(x-x*x)*(x-x*x) +
(z-z*z)*(z-z*z)*(y-y*y)*(y-y*y)) +
(1.-6.*x+6.*x*x)*(1.-6.*y+6.*y*y)*(z-z*z)*(z-z*z) +
(1.-6.*x+6.*x*x)*(1.-6.*z+6.*z*z)*(y-y*y)*(y-y*y) +
(1.-6.*y+6.*y*y)*(1.-6.*z+6.*z*z)*(x-x*x)*(x-x*x));
}
void assemble_biharmonic(EquationSystems& es,
const std::string& system_name)
{
#ifdef LIBMESH_ENABLE_AMR
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
libmesh_assert_equal_to (system_name, "Biharmonic");
PerfLog perf_log ("Matrix Assembly",false);
const MeshBase& mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Biharmonic");
const DofMap& dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(0);
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
AutoPtr<QBase> qrule(fe_type.default_quadrature_rule(dim));
fe->attach_quadrature_rule (qrule.get());
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
AutoPtr<QBase> qface(fe_type.default_quadrature_rule(dim-1));
fe_face->attach_quadrature_rule (qface.get());
const std::vector<Real>& JxW = fe->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<RealTensor> >& d2phi = fe->get_d2phi();
std::vector<Real> shape_laplacian;
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());
shape_laplacian.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<phi.size(); i++)
{
shape_laplacian[i] = d2phi[i][qp](0,0)+d2phi[i][qp](1,1);
if (dim == 3)
shape_laplacian[i] += d2phi[i][qp](2,2);
}
for (unsigned int i=0; i<phi.size(); i++)
for (unsigned int j=0; j<phi.size(); j++)
Ke(i,j) += JxW[qp]*
shape_laplacian[i]*shape_laplacian[j];
}
perf_log.pop ("Ke");
{
perf_log.push ("BCs");
const Real penalty = 1e10;
const Real penalty2 = 1e10;
for (unsigned int s=0; s<elem->n_sides(); s++)
if (elem->neighbor(s) == NULL)
{
const std::vector<std::vector<Real> >& phi_face =
fe_face->get_phi();
const std::vector<std::vector<RealGradient> >& dphi_face =
fe_face->get_dphi();
const std::vector<Real>& JxW_face = fe_face->get_JxW();
const std::vector<Point>& qface_point = fe_face->get_xyz();
const std::vector<Point>& face_normals =
fe_face->get_normals();
fe_face->reinit(elem, s);
for (unsigned int qp=0; qp<qface->n_points(); qp++)
{
Number value = exact_solution(qface_point[qp],
es.parameters, "null",
"void");
Gradient flux = exact_2D_derivative(qface_point[qp],
es.parameters,
"null", "void");
for (unsigned int i=0; i<phi_face.size(); i++)
for (unsigned int j=0; j<phi_face.size(); j++)
Ke(i,j) += JxW_face[qp] *
(penalty * phi_face[i][qp] *
phi_face[j][qp] + penalty2
* (dphi_face[i][qp] *
face_normals[qp]) *
(dphi_face[j][qp] *
face_normals[qp]));
for (unsigned int i=0; i<phi_face.size(); i++)
Fe(i) += JxW_face[qp] *
(penalty * value * phi_face[i][qp]
+ penalty2 *
(flux * face_normals[qp])
* (dphi_face[i][qp]
* face_normals[qp]));
}
}
perf_log.pop ("BCs");
}
for (unsigned int qp=0; qp<qrule->n_points(); qp++)
for (unsigned int i=0; i<phi.size(); i++)
Fe(i) += JxW[qp]*phi[i][qp]*forcing_function(q_point[qp]);
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");
}
#else
#endif // #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
#endif // #ifdef LIBMESH_ENABLE_AMR
}
The console output of the program:
***************************************************************
* Running Example adaptivity_ex4:
* 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, "Biharmonic"
Type "LinearImplicit"
Variables="u"
Finite Element Types="HERMITE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="THIRD", "THIRD"
n_dofs()=36
n_local_dofs()=24
n_constrained_dofs()=0
n_local_constrained_dofs()=0
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 17.3333
Average Off-Processor Bandwidth <= 6.22222
Maximum On-Processor Bandwidth <= 24
Maximum Off-Processor Bandwidth <= 12
DofMap Constraints
Number of DoF Constraints = 0
Number of Node Constraints = 0
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=25
n_local_nodes()=15
n_elem()=4
n_local_elem()=2
n_active_elem()=4
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
EquationSystems
n_systems()=1
System #0, "Biharmonic"
Type "LinearImplicit"
Variables="u"
Finite Element Types="HERMITE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="THIRD", "THIRD"
n_dofs()=36
n_local_dofs()=24
n_constrained_dofs()=0
n_local_constrained_dofs()=0
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 17.3333
Average Off-Processor Bandwidth <= 6.22222
Maximum On-Processor Bandwidth <= 24
Maximum Off-Processor Bandwidth <= 12
DofMap Constraints
Number of DoF Constraints = 0
Number of Node Constraints = 0
Beginning Solve 0
Linear solver converged at step: 7, final residual: 1.13351e-15
L2-Norm is: 0.384025
H1-Norm is: 1.98976
H2-Norm is: 14.3417
L2-Error is: 0.0335358
H1-Error is: 0.267039
H2-Error is: 3.51162
Refining the mesh...
Mean Error: 2.14746e-14
Error Variance: 1.78051e-29
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=41
n_local_nodes()=21
n_elem()=8
n_local_elem()=4
n_active_elem()=7
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
EquationSystems
n_systems()=1
System #0, "Biharmonic"
Type "LinearImplicit"
Variables="u"
Finite Element Types="HERMITE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="THIRD", "THIRD"
n_dofs()=56
n_local_dofs()=32
n_constrained_dofs()=8
n_local_constrained_dofs()=8
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 21.1429
Average Off-Processor Bandwidth <= 10.8571
Maximum On-Processor Bandwidth <= 32
Maximum Off-Processor Bandwidth <= 32
DofMap Constraints
Number of DoF Constraints = 8
Average DoF Constraint Length= 4
Number of Node Constraints = 9
Maximum Node Constraint Length= 5
Average Node Constraint Length= 3.22222
Beginning Solve 1
Linear solver converged at step: 7, final residual: 2.91791e-12
L2-Norm is: 0.385688
H1-Norm is: 1.99188
H2-Norm is: 14.3773
L2-Error is: 0.0312488
H1-Error is: 0.251755
H2-Error is: 3.36347
Refining the mesh...
Mean Error: 0.904773
Error Variance: 0.231378
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=69
n_local_nodes()=37
n_elem()=16
n_local_elem()=8
n_active_elem()=13
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
EquationSystems
n_systems()=1
System #0, "Biharmonic"
Type "LinearImplicit"
Variables="u"
Finite Element Types="HERMITE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="THIRD", "THIRD"
n_dofs()=88
n_local_dofs()=56
n_constrained_dofs()=8
n_local_constrained_dofs()=8
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 24.5455
Average Off-Processor Bandwidth <= 5.81818
Maximum On-Processor Bandwidth <= 52
Maximum Off-Processor Bandwidth <= 20
DofMap Constraints
Number of DoF Constraints = 8
Average DoF Constraint Length= 4
Number of Node Constraints = 9
Maximum Node Constraint Length= 5
Average Node Constraint Length= 3.22222
Beginning Solve 2
Linear solver converged at step: 14, final residual: 6.63466e-12
L2-Norm is: 0.39603
H1-Norm is: 2.01114
H2-Norm is: 14.5761
L2-Error is: 0.020273
H1-Error is: 0.166208
H2-Error is: 2.36714
Refining the mesh...
Mean Error: 1.07284
Error Variance: 1.26032
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=125
n_local_nodes()=68
n_elem()=32
n_local_elem()=16
n_active_elem()=25
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
EquationSystems
n_systems()=1
System #0, "Biharmonic"
Type "LinearImplicit"
Variables="u"
Finite Element Types="HERMITE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="THIRD", "THIRD"
n_dofs()=152
n_local_dofs()=92
n_constrained_dofs()=32
n_local_constrained_dofs()=20
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 31.6842
Average Off-Processor Bandwidth <= 5.26316
Maximum On-Processor Bandwidth <= 64
Maximum Off-Processor Bandwidth <= 28
DofMap Constraints
Number of DoF Constraints = 32
Average DoF Constraint Length= 4
Number of Node Constraints = 32
Maximum Node Constraint Length= 5
Average Node Constraint Length= 3.5
Beginning Solve 3
Linear solver converged at step: 23, final residual: 3.50547e-12
L2-Norm is: 0.405284
H1-Norm is: 2.03173
H2-Norm is: 14.7515
L2-Error is: 0.00185729
H1-Error is: 0.0267654
H2-Error is: 0.71977
Refining the mesh...
Mean Error: 0.210827
Error Variance: 0.0183284
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=249
n_local_nodes()=133
n_elem()=68
n_local_elem()=35
n_active_elem()=52
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
EquationSystems
n_systems()=1
System #0, "Biharmonic"
Type "LinearImplicit"
Variables="u"
Finite Element Types="HERMITE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="THIRD", "THIRD"
n_dofs()=292
n_local_dofs()=168
n_constrained_dofs()=72
n_local_constrained_dofs()=40
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 32.4384
Average Off-Processor Bandwidth <= 4.38356
Maximum On-Processor Bandwidth <= 60
Maximum Off-Processor Bandwidth <= 32
DofMap Constraints
Number of DoF Constraints = 72
Average DoF Constraint Length= 4
Number of Node Constraints = 75
Maximum Node Constraint Length= 9
Average Node Constraint Length= 3.4
Beginning Solve 4
Linear solver converged at step: 28, final residual: 2.55332e-12
L2-Norm is: 0.405907
H1-Norm is: 2.03286
H2-Norm is: 14.7613
L2-Error is: 0.00113017
H1-Error is: 0.016604
H2-Error is: 0.482316
Refining the mesh...
Mean Error: 0.108787
Error Variance: 0.00898017
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=485
n_local_nodes()=255
n_elem()=140
n_local_elem()=70
n_active_elem()=106
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
EquationSystems
n_systems()=1
System #0, "Biharmonic"
Type "LinearImplicit"
Variables="u"
Finite Element Types="HERMITE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="THIRD", "THIRD"
n_dofs()=548
n_local_dofs()=304
n_constrained_dofs()=120
n_local_constrained_dofs()=68
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 35.0949
Average Off-Processor Bandwidth <= 2.51095
Maximum On-Processor Bandwidth <= 68
Maximum Off-Processor Bandwidth <= 24
DofMap Constraints
Number of DoF Constraints = 120
Average DoF Constraint Length= 4
Number of Node Constraints = 120
Maximum Node Constraint Length= 9
Average Node Constraint Length= 3.5
Beginning Solve 5
Linear solver converged at step: 53, final residual: 1.28271e-11
L2-Norm is: 0.406259
H1-Norm is: 2.03204
H2-Norm is: 14.7673
L2-Error is: 0.000356962
H1-Error is: 0.00613803
H2-Error is: 0.227619
Refining the mesh...
Mean Error: 0.0264705
Error Variance: 0.00117171
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=957
n_local_nodes()=497
n_elem()=280
n_local_elem()=141
n_active_elem()=211
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
EquationSystems
n_systems()=1
System #0, "Biharmonic"
Type "LinearImplicit"
Variables="u"
Finite Element Types="HERMITE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="THIRD", "THIRD"
n_dofs()=1072
n_local_dofs()=576
n_constrained_dofs()=304
n_local_constrained_dofs()=160
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 37.0746
Average Off-Processor Bandwidth <= 2.26866
Maximum On-Processor Bandwidth <= 76
Maximum Off-Processor Bandwidth <= 28
DofMap Constraints
Number of DoF Constraints = 304
Average DoF Constraint Length= 4
Number of Node Constraints = 303
Maximum Node Constraint Length= 9
Average Node Constraint Length= 3.50825
Beginning Solve 6
Linear solver converged at step: 160, final residual: 1.868e-11
L2-Norm is: 0.406309
H1-Norm is: 2.03185
H2-Norm is: 14.7684
L2-Error is: 8.72883e-05
H1-Error is: 0.00241104
H2-Error is: 0.13224
Refining the mesh...
Mean Error: 0.0109873
Error Variance: 0.000156803
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=1885
n_local_nodes()=966
n_elem()=560
n_local_elem()=285
n_active_elem()=421
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
EquationSystems
n_systems()=1
System #0, "Biharmonic"
Type "LinearImplicit"
Variables="u"
Finite Element Types="HERMITE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="THIRD", "THIRD"
n_dofs()=2088
n_local_dofs()=1096
n_constrained_dofs()=616
n_local_constrained_dofs()=316
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 38.1916
Average Off-Processor Bandwidth <= 1.59387
Maximum On-Processor Bandwidth <= 80
Maximum Off-Processor Bandwidth <= 44
DofMap Constraints
Number of DoF Constraints = 616
Average DoF Constraint Length= 4
Number of Node Constraints = 609
Maximum Node Constraint Length= 9
Average Node Constraint Length= 3.52874
Beginning Solve 7
Linear solver converged at step: 257, final residual: 3.94494e-11
L2-Norm is: 0.406333
H1-Norm is: 2.03181
H2-Norm is: 14.7688
L2-Error is: 5.02694e-05
H1-Error is: 0.00140301
H2-Error is: 0.082078
************************************************************************************************************************
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***
************************************************************************************************************************
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------
/workspace/libmesh/examples/adaptivity/adaptivity_ex4/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:38:51 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): 7.459e+00 1.00000 7.459e+00
Objects: 5.210e+02 1.00000 5.210e+02
Flops: 2.424e+08 1.33814 2.117e+08 4.235e+08
Flops/sec: 3.249e+07 1.33814 2.839e+07 5.677e+07
MPI Messages: 6.990e+02 1.00000 6.990e+02 1.398e+03
MPI Message Lengths: 6.476e+05 1.00000 9.265e+02 1.295e+06
MPI Reductions: 1.721e+03 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: 7.4591e+00 100.0% 4.2347e+08 100.0% 1.398e+03 100.0% 9.265e+02 100.0% 1.720e+03 99.9%
------------------------------------------------------------------------------------------------------------------------
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 549 1.0 1.8737e-02 4.8 1.19e+07 1.1 0.0e+00 0.0e+00 5.5e+02 0 5 0 0 32 0 5 0 0 32 1196
VecNorm 579 1.0 2.1019e-03 1.6 8.29e+05 1.1 0.0e+00 0.0e+00 5.8e+02 0 0 0 0 34 0 0 0 0 34 744
VecScale 571 1.0 2.7800e-04 1.2 4.12e+05 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 2797
VecCopy 51 1.0 2.6941e-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
VecSet 676 1.0 2.6894e-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 44 1.0 6.1173e-03 1.0 5.72e+04 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 18
VecMAXPY 571 1.0 3.3858e-03 1.1 1.27e+07 1.1 0.0e+00 0.0e+00 0.0e+00 0 6 0 0 0 0 6 0 0 0 7069
VecAssemblyBegin 87 1.0 6.9883e-02 1.3 0.00e+00 0.0 3.4e+01 4.1e+02 2.0e+02 1 0 2 1 12 1 0 2 1 12 0
VecAssemblyEnd 87 1.0 7.1526e-05 1.4 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 607 1.0 7.4601e-04 1.0 0.00e+00 0.0 1.2e+03 8.4e+02 0.0e+00 0 0 87 78 0 0 0 87 78 0 0
VecScatterEnd 607 1.0 2.2985e-0236.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
VecNormalize 571 1.0 2.5153e-03 1.4 1.24e+06 1.1 0.0e+00 0.0e+00 5.7e+02 0 1 0 0 33 0 1 0 0 33 928
MatMult 571 1.0 3.3731e-02 2.6 3.18e+07 1.2 1.1e+03 8.3e+02 0.0e+00 0 14 82 73 0 0 14 82 73 0 1761
MatSolve 579 1.0 5.1167e-02 1.4 1.56e+08 1.4 0.0e+00 0.0e+00 0.0e+00 1 64 0 0 0 1 64 0 0 0 5287
MatLUFactorNum 8 1.0 1.0812e-02 1.7 2.85e+07 1.8 0.0e+00 0.0e+00 0.0e+00 0 11 0 0 0 0 11 0 0 0 4140
MatILUFactorSym 8 1.0 5.1319e-02 1.5 0.00e+00 0.0 0.0e+00 0.0e+00 2.4e+01 1 0 0 0 1 1 0 0 0 1 0
MatAssemblyBegin 16 1.0 6.1067e-02 2.7 0.00e+00 0.0 3.0e+01 7.9e+03 3.2e+01 1 0 2 18 2 1 0 2 18 2 0
MatAssemblyEnd 16 1.0 1.5497e-03 1.2 0.00e+00 0.0 3.2e+01 1.1e+02 6.4e+01 0 0 2 0 4 0 0 2 0 4 0
MatGetRowIJ 8 1.0 7.1526e-06 1.8 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
MatGetOrdering 8 1.0 2.9230e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 3.2e+01 0 0 0 0 2 0 0 0 0 2 0
MatZeroEntries 24 1.0 1.9360e-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 549 1.0 2.2039e-02 2.9 2.38e+07 1.1 0.0e+00 0.0e+00 5.5e+02 0 11 0 0 32 0 11 0 0 32 2035
KSPSetUp 16 1.0 1.4520e-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 8 1.0 1.4782e-01 1.0 2.42e+08 1.3 1.1e+03 8.3e+02 1.2e+03 2100 82 73 70 2100 82 73 70 2865
PCSetUp 16 1.0 6.3536e-02 1.5 2.85e+07 1.8 0.0e+00 0.0e+00 7.2e+01 1 11 0 0 4 1 11 0 0 4 704
PCSetUpOnBlocks 8 1.0 6.2851e-02 1.6 2.85e+07 1.8 0.0e+00 0.0e+00 5.6e+01 1 11 0 0 3 1 11 0 0 3 712
PCApply 579 1.0 5.4528e-02 1.3 1.56e+08 1.4 0.0e+00 0.0e+00 0.0e+00 1 64 0 0 0 1 64 0 0 0 4961
------------------------------------------------------------------------------------------------------------------------
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 332 332 1317344 0
Vector Scatter 30 30 31080 0
Index Set 86 86 88748 0
IS L to G Mapping 8 8 4512 0
Matrix 32 32 5856368 0
Krylov Solver 16 16 154880 0
Preconditioner 16 16 14272 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 1.19209e-07
Average time for MPI_Barrier(): 6.19888e-07
Average time for zero size MPI_Send(): 1.2517e-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:51 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=7.47175, Active time=7.33729 |
----------------------------------------------------------------------------------------------------------------
| 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() 8 0.0324 0.004045 0.0523 0.006536 0.44 0.71 |
| build_constraint_matrix() 415 0.0133 0.000032 0.0133 0.000032 0.18 0.18 |
| build_sparsity() 8 0.0747 0.009336 0.1267 0.015843 1.02 1.73 |
| cnstrn_elem_mat_vec() 415 0.0080 0.000019 0.0080 0.000019 0.11 0.11 |
| create_dof_constraints() 8 0.1273 0.015908 1.8322 0.229024 1.73 24.97 |
| distribute_dofs() 8 0.0197 0.002459 0.0523 0.006536 0.27 0.71 |
| dof_indices() 4471 0.5271 0.000118 0.5271 0.000118 7.18 7.18 |
| enforce_constraints_exactly() 7 0.0014 0.000202 0.0014 0.000202 0.02 0.02 |
| old_dof_indices() 830 0.0954 0.000115 0.0954 0.000115 1.30 1.30 |
| prepare_send_list() 8 0.0003 0.000038 0.0003 0.000038 0.00 0.00 |
| reinit() 8 0.0311 0.003889 0.0311 0.003889 0.42 0.42 |
| |
| EquationSystems |
| build_solution_vector() 1 0.1902 0.190189 0.2244 0.224413 2.59 3.06 |
| |
| ErrorVector |
| mean() 14 0.0000 0.000004 0.0000 0.000004 0.00 0.00 |
| variance() 7 0.0001 0.000007 0.0001 0.000007 0.00 0.00 |
| |
| ExodusII_IO |
| write_nodal_data() 1 0.0038 0.003758 0.0038 0.003758 0.05 0.05 |
| |
| FE |
| compute_shape_functions() 2757 1.6630 0.000603 1.6630 0.000603 22.67 22.67 |
| init_shape_functions() 2757 3.4243 0.001242 3.4243 0.001242 46.67 46.67 |
| inverse_map() 14656 0.1099 0.000008 0.1099 0.000008 1.50 1.50 |
| |
| FEMap |
| compute_affine_map() 3912 0.0461 0.000012 0.0461 0.000012 0.63 0.63 |
| compute_face_map() 1090 0.0256 0.000023 0.0611 0.000056 0.35 0.83 |
| init_face_shape_functions() 1090 0.0071 0.000006 0.0071 0.000006 0.10 0.10 |
| init_reference_to_physical_map() 2757 0.1267 0.000046 0.1267 0.000046 1.73 1.73 |
| |
| JumpErrorEstimator |
| estimate_error() 7 0.0322 0.004600 1.9852 0.283603 0.44 27.06 |
| |
| LocationMap |
| find() 2780 0.0109 0.000004 0.0109 0.000004 0.15 0.15 |
| init() 14 0.0052 0.000370 0.0052 0.000370 0.07 0.07 |
| |
| Mesh |
| all_second_order() 1 0.0004 0.000352 0.0004 0.000352 0.00 0.00 |
| contract() 7 0.0002 0.000034 0.0007 0.000099 0.00 0.01 |
| find_neighbors() 9 0.0129 0.001432 0.0143 0.001587 0.18 0.19 |
| renumber_nodes_and_elem() 25 0.0015 0.000059 0.0015 0.000059 0.02 0.02 |
| |
| MeshCommunication |
| compute_hilbert_indices() 10 0.0037 0.000374 0.0037 0.000374 0.05 0.05 |
| find_global_indices() 10 0.0020 0.000203 0.0082 0.000820 0.03 0.11 |
| parallel_sort() 10 0.0014 0.000137 0.0017 0.000169 0.02 0.02 |
| |
| MeshOutput |
| write_equation_systems() 1 0.0001 0.000116 0.2283 0.228336 0.00 3.11 |
| |
| MeshRefinement |
| _coarsen_elements() 14 0.0003 0.000022 0.0003 0.000025 0.00 0.00 |
| _refine_elements() 14 0.0123 0.000880 0.0347 0.002480 0.17 0.47 |
| add_point() 2780 0.0098 0.000004 0.0211 0.000008 0.13 0.29 |
| make_coarsening_compatible() 29 0.0063 0.000217 0.0063 0.000217 0.09 0.09 |
| make_refinement_compatible() 29 0.0005 0.000018 0.0006 0.000020 0.01 0.01 |
| |
| MeshTools::Generation |
| build_cube() 1 0.0004 0.000352 0.0004 0.000352 0.00 0.00 |
| |
| MetisPartitioner |
| partition() 9 0.0132 0.001470 0.0212 0.002350 0.18 0.29 |
| |
| Parallel |
| allgather() 74 0.0005 0.000007 0.0006 0.000009 0.01 0.01 |
| max(bool) 73 0.0012 0.000017 0.0012 0.000017 0.02 0.02 |
| max(scalar) 1539 0.1173 0.000076 0.1173 0.000076 1.60 1.60 |
| max(vector) 361 0.0023 0.000006 0.0049 0.000014 0.03 0.07 |
| min(bool) 1879 0.0043 0.000002 0.0043 0.000002 0.06 0.06 |
| min(scalar) 1470 0.1549 0.000105 0.1549 0.000105 2.11 2.11 |
| min(vector) 361 0.0024 0.000007 0.0053 0.000015 0.03 0.07 |
| probe() 86 0.0010 0.000011 0.0010 0.000011 0.01 0.01 |
| receive() 86 0.0004 0.000005 0.0014 0.000017 0.01 0.02 |
| send() 86 0.0002 0.000003 0.0002 0.000003 0.00 0.00 |
| send_receive() 106 0.0009 0.000008 0.0027 0.000026 0.01 0.04 |
| sum() 170 0.0008 0.000005 0.0097 0.000057 0.01 0.13 |
| |
| Parallel::Request |
| wait() 86 0.0001 0.000002 0.0001 0.000002 0.00 0.00 |
| |
| Partitioner |
| set_node_processor_ids() 9 0.0030 0.000337 0.0038 0.000418 0.04 0.05 |
| set_parent_processor_ids() 9 0.0012 0.000132 0.0012 0.000132 0.02 0.02 |
| |
| PetscLinearSolver |
| solve() 8 0.2118 0.026477 0.2118 0.026477 2.89 2.89 |
| |
| ProjectVector |
| operator() 7 0.0607 0.008669 1.1635 0.166217 0.83 15.86 |
| |
| System |
| assemble() 8 0.0754 0.009429 0.5747 0.071839 1.03 7.83 |
| project_vector() 7 0.0579 0.008278 1.2708 0.181547 0.79 17.32 |
----------------------------------------------------------------------------------------------------------------
| Totals: 47403 7.3373 100.00 |
----------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example adaptivity_ex4:
* 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 4 - Biharmonic Equation
This example solves the Biharmonic equation on a square or cube, using a Galerkin formulation with C1 elements approximating the H^2_0 function space. The initial mesh contains two TRI6, one QUAD9 or one HEX27 An input file named "ex15.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 "ex15.in", the variable "max_r_steps" controls the number of refinement steps, and "max_r_level" controls the maximum element refinement level.
LibMesh include files.