#include "mesh.h"
#include "equation_systems.h"
#include "linear_implicit_system.h"
#include "gmv_io.h"
#include "fe.h"
#include "quadrature.h"
#include "dense_matrix.h"
#include "dense_vector.h"
#include "sparse_matrix.h"
#include "mesh_generation.h"
#include "mesh_modification.h"
#include "mesh_refinement.h"
#include "error_vector.h"
#include "fourth_error_estimators.h"
#include "getpot.h"
#include "exact_solution.h"
#include "dof_map.h"
#include "numeric_vector.h"
#include "elem.h"
#include "tensor_value.h"
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.
libMesh::init (argc, argv);
#ifndef ENABLE_AMR
std::cerr << "ERROR: This example requires libMesh to be\n"
<< "compiled with AMR support!"
<< std::endl;
return 0;
#else
#ifndef ENABLE_SECOND_DERIVATIVES
std::cerr << "ERROR: This example requires the library to be "
<< "compiled with second derivatives support!"
<< std::endl;
return 0;
#else
{
Parse the input file
GetPot input_file("ex15.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);
We have only defined 2 and 3 dimensional problems
assert (dim == 2 || dim == 3);
Currently only the Hermite cubics give a 3D C^1 basis
assert (dim == 2 || approx_type == "HERMITE");
Create a dim-dimensional mesh.
Mesh mesh (dim);
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 gmv_file = output_file;
gmv_file += ".gmv";
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 gmv 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.
{
Creates 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
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 * 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);
Convenient reference to the system
LinearImplicitSystem& system =
equation_systems.get_system<LinearImplicitSystem>("Biharmonic");
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::cerr << "Mean Error: " << error.mean() <<
std::endl;
std::cerr << "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 ();
}
}
Write out the solution
After solving the system write the solution
to a GMV-formatted plot file.
GMVIO (mesh).write_equation_systems (gmv_file,
equation_systems);
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 libMesh::close ();
#endif // #ifndef ENABLE_SECOND_DERIVATIVES
#endif // #ifndef 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 ENABLE_AMR
#ifdef ENABLE_SECOND_DERIVATIVES
It is a good idea to make sure we are assembling
the proper system.
assert (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 Mesh& 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<unsigned int> 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.start_event("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.stop_event("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.start_event ("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.stop_event ("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.start_event ("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.stop_event ("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 PetscMatrix::add_matrix()
and \p PetscVector::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.start_event ("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.stop_event ("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 ENABLE_SECOND_DERIVATIVES
#endif // #ifdef ENABLE_AMR
}
The program without comments:
#include "mesh.h"
#include "equation_systems.h"
#include "linear_implicit_system.h"
#include "gmv_io.h"
#include "fe.h"
#include "quadrature.h"
#include "dense_matrix.h"
#include "dense_vector.h"
#include "sparse_matrix.h"
#include "mesh_generation.h"
#include "mesh_modification.h"
#include "mesh_refinement.h"
#include "error_vector.h"
#include "fourth_error_estimators.h"
#include "getpot.h"
#include "exact_solution.h"
#include "dof_map.h"
#include "numeric_vector.h"
#include "elem.h"
#include "tensor_value.h"
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)
{
libMesh::init (argc, argv);
#ifndef ENABLE_AMR
std::cerr << "ERROR: This example requires libMesh to be\n"
<< "compiled with AMR support!"
<< std::endl;
return 0;
#else
#ifndef ENABLE_SECOND_DERIVATIVES
std::cerr << "ERROR: This example requires the library to be "
<< "compiled with second derivatives support!"
<< std::endl;
return 0;
#else
{
GetPot input_file("ex15.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);
assert (dim == 2 || dim == 3);
assert (dim == 2 || approx_type == "HERMITE");
Mesh mesh (dim);
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 gmv_file = output_file;
gmv_file += ".gmv";
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
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 * 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);
LinearImplicitSystem& system =
equation_systems.get_system<LinearImplicitSystem>("Biharmonic");
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::cerr << "Mean Error: " << error.mean() <<
std::endl;
std::cerr << "Error Variance: " << error.variance() <<
std::endl;
mesh_refinement.refine_and_coarsen_elements();
}
else
{
mesh_refinement.uniformly_refine(1);
}
equation_systems.reinit ();
}
}
GMVIO (mesh).write_equation_systems (gmv_file,
equation_systems);
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 libMesh::close ();
#endif // #ifndef ENABLE_SECOND_DERIVATIVES
#endif // #ifndef 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 ENABLE_AMR
#ifdef ENABLE_SECOND_DERIVATIVES
assert (system_name == "Biharmonic");
PerfLog perf_log ("Matrix Assembly",false);
const Mesh& 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<unsigned int> 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.start_event("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.stop_event("elem init");
perf_log.start_event ("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.stop_event ("Ke");
{
perf_log.start_event ("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.stop_event ("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.start_event ("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.stop_event ("matrix insertion");
}
#else
#endif // #ifdef ENABLE_SECOND_DERIVATIVES
#endif // #ifdef ENABLE_AMR
}
The console output of the program:
***************************************************************
* Running Example ./ex15-devel
***************************************************************
EquationSystems
n_systems()=1
System "Biharmonic"
Type "LinearImplicit"
Variables="u"
Finite Element Types="HERMITE"
Approximation Orders="THIRD"
n_dofs()=36
n_local_dofs()=36
n_constrained_dofs()=0
n_vectors()=1
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=25
n_elem()=4
n_local_elem()=4
n_active_elem()=4
n_subdomains()=1
n_processors()=1
processor_id()=0
EquationSystems
n_systems()=1
System "Biharmonic"
Type "LinearImplicit"
Variables="u"
Finite Element Types="HERMITE"
Approximation Orders="THIRD"
n_dofs()=36
n_local_dofs()=36
n_constrained_dofs()=0
n_vectors()=1
Beginning Solve 0
Linear solver converged at step: 31, final residual: 6.00383e-19
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...
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=81
n_elem()=20
n_local_elem()=20
n_active_elem()=16
n_subdomains()=1
n_processors()=1
processor_id()=0
EquationSystems
n_systems()=1
System "Biharmonic"
Type "LinearImplicit"
Variables="u"
Finite Element Types="HERMITE"
Approximation Orders="THIRD"
n_dofs()=100
n_local_dofs()=100
n_constrained_dofs()=0
n_vectors()=1
Beginning Solve 1
Linear solver converged at step: 21, final residual: 1.78941e-17
L2-Norm is: 0.404988
H1-Norm is: 2.02995
H2-Norm is: 14.7459
L2-Error is: 0.0020746
H1-Error is: 0.0316727
H2-Error is: 0.822125
Refining the mesh...
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=289
n_elem()=84
n_local_elem()=84
n_active_elem()=64
n_subdomains()=1
n_processors()=1
processor_id()=0
EquationSystems
n_systems()=1
System "Biharmonic"
Type "LinearImplicit"
Variables="u"
Finite Element Types="HERMITE"
Approximation Orders="THIRD"
n_dofs()=324
n_local_dofs()=324
n_constrained_dofs()=0
n_vectors()=1
Beginning Solve 2
Linear solver converged at step: 22, final residual: 4.7073e-17
L2-Norm is: 0.406264
H1-Norm is: 2.03164
H2-Norm is: 14.7676
L2-Error is: 0.000129445
H1-Error is: 0.00390589
H2-Error is: 0.202531
Refining the mesh...
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=1089
n_elem()=340
n_local_elem()=340
n_active_elem()=256
n_subdomains()=1
n_processors()=1
processor_id()=0
EquationSystems
n_systems()=1
System "Biharmonic"
Type "LinearImplicit"
Variables="u"
Finite Element Types="HERMITE"
Approximation Orders="THIRD"
n_dofs()=1156
n_local_dofs()=1156
n_constrained_dofs()=0
n_vectors()=1
Beginning Solve 3
Linear solver converged at step: 48, final residual: 1.86903e-16
L2-Norm is: 0.406344
H1-Norm is: 2.03174
H2-Norm is: 14.7689
L2-Error is: 8.07721e-06
H1-Error is: 0.000486566
H2-Error is: 0.050454
***************************************************************
* Done Running Example ./ex15-devel
***************************************************************
Example 15 - 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.