Add a mixed-dimension-element example code

Along with some verification that it's giving the right result.

examples/fem_system/fem_system_ex4/Makefile.am

@@ -0,0 +1,21 @@
+example_name  = fem_system_ex4
+check_SCRIPTS = run.sh
+install_dir   = $(examples_install_path)/fem_system/ex1
+data          = fem_system_ex4.C fem_system_ex4.in heatsystem.C heatsystem.h run.sh
+sources       = $(data) run.sh
+
+CLEANFILES = out_*.e
+
+# also need to link files for VPATH builds
+if LIBMESH_VPATH_BUILD
+  BUILT_SOURCES = .linkstamp
+.linkstamp:
+	-rm -f fem_system_ex4.in && $(LN_S) -f $(srcdir)/fem_system_ex4.in .
+	$(AM_V_GEN)touch .linkstamp
+
+  CLEANFILES += fem_system_ex4.in .linkstamp
+endif
+
+##############################################
+# include common example environment
+include $(top_srcdir)/examples/Make.common

examples/fem_system/fem_system_ex4/fem_system_ex4.C

@@ -0,0 +1,302 @@
+/* The libMesh Finite Element Library. */
+/* Copyright (C) 2003  Benjamin S. Kirk */
+
+/* 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>
+//
+// 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/equation_systems.h"
+#include "libmesh/error_vector.h"
+#include "libmesh/exact_solution.h"
+#include "libmesh/getpot.h"
+#include "libmesh/gmv_io.h"
+#include "libmesh/exodusII_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
+
+  // 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
+  if (dim == 2)
+    {
+      MeshTools::Generation::build_square
+        (mesh,
+         coarsegridsize,
+         coarsegridsize*2/3, // int arithmetic best we can do here
+         0., 1.,
+         0., 2./3.,
+         QUAD9);
+      bcid = 2; // +y in 2D
+    }
+  else if (dim == 3)
+    {
+      MeshTools::Generation::build_cube
+        (mesh,
+         coarsegridsize,
+         coarsegridsize*2/3,
+         coarsegridsize,
+         0., 1.,
+         0., 2./3.,
+         0., 1.,
+         HEX27);
+    }
+
+  {
+    // 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();
+  }
+
+  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 = 0;
+  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();
+      std::cout << "Adaptive step " << a_step << ": " << std::endl;
+      if (global_tolerance != 0.)
+        std::cout << "Global_error = " << global_error
+                  << std::endl;
+      if (global_tolerance != 0.)
+        std::cout << "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();
+
+          std::cout << "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
+  // We want to write the file in the ExodusII format, but we don't
+  // yet support mixed-dimension meshes there?
+/*  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
+  std::cout << "L2-Error is: " << err << std::endl;
+
+  libmesh_assert_less(err, 2e-3);
+
+#endif // #ifdef LIBMESH_HAVE_FPARSER
+
+#endif // #ifndef LIBMESH_ENABLE_AMR
+
+  // All done.
+  return 0;
+}

examples/fem_system/fem_system_ex4/fem_system_ex4.in

@@ -0,0 +1,26 @@
+# 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 = 400
+
+# The coarse grid size from which to start adaptivity
+coarsegridsize = 20
+
+# The maximum number of adaptive steps per timestep
+max_adaptivesteps = 0
+
+# Turn this off to see the NewtonSolver chatter
+solver_quiet = false
+
+# Solve the 2D or 3D problem
+dimension = 2
+
+# Nonlinear solver tolerance
+relative_step_tolerance = 1e-9
+
+relative_residual_tolerance = 1.e-9