// 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

// <h1>FEMSystem Example 4 - Mixed-dimension heat transfer equation
// with FEMSystem</h1>
// \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 <iomanip>

// 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.hpp"
#include "libmesh/diff_solver.h"
#include "libmesh/euler_solver.h"
#include "libmesh/steady_solver.h"

// user defined literal for libmesh's Real
constexpr auto operator"" _r(long double n) { return libMesh::Real(n); }

// The main program.
int main(int argc, char** argv) {
	// Initialize libMesh.
	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() !=
	                             libMesh::EIGEN_SOLVERS,
	                         "--enable-petsc");

	// This doesn't converge without at least double precision
	libmesh_example_requires(sizeof(libMesh::Real) > 4,
	                         "--disable-singleprecision");

	// Parse the input file
	GetPot infile{"../fem_system_ex4.in"};

	// Read in parameters from the input file
	const auto& global_tolerance{infile("global_tolerance", 0._r)};
	const auto& nelem_target{infile("n_elements", 400u)};
	const auto& deltat{infile("deltat", 0.005_r)};
	const auto& coarsegridsize{infile("coarsegridsize", 20u)};
	const auto& coarserefinements{infile("coarserefinements", 0u)};
	const auto& max_adaptivesteps{infile("max_adaptivesteps", 10u)};
	const auto& dim{infile("dimension", 2u)};

	// 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.
	libMesh::Mesh mesh{init.comm()};

	// And an object to refine it
	libMesh::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.
	libMesh::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
	mesh.get_boundary_info().add_elements({bcid}, 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.
	std::for_each(mesh.elements_begin(), mesh.elements_end(), [&](auto elem) {
		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.
	libMesh::EquationSystems equation_systems{mesh};

	// Declare the system "Heat" and its variables.
	auto& system = equation_systems.add_system<HeatSystem>("Heat");

	// Solve this as a steady system
	system.time_solver = libMesh::UniquePtr<libMesh::TimeSolver>(
	    new libMesh::SteadySolver(system));

	// Initialize the system
	equation_systems.init();

	// Set the time stepping options
	system.deltat = deltat;

	// And the nonlinear solver options
	auto& 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", 15u);
	solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3_r);
	solver.relative_residual_tolerance =
	    infile("relative_residual_tolerance", 0.0_r);
	solver.absolute_residual_tolerance =
	    infile("absolute_residual_tolerance", 0.0_r);

	// And the linear solver options
	solver.max_linear_iterations = infile("max_linear_iterations", 50000u);
	solver.initial_linear_tolerance =
	    infile("initial_linear_tolerance", 1.e-3_r);

	// Print information about the system to the screen.
	equation_systems.print_info();

	// Adaptively solve the steady solution
	auto a_step = max_adaptivesteps;
	for (; a_step != max_adaptivesteps; ++a_step) {
		system.solve();

		system.postprocess();

		libMesh::UniquePtr<libMesh::ErrorEstimator> 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);

			auto u = new libMesh::UniformRefinementEstimator;

			// The lid-driven cavity problem isn't in H1, so
			// lets estimate L2 error
			u->error_norm = libMesh::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 libMesh::KellyErrorEstimator);
		}

		libMesh::ErrorVector error;
		error_estimator->estimate_error(system, error);

		// Print out status at each adaptive step.
		libMesh::out << "Adaptive step " << a_step << ": " << '\n';

		const auto& global_error = error.l2_norm();
		if (global_tolerance != 0.)
			libMesh::out << "Global_error = " << global_error << '\n';

		if (global_tolerance != 0.)
			libMesh::out << "Worst element error = " << error.maximum()
			             << ", mean = " << error.mean() << '\n';

		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.\n";
	}
	// Do one last solve if necessary
	if (a_step == max_adaptivesteps) {
		system.solve();

		system.postprocess();
	}

#ifdef LIBMESH_HAVE_EXODUS_API
	libMesh::ExodusII_IO(mesh).write_equation_systems("out.e",
	                                                  equation_systems);
#endif // #ifdef LIBMESH_HAVE_EXODUS_API

#ifdef LIBMESH_HAVE_GMV
	libMesh::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
	libMesh::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)";
	libMesh::ParsedFunction<libMesh::Number> exact_func(exact_str);
	exact_sol.attach_exact_value(0, &exact_func);
	exact_sol.compute_error("Heat", "T");

	libMesh::Number err{exact_sol.l2_error("Heat", "T")};

	// Print out the error value
	libMesh::out << "L2-Error is: " << err << '\n';

	libmesh_assert_less(libmesh_real(err), 2e-3);

#endif // #ifdef LIBMESH_HAVE_FPARSER

#endif // #ifndef LIBMESH_ENABLE_AMR

	// All done.
	return 0;
}