// The libMesh Finite Element Library.
// Copyright (C) 2002-2016 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
FEMSystem Example 4 - Mixed-dimension heat transfer equation
// with FEMSystem
// \author Roy Stogner
// \date 2015
//
// This example shows how elements of multiple dimensions can be
// linked and computed upon simultaneously within the
// DifferentiableSystem class framework
// C++ includes
#include
// Basic include files
#include "libmesh/elem.h"
#include "libmesh/equation_systems.h"
#include "libmesh/error_vector.h"
#include "libmesh/exact_solution.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/getpot.h"
#include "libmesh/gmv_io.h"
#include "libmesh/kelly_error_estimator.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/parsed_function.h"
#include "libmesh/uniform_refinement_estimator.h"
// The systems and solvers we may use
#include "heatsystem.h"
#include "libmesh/diff_solver.h"
#include "libmesh/euler_solver.h"
#include "libmesh/steady_solver.h"
// Bring in everything from the libMesh namespace
using namespace libMesh;
// The main program.
int main(int argc, char** argv) {
// Initialize libMesh.
LibMeshInit init(argc, argv);
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_requires(false, "--enable-amr");
#else
// This doesn't converge with Eigen BICGSTAB for some reason...
libmesh_example_requires(libMesh::default_solver_package() != EIGEN_SOLVERS,
"--enable-petsc");
// This doesn't converge without at least double precision
libmesh_example_requires(sizeof(Real) > 4, "--disable-singleprecision");
// Parse the input file
GetPot infile("../fem_system_ex4.in");
// Read in parameters from the input file
const Real global_tolerance = infile("global_tolerance", 0.);
const unsigned int nelem_target = infile("n_elements", 400);
const Real deltat = infile("deltat", 0.005);
const unsigned int coarsegridsize = infile("coarsegridsize", 20);
const unsigned int coarserefinements = infile("coarserefinements", 0);
const unsigned int max_adaptivesteps = infile("max_adaptivesteps", 10);
const unsigned int dim = infile("dimension", 2);
// Skip higher-dimensional examples on a lower-dimensional libMesh build
libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
// We have only defined 2 and 3 dimensional problems
libmesh_assert(dim == 2 || dim == 3);
// Create a mesh, with dimension to be overridden later, distributed
// across the default MPI communicator.
Mesh mesh(init.comm());
// And an object to refine it
MeshRefinement mesh_refinement(mesh);
mesh_refinement.coarsen_by_parents() = true;
mesh_refinement.absolute_global_tolerance() = global_tolerance;
mesh_refinement.nelem_target() = nelem_target;
mesh_refinement.refine_fraction() = 0.3;
mesh_refinement.coarsen_fraction() = 0.3;
mesh_refinement.coarsen_threshold() = 0.1;
// Use the MeshTools::Generation mesh generator to create a uniform
// grid on the square or cube. We crop the domain at y=2/3 to allow
// for a homogeneous Neumann BC in our benchmark there.
boundary_id_type bcid = 3; // +y in 3D
mesh.read("../meshes/bridge.xda");
{
// Add boundary elements corresponding to the +y boundary of our
// volume mesh
std::set bcids;
bcids.insert(bcid);
mesh.get_boundary_info().add_elements(bcids, mesh);
mesh.prepare_for_use();
}
// To work around ExodusII file format limitations, we need elements
// of different dimensionality to belong to different subdomains.
// Our interior elements defaulted to subdomain id 0, so we'll set
// boundary elements to subdomain 2.
{
const MeshBase::element_iterator end_el = mesh.elements_end();
for (MeshBase::element_iterator el = mesh.elements_begin();
el != end_el; ++el) {
Elem* elem = *el;
if (elem->dim() < dim) elem->subdomain_id() = 2;
}
}
mesh_refinement.uniformly_refine(coarserefinements);
// Print information about the mesh to the screen.
mesh.print_info();
// Create an equation systems object.
EquationSystems equation_systems(mesh);
// Declare the system "Heat" and its variables.
HeatSystem& system = equation_systems.add_system("Heat");
// Solve this as a steady system
system.time_solver = UniquePtr(new SteadySolver(system));
// Initialize the system
equation_systems.init();
// Set the time stepping options
system.deltat = deltat;
// And the nonlinear solver options
DiffSolver& solver = *(system.time_solver->diff_solver().get());
solver.quiet = infile("solver_quiet", true);
solver.verbose = !solver.quiet;
solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 15);
solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3);
solver.relative_residual_tolerance =
infile("relative_residual_tolerance", 0.0);
solver.absolute_residual_tolerance =
infile("absolute_residual_tolerance", 0.0);
// And the linear solver options
solver.max_linear_iterations = infile("max_linear_iterations", 50000);
solver.initial_linear_tolerance = infile("initial_linear_tolerance", 1.e-3);
// Print information about the system to the screen.
equation_systems.print_info();
// Adaptively solve the steady solution
unsigned int a_step = max_adaptivesteps;
for (; a_step != max_adaptivesteps; ++a_step) {
system.solve();
system.postprocess();
ErrorVector error;
UniquePtr error_estimator;
// To solve to a tolerance in this problem we
// need a better estimator than Kelly
if (global_tolerance != 0.) {
// We can't adapt to both a tolerance and a mesh
// size at once
libmesh_assert_equal_to(nelem_target, 0);
UniformRefinementEstimator* u = new UniformRefinementEstimator;
// The lid-driven cavity problem isn't in H1, so
// lets estimate L2 error
u->error_norm = L2;
error_estimator.reset(u);
} else {
// If we aren't adapting to a tolerance we need a
// target mesh size
libmesh_assert_greater(nelem_target, 0);
// Kelly is a lousy estimator to use for a problem
// not in H1 - if we were doing more than a few
// timesteps we'd need to turn off or limit the
// maximum level of our adaptivity eventually
error_estimator.reset(new KellyErrorEstimator);
}
error_estimator->estimate_error(system, error);
// Print out status at each adaptive step.
Real global_error = error.l2_norm();
libMesh::out << "Adaptive step " << a_step << ": " << std::endl;
if (global_tolerance != 0.)
libMesh::out << "Global_error = " << global_error << std::endl;
if (global_tolerance != 0.)
libMesh::out << "Worst element error = " << error.maximum()
<< ", mean = " << error.mean() << std::endl;
if (global_tolerance != 0.) {
// If we've reached our desired tolerance, we
// don't need any more adaptive steps
if (global_error < global_tolerance) break;
mesh_refinement.flag_elements_by_error_tolerance(error);
} else {
// If flag_elements_by_nelem_target returns true, this
// should be our last adaptive step.
if (mesh_refinement.flag_elements_by_nelem_target(error)) {
mesh_refinement.refine_and_coarsen_elements();
equation_systems.reinit();
a_step = max_adaptivesteps;
break;
}
}
// Carry out the adaptive mesh refinement/coarsening
mesh_refinement.refine_and_coarsen_elements();
equation_systems.reinit();
libMesh::out << "Refined mesh to " << mesh.n_active_elem()
<< " active elements and "
<< equation_systems.n_active_dofs() << " active dofs."
<< std::endl;
}
// Do one last solve if necessary
if (a_step == max_adaptivesteps) {
system.solve();
system.postprocess();
}
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems);
#endif // #ifdef LIBMESH_HAVE_EXODUS_API
#ifdef LIBMESH_HAVE_GMV
GMVIO(mesh).write_equation_systems("out.gmv", equation_systems);
#endif // #ifdef LIBMESH_HAVE_GMV
#ifdef LIBMESH_HAVE_FPARSER
// Check that we got close to the analytic solution
ExactSolution exact_sol(equation_systems);
const std::string exact_str =
(dim == 2) ? "sin(pi*x)*sin(pi*y)" : "sin(pi*x)*sin(pi*y)*sin(pi*z)";
ParsedFunction exact_func(exact_str);
exact_sol.attach_exact_value(0, &exact_func);
exact_sol.compute_error("Heat", "T");
Number err = exact_sol.l2_error("Heat", "T");
// Print out the error value
libMesh::out << "L2-Error is: " << err << std::endl;
libmesh_assert_less(libmesh_real(err), 2e-3);
#endif // #ifdef LIBMESH_HAVE_FPARSER
#endif // #ifndef LIBMESH_ENABLE_AMR
// All done.
return 0;
}