replaced existing simple example with a slightly modified version of libMesh FEM-example 4

libmesh/.clang-format

+Language: Cpp
+DerivePointerAlignment: false
+PointerAlignment: Left
+IndentWidth: 4
+TabWidth: 4
+UseTab: ForIndentation
+ColumnLimit: 80
+AllowShortFunctionsOnASingleLine: true
+AllowShortIfStatementsOnASingleLine: true
+AllowShortLoopsOnASingleLine: true
+AllowShortBlocksOnASingleLine: true
\ No newline at end of file

libmesh/CMakeLists.txt

@@ -2,6 +2,7 @@
 cmake_minimum_required(VERSION 2.8.10 FATAL_ERROR)
 
 find_package(MPI)      # Built-in in CMake
+find_package(Threads)
 
 include(${CMAKE_CURRENT_SOURCE_DIR}/../../cmake/common.cmake)
 
@@ -22,15 +23,15 @@ if (MPI_FOUND)
 	enable_language(CXX)
     include_directories(${MPI_INCLUDE_PATH} ${LIBMESH_INCLUDE_PATH})
     link_directories(${LIBMESH_LIBRARY_PATH}) 
-    add_executable(${NAME} main.cpp)
+    add_executable(${NAME} fem_system_ex4.C heatsystem.C)
     set(CMAKE_BUILD_TYPE RelWithDebInfo)
     if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
     	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native")
     elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
-    	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -xHost -std=c++14")
+    	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -xHost -std=c++11")
     endif()
-    set_target_properties(${NAME} PROPERTIES CXX_STANDARD 14 CXX_STANDARD_REQUIRED YES)
-    target_link_libraries(${NAME} ${LIBMESH_LIBRARY} ${MPI_LIBRARIES})
+    set_target_properties(${NAME} PROPERTIES CXX_STANDARD 11 CXX_STANDARD_REQUIRED YES)
+    target_link_libraries(${NAME} ${LIBMESH_LIBRARY} ${MPI_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT})
     install(TARGETS ${NAME} DESTINATION bin)
     message("** Enabling '${NAME}': with MPI")
 else()

libmesh/README.md

+# README - libMesh Example
+
+## Description
+
+The libMesh library provides a framework for the numerical simulation of partial differential equations using arbitrary unstructured discretizations on serial and parallel platforms. (http://libmesh.github.io/index.html)
+
+This example is a slightly modified version of the official libMesh example "fem_system_ex4" (https://github.com/libMesh/libmesh/tree/master/examples/fem_system/fem_system_ex4). In this example, a heat transfer equation is solved with an FEM System.
+
+The example demonstrates:
+ * Initialization of libMesh
+ * Definition of a equation system
+ * Usage of adaptive mesh refinement
+ * IO: reading and exporting meshes in different file formats
+
+The example is structured as followed:
+ * `fem_system_ex4.C`: The main program, Initialization, IO, refinement
+ * `heatsystem.C|h`: Heat system with laplace heat equation
+ * `meshes/bridge.e|xda`: FEM-Mesh of a bridge
+ 
+ Screenshot of the result:
+ ![Screenshot of Result](hpc_kernel_samples/unstructured_grids/libmesh/meshes/bridge_screenshot.PNG)
+ 
+ ## Release Date
+ 
+ 2016-09-08
+ 
+ ## Version History
+
+ *  2016-09-08 Initial Release on PRACE CodeVault repository
+ 
+ ## Contributors
+
+ * Thomas Steinreiter - [thomas.steinreiter@risc-software.at](mailto:thomas.steinreiter@risc-software.at)
+
+ ## Copyright
+
+This code is available under LGPL, Version 2.1 - see also the license file in the CodeVault root directory.
+
+## Languages
+
+This sample is written in C++ 11.
+
+
+## Parallelisation
+
+This sample uses PETSCs internal parallelisation.
+
+
+## Level of the code sample complexity
+
+Advanced
+
+
+## Compiling
+
+Follow the compilation instructions given in the main directory of the kernel samples directory (`/hpc_kernel_samples/README.md`).
+
+Note: libMesh needs to be built against the correct version of MPI. PETSC is needed for parallelism.
+
+## Running
+
+To run the program, use something similar to
+
+	mpiexec -n [nprocs] ./8_unstructured_libmesh
+
+either on the command line or in your batch script, where `nprocs` specifies the number of processes used. Node: if `nprocs` > 1, libMesh must be built with PETSC enabled.
+ 
+### Example
+ 
+If you run
+	
+	mpiexec -n 8 ./8_unstructured_libmesh
+
+the output should look similar to
+
+	 Mesh Information:
+	  elem_dimensions()={2, 3}
+	  spatial_dimension()=3
+	  n_nodes()=25227
+	    n_local_nodes()=3322
+	  n_elem()=127818
+	    n_local_elem()=15972
+	    n_active_elem()=127818
+	  n_subdomains()=2
+	  n_partitions()=8
+	  n_processors()=8
+	  n_threads()=1
+	  processor_id()=0
+	
+	 EquationSystems
+	  n_systems()=1
+	   System #0, "Heat"
+	    Type "Implicit"
+	    Variables="T"
+	    Finite Element Types="LAGRANGE"
+	    Approximation Orders="FIRST"
+	    n_dofs()=25227
+	    n_local_dofs()=3322
+	    n_constrained_dofs()=1887
+	    n_local_constrained_dofs()=122
+	    n_vectors()=1
+	    n_matrices()=1
+	    DofMap Sparsity
+	      Average  On-Processor Bandwidth <= 13.6686
+	      Average Off-Processor Bandwidth <= 0.533635
+	      Maximum  On-Processor Bandwidth <= 27
+	      Maximum Off-Processor Bandwidth <= 15
+	    DofMap Constraints
+	      Number of DoF Constraints = 1887
+	      Average DoF Constraint Length= 0
+	
+	Assembling the System
+	*** Warning, This code is deprecated, and likely to be removed in future library versions! /usr/local/include/libmesh/libmesh_common.h, line 497, compiled Sep  1 2016 at 15:28:09 ***
+	Nonlinear Residual: 8495.66
+	Linear solve starting, tolerance 0.001
+	Linear solve finished, step 55, residual 3.10027
+	Trying full Newton step
+	  Current Residual: 14.2871
+	  Nonlinear step: |du|/|u| = 1, |du| = 510394
+	Assembling the System
+	Nonlinear Residual: 14.2871
+	Linear solve starting, tolerance 0.001
+	Linear solve finished, step 58, residual 0.00263574
+	Trying full Newton step
+	  Current Residual: 0.012002
+	  Nonlinear step: |du|/|u| = 0.000564056, |du| = 287.935
+	Assembling the System
+	Nonlinear Residual: 0.012002
+	Linear solve starting, tolerance 1.2002e-05
+	Linear solve finished, step 94, residual 2.75721e-08
+	Trying full Newton step
+	  Current Residual: 1.29858e-07
+	  Nonlinear solver converged, step 2, residual reduction 1.52852e-11 < 1e-07
+	  Nonlinear solver relative step size 5.66919e-07 > 1e-07
+	L2-Error is: 59155
\ No newline at end of file

libmesh/fem_system_ex4.C

+// 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.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<boundary_id_type> 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<HeatSystem>("Heat");
+
+	// Solve this as a steady system
+	system.time_solver = UniquePtr<TimeSolver>(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<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);
+
+			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<Number> 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;
+}

libmesh/fem_system_ex4.in

+# The global FEM error tolerance at each timestep
+# Make this nonzero to solve to a specified tolerance
+# This will probably break with KellyErrorIndicator
+# const Real global_tolerance = 1.e-3;
+global_tolerance = 0.0
+
+# The desired number of active mesh elements
+# Make this nonzero to solve to a specified mesh size
+n_elements = 1200
+
+# The coarse grid size from which to start adaptivity
+coarsegridsize = 20
+
+# The maximum number of adaptive steps per timestep
+max_adaptivesteps = 2
+
+# Turn this off to see the NewtonSolver chatter
+solver_quiet = false
+
+# Solve the 2D or 3D problem
+dimension = 3
+
+# Nonlinear solver tolerance
+relative_step_tolerance = 1e-7
+
+relative_residual_tolerance = 1.e-7

libmesh/heatsystem.C

+#include "heatsystem.h"
+
+#include "libmesh/dirichlet_boundaries.h"
+#include "libmesh/dof_map.h"
+#include "libmesh/elem.h"
+#include "libmesh/fe_base.h"
+#include "libmesh/fe_interface.h"
+#include "libmesh/fem_context.h"
+#include "libmesh/getpot.h"
+#include "libmesh/mesh.h"
+#include "libmesh/quadrature.h"
+#include "libmesh/string_to_enum.h"
+#include "libmesh/zero_function.h"
+
+using namespace libMesh;
+
+void HeatSystem::init_data() {
+	T_var = this->add_variable("T", static_cast<Order>(_fe_order),
+	                           Utility::string_to_enum<FEFamily>(_fe_family));
+
+	std::set<boundary_id_type> bdys{0, 1, 2};
+	std::vector<unsigned int> T_only(1, T_var);
+	ZeroFunction<Number> zero;
+
+	this->get_dof_map().add_dirichlet_boundary(
+	    DirichletBoundary(bdys, T_only, &zero));
+
+	// Do the parent's initialization after variables are defined
+	FEMSystem::init_data();
+
+	this->time_evolving(0);
+}
+
+void HeatSystem::init_context(DiffContext& context) {
+	FEMContext& c = libmesh_cast_ref<FEMContext&>(context);
+
+	const std::set<unsigned char>& elem_dims = c.elem_dimensions();
+
+	for (std::set<unsigned char>::const_iterator dim_it = elem_dims.begin();
+	     dim_it != elem_dims.end(); ++dim_it) {
+		const unsigned char dim = *dim_it;
+
+		FEBase* fe = libmesh_nullptr;
+
+		c.get_element_fe(T_var, fe, dim);
+
+		fe->get_JxW();  // For integration
+		fe->get_dphi(); // For bilinear form
+		fe->get_xyz();  // For forcing
+		fe->get_phi();  // For forcing
+	}
+
+	FEMSystem::init_context(context);
+}
+
+bool HeatSystem::element_time_derivative(bool request_jacobian,
+                                         DiffContext& context) {
+	FEMContext& c = libmesh_cast_ref<FEMContext&>(context);
+
+	const unsigned int mesh_dim = c.get_system().get_mesh().mesh_dimension();
+
+	// First we get some references to cell-specific data that
+	// will be used to assemble the linear system.
+	const unsigned int dim = c.get_elem().dim();
+	FEBase* fe = libmesh_nullptr;
+	c.get_element_fe(T_var, fe, dim);
+
+	// Element Jacobian * quadrature weights for interior integration
+	const std::vector<Real>& JxW = fe->get_JxW();
+
+	const std::vector<Point>& xyz = fe->get_xyz();
+
+	const std::vector<std::vector<Real>>& phi = fe->get_phi();
+
+	const std::vector<std::vector<RealGradient>>& dphi = fe->get_dphi();
+
+	// The number of local degrees of freedom in each variable
+	const unsigned int n_T_dofs = c.get_dof_indices(T_var).size();
+
+	// The subvectors and submatrices we need to fill:
+	DenseSubMatrix<Number>& K = c.get_elem_jacobian(T_var, T_var);
+	DenseSubVector<Number>& F = c.get_elem_residual(T_var);
+
+	// Now we will build the element Jacobian and residual.
+	// Constructing the residual requires the solution and its
+	// gradient from the previous timestep.  This must be
+	// calculated at each quadrature point by summing the
+	// solution degree-of-freedom values by the appropriate
+	// weight functions.
+	unsigned int n_qpoints = c.get_element_qrule().n_points();
+
+	for (unsigned int qp = 0; qp != n_qpoints; qp++) {
+		// Compute the solution gradient at the Newton iterate
+		Gradient grad_T = c.interior_gradient(T_var, qp);
+
+		const Number k = _k[dim];
+
+		const Point& p = xyz[qp];
+
+		const Number u_exact = +1.0;
+
+		// Only apply forcing to interior elements
+		const Number forcing =
+		    (dim == mesh_dim) ? -k * u_exact * (dim * libMesh::pi * libMesh::pi)
+		                      : 0;
+
+		const Number JxWxNK = JxW[qp] * -k;
+
+		for (unsigned int i = 0; i != n_T_dofs; i++)
+			F(i) += JxWxNK * (grad_T * dphi[i][qp] + forcing * phi[i][qp]);
+		if (request_jacobian) {
+			const Number JxWxNKxD =
+			    JxWxNK * context.get_elem_solution_derivative();
+
+			for (unsigned int i = 0; i != n_T_dofs; i++)
+				for (unsigned int j = 0; j != n_T_dofs; ++j)
+					K(i, j) += JxWxNKxD * (dphi[i][qp] * dphi[j][qp]);
+		}
+	} // end of the quadrature point qp-loop
+
+	return request_jacobian;
+}

libmesh/heatsystem.h

+#include "libmesh/enum_fe_family.h"
+#include "libmesh/fem_system.h"
+
+// FEMSystem, TimeSolver and  NewtonSolver will handle most tasks,
+// but we must specify element residuals
+class HeatSystem : public libMesh::FEMSystem {
+  public:
+	// Constructor
+	HeatSystem(libMesh::EquationSystems& es, const std::string& name,
+	           const unsigned int number)
+	    : libMesh::FEMSystem(es, name, number), _fe_family("LAGRANGE"),
+	      _fe_order(1) {
+		// Get the conductivity ratios right for both 2D and 3D
+		// benchmarks
+		_k[1] = 1 / libMesh::pi / std::sqrt(libMesh::Real(3));
+		_k[2] = 1;