Some mods for compiling

Const parameter for mesh copy assignment

.gitlab/lassen.yml

@@ -30,5 +30,5 @@
 
 opt_mpi_cuda_xl_16_1_1_8:
   variables:
-    SPEC: "%xl@16.1.1.8 +mpi +cuda cuda_arch=sm_70"
+    SPEC: "%xl@16.1.1.8 +mpi +cuda cuda_arch=70"
   extends: .build_and_test_on_lassen

CHANGELOG

@@ -10,6 +10,7 @@
 
 Version 4.3.1 (development)
 ===========================
+- Added support for hr-adaptivity using TMOP-based error estimator.
 
 
 Version 4.3, released on July 29, 2021

CMakeLists.txt

@@ -267,15 +267,6 @@ if (MFEM_USE_SUNDIALS)
   find_package(SUNDIALS REQUIRED ${SUNDIALS_COMPONENTS})
 endif()
 
-# EPIC
-if (MFEM_USE_EPIC)
-  if (NOT (MFEM_USE_MPI AND MFEM_USE_SUNDIALS AND MFEM_USE_LAPACK) )
-    message(FATAL_ERROR " *** EPIC requires that MPI, SUNDIALS and LAPACK be enabled.")
-  else()
-    find_package(EPIC REQUIRED SUNDIALS NVector_Serial NVector_Parallel BLAS LAPACK)
-  endif()
-endif()
-
 # Mesquite
 if (MFEM_USE_MESQUITE)
   find_package(Mesquite REQUIRED)
@@ -437,7 +428,7 @@ endif()
 #    With newer versions of SuiteSparse which include METIS header using 64-bit
 #    integers, the METIS header (with 32-bit indices, as used by mfem) needs to
 #    be before SuiteSparse.
-set(MFEM_TPLS MPI_CXX OPENMP HYPRE BLAS LAPACK SuperLUDist METIS SuiteSparse SUNDIALS EPIC PETSC
+set(MFEM_TPLS MPI_CXX OPENMP HYPRE BLAS LAPACK SuperLUDist METIS SuiteSparse SUNDIALS PETSC
     SLEPC MESQUITE MUMPS STRUMPACK AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB NETCDF
     MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE ADIOS2
     CUSPARSE MKL_CPARDISO AMGX CALIPER)

config/cmake/MFEMConfig.cmake.in

@@ -29,7 +29,6 @@ set(MFEM_USE_LEGACY_OPENMP @MFEM_USE_LEGACY_OPENMP@)
 set(MFEM_USE_MEMALLOC @MFEM_USE_MEMALLOC@)
 set(MFEM_TIMER_TYPE @MFEM_TIMER_TYPE@)
 set(MFEM_USE_SUNDIALS @MFEM_USE_SUNDIALS@)
-set(MFEM_USE_EPIC @MFEM_USE_EPIC@)
 set(MFEM_USE_MESQUITE @MFEM_USE_MESQUITE@)
 set(MFEM_USE_SUITESPARSE @MFEM_USE_SUITESPARSE@)
 set(MFEM_USE_SUPERLU @MFEM_USE_SUPERLU@)

config/cmake/config.hpp.in

@@ -165,9 +165,6 @@
 // Enable MFEM functionality based on the SUNDIALS libraries.
 #cmakedefine MFEM_USE_SUNDIALS
 
-// Enable MFEM functionality based on the EPIC libraries.
-#cmakedefine MFEM_USE_EPIC
-
 // Version of HYPRE used for building MFEM.
 #cmakedefine MFEM_HYPRE_VERSION @MFEM_HYPRE_VERSION@
 

config/cmake/modules/FindEPIC.cmake

@@ -1,21 +0,0 @@
-# Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
-# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-# LICENSE and NOTICE for details. LLNL-CODE-806117.
-#
-# This file is part of the MFEM library. For more information and source code
-# availability visit https://mfem.org.
-#
-# MFEM is free software; you can redistribute it and/or modify it under the
-# terms of the BSD-3 license. We welcome feedback and contributions, see file
-# CONTRIBUTING.md for details.
-
-# Defines the following variables:
-#   - EPIC_FOUND
-#   - EPIC_LIBRARIES
-#   - EPIC_INCLUDE_DIRS
-
-include(MfemCmakeUtilities)
-mfem_find_package(EPIC EPIC EPIC_DIR
-       "include" Epic.h "lib" epic1.0.0
-  "Paths to headers required by EPIC." "Libraries required by EPIC.")
-

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -759,7 +759,7 @@ function(mfem_export_mk_files)
   set(CONFIG_MK_BOOL_VARS MFEM_USE_MPI MFEM_USE_METIS MFEM_USE_METIS_5
       MFEM_DEBUG MFEM_USE_EXCEPTIONS MFEM_USE_ZLIB MFEM_USE_LIBUNWIND
       MFEM_USE_LAPACK MFEM_THREAD_SAFE MFEM_USE_OPENMP MFEM_USE_LEGACY_OPENMP
-      MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS MFEM_USE_EPIC MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE
+      MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE
       MFEM_USE_SUPERLU MFEM_USE_STRUMPACK MFEM_USE_GINKGO MFEM_USE_AMGX
       MFEM_USE_GNUTLS MFEM_USE_GSLIB MFEM_USE_NETCDF MFEM_USE_PETSC
       MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_CONDUIT MFEM_USE_PUMI

config/config.hpp.in

@@ -85,9 +85,6 @@
 // Enable MFEM functionality based on the SUNDIALS libraries.
 // #define MFEM_USE_SUNDIALS
 
-// Enable MFEM functionality based on the EPIC libraries.
-// #define MFEM_USE_EPIC
-
 // Enable MFEM functionality based on the Mesquite library.
 // #define MFEM_USE_MESQUITE
 

config/config.mk.in

@@ -29,7 +29,6 @@ MFEM_USE_OPENMP        = @MFEM_USE_OPENMP@
 MFEM_USE_MEMALLOC      = @MFEM_USE_MEMALLOC@
 MFEM_TIMER_TYPE        = @MFEM_TIMER_TYPE@
 MFEM_USE_SUNDIALS      = @MFEM_USE_SUNDIALS@
-MFEM_USE_EPIC          = @MFEM_USE_EPIC@
 MFEM_USE_MESQUITE      = @MFEM_USE_MESQUITE@
 MFEM_USE_SUITESPARSE   = @MFEM_USE_SUITESPARSE@
 MFEM_USE_SUPERLU       = @MFEM_USE_SUPERLU@

config/defaults.cmake

@@ -30,7 +30,6 @@ option(MFEM_USE_OPENMP "Enable the OpenMP backend" OFF)
 option(MFEM_USE_LEGACY_OPENMP "Enable legacy OpenMP usage" OFF)
 option(MFEM_USE_MEMALLOC "Enable the internal MEMALLOC option." ON)
 option(MFEM_USE_SUNDIALS "Enable SUNDIALS usage" OFF)
-option(MFEM_USE_EPIC "Enable EPIC usage" OFF)
 option(MFEM_USE_MESQUITE "Enable MESQUITE usage" OFF)
 option(MFEM_USE_SUITESPARSE "Enable SuiteSparse usage" OFF)
 option(MFEM_USE_SUPERLU "Enable SuperLU_DIST usage" OFF)
@@ -116,9 +115,6 @@ set(SUNDIALS_DIR "${MFEM_DIR}/../sundials-5.0.0/instdir" CACHE PATH
 # set(SUNDIALS_REQUIRED_PACKAGES "SuiteSparse/KLU/AMD/BTF/COLAMD/config"
 #     CACHE STRING "Additional packages required by SUNDIALS.")
 
-set(EPIC_DIR "${MFEM_DIR}/../epic-cpp/instdir" CACHE PATH
-    "Path to the EPIC library.")
-
 set(MESQUITE_DIR "${MFEM_DIR}/../mesquite-2.99" CACHE PATH
     "Path to the Mesquite library.")
 

config/defaults.mk

@@ -122,7 +122,6 @@ MFEM_USE_LEGACY_OPENMP = NO
 MFEM_USE_MEMALLOC      = YES
 MFEM_TIMER_TYPE        = $(if $(NOTMAC),2,4)
 MFEM_USE_SUNDIALS      = NO
-MFEM_USE_EPIC          = NO
 MFEM_USE_MESQUITE      = NO
 MFEM_USE_SUITESPARSE   = NO
 MFEM_USE_SUPERLU       = NO
@@ -232,11 +231,6 @@ endif
 # If SUNDIALS was built with KLU:
 # MFEM_USE_SUITESPARSE = YES
 
-# EPIC library configuration
-MESQUITE_DIR = @MFEM_DIR@/../epic-cpp/instdir
-MESQUITE_OPT = -I$(EPIC_DIR)/include
-MESQUITE_LIB = -L$(EPIC_DIR)/lib -lepic1.0.0
-
 # MESQUITE library configuration
 MESQUITE_DIR = @MFEM_DIR@/../mesquite-2.99
 MESQUITE_OPT = -I$(MESQUITE_DIR)/include

examples/CMakeLists.txt

@@ -159,11 +159,6 @@ if (MFEM_USE_AMGX)
   add_subdirectory(amgx)
 endif()
 
-# Include the examples/epic directory if EPIC is enabled.
-if (MFEM_USE_EPIC)
-  add_subdirectory(epic)
-endif()
-
 # Include the examples/ginkgo directory if GINKGO is enabled.
 if (MFEM_USE_GINKGO)
   add_subdirectory(ginkgo)

examples/epic/CMakeLists.txt

@@ -1,64 +0,0 @@
-# Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
-# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-# LICENSE and NOTICE for details. LLNL-CODE-806117.
-#
-# This file is part of the MFEM library. For more information and source code
-# availability visit https://mfem.org.
-#
-# MFEM is free software; you can redistribute it and/or modify it under the
-# terms of the BSD-3 license. We welcome feedback and contributions, see file
-# CONTRIBUTING.md for details.
-
-set(EPIC_EXAMPLES_SRCS)
-list(APPEND EPIC_EXAMPLES_SRCS
-  ex16.cpp
-  )
-
-if (MFEM_USE_MPI)
-list(APPEND EPIC_EXAMPLES_SRCS
-    ex16p.cpp
-    )
-endif()
-
-# Include the source directory where mfem.hpp and mfem-performance.hpp are.
-include_directories(BEFORE ${PROJECT_BINARY_DIR})
-
-# Add "test_epic" target, see below.
-add_custom_target(test_epic
-  ${CMAKE_CTEST_COMMAND} -R epic USES_TERMINAL)
-
-# Add one executable per cpp file, adding "epic_" as prefix. Sets
-# "test_epic" as a target that depends on the given examples.
-set(PFX epic_)
-add_mfem_examples(EPIC_EXAMPLES_SRCS ${PFX} "" test_epic)
-
-# Testing.
-# The EPIC tests can be run separately using the target "test_epic"
-# which builds the examples and runs:
-#   ctest -R epic
-
-# Example 16: use the default options
-
-# Add the tests: one test per source file.
-foreach(SRC_FILE ${EPIC_EXAMPLES_SRCS})
-  get_filename_component(SRC_FILENAME ${SRC_FILE} NAME)
-  string(REPLACE ".cpp" "" TEST_NAME ${SRC_FILENAME})
-  string(TOUPPER ${TEST_NAME} UP_TEST_NAME)
-  set(TEST_NAME ${PFX}${TEST_NAME})
-
-  set(THIS_TEST_OPTIONS "-no-vis")
-  list(APPEND THIS_TEST_OPTIONS ${${UP_TEST_NAME}_TEST_OPTS})
-  # message(STATUS "Test ${TEST_NAME} options: ${THIS_TEST_OPTIONS}")
-
-  if (NOT (${TEST_NAME} MATCHES ".*p$"))
-    add_test(NAME ${TEST_NAME}_ser
-      COMMAND ${TEST_NAME} ${THIS_TEST_OPTIONS})
-  else()
-    add_test(NAME ${TEST_NAME}_np=4
-      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
-      ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:${TEST_NAME}> ${THIS_TEST_OPTIONS}
-      ${MPIEXEC_POSTFLAGS})
-  endif()
-endforeach()
-

examples/epic/README

@@ -1,17 +0,0 @@
-                    Finite Element Discretization Library
-                                   __
-                       _ __ ___   / _|  ___  _ __ ___
-                      | '_ ` _ \ | |_  / _ \| '_ ` _ \
-                      | | | | | ||  _||  __/| | | | | |
-                      |_| |_| |_||_|   \___||_| |_| |_|
-
-                               http://mfem.org
-
-This directory contains modifications of the example codes that illustrate the
-use of MFEM features based on the EPIC suite of time integration.
-
-To build these examples, make sure that MFEM is configured with the option
-"MFEM_USE_EPIC = YES".
-
-We recommend comparing the original example codes with the corresponding files
-in the current directory.

examples/epic/ex16.cpp

@@ -1,610 +0,0 @@
-//                                MFEM Example 16
-//                               EPIC Modification
-//
-// Compile with: make ex16
-//
-// Sample runs:  ex16
-//               ex16 -m ../../data/inline-tri.mesh
-//               ex16 -m ../../data/disc-nurbs.mesh -tf 2
-//               ex16 -s 8 -a 1.0 -k 0.0 -dt 1e-4 -tf 5e-2 -vs 25
-//               ex16 -m ../../data/fichera-q2.mesh
-//               ex16 -m ../../data/escher.mesh
-//               ex16 -m ../../data/beam-tet.mesh -tf 10 -dt 0.1
-//               ex16 -m ../../data/amr-quad.mesh -o 4 -r 0
-//               ex16 -m ../../data/amr-hex.mesh -o 2 -r 0
-//
-// Description:  This example solves a time dependent nonlinear heat equation
-//               problem of the form du/dt = C(u), with a non-linear diffusion
-//               operator C(u) = \nabla \cdot (\kappa + \alpha u) \nabla u.
-//
-//               We recommend viewing examples 2, 9 and 10 before viewing this
-//               example.
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-class ImplicitSolveOperator;
-class JacobianOperator;
-
-/** After spatial discretization, the conduction model can be written as:
- *
- *     du/dt = M^{-1}(-K(u) u)
- *
- *  where u is the vector representing the temperature, M is the mass matrix,
- *  and K is the diffusion operator with diffusivity depending on u:
- *  (\kappa + \alpha u).
- *
- *  Class ConductionOperator represents the right-hand side of the above ODE.
- */
-class ConductionOperator : public TimeDependentOperator
-{
-protected:
-   FiniteElementSpace &fespace;
-   Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
-
-   BilinearForm *M;
-   mutable BilinearForm *K;
-   mutable BilinearForm *dK;
-   mutable BilinearForm *J_K;
-
-   SparseMatrix Mmat;
-   mutable SparseMatrix J_K_mat;
-
-   mutable CGSolver M_solver; // Krylov solver for inverting the mass matrix M
-   DSmoother M_prec;  // Preconditioner for the mass matrix M
-
-   CGSolver Jg_solver; // Krylov solver for inverting the Jacobian in the nonlinear solve
-   DSmoother Jg_prec;  // Preconditioner for the Jacobian Jg
-
-   NewtonSolver newton_solver;
-   mutable JacobianOperator *jac;
-
-   double alpha, kappa;
-
-   mutable Vector z; // auxiliary vector
-
-   mutable int nRhsMult, nSetJac, nJacMult, nImpSolve, nImpIter, nImpMult, nImpSet;
-
-public:
-   Vector u0;
-
-   ConductionOperator(FiniteElementSpace &f, double alpha, double kappa, const Vector &u);
-
-   void UpdateStats();
-   void PrintStats(ostream& out);
-
-   void ExtractJacobians(const Vector& x, std::ostream &out, std::ostream &out2);
-
-   BilinearForm& GetKLambda(const Vector& u) const;
-   BilinearForm& GetdKLambda(const Vector& u) const;
-
-   virtual void Mult(const Vector &u, Vector &du_dt) const;
-   virtual Operator& GetGradient(const Vector &k) const;
-
-   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
-
-   virtual ~ConductionOperator();
-};
-
-class ImplicitSolveOperator : public Operator
-{
-private:
-   double dt;
-   const Vector* x;
-   ConductionOperator* oper;
-
-   const SparseMatrix* M;
-   mutable SparseMatrix* Jg;
-
-   mutable Vector u, z;
-   mutable int nMult, nSet;
-
-public:
-   ImplicitSolveOperator(ConductionOperator* oper, const SparseMatrix* M, double dt, const Vector* x);
-
-   int GetnMult() { return nMult; }
-   int GetnSet() { return nSet; }
-   virtual void Mult(const Vector &k, Vector &gk) const;
-   virtual Operator &GetGradient(const Vector &k) const;
-};
-
-class JacobianOperator : public Operator
-{
-private:
-   Operator* J;
-   Operator* M_solver;
-
-   mutable int nMult;
-   mutable Vector z;
-public:
-   JacobianOperator(Operator* J, Operator* M_solver);
-
-   int GetnMult() { return nMult; }
-
-   void ExtractJacobian(const Vector& x, std::ostream &out);
-   virtual void Mult(const Vector &k, Vector &gk) const;
-};
-
-double InitialTemperature(const Vector &x);
-
-int main(int argc, char *argv[])
-{
-   // 1. Parse command-line options.
-   const char *mesh_file = "../../data/star.mesh";
-   int ref_levels = 2;
-   int order = 2;
-   int ode_solver_type = 8; // Exponential Euler
-   double t_final = 0.5;
-   double dt = 1.0e-2;
-   double alpha = 1.0e-2;
-   double kappa = 0.5;
-   bool visualization = true;
-   bool visit = false;
-   int vis_steps = 5;
-
-   int precision = 8;
-   cout.precision(precision);
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ref_levels, "-r", "--refine",
-                  "Number of times to refine the mesh uniformly.");
-   args.AddOption(&order, "-o", "--order",
-                  "Order (degree) of the finite elements.");
-   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
-                  "ODE solver:\n\t"
-                  "1  - Forward Euler,\n\t"
-                  "2  - RK2,\n\t"
-                  "3  - RK3 SSP,\n\t"
-                  "4  - RK4,\n\t"
-                  "5  - Backward Euler,\n\t"
-                  "6  - SDIRK 2,\n\t"
-                  "7  - SDIRK 3,\n\t"
-                  "8  - EPIC (exponential euler)\n\t");
-   args.AddOption(&t_final, "-tf", "--t-final",
-                  "Final time; start time is 0.");
-   args.AddOption(&dt, "-dt", "--time-step",
-                  "Time step.");
-   args.AddOption(&alpha, "-a", "--alpha",
-                  "Alpha coefficient.");
-   args.AddOption(&kappa, "-k", "--kappa",
-                  "Kappa coefficient offset.");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
-                  "--no-visit-datafiles",
-                  "Save data files for VisIt (visit.llnl.gov) visualization.");
-   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
-                  "Visualize every n-th timestep.");
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(cout);
-      return 1;
-   }
-   if (ode_solver_type < 1 || ode_solver_type > 9)
-   {
-      cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
-      return 3;
-   }
-   args.PrintOptions(cout);
-
-   // 2. Read the mesh from the given mesh file. We can handle triangular,
-   //    quadrilateral, tetrahedral and hexahedral meshes with the same code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   int dim = mesh->Dimension();
-
-   // 3. Refine the mesh to increase the resolution. In this example we do
-   //    'ref_levels' of uniform refinement, where 'ref_levels' is a
-   //    command-line parameter.
-   for (int lev = 0; lev < ref_levels; lev++)
-   {
-      mesh->UniformRefinement();
-   }
-
-   // 4. Define the vector finite element space representing the current and the
-   //    initial temperature, u_ref.
-   H1_FECollection fe_coll(order, dim);
-   FiniteElementSpace fespace(mesh, &fe_coll);
-
-   int fe_size = fespace.GetTrueVSize();
-   cout << "Number of temperature unknowns: " << fe_size << endl;
-
-   GridFunction u_gf(&fespace);
-
-   // 5. Set the initial conditions for u. All boundaries are considered
-   //    natural.
-   FunctionCoefficient u_0(InitialTemperature);
-   u_gf.ProjectCoefficient(u_0);
-   Vector u;
-   u_gf.GetTrueDofs(u);
-
-   // 6. Initialize the conduction operator and the visualization.
-   ConductionOperator oper(fespace, alpha, kappa, u);
-
-   u_gf.SetFromTrueDofs(u);
-   {
-      ofstream omesh("ex16.mesh");
-      omesh.precision(precision);
-      mesh->Print(omesh);
-      ofstream osol("ex16-init.gf");
-      osol.precision(precision);
-      u_gf.Save(osol);
-   }
-
-   VisItDataCollection visit_dc("Example16", mesh);
-   visit_dc.RegisterField("temperature", &u_gf);
-   if (visit)
-   {
-      visit_dc.SetCycle(0);
-      visit_dc.SetTime(0.0);
-      visit_dc.Save();
-   }
-
-   socketstream sout;
-   if (visualization)
-   {
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      sout.open(vishost, visport);
-      if (!sout)
-      {
-         cout << "Unable to connect to GLVis server at "
-              << vishost << ':' << visport << endl;
-         visualization = false;
-         cout << "GLVis visualization disabled.\n";
-      }
-      else
-      {
-         sout.precision(precision);
-         sout << "solution\n" << *mesh << u_gf;
-         sout << "pause\n";
-         sout << flush;
-         cout << "GLVis visualization paused."
-              << " Press space (in the GLVis window) to resume it.\n";
-      }
-   }
-
-   // 7. Define the ODE solver used for time integration.
-   double t = 0.0;
-   ODESolver *ode_solver = NULL;
-   switch (ode_solver_type)
-   {
-      // MFEM explicit methods
-      case 1: ode_solver = new ForwardEulerSolver; break;
-      case 2: ode_solver = new RK2Solver(0.5); break; // midpoint method
-      case 3: ode_solver = new RK3SSPSolver; break;
-      case 4: ode_solver = new RK4Solver; break;
-      // MFEM implicit L-stable methods
-      case 5: ode_solver = new BackwardEulerSolver; break;
-      case 6: ode_solver = new SDIRK23Solver(2); break;
-      case 7: ode_solver = new SDIRK33Solver; break;
-      // EPIC 
-      case 8: ode_solver = new EPI2();break;
-      case 9: ode_solver = new EPIRK4(); break;
-   }
-
-   // Initialize integrators
-   ode_solver->Init(oper);
-
-   // 8. Perform time-integration (looping over the time iterations, ti, with a
-   //    time-step dt).
-   cout << "Integrating the ODE ..." << endl;
-   tic_toc.Clear();
-   tic_toc.Start();
-
-   /*ofstream out_jac_an("jacobian_an.txt");
-   ofstream out_jac_fd("jacobian_fd.txt");
-   oper.ExtractJacobians(u, out_jac_fd, out_jac_an);*/
-
-   bool last_step = false;
-   int ti;
-   for (ti = 1; !last_step; ti++)
-   {
-      double dt_real = min(dt, t_final - t);
-
-      // Note that since we are using the "one-step" mode of the SUNDIALS
-      // solvers, they will, generally, step over the final time and will not
-      // explicitly perform the interpolation to t_final as they do in the
-      // "normal" step mode.
-      ode_solver->Step(u, t, dt_real);
-
-      oper.UpdateStats();
-
-      last_step = (t >= t_final - 1e-8*dt);
-
-      if (last_step || (ti % vis_steps) == 0) {
-         cout << "step " << ti << ", t = " << t << endl;
-
-         u_gf.SetFromTrueDofs(u);
-         if (visualization) {
-            sout << "solution\n" << *mesh << u_gf << flush;
-         }
-
-         if (visit) {
-            visit_dc.SetCycle(ti);
-            visit_dc.SetTime(t);
-            visit_dc.Save();
-         }
-      }
-   }
-   tic_toc.Stop();
-   double comp_time = tic_toc.RealTime();
-   cout << "Done, " << comp_time << "s." << endl;
-
-   // 9. Save the final solution. This output can be viewed later using GLVis:
-   //    "glvis -m ex16.mesh -g ex16-final.gf".
-   {
-      ofstream osol("ex16-final.gf");
-      osol.precision(precision);
-      u_gf.Save(osol);
-
-      ofstream ostats("ex16-stats.txt");
-      ostats << "time " << comp_time << endl;
-      oper.PrintStats(ostats);
-   }
-
-   // 10. Free the used memory.
-   delete ode_solver;
-   delete mesh;
-
-   return 0;
-}
-
-ConductionOperator::ConductionOperator(FiniteElementSpace &f, double al, double kap, const Vector &u)
-   : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL), dK(NULL), J_K(NULL), jac(NULL), z(height), u0(height),
-     nRhsMult(0), nSetJac(0), nJacMult(0), nImpSolve(0), nImpIter(0), nImpMult(0), nImpSet(0)
-{
-   const double rel_tol = 1e-8;
-
-   M = new BilinearForm(&fespace);
-   M->AddDomainIntegrator(new MassIntegrator());
-   M->Assemble();
-   M->FormSystemMatrix(ess_tdof_list, Mmat);
-
-   M_solver.iterative_mode = false;
-   M_solver.SetRelTol(rel_tol);
-   M_solver.SetAbsTol(0.0);
-   M_solver.SetMaxIter(50);
-   M_solver.SetPrintLevel(0);
-   M_solver.SetPreconditioner(M_prec);
-   M_solver.SetOperator(Mmat);
-
-   Jg_solver.SetRelTol(rel_tol);
-   Jg_solver.SetAbsTol(0.0);
-   Jg_solver.SetMaxIter(50);
-   Jg_solver.SetPrintLevel(0);
-   Jg_solver.SetPreconditioner(Jg_prec);
-
-   newton_solver.SetMaxIter(10);
-   newton_solver.SetRelTol(rel_tol);
-   newton_solver.SetPrintLevel(-1);
-   newton_solver.SetSolver(Jg_solver);
-   newton_solver.SetMaxIter(100);
-   newton_solver.iterative_mode = false;
-
-   alpha = al;
-   kappa = kap;
-}
-
-void ConductionOperator::UpdateStats()
-{
-   if (jac)
-   {
-      nJacMult += jac->GetnMult();
-   }
-}
-
-void ConductionOperator::PrintStats(ostream &out)
-{
-   out  << "nRhsMult " << nRhsMult << endl
-        << "nSetJac " << nSetJac << endl
-        << "nJacMult " << nJacMult  << endl
-        << "nImplicitSolve " << nImpSolve  << endl
-        << "nImplicitIter " << nImpIter << endl
-        << "nImplicitMult " << nImpMult << endl
-        << "nImplicitSet " << nImpSet << endl;
-}
-
-BilinearForm& ConductionOperator::GetKLambda(const Vector &u) const
-{
-   GridFunction conductivity_gf(&fespace);
-   conductivity_gf.SetFromTrueDofs(u);
-   for (int i = 0; i < conductivity_gf.Size(); i++)
-   {
-      conductivity_gf(i) = kappa + alpha*conductivity_gf(i);
-   }
-
-   GridFunctionCoefficient conductivity_coeff(&conductivity_gf);
-
-   delete K;
-   K = new BilinearForm(&fespace);
-   K->AddDomainIntegrator(new DiffusionIntegrator(conductivity_coeff));
-   K->Assemble();
-
-   return *K;
-}
-
-BilinearForm& ConductionOperator::GetdKLambda(const Vector &u) const
-{
-   GridFunction conductivity_gf(&fespace);
-   conductivity_gf.SetFromTrueDofs(u);
-   for (int i = 0; i < conductivity_gf.Size(); i++)
-   {
-      conductivity_gf(i) = kappa + alpha*conductivity_gf(i);
-   }
-
-   // Define diffusion form with conductivity = kappa(u0)
-   GridFunctionCoefficient conductivity_coeff(&conductivity_gf);
-
-   // Define advection form with velocity = grad kappa(u0)
-   GridFunction neg_cond_gf(conductivity_gf);
-   neg_cond_gf.Neg();
-   GradientGridFunctionCoefficient velocity_coeff(&neg_cond_gf);
-
-   delete dK;
-   dK = new BilinearForm(&fespace);
-
-   dK->AddDomainIntegrator(new DiffusionIntegrator(conductivity_coeff));
-   dK->AddDomainIntegrator(new MixedScalarWeakDivergenceIntegrator(velocity_coeff));
-   dK->Assemble();
-
-   return *dK;
-}
-
-void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
-{
-   // Compute:
-   //    du_dt = M^{-1}*-K(u)
-   // for du_dt
-   GetKLambda(u);
-   K->Mult(u, z);
-   z.Neg(); // z = -z
-   M_solver.Mult(z, du_dt);
-   nRhsMult++;
-}
-
-void ConductionOperator::ImplicitSolve(const double dt, const Vector &x, Vector &k)
-{
-   ImplicitSolveOperator imp_oper(this, &this->Mmat, dt, &x);
-   newton_solver.SetOperator(imp_oper);
-
-   Vector zero; // empty vector is interpreted as zero r.h.s. by NewtonSolver
-   newton_solver.Mult(zero, k);
-   MFEM_VERIFY(newton_solver.GetConverged(), "Newton solver did not converge.");
-
-   nImpSolve++;
-   nImpMult += imp_oper.GetnMult();
-   nImpSet  += imp_oper.GetnSet();
-   nImpIter += newton_solver.GetNumIterations();
-}
-
-Operator &ConductionOperator::GetGradient(const Vector &u) const
-{
-   delete jac;
-   GetdKLambda(u);
-   jac = new JacobianOperator(dK, &M_solver);
-
-   nSetJac++;
-
-   return *jac;
-}
-
-ConductionOperator::~ConductionOperator()
-{
-   delete M;
-   delete K;
-   delete dK;
-   delete J_K;
-   delete jac;
-}
-
-ImplicitSolveOperator::ImplicitSolveOperator(ConductionOperator *oper_, const SparseMatrix* M_, double dt_, const Vector* x_):
-   Operator(oper_->Height()), oper(oper_), M(M_), dt(dt_), x(x_), u(height), z(height), Jg(NULL), nMult(0), nSet(0)
-{ }
-
-
-void ImplicitSolveOperator::Mult(const Vector& y, Vector& gy) const
-{
-   // Compute gy = g(y) = My + dt K(lambda(u)) u
-   // with u = x + dt y
-   add(*x, dt, y, u);
-   BilinearForm& K = oper->GetKLambda(u);
-   K.Mult(u, gy);
-
-   M->AddMult(y, gy);
-
-   nMult++;
-}
-
-Operator& ImplicitSolveOperator::GetGradient(const Vector &k) const
-{
-   add(*x, dt, k, u);
-
-   BilinearForm& dK = oper->GetdKLambda(u);
-   Array<int> ess_tdof_list;
-   SparseMatrix dK_mat;
-   dK.FormSystemMatrix(ess_tdof_list, dK_mat);
-
-   delete Jg;
-   Jg = Add(1.0, *M, dt, dK_mat);
-
-   nSet++;
-   return *Jg;
-}
-
-JacobianOperator::JacobianOperator(Operator* J_, Operator* M_solver_):
-   Operator(M_solver_->Height()), J(J_), M_solver(M_solver_), z(height), nMult(0)
-{ }
-
-void JacobianOperator::Mult(const Vector &v, Vector &Jv) const
-{
-   Vector temp(v);
-   J->Mult(v, z);
-   z.Neg(); // z = -z
-   M_solver->Mult(z, Jv);
-   nMult++;
-}
-
-
-void ConductionOperator::ExtractJacobians(const Vector& x, std::ostream &out, std::ostream &out2)
-{
-   int n = x.Size();
-
-   Vector e(n);
-   e = 0.0;
-
-   double eps = 1e-8;
-   Vector fx(n), fx_eps(n), x_eps(n);
-   Mult(x, fx);
-
-   DenseMatrix J(n);
-
-   for (int i = 0; i < n; i++)
-   {
-      e[i] = 1.0;
-      add(x, eps, e, x_eps);
-      Mult(x_eps, fx_eps);
-      fx_eps -= fx;
-      fx_eps /= eps;
-      J.SetCol(i, fx_eps);
-      e[i] = 0.0;
-   }
-
-   J.PrintMatlab(out);
-   GetGradient(x);
-   jac->ExtractJacobian(x, out2);
-}
-
-void JacobianOperator::ExtractJacobian(const Vector& x, std::ostream &out)
-{
-   int n = z.Size();
-
-   Vector e(n);
-   e= 0.0;
-
-   Vector J_i(n);
-   DenseMatrix J(n);
-
-   for (int i = 0; i < n; i++)
-   {
-      e[i] = 1.0;
-      Mult(e, J_i);
-      J.SetCol(i, J_i);
-      e[i] = 0.0;
-   }
-
-   J.PrintMatlab(out);
-}
-
-double InitialTemperature(const Vector &x)
-{
-   if (x.Norml2() < 0.5) { return 2.0; }
-   else                  { return 1.0; }
-}

examples/epic/ex16p.cpp

@@ -1,494 +0,0 @@
-//                         MFEM Example 16 - Parallel Version
-//                             SUNDIALS Modification
-//
-// Compile with: make ex16p
-//
-// Sample runs:
-//     mpirun -np 4 ex16p
-//     mpirun -np 4 ex16p -m ../../data/inline-tri.mesh
-//     mpirun -np 4 ex16p -m ../../data/disc-nurbs.mesh -tf 2
-//     mpirun -np 4 ex16p -s 12 -a 0.0 -k 1.0
-//     mpirun -np 4 ex16p -s 8 -a 1.0 -k 0.0 -dt 4e-6 -tf 2e-2 -vs 50
-//     mpirun -np 8 ex16p -s 9 -a 0.5 -k 0.5 -o 4 -dt 8e-6 -tf 2e-2 -vs 50
-//     mpirun -np 4 ex16p -s 10 -dt 2.0e-4 -tf 4.0e-2
-//     mpirun -np 16 ex16p -m ../../data/fichera-q2.mesh
-//     mpirun -np 16 ex16p -m ../../data/escher-p2.mesh
-//     mpirun -np 8 ex16p -m ../../data/beam-tet.mesh -tf 10 -dt 0.1
-//     mpirun -np 4 ex16p -m ../../data/amr-quad.mesh -o 4 -rs 0 -rp 0
-//     mpirun -np 4 ex16p -m ../../data/amr-hex.mesh -o 2 -rs 0 -rp 0
-//
-// Description:  This example solves a time dependent nonlinear heat equation
-//               problem of the form du/dt = C(u), with a non-linear diffusion
-//               operator C(u) = \nabla \cdot (\kappa + \alpha u) \nabla u.
-//
-//               We recommend viewing examples 2, 9 and 10 before viewing this
-//               example.
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-/** After spatial discretization, the conduction model can be written as:
- *
- *     du/dt = M^{-1}(-Ku)
- *
- *  where u is the vector representing the temperature, M is the mass matrix,
- *  and K is the diffusion operator with diffusivity depending on u:
- *  (\kappa + \alpha u).
- *
- *  Class ConductionOperator represents the right-hand side of the above ODE.
- */
-class ConductionOperator : public TimeDependentOperator
-{
-protected:
-   ParFiniteElementSpace &fespace;
-   Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
-
-   ParBilinearForm *M;
-   ParBilinearForm *K;
-
-   HypreParMatrix Mmat;
-   HypreParMatrix Kmat;
-   HypreParMatrix *T; // T = M + dt K
-   double current_dt;
-
-   CGSolver M_solver;    // Krylov solver for inverting the mass matrix M
-   HypreSmoother M_prec; // Preconditioner for the mass matrix M
-
-   CGSolver T_solver;    // Implicit solver for T = M + dt K
-   HypreSmoother T_prec; // Preconditioner for the implicit solver
-
-   double alpha, kappa;
-
-   mutable Vector z; // auxiliary vector
-
-public:
-   ConductionOperator(ParFiniteElementSpace &f, double alpha, double kappa,
-                      const Vector &u);
-
-   virtual void Mult(const Vector &u, Vector &du_dt) const;
-
-   /** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k.
-       This is the only requirement for high-order SDIRK implicit integration.*/
-   virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k);
-
-   /** Setup the system (M + dt K) x = M b. This method is used by the implicit
-       SUNDIALS solvers. */
-   virtual int SUNImplicitSetup(const Vector &x, const Vector &fx,
-                                int jok, int *jcur, double gamma);
-
-   /** Solve the system (M + dt K) x = M b. This method is used by the implicit
-       SUNDIALS solvers. */
-   virtual int SUNImplicitSolve(const Vector &b, Vector &x, double tol);
-
-   /// Update the diffusion BilinearForm K using the given true-dof vector `u`.
-   void SetParameters(const Vector &u);
-
-   virtual ~ConductionOperator();
-};
-
-double InitialTemperature(const Vector &x);
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI.
-   int num_procs, myid;
-   MPI_Init(&argc, &argv);
-   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
-   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
-
-   // 2. Parse command-line options.
-   const char *mesh_file = "../../data/star.mesh";
-   int ser_ref_levels = 2;
-   int par_ref_levels = 1;
-   int order = 2;
-   int ode_solver_type = 8; // Exponential Euler 
-   double t_final = 0.5;
-   double dt = 1.0e-2;
-   double alpha = 1.0e-2;
-   double kappa = 0.5;
-   bool visualization = true;
-   bool visit = false;
-   int vis_steps = 5;
-
-   int precision = 8;
-   cout.precision(precision);
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
-                  "Number of times to refine the mesh uniformly in serial.");
-   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
-                  "Number of times to refine the mesh uniformly in parallel.");
-   args.AddOption(&order, "-o", "--order",
-                  "Order (degree) of the finite elements.");
-   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
-                  "ODE solver:\n\t"
-                  "1  - Forward Euler,\n\t"
-                  "2  - RK2,\n\t"
-                  "3  - RK3 SSP,\n\t"
-                  "4  - RK4,\n\t"
-                  "5  - Backward Euler,\n\t"
-                  "6  - SDIRK 2,\n\t"
-                  "7  - SDIRK 3,\n\t"
-                  "8  - Exponential Euler,\n\t");
-   args.AddOption(&t_final, "-tf", "--t-final",
-                  "Final time; start time is 0.");
-   args.AddOption(&dt, "-dt", "--time-step",
-                  "Time step.");
-   args.AddOption(&alpha, "-a", "--alpha",
-                  "Alpha coefficient.");
-   args.AddOption(&kappa, "-k", "--kappa",
-                  "Kappa coefficient offset.");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
-                  "--no-visit-datafiles",
-                  "Save data files for VisIt (visit.llnl.gov) visualization.");
-   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
-                  "Visualize every n-th timestep.");
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(cout);
-      MPI_Finalize();
-      return 1;
-   }
-
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
-
-   // check for vaild ODE solver option
-   if (ode_solver_type < 1 || ode_solver_type > 8)
-   {
-      if (myid == 0)
-      {
-         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
-      }
-      MPI_Finalize();
-      return 1;
-   }
-
-   // 3. Read the serial mesh from the given mesh file on all processors. We can
-   //    handle triangular, quadrilateral, tetrahedral and hexahedral meshes
-   //    with the same code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   int dim = mesh->Dimension();
-
-   // 4. Refine the mesh in serial to increase the resolution. In this example
-   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
-   //    a command-line parameter.
-   for (int lev = 0; lev < ser_ref_levels; lev++)
-   {
-      mesh->UniformRefinement();
-   }
-
-   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh further in parallel to increase the resolution. Once the
-   //    parallel mesh is defined, the serial mesh can be deleted.
-   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
-   delete mesh;
-   for (int lev = 0; lev < par_ref_levels; lev++)
-   {
-      pmesh->UniformRefinement();
-   }
-
-   // 6. Define the vector finite element space representing the current and the
-   //    initial temperature, u_ref.
-   H1_FECollection fe_coll(order, dim);
-   ParFiniteElementSpace fespace(pmesh, &fe_coll);
-
-   int fe_size = fespace.GlobalTrueVSize();
-   if (myid == 0)
-   {
-      cout << "Number of temperature unknowns: " << fe_size << endl;
-   }
-
-   ParGridFunction u_gf(&fespace);
-
-   // 7. Set the initial conditions for u. All boundaries are considered
-   //    natural.
-   FunctionCoefficient u_0(InitialTemperature);
-   u_gf.ProjectCoefficient(u_0);
-   Vector u;
-   u_gf.GetTrueDofs(u);
-
-   // 8. Initialize the conduction operator and the VisIt visualization.
-   ConductionOperator oper(fespace, alpha, kappa, u);
-
-   u_gf.SetFromTrueDofs(u);
-   {
-      ostringstream mesh_name, sol_name;
-      mesh_name << "ex16-mesh." << setfill('0') << setw(6) << myid;
-      sol_name << "ex16-init." << setfill('0') << setw(6) << myid;
-      ofstream omesh(mesh_name.str().c_str());
-      omesh.precision(precision);
-      pmesh->Print(omesh);
-      ofstream osol(sol_name.str().c_str());
-      osol.precision(precision);
-      u_gf.Save(osol);
-   }
-
-   VisItDataCollection visit_dc("Example16-Parallel", pmesh);
-   visit_dc.RegisterField("temperature", &u_gf);
-   if (visit)
-   {
-      visit_dc.SetCycle(0);
-      visit_dc.SetTime(0.0);
-      visit_dc.Save();
-   }
-
-   socketstream sout;
-   if (visualization)
-   {
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      sout.open(vishost, visport);
-      sout << "parallel " << num_procs << " " << myid << endl;
-      int good = sout.good(), all_good;
-      MPI_Allreduce(&good, &all_good, 1, MPI_INT, MPI_MIN, pmesh->GetComm());
-      if (!all_good)
-      {
-         sout.close();
-         visualization = false;
-         if (myid == 0)
-         {
-            cout << "Unable to connect to GLVis server at "
-                 << vishost << ':' << visport << endl;
-            cout << "GLVis visualization disabled.\n";
-         }
-      }
-      else
-      {
-         sout.precision(precision);
-         sout << "solution\n" << *pmesh << u_gf;
-         sout << "pause\n";
-         sout << flush;
-         if (myid == 0)
-         {
-            cout << "GLVis visualization paused."
-                 << " Press space (in the GLVis window) to resume it.\n";
-         }
-      }
-   }
-
-   // 9. Define the ODE solver used for time integration.
-   double t = 0.0;
-   ODESolver *ode_solver = NULL;
-   EPICSolver *epic_solver = NULL;
-   switch (ode_solver_type)
-   {
-      // MFEM explicit methods
-      case 1: ode_solver = new ForwardEulerSolver; break;
-      case 2: ode_solver = new RK2Solver(0.5); break; // midpoint method
-      case 3: ode_solver = new RK3SSPSolver; break;
-      case 4: ode_solver = new RK4Solver; break;
-      // MFEM implicit L-stable methods
-      case 5: ode_solver = new BackwardEulerSolver; break;
-      case 6: ode_solver = new SDIRK23Solver(2); break;
-      case 7: ode_solver = new SDIRK33Solver; break;
-      // EPIC
-      case 8:
-         epic_solver = new EPICSolver();
-	 epic_solver->Init(oper);
-	 ode_solver = epic_solver;
-	 break;
-   }
-
-   // Initialize MFEM integrators
-   ode_solver->Init(oper);
-
-   // 10. Perform time-integration (looping over the time iterations, ti, with a
-   //     time-step dt).
-   if (myid == 0)
-   {
-      cout << "Integrating the ODE ..." << endl;
-   }
-   tic_toc.Clear();
-   tic_toc.Start();
-
-   bool last_step = false;
-   for (int ti = 1; !last_step; ti++)
-   {
-      double dt_real = min(dt, t_final - t);
-
-      // Note that since we are using the "one-step" mode of the SUNDIALS
-      // solvers, they will, generally, step over the final time and will not
-      // explicitly perform the interpolation to t_final as they do in the
-      // "normal" step mode.
-
-      ode_solver->Step(u, t, dt_real);
-
-      last_step = (t >= t_final - 1e-8*dt);
-
-      if (last_step || (ti % vis_steps) == 0)
-      {
-         if (myid == 0)
-         {
-            cout << "step " << ti << ", t = " << t << endl;
-         }
-
-         u_gf.SetFromTrueDofs(u);
-         if (visualization)
-         {
-            sout << "parallel " << num_procs << " " << myid << "\n";
-            sout << "solution\n" << *pmesh << u_gf << flush;
-         }
-
-         if (visit)
-         {
-            visit_dc.SetCycle(ti);
-            visit_dc.SetTime(t);
-            visit_dc.Save();
-         }
-      }
-      oper.SetParameters(u);
-   }
-   tic_toc.Stop();
-   if (myid == 0)
-   {
-      cout << "Done, " << tic_toc.RealTime() << "s." << endl;
-   }
-
-   // 11. Save the final solution in parallel. This output can be viewed later
-   //     using GLVis: "glvis -np <np> -m ex16-mesh -g ex16-final".
-   {
-      ostringstream sol_name;
-      sol_name << "ex16-final." << setfill('0') << setw(6) << myid;
-      ofstream osol(sol_name.str().c_str());
-      osol.precision(precision);
-      u_gf.Save(osol);
-   }
-
-   // 12. Free the used memory.
-   delete ode_solver;
-   delete pmesh;
-
-   MPI_Finalize();
-
-   return 0;
-}
-
-ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
-                                       double kap, const Vector &u)
-   : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
-     T(NULL),
-     M_solver(f.GetComm()), T_solver(f.GetComm()), z(height)
-{
-   const double rel_tol = 1e-8;
-
-   M = new ParBilinearForm(&fespace);
-   M->AddDomainIntegrator(new MassIntegrator());
-   M->Assemble(0); // keep sparsity pattern of M and K the same
-   M->FormSystemMatrix(ess_tdof_list, Mmat);
-
-   M_solver.iterative_mode = false;
-   M_solver.SetRelTol(rel_tol);
-   M_solver.SetAbsTol(0.0);
-   M_solver.SetMaxIter(100);
-   M_solver.SetPrintLevel(0);
-   M_prec.SetType(HypreSmoother::Jacobi);
-   M_solver.SetPreconditioner(M_prec);
-   M_solver.SetOperator(Mmat);
-
-   alpha = al;
-   kappa = kap;
-
-   T_solver.iterative_mode = false;
-   T_solver.SetRelTol(rel_tol);
-   T_solver.SetAbsTol(0.0);
-   T_solver.SetMaxIter(100);
-   T_solver.SetPrintLevel(0);
-   T_solver.SetPreconditioner(T_prec);
-
-   SetParameters(u);
-}
-
-void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
-{
-   // Compute:
-   //    du_dt = M^{-1}*-K(u)
-   // for du_dt
-   Kmat.Mult(u, z);
-   z.Neg(); // z = -z
-   M_solver.Mult(z, du_dt);
-}
-
-void ConductionOperator::ImplicitSolve(const double dt,
-                                       const Vector &u, Vector &du_dt)
-{
-   // Solve the equation:
-   //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
-   // for du_dt
-   if (T) { delete T; }
-   T = Add(1.0, Mmat, dt, Kmat);
-   T_solver.SetOperator(*T);
-   Kmat.Mult(u, z);
-   z.Neg();
-   T_solver.Mult(z, du_dt);
-}
-
-int ConductionOperator::SUNImplicitSetup(const Vector &x,
-                                         const Vector &fx, int jok, int *jcur,
-                                         double gamma)
-{
-   // Setup the ODE Jacobian T = M + gamma K.
-   if (T) { delete T; }
-   T = Add(1.0, Mmat, gamma, Kmat);
-   T_solver.SetOperator(*T);
-   *jcur = 1;
-   return (0);
-}
-
-int ConductionOperator::SUNImplicitSolve(const Vector &b, Vector &x, double tol)
-{
-   // Solve the system A x = z => (M - gamma K) x = M b.
-   Mmat.Mult(b, z);
-   T_solver.Mult(z, x);
-   return (0);
-}
-
-void ConductionOperator::SetParameters(const Vector &u)
-{
-   ParGridFunction u_alpha_gf(&fespace);
-   u_alpha_gf.SetFromTrueDofs(u);
-   for (int i = 0; i < u_alpha_gf.Size(); i++)
-   {
-      u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i);
-   }
-
-   delete K;
-   K = new ParBilinearForm(&fespace);
-
-   GridFunctionCoefficient u_coeff(&u_alpha_gf);
-
-   K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff));
-   K->Assemble(0); // keep sparsity pattern of M and K the same
-   K->FormSystemMatrix(ess_tdof_list, Kmat);
-}
-
-ConductionOperator::~ConductionOperator()
-{
-   delete T;
-   delete M;
-   delete K;
-}
-
-double InitialTemperature(const Vector &x)
-{
-   if (x.Norml2() < 0.5)
-   {
-      return 2.0;
-   }
-   else
-   {
-      return 1.0;
-   }
-}

examples/epic/makefile

@@ -1,76 +0,0 @@
-# Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
-# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-# LICENSE and NOTICE for details. LLNL-CODE-806117.
-#
-# This file is part of the MFEM library. For more information and source code
-# availability visit https://mfem.org.
-#
-# MFEM is free software; you can redistribute it and/or modify it under the
-# terms of the BSD-3 license. We welcome feedback and contributions, see file
-# CONTRIBUTING.md for details.
-
-# Use the MFEM build directory
-MFEM_DIR ?= ../..
-MFEM_BUILD_DIR ?= ../..
-SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/examples/epic/,)
-CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
-# Use the MFEM install directory
-# MFEM_INSTALL_DIR = ../../mfem
-# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
-
-MFEM_LIB_FILE = mfem_is_not_built
--include $(CONFIG_MK)
-
-SEQ_EXAMPLES = ex16
-PAR_EXAMPLES = ex16p
-ifeq ($(MFEM_USE_MPI),NO)
-   EXAMPLES = $(SEQ_EXAMPLES)
-else
-   EXAMPLES = $(PAR_EXAMPLES) $(SEQ_EXAMPLES)
-endif
-
-.SUFFIXES:
-.SUFFIXES: .o .cpp .mk
-.PHONY: all clean clean-build clean-exec
-
-# Remove built-in rule
-%: %.cpp
-
-# Replace the default implicit rule for *.cpp files
-%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
-	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
-
-all: $(EXAMPLES)
-
-ifeq ($(MFEM_USE_EPIC),NO)
-$(EXAMPLES):
-	$(error MFEM is not configured with EPIC)
-endif
-
-MFEM_TESTS = EXAMPLES
-include $(MFEM_TEST_MK)
-
-# Testing: Parallel vs. serial runs
-RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
-SERIAL_NAME := Serial EPIC example
-PARALLEL_NAME := Parallel EPIC example
-%-test-par: %
-	@$(call mfem-test,$<, $(RUN_MPI), $(PARALLEL_NAME))
-%-test-seq: %
-	@$(call mfem-test,$<,, $(SERIAL_NAME))
-
-# Testing: "test" target and mfem-test* variables are defined in config/test.mk
-
-# Generate an error message if the MFEM library is not built and exit
-$(MFEM_LIB_FILE):
-	$(error The MFEM library is not built)
-
-clean: clean-build clean-exec
-
-clean-build:
-	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
-	rm -rf *.dSYM *.TVD.*breakpoints
-
-clean-exec:
-	@rm -f deformed.* velocity.* elastic_energy.*
-	@rm -f ex16.mesh ex16-mesh.* ex16-init.* ex16-final.* Example16*

examples/makefile

@@ -45,9 +45,6 @@ endif
 ifeq ($(MFEM_USE_HIOP),YES)
    SUBDIRS += hiop
 endif
-ifeq ($(MFEM_USE_EPIC),YES)
-   SUBDIRS += epic
-endif
 ifeq ($(MFEM_USE_PETSC),YES)
    SUBDIRS += petsc
 endif

fem/CMakeLists.txt

@@ -105,6 +105,7 @@ set(SRCS
   tmop/tmop_pa_w3.cpp
   tmop/tmop_pa_w3_c0.cpp
   tmop_tools.cpp
+  tmop_amr.cpp
   gslib.cpp
   transfer.cpp
   lor.cpp
@@ -164,6 +165,7 @@ set(HDRS
   tmop.hpp
   tmop/tmop_pa.hpp
   tmop_tools.hpp
+  tmop_amr.hpp
   gslib.hpp
   transfer.hpp
   lor.hpp

fem/fem.hpp

@@ -34,6 +34,7 @@
 #include "staticcond.hpp"
 #include "tmop.hpp"
 #include "tmop_tools.hpp"
+#include "tmop_amr.hpp"
 #include "gslib.hpp"
 #include "restriction.hpp"
 #include "quadinterpolator.hpp"

fem/gridfunc.hpp

@@ -95,6 +95,12 @@ public:
       : Vector(data, f->GetVSize())
    { fes = f; fec = NULL; fes_sequence = f->GetSequence(); UseDevice(true); }
 
+   /** @brief Construct a GridFunction using previously allocated Vector @a base
+       starting at the given offset, @a base_offset. */
+   GridFunction(FiniteElementSpace *f, Vector &base, int base_offset = 0)
+      : Vector(base, base_offset, f->GetVSize())
+   { fes = f; fec = NULL; fes_sequence = f->GetSequence(); UseDevice(true); }
+
    /// Construct a GridFunction on the given Mesh, using the data from @a input.
    /** The content of @a input should be in the format created by the method
        Save(). The reconstructed FiniteElementSpace and FiniteElementCollection

fem/gslib.cpp

@@ -610,7 +610,12 @@ void FindPointsGSLIB::Interpolate(const GridFunction &field_in,
       {
          if (gsl_code[i] == 1) { indl2.Append(i); }
       }
-      if (indl2.Size() == 0) { return; } // no points on element borders
+      int borderPts = indl2.Size();
+#ifdef MFEM_USE_MPI
+      MPI_Allreduce(MPI_IN_PLACE, &borderPts, 1, MPI_INT, MPI_SUM, gsl_comm->c);
+#endif
+      if (borderPts == 0) { return; } // no points on element borders
+
 
       Vector field_out_l2(field_out.Size());
       VectorGridFunctionCoefficient field_in_dg(&field_in);

fem/pbilinearform.cpp

@@ -317,10 +317,6 @@ ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
 void ParBilinearForm::TrueAddMult(const Vector &x, Vector &y, const double a)
 const
 {
-   MFEM_VERIFY(interior_face_integs.Size() == 0,
-               "the case of interior face integrators is not"
-               " implemented");
-
    if (X.ParFESpace() != pfes)
    {
       X.SetSpace(pfes);
@@ -334,6 +330,9 @@ const
    }
    else
    {
+      MFEM_VERIFY(interior_face_integs.Size() == 0,
+                  "the case of interior face integrators is not"
+                  " implemented");
       mat->Mult(X, Y);
    }
    pfes->Dof_TrueDof_Matrix()->MultTranspose(a, Y, 1.0, y);

fem/pgridfunc.hpp

@@ -65,6 +65,11 @@ public:
    ParGridFunction(ParFiniteElementSpace *pf, double *data) :
       GridFunction(pf, data), pfes(pf) { }
 
+   /** @brief Construct a ParGridFunction using previously allocated Vector
+       @a base starting at the given offset, @a base_offset. */
+   ParGridFunction(ParFiniteElementSpace *pf, Vector &base, int base_offset = 0)
+      : GridFunction(pf, base, base_offset), pfes(pf) { }
+
    /// Construct a ParGridFunction using a GridFunction as external data.
    /** The parallel space @a *pf and the space used by @a *gf should match. The
        data from @a *gf is used as the local data of the ParGridFunction on each

fem/prestriction.cpp

@@ -33,6 +33,10 @@ ParL2FaceRestriction::ParL2FaceRestriction(const ParFiniteElementSpace &fes,
    // If fespace == L2
    const ParFiniteElementSpace &pfes =
       static_cast<const ParFiniteElementSpace&>(this->fes);
+
+   // Ensure the face neighbor data is constructed
+   pfes.GetParMesh()->ExchangeFaceNbrData();
+
    const FiniteElement *fe = pfes.GetFE(0);
    const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement*>(fe);
    MFEM_VERIFY(tfe != NULL &&

fem/restriction.cpp

@@ -15,6 +15,12 @@
 #include "../general/forall.hpp"
 #include <climits>
 
+#ifdef MFEM_USE_MPI
+
+#include "pfespace.hpp"
+
+#endif
+
 namespace mfem
 {
 
@@ -675,6 +681,19 @@ H1FaceRestriction::H1FaceRestriction(const FiniteElementSpace &fes,
      gather_indices(nf*dof)
 {
    if (nf==0) { return; }
+
+#ifdef MFEM_USE_MPI
+
+   // If the underlying finite element space is parallel, ensure the face
+   // neighbor information is generated.
+   if (const ParFiniteElementSpace *pfes
+       = dynamic_cast<const ParFiniteElementSpace*>(&fes))
+   {
+      pfes->GetParMesh()->ExchangeFaceNbrData();
+   }
+
+#endif
+
    // If fespace == H1
    const FiniteElement *fe = fes.GetFE(0);
    const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement*>(fe);

fem/tmop.cpp

@@ -1314,33 +1314,61 @@ static inline void device_copy(double *d_dest, const double *d_src, int size)
 } // namespace internal
 
 #ifdef MFEM_USE_MPI
-void DiscreteAdaptTC::FinalizeParDiscreteTargetSpec(const ParGridFunction
-                                                    &tspec_)
+void DiscreteAdaptTC::FinalizeParDiscreteTargetSpec(const ParGridFunction &t)
 {
    MFEM_VERIFY(adapt_eval, "SetAdaptivityEvaluator() has not been called!")
    MFEM_VERIFY(ncomp > 0, "No target specifications have been set!");
 
-   ParFiniteElementSpace *ptspec_fes = tspec_.ParFESpace();
+   ParFiniteElementSpace *ptspec_fes = t.ParFESpace();
 
    adapt_eval->SetParMetaInfo(*ptspec_fes->GetParMesh(),
                               *ptspec_fes->FEColl(), ncomp);
-   adapt_eval->SetInitialField(*tspec_fes->GetMesh()->GetNodes(), tspec);
+   adapt_eval->SetInitialField(*ptspec_fes->GetMesh()->GetNodes(), tspec);
 
    tspec_sav = tspec;
 
    delete tspec_fesv;
-   tspec_fesv = new FiniteElementSpace(tspec_fes->GetMesh(),
-                                       tspec_fes->FEColl(), ncomp);
+   tspec_fesv = new FiniteElementSpace(ptspec_fes->GetMesh(),
+                                       ptspec_fes->FEColl(), ncomp);
+
+   delete ptspec_fesv;
+   ptspec_fesv = new ParFiniteElementSpace(ptspec_fes->GetParMesh(),
+                                           ptspec_fes->FEColl(), ncomp);
+
+   delete tspec_pgf;
+   tspec_pgf = new ParGridFunction(ptspec_fesv, tspec);
+   tspec_gf = tspec_pgf;
+}
+
+void DiscreteAdaptTC::ParUpdateAfterMeshTopologyChange()
+{
+   ptspec_fesv->Update();
+   if (tspec_fesv)
+   {
+      delete tspec_fesv;
+      tspec_fesv = new FiniteElementSpace(ptspec_fesv->GetMesh(),
+                                          ptspec_fesv->FEColl(), ncomp);
+   }
+   tspec_pgf->Update();
+   tspec_gf = tspec_pgf;
+   tspec.SetDataAndSize(tspec_pgf->GetData(), tspec_pgf->Size());
+   tspec_sav = tspec;
+
+   adapt_eval->SetParMetaInfo(*ptspec_fesv->GetParMesh(),
+                              *ptspec_fesv->FEColl(), ncomp);
+   adapt_eval->SetInitialField(*ptspec_fesv->GetMesh()->GetNodes(), tspec);
 }
 
 void DiscreteAdaptTC::SetTspecAtIndex(int idx, const ParGridFunction &tspec_)
 {
-   const int vdim     = tspec_.FESpace()->GetVDim(),
-             dof_cnt  = tspec_.Size()/vdim;
+   const int vdim = tspec_.FESpace()->GetVDim(),
+             ndof = tspec_.FESpace()->GetNDofs();
+   MFEM_VERIFY(ndof == tspec.Size()/ncomp, "Inconsistency in SetTspecAtIndex.");
+
    const auto tspec__d = tspec_.Read();
    auto tspec_d = tspec.ReadWrite();
-   const int offset = idx*dof_cnt;
-   internal::device_copy(tspec_d + offset, tspec__d, dof_cnt*vdim);
+   const int offset = idx*ndof;
+   internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
    FinalizeParDiscreteTargetSpec(tspec_);
 }
 
@@ -1360,78 +1388,71 @@ void DiscreteAdaptTC::SetParDiscreteTargetSkew(const ParGridFunction &tspec_)
    FinalizeParDiscreteTargetSpec(tspec_);
 }
 
-void DiscreteAdaptTC::SetParDiscreteTargetAspectRatio(const ParGridFunction
-                                                      &tspec_)
+void DiscreteAdaptTC::SetParDiscreteTargetAspectRatio(const ParGridFunction &ar)
 {
-   if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, tspec_); return; }
+   if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, ar); return; }
    aspectratioidx = ncomp;
-   SetDiscreteTargetBase(tspec_);
-   FinalizeParDiscreteTargetSpec(tspec_);
+   SetDiscreteTargetBase(ar);
+   FinalizeParDiscreteTargetSpec(ar);
 }
 
-void DiscreteAdaptTC::SetParDiscreteTargetOrientation(const ParGridFunction
-                                                      &tspec_)
+void DiscreteAdaptTC::SetParDiscreteTargetOrientation(const ParGridFunction &o)
 {
-   if (orientationidx > -1) { SetTspecAtIndex(orientationidx, tspec_); return; }
+   if (orientationidx > -1) { SetTspecAtIndex(orientationidx, o); return; }
    orientationidx = ncomp;
-   SetDiscreteTargetBase(tspec_);
-   FinalizeParDiscreteTargetSpec(tspec_);
+   SetDiscreteTargetBase(o);
+   FinalizeParDiscreteTargetSpec(o);
 }
 
 void DiscreteAdaptTC::SetParDiscreteTargetSpec(const ParGridFunction &tspec_)
 {
    SetParDiscreteTargetSize(tspec_);
-   FinalizeParDiscreteTargetSpec(tspec_);
 }
 #endif // MFEM_USE_MPI
 
 void DiscreteAdaptTC::SetDiscreteTargetBase(const GridFunction &tspec_)
 {
-   const int vdim     = tspec_.FESpace()->GetVDim(),
-             dof_cnt  = tspec_.Size()/vdim;
+   const int vdim = tspec_.FESpace()->GetVDim(),
+             ndof = tspec_.FESpace()->GetNDofs();
 
    ncomp += vdim;
 
-   delete tspec_fes;
-   tspec_fes = new FiniteElementSpace(tspec_.FESpace()->GetMesh(),
-                                      tspec_.FESpace()->FEColl(), 1);
-
    // need to append data to tspec
    // make a copy of tspec->tspec_temp, increase its size, and
    // copy data from tspec_temp -> tspec, then add new entries
    Vector tspec_temp = tspec;
    tspec.UseDevice(true);
    tspec_sav.UseDevice(true);
-   tspec.SetSize(ncomp*dof_cnt);
+   tspec.SetSize(ncomp*ndof);
 
    const auto tspec_temp_d = tspec_temp.Read();
    auto tspec_d = tspec.ReadWrite();
    internal::device_copy(tspec_d, tspec_temp_d, tspec_temp.Size());
 
    const auto tspec__d = tspec_.Read();
-   const int offset = (ncomp-vdim)*dof_cnt;
-   internal::device_copy(tspec_d + offset, tspec__d, dof_cnt*vdim);
+   const int offset = (ncomp-vdim)*ndof;
+   internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
 }
 
 void DiscreteAdaptTC::SetTspecAtIndex(int idx, const GridFunction &tspec_)
 {
-   const int vdim     = tspec_.FESpace()->GetVDim(),
-             dof_cnt  = tspec_.Size()/vdim;
+   const int vdim = tspec_.FESpace()->GetVDim(),
+             ndof = tspec_.FESpace()->GetNDofs();
+   MFEM_VERIFY(ndof == tspec.Size()/ncomp, "Inconsistency in SetTargetSpec.");
 
    const auto tspec__d = tspec_.Read();
    auto tspec_d = tspec.ReadWrite();
-   const int offset = idx*dof_cnt;
-   internal::device_copy(tspec_d + offset, tspec__d, dof_cnt*vdim);
-   FinalizeSerialDiscreteTargetSpec();
+   const int offset = idx*ndof;
+   internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
+   FinalizeSerialDiscreteTargetSpec(tspec_);
 }
 
 void DiscreteAdaptTC::SetSerialDiscreteTargetSize(const GridFunction &tspec_)
 {
-
    if (sizeidx > -1) { SetTspecAtIndex(sizeidx, tspec_); return; }
    sizeidx = ncomp;
    SetDiscreteTargetBase(tspec_);
-   FinalizeSerialDiscreteTargetSpec();
+   FinalizeSerialDiscreteTargetSpec(tspec_);
 }
 
 void DiscreteAdaptTC::SetSerialDiscreteTargetSkew(const GridFunction &tspec_)
@@ -1439,32 +1460,31 @@ void DiscreteAdaptTC::SetSerialDiscreteTargetSkew(const GridFunction &tspec_)
    if (skewidx > -1) { SetTspecAtIndex(skewidx, tspec_); return; }
    skewidx = ncomp;
    SetDiscreteTargetBase(tspec_);
-   FinalizeSerialDiscreteTargetSpec();
+   FinalizeSerialDiscreteTargetSpec(tspec_);
 }
 
-void DiscreteAdaptTC::SetSerialDiscreteTargetAspectRatio(
-   const GridFunction &tspec_)
+void DiscreteAdaptTC::SetSerialDiscreteTargetAspectRatio(const GridFunction &ar)
 {
-   if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, tspec_); return; }
+   if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, ar); return; }
    aspectratioidx = ncomp;
-   SetDiscreteTargetBase(tspec_);
-   FinalizeSerialDiscreteTargetSpec();
+   SetDiscreteTargetBase(ar);
+   FinalizeSerialDiscreteTargetSpec(ar);
 }
 
-void DiscreteAdaptTC::SetSerialDiscreteTargetOrientation(
-   const GridFunction &tspec_)
+void DiscreteAdaptTC::SetSerialDiscreteTargetOrientation(const GridFunction &o)
 {
-   if (orientationidx > -1) { SetTspecAtIndex(orientationidx, tspec_); return; }
+   if (orientationidx > -1) { SetTspecAtIndex(orientationidx, o); return; }
    orientationidx = ncomp;
-   SetDiscreteTargetBase(tspec_);
-   FinalizeSerialDiscreteTargetSpec();
+   SetDiscreteTargetBase(o);
+   FinalizeSerialDiscreteTargetSpec(o);
 }
 
-void DiscreteAdaptTC::FinalizeSerialDiscreteTargetSpec()
+void DiscreteAdaptTC::FinalizeSerialDiscreteTargetSpec(const GridFunction &t)
 {
    MFEM_VERIFY(adapt_eval, "SetAdaptivityEvaluator() has not been called!")
    MFEM_VERIFY(ncomp > 0, "No target specifications have been set!");
 
+   const FiniteElementSpace *tspec_fes = t.FESpace();
    adapt_eval->SetSerialMetaInfo(*tspec_fes->GetMesh(),
                                  *tspec_fes->FEColl(), ncomp);
    adapt_eval->SetInitialField(*tspec_fes->GetMesh()->GetNodes(), tspec);
@@ -1474,12 +1494,40 @@ void DiscreteAdaptTC::FinalizeSerialDiscreteTargetSpec()
    delete tspec_fesv;
    tspec_fesv = new FiniteElementSpace(tspec_fes->GetMesh(),
                                        tspec_fes->FEColl(), ncomp);
+
+   delete tspec_gf;
+   tspec_gf = new GridFunction(tspec_fesv, tspec);
+}
+
+void DiscreteAdaptTC::GetDiscreteTargetSpec(GridFunction &tspec_, int idx)
+{
+   if (idx < 0) { return; }
+   const int ndof = tspec_.FESpace()->GetNDofs(),
+             vdim = tspec_.FESpace()->GetVDim();
+   MFEM_VERIFY(ndof == tspec.Size()/ncomp,
+               "Inconsistency in GetSerialDiscreteTargetSpec.");
+
+   for (int i = 0; i < ndof*vdim; i++)
+   {
+      tspec_(i) = tspec(i + idx*ndof);
+   }
+}
+
+void DiscreteAdaptTC::UpdateAfterMeshTopologyChange()
+{
+   tspec_fesv->Update();
+   tspec_gf->Update();
+   tspec.SetDataAndSize(tspec_gf->GetData(), tspec_gf->Size());
+   tspec_sav = tspec;
+
+   adapt_eval->SetSerialMetaInfo(*tspec_fesv->GetMesh(),
+                                 *tspec_fesv->FEColl(), ncomp);
+   adapt_eval->SetInitialField(*tspec_fesv->GetMesh()->GetNodes(), tspec);
 }
 
 void DiscreteAdaptTC::SetSerialDiscreteTargetSpec(const GridFunction &tspec_)
 {
    SetSerialDiscreteTargetSize(tspec_);
-   FinalizeSerialDiscreteTargetSpec();
 }
 
 
@@ -1509,7 +1557,7 @@ void DiscreteAdaptTC::UpdateTargetSpecificationAtNode(const FiniteElement &el,
    MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
 
    Array<int> dofs;
-   tspec_fes->GetElementDofs(T.ElementNo, dofs);
+   tspec_fesv->GetElementDofs(T.ElementNo, dofs);
    const int cnt = tspec.Size()/ncomp; // dofs per scalar-field
 
    for (int i = 0; i < ncomp; i++)
@@ -1524,7 +1572,7 @@ void DiscreteAdaptTC::RestoreTargetSpecificationAtNode(ElementTransformation &T,
    MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
 
    Array<int> dofs;
-   tspec_fes->GetElementDofs(T.ElementNo, dofs);
+   tspec_fesv->GetElementDofs(T.ElementNo, dofs);
    const int cnt = tspec.Size()/ncomp;
    for (int i = 0; i < ncomp; i++)
    {
@@ -1532,6 +1580,40 @@ void DiscreteAdaptTC::RestoreTargetSpecificationAtNode(ElementTransformation &T,
    }
 }
 
+void DiscreteAdaptTC::SetTspecFromIntRule(int e_id,
+                                          const IntegrationRule &intrule)
+{
+   switch (target_type)
+   {
+      case IDEAL_SHAPE_GIVEN_SIZE:
+      case GIVEN_SHAPE_AND_SIZE:
+      {
+         const int ndofs = tspec_fesv->GetFE(e_id)->GetDof(),
+                   ntspec_dofs = ndofs*ncomp;
+
+         Vector tspec_vals(ntspec_dofs);
+
+         Array<int> dofs;
+         tspec_fesv->GetElementVDofs(e_id, dofs);
+         tspec.GetSubVector(dofs, tspec_vals);
+         DenseMatrix tr;
+         tspec_gf->GetVectorValues(e_id, intrule, tspec_refine, tr);
+         tspec_refine.Transpose();
+         break;
+      }
+      default:
+         MFEM_ABORT("Incompatible target type for discrete adaptation!");
+   }
+}
+
+void DiscreteAdaptTC::SetTspecDataForDerefinement(FiniteElementSpace *fes)
+{
+   coarse_tspec_fesv = fes;
+   const Operator *c_op = fes->GetUpdateOperator();
+   tspec_derefine.SetSize(c_op->Height());
+   c_op->Mult(tspec, tspec_derefine);
+}
+
 void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
                                             const IntegrationRule &ir,
                                             const Vector &elfun,
@@ -1542,6 +1624,8 @@ void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
              nqp = ir.GetNPoints();
    Jtrcomp.SetSize(dim, dim, 4*nqp);
 
+   FiniteElementSpace *src_fes = tspec_fesv;
+
    switch (target_type)
    {
       case IDEAL_SHAPE_GIVEN_SIZE:
@@ -1550,7 +1634,7 @@ void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
          const DenseMatrix &Wideal =
             Geometries.GetGeomToPerfGeomJac(fe.GetGeomType());
          const int dim = Wideal.Height(),
-                   ndofs = tspec_fes->GetFE(e_id)->GetDof(),
+                   ndofs = tspec_fesv->GetFE(e_id)->GetDof(),
                    ntspec_dofs = ndofs*ncomp;
 
          Vector shape(ndofs), tspec_vals(ntspec_dofs), par_vals,
@@ -1561,11 +1645,29 @@ void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
          tspec_fesv->GetElementVDofs(e_id, dofs);
          tspec.UseDevice(true);
          tspec.GetSubVector(dofs, tspec_vals);
+         if (tspec_refine.NumCols() > 0) // Refinement
+         {
+            MFEM_VERIFY(amr_el >= 0, " Target being constructed for an AMR element.");
+            for (int i = 0; i < ncomp; i++)
+            {
+               for (int j = 0; j < ndofs; j++)
+               {
+                  tspec_vals(j + i*ndofs) = tspec_refine(j + amr_el*ndofs, i);
+               }
+            }
+         }
+         else if (tspec_derefine.Size() > 0) // Derefinement
+         {
+            dofs.SetSize(0);
+            coarse_tspec_fesv->GetElementVDofs(e_id, dofs);
+            tspec_derefine.GetSubVector(dofs, tspec_vals);
+            src_fes = coarse_tspec_fesv;
+         }
 
          for (int q = 0; q < nqp; q++)
          {
             const IntegrationPoint &ip = ir.IntPoint(q);
-            tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+            src_fes->GetFE(e_id)->CalcShape(ip, shape);
             Jtr(q) = Wideal; // Initialize to identity
             for (int d = 0; d < 4; d++)
             {
@@ -1576,9 +1678,16 @@ void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
             if (sizeidx != -1) // Set size
             {
                par_vals.SetDataAndSize(tspec_vals.GetData()+sizeidx*ndofs, ndofs);
-               const double min_size = par_vals.Min();
-               MFEM_VERIFY(min_size > 0.0,
-                           "Non-positive size propagated in the target definition.");
+               double min_size = par_vals.Min();//0.001; //
+               if (lim_min_size > 0.)
+               {
+                  min_size = lim_min_size;
+               }
+               else
+               {
+                  MFEM_VERIFY(min_size > 0.0,
+                              "Non-positive size propagated in the target definition.");
+               }
                const double size = std::max(shape * par_vals, min_size);
                Jtr(q).Set(std::pow(size, 1.0/dim), Jtr(q));
                DenseMatrix Jtrcomp_q(Jtrcomp.GetData(0 + 4*q), dim, dim);
@@ -1593,6 +1702,9 @@ void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
                {
                   par_vals.SetDataAndSize(tspec_vals.GetData()+
                                           aspectratioidx*ndofs, ndofs);
+                  const double min_size = par_vals.Min();
+                  MFEM_VERIFY(min_size > 0.0,
+                              "Non-positive aspect-ratio propagated in the target definition.");
 
                   const double aspectratio = shape * par_vals;
                   D_rho = 0.;
@@ -1777,7 +1889,7 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
 
                grad_phys.Mult(par_vals, grad_ptr_c1);
                Vector grad_q(dim);
-               tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+               tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
                grad_e_c1.MultTranspose(shape, grad_q);
 
                const double min_size = par_vals.Min();
@@ -1810,7 +1922,7 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
 
                   grad_phys.Mult(par_vals, grad_ptr_c1);
                   Vector grad_q(dim);
-                  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+                  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
                   grad_e_c1.MultTranspose(shape, grad_q);
 
                   const double aspectratio = shape * par_vals;
@@ -1841,7 +1953,7 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
                   grad_phys.Mult(par_vals_c2, grad_ptr_c2);
                   grad_phys.Mult(par_vals_c3, grad_ptr_c3);
                   Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
-                  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+                  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
                   grad_e_c1.MultTranspose(shape, grad_q1);
                   grad_e_c2.MultTranspose(shape, grad_q2);
                   grad_e_c3.MultTranspose(shape, grad_q3);
@@ -1880,7 +1992,7 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
 
                   grad_phys.Mult(par_vals, grad_ptr_c1);
                   Vector grad_q(dim);
-                  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+                  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
                   grad_e_c1.MultTranspose(shape, grad_q);
 
                   const double skew = shape * par_vals;
@@ -1913,7 +2025,7 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
                   grad_phys.Mult(par_vals_c2, grad_ptr_c2);
                   grad_phys.Mult(par_vals_c3, grad_ptr_c3);
                   Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
-                  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+                  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
                   grad_e_c1.MultTranspose(shape, grad_q1);
                   grad_e_c2.MultTranspose(shape, grad_q2);
                   grad_e_c3.MultTranspose(shape, grad_q3);
@@ -1960,7 +2072,7 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
 
                   grad_phys.Mult(par_vals, grad_ptr_c1);
                   Vector grad_q(dim);
-                  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+                  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
                   grad_e_c1.MultTranspose(shape, grad_q);
 
                   const double theta = shape * par_vals;
@@ -1991,7 +2103,7 @@ void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
                   grad_phys.Mult(par_vals_c2, grad_ptr_c2);
                   grad_phys.Mult(par_vals_c3, grad_ptr_c3);
                   Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
-                  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
+                  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
                   grad_e_c1.MultTranspose(shape, grad_q1);
                   grad_e_c2.MultTranspose(shape, grad_q2);
                   grad_e_c3.MultTranspose(shape, grad_q3);
@@ -2071,7 +2183,7 @@ void DiscreteAdaptTC::UpdateGradientTargetSpecification(const Vector &x,
 {
    if (use_flag && good_tspec_grad) { return; }
 
-   const int dim = tspec_fes->GetFE(0)->GetDim(),
+   const int dim = tspec_fesv->GetFE(0)->GetDim(),
              cnt = x.Size()/dim;
 
    tspec_pert1h.SetSize(x.Size()*ncomp);
@@ -2097,7 +2209,7 @@ void DiscreteAdaptTC::UpdateHessianTargetSpecification(const Vector &x,
 
    if (use_flag && good_tspec_hess) { return; }
 
-   const int dim    = tspec_fes->GetFE(0)->GetDim(),
+   const int dim    = tspec_fesv->GetFE(0)->GetDim(),
              cnt    = x.Size()/dim,
              totmix = 1+2*(dim-2);
 
@@ -2145,6 +2257,16 @@ void DiscreteAdaptTC::UpdateHessianTargetSpecification(const Vector &x,
    good_tspec_hess = use_flag;
 }
 
+DiscreteAdaptTC::~DiscreteAdaptTC()
+{
+   delete tspec_gf;
+   delete adapt_eval;
+   delete tspec_fesv;
+#ifdef MFEM_USE_MPI
+   delete ptspec_fesv;
+#endif
+}
+
 void AdaptivityEvaluator::SetSerialMetaInfo(const Mesh &m,
                                             const FiniteElementCollection &fec,
                                             int num_comp)
@@ -2258,6 +2380,7 @@ void TMOP_Integrator::EnableAdaptiveLimiting(const ParGridFunction &z0,
                                              AdaptivityEvaluator &ae)
 {
    zeta_0 = &z0;
+   pzeta_0 = &z0;
    delete zeta;
    zeta   = new GridFunction(z0);
    coeff_zeta = &coeff;
@@ -2270,6 +2393,33 @@ void TMOP_Integrator::EnableAdaptiveLimiting(const ParGridFunction &z0,
 }
 #endif
 
+void TMOP_Integrator::UpdateAfterMeshTopologyChange()
+{
+   if (zeta)
+   {
+      zeta->Update();
+      adapt_eval->SetSerialMetaInfo(*zeta->FESpace()->GetMesh(),
+                                    *zeta->FESpace()->FEColl(), 1);
+      adapt_eval->SetInitialField
+      (*zeta->FESpace()->GetMesh()->GetNodes(), *zeta);
+   }
+}
+
+#ifdef MFEM_USE_MPI
+void TMOP_Integrator::ParUpdateAfterMeshTopologyChange()
+{
+   if (zeta)
+   {
+      zeta->Update();
+      adapt_eval->SetParMetaInfo(*pzeta_0->ParFESpace()->GetParMesh(),
+                                 *pzeta_0->ParFESpace()->FEColl(), 1);
+      adapt_eval->SetInitialField
+      (*zeta->FESpace()->GetMesh()->GetNodes(), *zeta);
+   }
+}
+#endif
+
+
 double TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
                                          ElementTransformation &T,
                                          const Vector &elfun)
@@ -2378,6 +2528,145 @@ double TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
 
    return energy;
 }
+
+double TMOP_Integrator::GetRefinementElementEnergy(const FiniteElement &el,
+                                                   ElementTransformation &T,
+                                                   const Vector &elfun,
+                                                   const IntegrationRule &irule)
+{
+   int dof = el.GetDof(), dim = el.GetDim(),
+       NEsplit = elfun.Size() / (dof*dim), el_id = T.ElementNo;
+   double energy = 0.;
+
+   TargetConstructor *tc = const_cast<TargetConstructor *>(targetC);
+   DiscreteAdaptTC *dtc = dynamic_cast<DiscreteAdaptTC *>(tc);
+   // For DiscreteAdaptTC the GridFunctions used to set the targets must be
+   // mapped onto the fine elements.
+   if (dtc) { dtc->SetTspecFromIntRule(el_id, irule); }
+
+   for (int e = 0; e < NEsplit; e++)
+   {
+      DSh.SetSize(dof, dim);
+      Jrt.SetSize(dim);
+      Jpr.SetSize(dim);
+      Jpt.SetSize(dim);
+      Vector elfun_child(dof*dim);
+      for (int i = 0; i < dof; i++)
+      {
+         for (int d = 0; d < dim; d++)
+         {
+            // elfun is (xe1,xe2,...xen,ye1,ye2...yen) and has nodal coordinates
+            // for all the children element of the parent element being considered.
+            // So we must index and get (xek, yek) i.e. nodal coordinates for
+            // the fine element being considered.
+            elfun_child(i + d*dof) = elfun(i + e*dof + d*dof*NEsplit);
+         }
+      }
+      PMatI.UseExternalData(elfun_child.GetData(), dof, dim);
+
+      const IntegrationRule &ir = EnergyIntegrationRule(el);
+
+      double el_energy = 0;
+      DenseTensor Jtr(dim, dim, ir.GetNPoints());
+      if (dtc)
+      {
+         // This is used to index into the tspec vector inside DiscreteAdaptTC.
+         dtc->SetRefinementSubElement(e);
+      }
+      targetC->ComputeElementTargets(el_id, el, ir, elfun_child, Jtr);
+
+      // Define ref->physical transformation, wn a Coefficient is specified.
+      IsoparametricTransformation *Tpr = NULL;
+      if (coeff1 || coeff0)
+      {
+         Tpr = new IsoparametricTransformation;
+         Tpr->SetFE(&el);
+         Tpr->ElementNo = T.ElementNo;
+         Tpr->ElementType = ElementTransformation::ELEMENT;
+         Tpr->Attribute = T.Attribute;
+         Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
+      }
+
+      for (int i = 0; i < ir.GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir.IntPoint(i);
+         const DenseMatrix &Jtr_i = Jtr(i);
+         h_metric->SetTargetJacobian(Jtr_i);
+         CalcInverse(Jtr_i, Jrt);
+         const double weight = ip.weight * Jtr_i.Det();
+
+         el.CalcDShape(ip, DSh);
+         MultAtB(PMatI, DSh, Jpr);
+         Mult(Jpr, Jrt, Jpt);
+
+         double val = metric_normal * h_metric->EvalW(Jpt);
+         if (coeff1) { val *= coeff1->Eval(*Tpr, ip); }
+
+         el_energy += weight * val;
+         delete Tpr;
+      }
+      energy += el_energy;
+   }
+   energy /= NEsplit;
+
+   if (dtc) { dtc->ResetRefinementTspecData(); }
+
+   return energy;
+}
+
+double TMOP_Integrator::GetDerefinementElementEnergy(const FiniteElement &el,
+                                                     ElementTransformation &T,
+                                                     const Vector &elfun)
+{
+   int dof = el.GetDof(), dim = el.GetDim();
+   double energy = 0.;
+
+   DSh.SetSize(dof, dim);
+   Jrt.SetSize(dim);
+   Jpr.SetSize(dim);
+   Jpt.SetSize(dim);
+   PMatI.UseExternalData(elfun.GetData(), dof, dim);
+
+   const IntegrationRule &ir = EnergyIntegrationRule(el);
+
+   energy = 0.0;
+   DenseTensor Jtr(dim, dim, ir.GetNPoints());
+   targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
+
+   // Define ref->physical transformation, wn a Coefficient is specified.
+   IsoparametricTransformation *Tpr = NULL;
+   if (coeff1)
+   {
+      Tpr = new IsoparametricTransformation;
+      Tpr->SetFE(&el);
+      Tpr->ElementNo = T.ElementNo;
+      Tpr->ElementType = ElementTransformation::ELEMENT;
+      Tpr->Attribute = T.Attribute;
+      Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
+   }
+
+   for (int i = 0; i < ir.GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir.IntPoint(i);
+      const DenseMatrix &Jtr_i = Jtr(i);
+      h_metric->SetTargetJacobian(Jtr_i);
+      CalcInverse(Jtr_i, Jrt);
+      const double weight = ip.weight * Jtr_i.Det();
+
+      el.CalcDShape(ip, DSh);
+      MultAtB(PMatI, DSh, Jpr);
+      Mult(Jpr, Jrt, Jpt);
+
+      double val = metric_normal * h_metric->EvalW(Jpt);
+      if (coeff1) { val *= coeff1->Eval(*Tpr, ip); }
+
+      energy += weight * val;
+   }
+
+   delete Tpr;
+   return energy;
+}
+
 void TMOP_Integrator::AssembleElementVector(const FiniteElement &el,
                                             ElementTransformation &T,
                                             const Vector &elfun, Vector &elvect)
@@ -3039,7 +3328,7 @@ void TMOP_Integrator::ComputeMinJac(const Vector &x,
    dx = detv_avg_min / dxscale;
 }
 
-void TMOP_Integrator::UpdateAfterMeshChange(const Vector &new_x)
+void TMOP_Integrator::UpdateAfterMeshPositionChange(const Vector &new_x)
 {
    if (discr_tc)
    {
@@ -3168,6 +3457,32 @@ void TMOPComboIntegrator::AssembleElementGrad(const FiniteElement &el,
    }
 }
 
+double TMOPComboIntegrator::GetRefinementElementEnergy(const FiniteElement &el,
+                                                       ElementTransformation &T,
+                                                       const Vector &elfun,
+                                                       const IntegrationRule &irule)
+{
+   double energy= 0.0;
+   for (int i = 0; i < tmopi.Size(); i++)
+   {
+      energy += tmopi[i]->GetRefinementElementEnergy(el, T, elfun, irule);
+   }
+   return energy;
+}
+
+double TMOPComboIntegrator::GetDerefinementElementEnergy(
+   const FiniteElement &el,
+   ElementTransformation &T,
+   const Vector &elfun)
+{
+   double energy= 0.0;
+   for (int i = 0; i < tmopi.Size(); i++)
+   {
+      energy += tmopi[i]->GetDerefinementElementEnergy(el, T, elfun);
+   }
+   return energy;
+}
+
 void TMOPComboIntegrator::EnableNormalization(const GridFunction &x)
 {
    const int cnt = tmopi.Size();

fem/tmop.hpp

@@ -1057,14 +1057,31 @@ protected:
    // eta1(x+h,y), eta2(x+h,y) ... etan(x+h,y), eta1(x,y+h), eta2(x,y+h) ...
    // same for tspec_pert2h and tspec_pertmix.
 
+   // DenseMatrix to hold target_spec values for the (children of the)
+   // element being refined to consider for h-refinement.
+   DenseMatrix tspec_refine;
+   // Vector to hold the target_spec values for the coarse version of the
+   // current mesh. Used for derefinement decision with hr-adaptivity.
+   Vector tspec_derefine;
+
    // Components of Target Jacobian at each quadrature point of an element. This
    // is required for computation of the derivative using chain rule.
    mutable DenseTensor Jtrcomp;
 
    // Note: do not use the Nodes of this space as they may not be on the
    // positions corresponding to the values of tspec.
-   const FiniteElementSpace *tspec_fes;
-   const FiniteElementSpace *tspec_fesv;
+   FiniteElementSpace *tspec_fesv;         //owned
+   FiniteElementSpace *coarse_tspec_fesv;  //not owned, derefinement FESpace
+   GridFunction *tspec_gf;                 //owned, uses tspec and tspec_fes
+   // discrete adaptivity
+#ifdef MFEM_USE_MPI
+   ParFiniteElementSpace *ptspec_fesv;     //owned, needed for derefinement to
+   // get update operator.
+   ParGridFunction *tspec_pgf;             // similar to tspec_gf
+#endif
+
+   int amr_el;
+   double lim_min_size;
 
    // These flags can be used by outside functions to avoid recomputing the
    // tspec and tspec_perth fields again on the same mesh.
@@ -1076,7 +1093,7 @@ protected:
 
    void SetDiscreteTargetBase(const GridFunction &tspec_);
    void SetTspecAtIndex(int idx, const GridFunction &tspec_);
-   void FinalizeSerialDiscreteTargetSpec();
+   void FinalizeSerialDiscreteTargetSpec(const GridFunction &tspec_);
 #ifdef MFEM_USE_MPI
    void SetTspecAtIndex(int idx, const ParGridFunction &tspec_);
    void FinalizeParDiscreteTargetSpec(const ParGridFunction &tspec_);
@@ -1088,16 +1105,16 @@ public:
         ncomp(0),
         sizeidx(-1), skewidx(-1), aspectratioidx(-1), orientationidx(-1),
         tspec(), tspec_sav(), tspec_pert1h(), tspec_pert2h(), tspec_pertmix(),
-        tspec_fes(NULL), tspec_fesv(NULL),
+        tspec_refine(), tspec_derefine(),
+        tspec_fesv(NULL), coarse_tspec_fesv(NULL), tspec_gf(NULL),
+#ifdef MFEM_USE_MPI
+        ptspec_fesv(NULL), tspec_pgf(NULL),
+#endif
+        amr_el(-1), lim_min_size(-0.1),
         good_tspec(false), good_tspec_grad(false), good_tspec_hess(false),
         adapt_eval(NULL) { }
 
-   virtual ~DiscreteAdaptTC()
-   {
-      delete adapt_eval;
-      delete tspec_fes;
-      delete tspec_fesv;
-   }
+   virtual ~DiscreteAdaptTC();
 
    /** @name Target specification methods.
        The following methods are used to specify geometric parameters of the
@@ -1128,6 +1145,20 @@ public:
    void ResetUpdateFlags()
    { good_tspec = good_tspec_grad = good_tspec_hess = false; }
 
+   /// Get one of the discrete fields from tspec.
+   void GetDiscreteTargetSpec(GridFunction &tspec_, int idx);
+   /// Get the FESpace associated with tspec.
+   FiniteElementSpace *GetTSpecFESpace() { return tspec_fesv; }
+   /// Get the entire tspec.
+   GridFunction *GetTSpecData() { return tspec_gf; }
+   /// Update all discrete fields based on tspec and update for AMR
+   void UpdateAfterMeshTopologyChange();
+
+#ifdef MFEM_USE_MPI
+   ParFiniteElementSpace *GetTSpecParFESpace() { return ptspec_fesv; }
+   void ParUpdateAfterMeshTopologyChange();
+#endif
+
    /** Used to update the target specification after the mesh has changed. The
        new mesh positions are given by new_x. If @a use_flags is true, repeated
        calls won't do anything until ResetUpdateFlags() is called. */
@@ -1184,6 +1215,36 @@ public:
                                               const Vector &elfun,
                                               IsoparametricTransformation &Tpr,
                                               DenseTensor &dJtr) const;
+
+   // Generates tspec_vals for target construction using intrule
+   // Used for the refinement component in hr-adaptivity.
+   void SetTspecFromIntRule(int e_id, const IntegrationRule &intrule);
+
+   // Targets based on discrete functions can result in invalid (negative)
+   // size at the quadrature points. This method can be used to set a
+   // minimum target size.
+   void SetMinSizeForTargets(double min_size_) { lim_min_size = min_size_; }
+
+   /// Computes target specification data with respect to the coarse FE space.
+   void SetTspecDataForDerefinement(FiniteElementSpace *fes);
+
+   // Reset refinement data associated with h-adaptivity component.
+   void ResetRefinementTspecData()
+   {
+      tspec_refine.Clear();
+      amr_el = -1;
+   }
+
+   // Reset derefinement data associated with h-adaptivity component.
+   void ResetDerefinementTspecData()
+   {
+      tspec_derefine.Destroy();
+      coarse_tspec_fesv = NULL;
+   }
+
+   // Used to specify the fine element for determining energy of children of a
+   // parent element.
+   void SetRefinementSubElement(int amr_el_) { amr_el = amr_el_; }
 };
 
 class TMOPNewtonSolver;
@@ -1201,6 +1262,7 @@ protected:
    friend class TMOPNewtonSolver;
    friend class TMOPComboIntegrator;
 
+   TMOP_QualityMetric *h_metric;
    TMOP_QualityMetric *metric;        // not owned
    const TargetConstructor *targetC;  // not owned
 
@@ -1227,6 +1289,9 @@ protected:
 
    // Adaptive limiting.
    const GridFunction *zeta_0;       // Not owned.
+#ifdef MFEM_USE_MPI
+   const ParGridFunction *pzeta_0;
+#endif
    GridFunction *zeta;               // Owned. Updated by adapt_eval.
    Coefficient *coeff_zeta;          // Not owned.
    AdaptivityEvaluator *adapt_eval;  // Not owned.
@@ -1337,7 +1402,7 @@ protected:
 #endif
    void ComputeMinJac(const Vector &x, const FiniteElementSpace &fes);
 
-   void UpdateAfterMeshChange(const Vector &new_x);
+   void UpdateAfterMeshPositionChange(const Vector &new_x);
 
    void DisableLimiting()
    {
@@ -1395,11 +1460,13 @@ protected:
    void ComputeAllElementTargets(const Vector &xe = Vector()) const;
 
 public:
-   /** @param[in] m  TMOP_QualityMetric that will be integrated (not owned).
-       @param[in] tc Target-matrix construction algorithm to use (not owned). */
-   TMOP_Integrator(TMOP_QualityMetric *m, TargetConstructor *tc)
-      : metric(m), targetC(tc), IntegRules(NULL), integ_order(-1),
-        coeff1(NULL), metric_normal(1.0),
+   /** @param[in] m    TMOP_QualityMetric for r-adaptivity (not owned).
+       @param[in] tc   Target-matrix construction algorithm to use (not owned).
+       @param[in] hm   TMOP_QualityMetric for h-adaptivity (not owned). */
+   TMOP_Integrator(TMOP_QualityMetric *m, TargetConstructor *tc,
+                   TMOP_QualityMetric *hm)
+      : h_metric(hm), metric(m), targetC(tc), IntegRules(NULL),
+        integ_order(-1), coeff1(NULL), metric_normal(1.0),
         nodes0(NULL), coeff0(NULL),
         lim_dist(NULL), lim_func(NULL), lim_normal(1.0),
         zeta_0(NULL), zeta(NULL), coeff_zeta(NULL), adapt_eval(NULL),
@@ -1407,6 +1474,9 @@ public:
         fdflag(false), dxscale(1.0e3), fd_call_flag(false), exact_action(false)
    { PA.enabled = false; }
 
+   TMOP_Integrator(TMOP_QualityMetric *m, TargetConstructor *tc)
+      : TMOP_Integrator(m, tc, m) { }
+
    ~TMOP_Integrator();
 
    /// Release the device memory of large PA allocations. This will copy device
@@ -1478,6 +1548,22 @@ public:
                                    ElementTransformation &T,
                                    const Vector &elfun);
 
+   /** @brief Computes the mean of the energies of the given element's children.
+
+       In addition to the inputs for GetElementEnergy, this function requires an
+       IntegrationRule to be specified that will give the decomposition of the
+       given element based on the refinement type being considered. */
+   virtual double GetRefinementElementEnergy(const FiniteElement &el,
+                                             ElementTransformation &T,
+                                             const Vector &elfun,
+                                             const IntegrationRule &irule);
+
+   /// This function is similar to GetElementEnergy, but ignores components
+   /// such as limiting etc. to compute the element energy.
+   virtual double GetDerefinementElementEnergy(const FiniteElement &el,
+                                               ElementTransformation &T,
+                                               const Vector &elfun);
+
    virtual void AssembleElementVector(const FiniteElement &el,
                                       ElementTransformation &T,
                                       const Vector &elfun, Vector &elvect);
@@ -1486,6 +1572,13 @@ public:
                                     ElementTransformation &T,
                                     const Vector &elfun, DenseMatrix &elmat);
 
+   TMOP_QualityMetric &GetAMRQualityMetric() { return *h_metric; }
+
+   void UpdateAfterMeshTopologyChange();
+#ifdef MFEM_USE_MPI
+   void ParUpdateAfterMeshTopologyChange();
+#endif
+
    // PA extension
    using NonlinearFormIntegrator::AssemblePA;
    virtual void AssemblePA(const FiniteElementSpace&);
@@ -1564,6 +1657,15 @@ public:
                                     ElementTransformation &T,
                                     const Vector &elfun, DenseMatrix &elmat);
 
+   virtual double GetRefinementElementEnergy(const FiniteElement &el,
+                                             ElementTransformation &T,
+                                             const Vector &elfun,
+                                             const IntegrationRule &irule);
+
+   virtual double GetDerefinementElementEnergy(const FiniteElement &el,
+                                               ElementTransformation &T,
+                                               const Vector &elfun);
+
    /// Normalization factor that considers all integrators in the combination.
    void EnableNormalization(const GridFunction &x);
 #ifdef MFEM_USE_MPI

fem/tmop_amr.cpp

@@ -0,0 +1,896 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "tmop_amr.hpp"
+
+namespace mfem
+{
+
+using namespace mfem;
+
+void TMOPRefinerEstimator::ComputeEstimates()
+{
+   bool iso = false;
+   bool aniso = false;
+   if (amrmetric == 1 || amrmetric == 2 || amrmetric == 58)
+   {
+      aniso = true;
+   }
+   if (amrmetric == 55 || amrmetric == 56 || amrmetric == 77 ||
+       amrmetric == 315 || amrmetric == 316 || amrmetric == 321)
+   {
+      iso = true;
+   }
+   if (amrmetric == 7 || amrmetric == 9)
+   {
+      iso = true; aniso = true;
+   }
+
+   MFEM_VERIFY(iso || aniso, "Metric type not supported in hr-adaptivity.");
+
+   const int dim = mesh->Dimension();
+   const int num_ref_types = 3 + 4*(dim-2);
+   const int NEorig = mesh->GetNE();
+
+   aniso_flags.SetSize(NEorig);
+   error_estimates.SetSize(NEorig);
+   Vector amr_base_energy(NEorig), amr_temp_energy(NEorig);
+   error_estimates = 1.*std::numeric_limits<float>::max();
+   aniso_flags = -1;
+   GetTMOPRefinementEnergy(0, amr_base_energy);
+
+   for (int i = 1; i < num_ref_types+1; i++)
+   {
+      if ( dim == 2 && i < 3  && aniso != true ) { continue; }
+      if ( dim == 2 && i == 3 && iso   != true ) { continue; }
+      if ( dim == 3 && i < 7  && aniso != true ) { continue; }
+      if ( dim == 3 && i == 7 && iso   != true ) { continue; }
+
+      GetTMOPRefinementEnergy(i, amr_temp_energy);
+
+      for (int e = 0; e < NEorig; e++)
+      {
+         if ( amr_temp_energy(e) < error_estimates(e) )
+         {
+            error_estimates(e) = amr_temp_energy(e);
+            aniso_flags[e] = i;
+         }
+      }
+   }
+   error_estimates *= energy_scaling_factor;
+
+   if (spat_gf)
+   {
+      L2_FECollection avg_fec(0, mesh->Dimension());
+      FiniteElementSpace avg_fes(spat_gf->FESpace()->GetMesh(), &avg_fec);
+      GridFunction elem_avg(&avg_fes);
+      spat_gf->GetElementAverages(elem_avg);
+      for (int i = 0; i < amr_base_energy.Size(); i++)
+      {
+         if (elem_avg(i) < spat_gf_critical) { amr_base_energy(i) = 0.; }
+      }
+   }
+
+   error_estimates -= amr_base_energy;
+   error_estimates *= -1; // error = E(parent) - scaling_factor*mean(E(children))
+   current_sequence = mesh->GetSequence();
+}
+
+void TMOPRefinerEstimator::GetTMOPRefinementEnergy(int reftype,
+                                                   Vector &el_energy_vec)
+{
+   const FiniteElementSpace *fes = mesh->GetNodalFESpace();
+   const int NE = fes->GetNE();
+   GridFunction *xdof = mesh->GetNodes();
+   xdof->SetTrueVector();
+   xdof->SetFromTrueVector();
+
+   el_energy_vec.SetSize(NE);
+   el_energy_vec = std::numeric_limits<float>::max();
+
+   for (int e = 0; e < NE; e++)
+   {
+      Geometry::Type gtype = fes->GetFE(e)->GetGeomType();
+      DenseMatrix tr, xsplit;
+      IntegrationRule *irule = NULL;
+
+      if ( (gtype == Geometry::TRIANGLE && reftype > 0 && reftype < 3) ||
+           (gtype == Geometry::CUBE && reftype > 0 && reftype < 7) ||
+           (gtype == Geometry::TETRAHEDRON && reftype > 0 && reftype < 7) )
+      {
+         continue;
+      }
+
+      switch (gtype)
+      {
+         case Geometry::TRIANGLE:
+         {
+            int ref_access = reftype == 0 ? 0 : 1;
+            xdof->GetVectorValues(e, *TriIntRule[ref_access], xsplit, tr);
+            irule = TriIntRule[ref_access];
+            break;
+         }
+         case Geometry::TETRAHEDRON:
+         {
+            int ref_access = reftype == 0 ? 0 : 1;
+            xdof->GetVectorValues(e, *TetIntRule[ref_access], xsplit, tr);
+            irule = TetIntRule[ref_access];
+            break;
+         }
+         case Geometry::SQUARE:
+         {
+            MFEM_VERIFY(QuadIntRule[reftype], " Integration rule does not exist.");
+            xdof->GetVectorValues(e, *QuadIntRule[reftype], xsplit, tr);
+            irule = QuadIntRule[reftype];
+            break;
+         }
+         case Geometry::CUBE:
+         {
+            int ref_access = reftype == 0 ? 0 : 1;
+            xdof->GetVectorValues(e, *HexIntRule[ref_access], xsplit, tr);
+            irule = HexIntRule[ref_access];
+            break;
+         }
+         default:
+            MFEM_ABORT("Incompatible geometry type!");
+      }
+      xsplit.Transpose();
+
+      el_energy_vec(e) = 0.; // Re-set to 0
+
+      // The data format is xe1,xe2,..xen,ye1,ye2..yen.
+      // We will reformat it inside GetRefinementElementEnergy
+      Vector elfun(xsplit.GetData(), xsplit.NumCols()*xsplit.NumRows());
+
+      Array<NonlinearFormIntegrator*> &integs = *(nlf->GetDNFI());
+      TMOP_Integrator *ti  = NULL;
+      TMOPComboIntegrator *co = NULL;
+      for (int i = 0; i < integs.Size(); i++)
+      {
+         ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
+         if (ti)
+         {
+            el_energy_vec(e) = ti->GetRefinementElementEnergy(*fes->GetFE(e),
+                                                              *mesh->GetElementTransformation(e),
+                                                              elfun,
+                                                              *irule);
+         }
+         co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
+         if (co)
+         {
+            Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
+            for (int j = 0; j < ati.Size(); j++)
+            {
+               el_energy_vec(e) += ati[j]->GetRefinementElementEnergy(*fes->GetFE(e),
+                                                                      *mesh->GetElementTransformation(e),
+                                                                      elfun,
+                                                                      *irule);
+            }
+         }
+      }
+   }
+}
+
+void TMOPRefinerEstimator::SetHexIntRules()
+{
+   HexIntRule.SetSize(1+1);
+   // Reftype = 0 -> original element
+   Mesh meshsplit = Mesh::MakeCartesian3D(1, 1, 1, Element::HEXAHEDRON);
+   Mesh base_mesh_copy(meshsplit);
+   HexIntRule[0] = SetIntRulesFromMesh(meshsplit);
+   meshsplit.Clear();
+
+   // Reftype = 7
+   for (int i = 7; i < 8; i++)
+   {
+      Array<Refinement> marked_elements;
+      Mesh mesh_ref(base_mesh_copy);
+      for (int e = 0; e < mesh_ref.GetNE(); e++)
+      {
+         marked_elements.Append(Refinement(e, i));
+      }
+      mesh_ref.GeneralRefinement(marked_elements, 1, 0);
+      HexIntRule[1] = SetIntRulesFromMesh(mesh_ref);
+      mesh_ref.Clear();
+   }
+}
+
+void TMOPRefinerEstimator::SetQuadIntRules()
+{
+   QuadIntRule.SetSize(3+1);
+
+   // Reftype = 0 -> original element
+   Mesh meshsplit = Mesh::MakeCartesian2D(1, 1, Element::QUADRILATERAL);
+   Mesh base_mesh_copy(meshsplit);
+   QuadIntRule[0] = SetIntRulesFromMesh(meshsplit);
+   meshsplit.Clear();
+
+   // Reftype = 1-3
+   for (int i = 1; i < 4; i++)
+   {
+      Array<Refinement> marked_elements;
+      Mesh mesh_ref(base_mesh_copy);
+      for (int e = 0; e < mesh_ref.GetNE(); e++)
+      {
+         marked_elements.Append(Refinement(e, i));
+      }
+      mesh_ref.GeneralRefinement(marked_elements, 1, 0);
+      QuadIntRule[i] = SetIntRulesFromMesh(mesh_ref);
+      mesh_ref.Clear();
+   }
+}
+
+void TMOPRefinerEstimator::SetTriIntRules()
+{
+   TriIntRule.SetSize(1+1);
+
+   // Reftype = 0 // original element
+   const int Nvert = 3, NEsplit = 1;
+   Mesh meshsplit(2, Nvert, NEsplit, 0 ,2);
+   const double tri_v[3][2] =
+   {
+      {0, 0}, {1, 0}, {0, 1}
+   };
+   const int tri_e[1][3] =
+   {
+      {0, 1, 2}
+   };
+
+   for (int j = 0; j < Nvert; j++)
+   {
+      meshsplit.AddVertex(tri_v[j]);
+   }
+   meshsplit.AddTriangle(tri_e[0], 1);
+   meshsplit.FinalizeTriMesh(1, 1, true);
+
+   Mesh base_mesh_copy(meshsplit);
+   TriIntRule[0] = SetIntRulesFromMesh(meshsplit);
+   meshsplit.Clear();
+
+   // no anisotropic refinements for triangle
+   // Reftype = 3
+   for (int i = 1; i < 2; i++)
+   {
+      Array<Refinement> marked_elements;
+      Mesh mesh_ref(base_mesh_copy);
+      for (int e = 0; e < mesh_ref.GetNE(); e++)
+      {
+         marked_elements.Append(Refinement(e, i));
+      }
+      mesh_ref.GeneralRefinement(marked_elements, 1, 0);
+      TriIntRule[i] = SetIntRulesFromMesh(mesh_ref);
+      mesh_ref.Clear();
+   }
+}
+
+void TMOPRefinerEstimator::SetTetIntRules()
+{
+   TetIntRule.SetSize(1+1);
+
+   // Reftype = 0 // original element
+   const int Nvert = 4, NEsplit = 1;
+   Mesh meshsplit(3, Nvert, NEsplit, 0, 3);
+   const double tet_v[4][3] =
+   {
+      {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}
+   };
+   const int tet_e[1][4] =
+   {
+      {0, 1, 2, 3}
+   };
+
+   for (int j = 0; j < Nvert; j++)
+   {
+      meshsplit.AddVertex(tet_v[j]);
+   }
+   meshsplit.AddTet(tet_e[0], 1);
+   meshsplit.FinalizeTetMesh(1, 1, true);
+
+   Mesh base_mesh_copy(meshsplit);
+   TetIntRule[0] = SetIntRulesFromMesh(meshsplit);
+   meshsplit.Clear();
+
+   // no anisotropic refinements for triangle
+   // Reftype = 7
+   for (int i = 1; i < 2; i++)
+   {
+      Array<Refinement> marked_elements;
+      Mesh mesh_ref(base_mesh_copy);
+      for (int e = 0; e < mesh_ref.GetNE(); e++)
+      {
+         marked_elements.Append(Refinement(e, i)); //ref_type will default to 7
+      }
+      mesh_ref.GeneralRefinement(marked_elements, 1, 0);
+      TetIntRule[i] = SetIntRulesFromMesh(mesh_ref);
+      mesh_ref.Clear();
+   }
+}
+
+IntegrationRule* TMOPRefinerEstimator::SetIntRulesFromMesh(Mesh &meshsplit)
+{
+   const int dim = meshsplit.Dimension();
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace nodal_fes(&meshsplit, &fec, dim);
+   meshsplit.SetNodalFESpace(&nodal_fes);
+
+   const int NEsplit = meshsplit.GetNE();
+   const int dof_cnt = nodal_fes.GetFE(0)->GetDof(),
+             pts_cnt = NEsplit * dof_cnt;
+
+   DenseMatrix pos(dof_cnt, dim);
+   Vector posV(pos.Data(), dof_cnt * dim);
+   Array<int> xdofs(dof_cnt * dim);
+
+   // Create an IntegrationRule on the nodes of the reference submesh.
+   IntegrationRule *irule = new IntegrationRule(pts_cnt);
+   GridFunction *nodesplit = meshsplit.GetNodes();
+
+   int pt_id = 0;
+   for (int i = 0; i < NEsplit; i++)
+   {
+      nodal_fes.GetElementVDofs(i, xdofs);
+      nodesplit->GetSubVector(xdofs, posV);
+      for (int j = 0; j < dof_cnt; j++)
+      {
+         if (dim == 2)
+         {
+            irule->IntPoint(pt_id).Set2(pos(j, 0), pos(j, 1));
+         }
+         else if (dim == 3)
+         {
+            irule->IntPoint(pt_id).Set3(pos(j, 0), pos(j, 1), pos(j, 2));
+         }
+         pt_id++;
+      }
+   }
+   return irule;
+}
+
+bool TMOPDeRefinerEstimator::GetDerefineEnergyForIntegrator(
+   TMOP_Integrator &tmopi,
+   Vector &fine_energy)
+{
+   DiscreteAdaptTC *tcd = tmopi.GetDiscreteAdaptTC();
+   fine_energy.SetSize(mesh->GetNE());
+
+   if (serial)
+   {
+      Mesh meshcopy(*mesh);
+      FiniteElementSpace *tcdfes = NULL;
+      if (tcd)
+      {
+         tcdfes = new FiniteElementSpace(*tcd->GetTSpecFESpace(), &meshcopy);
+      }
+
+      Vector local_err(meshcopy.GetNE());
+      local_err = 0.;
+      double threshold = std::numeric_limits<float>::max();
+      meshcopy.DerefineByError(local_err, threshold, 0, 1);
+
+      if (meshcopy.GetGlobalNE() == mesh->GetGlobalNE())
+      {
+         delete tcdfes;
+         return false;
+      }
+
+      if (tcd)
+      {
+         tcdfes->Update();
+         tcd->SetTspecDataForDerefinement(tcdfes);
+      }
+
+      Vector coarse_energy(meshcopy.GetNE());
+      GetTMOPDerefinementEnergy(meshcopy, tmopi, coarse_energy);
+      if (tcd) { tcd->ResetDerefinementTspecData(); }
+      GetTMOPDerefinementEnergy(*mesh, tmopi, fine_energy);
+
+      const CoarseFineTransformations &dtrans =
+         meshcopy.ncmesh->GetDerefinementTransforms();
+      Table coarse_to_fine;
+      dtrans.GetCoarseToFineMap(meshcopy, coarse_to_fine);
+
+      for (int pe = 0; pe < coarse_to_fine.Size(); pe++)
+      {
+         Array<int> tabrow;
+         coarse_to_fine.GetRow(pe, tabrow);
+         int nchild = tabrow.Size();
+         double parent_energy = coarse_energy(pe);
+         for (int fe = 0; fe < nchild; fe++)
+         {
+            int child = tabrow[fe];
+            MFEM_VERIFY(child < mesh->GetNE(), " invalid coarse to fine mapping");
+            fine_energy(child) -= parent_energy;
+         }
+      }
+      delete tcdfes;
+   }
+   else
+   {
+#ifdef MFEM_USE_MPI
+      ParMesh meshcopy(*pmesh);
+      ParFiniteElementSpace *tcdfes = NULL;
+      if (tcd)
+      {
+         tcdfes = new ParFiniteElementSpace(*tcd->GetTSpecParFESpace(), meshcopy);
+      }
+
+      Vector local_err(meshcopy.GetNE());
+      local_err = 0.;
+      double threshold = std::numeric_limits<float>::max();
+      meshcopy.DerefineByError(local_err, threshold, 0, 1);
+
+      if (meshcopy.GetGlobalNE() == pmesh->GetGlobalNE())
+      {
+         delete tcdfes;
+         return false;
+      }
+
+      if (tcd)
+      {
+         tcdfes->Update();
+         tcd->SetTspecDataForDerefinement(tcdfes);
+      }
+
+      Vector coarse_energy(meshcopy.GetNE());
+      GetTMOPDerefinementEnergy(meshcopy, tmopi, coarse_energy);
+      if (tcd) { tcd->ResetDerefinementTspecData(); }
+      GetTMOPDerefinementEnergy(*pmesh, tmopi, fine_energy);
+
+      const CoarseFineTransformations &dtrans =
+         meshcopy.pncmesh->GetDerefinementTransforms();
+      Table coarse_to_fine;
+      dtrans.GetCoarseToFineMap(meshcopy, coarse_to_fine);
+
+      for (int pe = 0; pe < meshcopy.GetNE(); pe++)
+      {
+         Array<int> tabrow;
+         coarse_to_fine.GetRow(pe, tabrow);
+         int nchild = tabrow.Size();
+         double parent_energy = coarse_energy(pe);
+         for (int fe = 0; fe < nchild; fe++)
+         {
+            int child = tabrow[fe];
+            MFEM_VERIFY(child < pmesh->GetNE(), " invalid coarse to fine mapping");
+            fine_energy(child) -= parent_energy;
+         }
+      }
+      delete tcdfes;
+#endif
+   }
+
+   // error_estimate(e) = energy(parent_of_e)-energy(e)
+   // Negative energy means derefinement is desirable.
+   fine_energy *= -1;
+   return true;
+}
+
+void TMOPDeRefinerEstimator::ComputeEstimates()
+{
+   Array<NonlinearFormIntegrator*> &integs = *(nlf->GetDNFI());
+   TMOP_Integrator *ti  = NULL;
+   TMOPComboIntegrator *co = NULL;
+   error_estimates.SetSize(mesh->GetNE());
+   error_estimates = 0.;
+   Vector fine_energy(mesh->GetNE());
+
+   for (int i = 0; i < integs.Size(); i++)
+   {
+      ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
+      if (ti)
+      {
+         bool deref = GetDerefineEnergyForIntegrator(*ti, fine_energy);
+         if (!deref) { error_estimates = 1; return; }
+         error_estimates += fine_energy;
+      }
+      co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
+      if (co)
+      {
+         Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
+         for (int j = 0; j < ati.Size(); j++)
+         {
+            bool deref = GetDerefineEnergyForIntegrator(*ati[j], fine_energy);
+            if (!deref) { error_estimates = 1; return; }
+            error_estimates += fine_energy;
+         }
+      }
+   }
+}
+
+void TMOPDeRefinerEstimator::GetTMOPDerefinementEnergy(Mesh &cmesh,
+                                                       TMOP_Integrator &tmopi,
+                                                       Vector &el_energy_vec)
+{
+   const int cNE = cmesh.GetNE();
+   el_energy_vec.SetSize(cNE);
+   const FiniteElementSpace *fespace = cmesh.GetNodalFESpace();
+
+   GridFunction *cxdof = cmesh.GetNodes();
+
+   Array<int> vdofs;
+   Vector el_x;
+   const FiniteElement *fe;
+   ElementTransformation *T;
+
+   for (int j = 0; j < cNE; j++)
+   {
+      fe = fespace->GetFE(j);
+      fespace->GetElementVDofs(j, vdofs);
+      T = cmesh.GetElementTransformation(j);
+      cxdof->GetSubVector(vdofs, el_x);
+      el_energy_vec(j) = tmopi.GetDerefinementElementEnergy(*fe, *T, el_x);
+   }
+}
+
+
+TMOPHRSolver::TMOPHRSolver(Mesh &mesh_, NonlinearForm &nlf_,
+                           TMOPNewtonSolver &tmopns_, GridFunction &x_,
+                           bool move_bnd_, bool hradaptivity_,
+                           int mesh_poly_deg_, int amr_metric_id_,
+                           int hr_iter_, int h_per_r_iter_) :
+   mesh(&mesh_), nlf(&nlf_), tmopns(&tmopns_), x(&x_),
+   gridfuncarr(), fespacearr(),
+   move_bnd(move_bnd_), hradaptivity(hradaptivity_),
+   mesh_poly_deg(mesh_poly_deg_), amr_metric_id(amr_metric_id_),
+   serial(true), hr_iter(hr_iter_), h_per_r_iter(h_per_r_iter_)
+{
+   if (!hradaptivity) { return; }
+   tmop_r_est = new TMOPRefinerEstimator(*mesh, *nlf, mesh_poly_deg,
+                                         amr_metric_id);
+   tmop_r = new ThresholdRefiner(*tmop_r_est);
+   tmop_r->SetTotalErrorFraction(0.0);
+   tmop_r_est->SetEnergyScalingFactor(1.);
+   tmop_dr_est= new TMOPDeRefinerEstimator(*mesh, *nlf);
+   tmop_dr = new ThresholdDerefiner(*tmop_dr_est);
+   AddGridFunctionForUpdate(x);
+}
+
+#ifdef MFEM_USE_MPI
+TMOPHRSolver::TMOPHRSolver(ParMesh &pmesh_, ParNonlinearForm &pnlf_,
+                           TMOPNewtonSolver &tmopns_, ParGridFunction &px_,
+                           bool move_bnd_, bool hradaptivity_,
+                           int mesh_poly_deg_, int amr_metric_id_,
+                           int hr_iter_, int h_per_r_iter_) :
+   mesh(&pmesh_), nlf(&pnlf_), tmopns(&tmopns_), x(&px_),
+   gridfuncarr(), fespacearr(),
+   move_bnd(move_bnd_), hradaptivity(hradaptivity_),
+   mesh_poly_deg(mesh_poly_deg_), amr_metric_id(amr_metric_id_),
+   pmesh(&pmesh_), pnlf(&pnlf_), pgridfuncarr(), pfespacearr(),
+   serial(false), hr_iter(hr_iter_), h_per_r_iter(h_per_r_iter_)
+{
+   if (!hradaptivity) { return; }
+   tmop_r_est = new TMOPRefinerEstimator(*pmesh, *pnlf, mesh_poly_deg,
+                                         amr_metric_id);
+   tmop_r = new ThresholdRefiner(*tmop_r_est);
+   tmop_r->SetTotalErrorFraction(0.0);
+   tmop_r_est->SetEnergyScalingFactor(1.);
+   tmop_dr_est= new TMOPDeRefinerEstimator(*pmesh, *pnlf);
+   tmop_dr = new ThresholdDerefiner(*tmop_dr_est);
+   AddGridFunctionForUpdate(&px_);
+}
+#endif
+
+void TMOPHRSolver::Mult()
+{
+   Vector b(0);
+   int myid = 0;
+   if (serial)
+   {
+      tmopns->SetOperator(*nlf);
+   }
+   else
+   {
+#ifdef MFEM_USE_MPI
+      myid = pnlf->ParFESpace()->GetMyRank();
+      tmopns->SetOperator(*pnlf);
+#endif
+   }
+   if (!hradaptivity)
+   {
+      tmopns->Mult(b, x->GetTrueVector());
+      if (tmopns->GetConverged() == false)
+      {
+         if (myid == 0) { mfem::out << "Nonlinear solver: rtol not achieved.\n"; }
+      }
+      x->SetFromTrueVector();
+      return;
+   }
+
+   bool radaptivity = true;
+
+   tmop_dr->Reset();
+   tmop_r->Reset();
+
+   if (serial)
+   {
+      for (int i_hr = 0; i_hr < hr_iter; i_hr++)
+      {
+         if (!radaptivity)
+         {
+            break;
+         }
+         mfem::out << i_hr << " r-adaptivity iteration.\n";
+
+         tmopns->SetOperator(*nlf);
+         tmopns->Mult(b, x->GetTrueVector());
+         x->SetFromTrueVector();
+
+         mfem::out << "TMOP energy after r-adaptivity: " <<
+                   nlf->GetGridFunctionEnergy(*x)/mesh->GetNE() <<
+                   ", Elements: " << mesh->GetNE() << std::endl;
+
+         for (int i_h = 0; i_h < h_per_r_iter; i_h++)
+         {
+            // Derefinement step.
+            if (mesh->ncmesh)
+            {
+               tmop_dr->Apply(*mesh);
+               Update();
+            }
+            mfem::out << "TMOP energy after derefinement: " <<
+                      nlf->GetGridFunctionEnergy(*x)/mesh->GetNE() <<
+                      ", Elements: " << mesh->GetNE() << std::endl;
+
+            // Refinement step.
+            tmop_r->Apply(*mesh);
+            Update();
+            mfem::out << "TMOP energy after   refinement: " <<
+                      nlf->GetGridFunctionEnergy(*x)/mesh->GetNE() <<
+                      ", Elements: " << mesh->GetNE() << std::endl;
+
+            if (!tmop_dr->Derefined() && tmop_r->Stop())
+            {
+               radaptivity = false;
+               mfem::out << "AMR stopping criterion satisfied. Stop.\n";
+               break;
+            }
+         } //n_h
+      } //n_hr
+   }
+   else
+   {
+#ifdef MFEM_USE_MPI
+      int NEGlob;
+      double tmopenergy;
+      for (int i_hr = 0; i_hr < hr_iter; i_hr++)
+      {
+         if (!radaptivity)
+         {
+            break;
+         }
+         if (myid == 0) { mfem::out << i_hr << " r-adaptivity iteration.\n"; }
+         tmopns->SetOperator(*pnlf);
+         tmopns->Mult(b, x->GetTrueVector());
+         x->SetFromTrueVector();
+
+         NEGlob = pmesh->GetGlobalNE();
+         tmopenergy = pnlf->GetParGridFunctionEnergy(*x) / NEGlob;
+         if (myid == 0)
+         {
+            mfem::out << "TMOP energy after r-adaptivity: " << tmopenergy <<
+                      ", Elements: " << NEGlob << std::endl;
+         }
+
+         for (int i_h = 0; i_h < h_per_r_iter; i_h++)
+         {
+            // Derefinement step.
+            if (pmesh->pncmesh)
+            {
+               RebalanceParNCMesh();
+               ParUpdate();
+
+               tmop_dr->Apply(*pmesh);
+               ParUpdate();
+            }
+            NEGlob = pmesh->GetGlobalNE();
+            tmopenergy = pnlf->GetParGridFunctionEnergy(*x) / NEGlob;
+            if (myid == 0)
+            {
+               mfem::out << "TMOP energy after derefinement: " << tmopenergy <<
+                         ", Elements: " << NEGlob << std::endl;
+            }
+
+            // Refinement step.
+            tmop_r->Apply(*pmesh);
+            ParUpdate();
+            NEGlob = pmesh->GetGlobalNE();
+            tmopenergy = pnlf->GetParGridFunctionEnergy(*x) / NEGlob;
+            if (myid == 0)
+            {
+               mfem::out << "TMOP energy after   refinement: " << tmopenergy <<
+                         ", Elements: " << NEGlob << std::endl;
+            }
+
+            if (!tmop_dr->Derefined() && tmop_r->Stop())
+            {
+               radaptivity = false;
+               if (myid == 0)
+               {
+                  mfem::out << "AMR stopping criterion satisfied. Stop.\n";
+               }
+               break;
+            }
+         } //n_r limit
+      } //n_hr
+#endif
+   }
+}
+
+#ifdef MFEM_USE_MPI
+void TMOPHRSolver::RebalanceParNCMesh()
+{
+   ParNCMesh *pncmesh = pmesh->pncmesh;
+   if (pncmesh)
+   {
+      const Table &dreftable = pncmesh->GetDerefinementTable();
+      Array<int> drefs, new_ranks;
+      for (int i = 0; i < dreftable.Size(); i++)
+      {
+         drefs.Append(i);
+      }
+      pncmesh->GetFineToCoarsePartitioning(drefs, new_ranks);
+      pmesh->Rebalance(new_ranks);
+   }
+}
+#endif
+
+void TMOPHRSolver::Update()
+{
+   // Update FESpace
+   for (int i = 0; i < fespacearr.Size(); i++)
+   {
+      fespacearr[i]->Update();
+   }
+   // Update nodal GF
+   for (int i = 0; i < gridfuncarr.Size(); i++)
+   {
+      gridfuncarr[i]->Update();
+      gridfuncarr[i]->SetTrueVector();
+      gridfuncarr[i]->SetFromTrueVector();
+   }
+
+   // Update Discrete Indicator for all the TMOP_Integrators in NonLinearForm
+   Array<NonlinearFormIntegrator*> &integs = *(nlf->GetDNFI());
+   TMOP_Integrator *ti  = NULL;
+   TMOPComboIntegrator *co = NULL;
+   DiscreteAdaptTC *dtc = NULL;
+   for (int i = 0; i < integs.Size(); i++)
+   {
+      ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
+      if (ti)
+      {
+         ti->UpdateAfterMeshTopologyChange();
+         dtc = ti->GetDiscreteAdaptTC();
+         if (dtc) { dtc->UpdateAfterMeshTopologyChange(); }
+      }
+      co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
+      if (co)
+      {
+         Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
+         for (int j = 0; j < ati.Size(); j++)
+         {
+            ati[j]->UpdateAfterMeshTopologyChange();
+            dtc = ati[j]->GetDiscreteAdaptTC();
+            if (dtc) { dtc->UpdateAfterMeshTopologyChange(); }
+         }
+      }
+   }
+
+   // Update the Nonlinear form and set Essential BC.
+   UpdateNonlinearFormAndBC(mesh, nlf);
+}
+
+#ifdef MFEM_USE_MPI
+void TMOPHRSolver::ParUpdate()
+{
+   // Update FESpace
+   for (int i = 0; i < pfespacearr.Size(); i++)
+   {
+      pfespacearr[i]->Update();
+   }
+   // Update nodal GF
+   for (int i = 0; i < pgridfuncarr.Size(); i++)
+   {
+      pgridfuncarr[i]->Update();
+      pgridfuncarr[i]->SetTrueVector();
+      pgridfuncarr[i]->SetFromTrueVector();
+   }
+
+   // Update Discrete Indicator
+   Array<NonlinearFormIntegrator*> &integs = *(nlf->GetDNFI());
+   TMOP_Integrator *ti  = NULL;
+   TMOPComboIntegrator *co = NULL;
+   DiscreteAdaptTC *dtc = NULL;
+   for (int i = 0; i < integs.Size(); i++)
+   {
+      ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
+      if (ti)
+      {
+         ti->ParUpdateAfterMeshTopologyChange();
+         dtc = ti->GetDiscreteAdaptTC();
+         if (dtc) { dtc->ParUpdateAfterMeshTopologyChange(); }
+      }
+      co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
+      if (co)
+      {
+         Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
+         for (int j = 0; j < ati.Size(); j++)
+         {
+            ati[j]->ParUpdateAfterMeshTopologyChange();
+            dtc = ati[j]->GetDiscreteAdaptTC();
+            if (dtc) { dtc->ParUpdateAfterMeshTopologyChange(); }
+         }
+      }
+   }
+
+   // Update the Nonlinear form and set Essential BC.
+   UpdateNonlinearFormAndBC(pmesh, pnlf);
+}
+#endif
+
+void TMOPHRSolver::UpdateNonlinearFormAndBC(Mesh *mesh, NonlinearForm *nlf)
+{
+   const FiniteElementSpace &fes = *mesh->GetNodalFESpace();
+
+   // Update Nonlinear form and Set Essential BC
+   nlf->Update();
+   const int dim = fes.GetFE(0)->GetDim();
+   if (move_bnd == false)
+   {
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      nlf->SetEssentialBC(ess_bdr);
+   }
+   else
+   {
+      const int nd  = fes.GetBE(0)->GetDof();
+      int n = 0;
+      for (int i = 0; i < mesh->GetNBE(); i++)
+      {
+         const int attr = mesh->GetBdrElement(i)->GetAttribute();
+         MFEM_VERIFY(!(dim == 2 && attr == 3),
+                     "Boundary attribute 3 must be used only for 3D meshes. "
+                     "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
+                     "components, rest for free nodes), or use -fix-bnd.");
+         if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
+         if (attr == 4) { n += nd * dim; }
+      }
+      Array<int> ess_vdofs(n), vdofs;
+      n = 0;
+      for (int i = 0; i < mesh->GetNBE(); i++)
+      {
+         const int attr = mesh->GetBdrElement(i)->GetAttribute();
+         fes.GetBdrElementVDofs(i, vdofs);
+         if (attr == 1) // Fix x components.
+         {
+            for (int j = 0; j < nd; j++)
+            { ess_vdofs[n++] = vdofs[j]; }
+         }
+         else if (attr == 2) // Fix y components.
+         {
+            for (int j = 0; j < nd; j++)
+            { ess_vdofs[n++] = vdofs[j+nd]; }
+         }
+         else if (attr == 3) // Fix z components.
+         {
+            for (int j = 0; j < nd; j++)
+            { ess_vdofs[n++] = vdofs[j+2*nd]; }
+         }
+         else if (attr == 4) // Fix all components.
+         {
+            for (int j = 0; j < vdofs.Size(); j++)
+            { ess_vdofs[n++] = vdofs[j]; }
+         }
+      }
+      nlf->SetEssentialVDofs(ess_vdofs);
+   }
+}
+
+}

fem/tmop_amr.hpp

@@ -0,0 +1,284 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_TMOP_AMR_HPP
+#define MFEM_TMOP_AMR_HPP
+
+#include "tmop_tools.hpp"
+#include "nonlinearform.hpp"
+#include "pnonlinearform.hpp"
+#include "estimators.hpp"
+#include "../mesh/mesh_operators.hpp"
+
+namespace mfem
+{
+
+class TMOPRefinerEstimator : public AnisotropicErrorEstimator
+{
+protected:
+   Mesh *mesh; // not owned
+   NonlinearForm *nlf; // not owned
+   int order;
+   int amrmetric;
+   Array<IntegrationRule *> TriIntRule, QuadIntRule, TetIntRule, HexIntRule;
+   long current_sequence;
+   Vector error_estimates;
+   Array<int> aniso_flags;
+   // An element is refined only if
+   // [mean TMOPEnergy(children)]*energy_scaling_factor < TMOPEnergy(parent)
+   double energy_scaling_factor;
+   GridFunction *spat_gf;   // If specified, can be used to specify the
+   double spat_gf_critical; // region where hr-adaptivity is done.
+
+   /// Check if the mesh of the solution was modified.
+   bool MeshIsModified()
+   {
+      long mesh_sequence = mesh->GetSequence();
+      MFEM_ASSERT(mesh_sequence >= current_sequence, "");
+      return (mesh_sequence > current_sequence);
+   }
+
+   /// Compute the element error estimates. For an element E in the mesh,
+   /// error(E) = TMOPEnergy(E)*energy_scaling_factor-Mean(TMOPEnergy(ChildofE)),
+   /// where TMOPEnergy of Children of E is obtained by assuming the element E
+   /// is refined using the refinement type being considered based on the TMOP
+   /// mesh quality metric.
+   void ComputeEstimates();
+
+   /// Construct the integration rules to model how each element type is split
+   /// using different refinement types. ref_type = 0 is the original element
+   /// and reftype \ in [1, 7] represent different refinement type based on
+   /// NCMesh class.
+   void SetQuadIntRules(); // supports ref_type = 1 to 3.
+   void SetTriIntRules(); // currently supports only isotropic refinement.
+   void SetHexIntRules(); // currently supports only isotropic refinement.
+   void SetTetIntRules(); // currently supports only isotropic refinement.
+
+   /// Get TMOP energy for each element corresponding to the refinement type
+   /// specified.
+   void GetTMOPRefinementEnergy(int reftype, Vector &el_energy_vec);
+
+   /// Use a mesh to setup an integration rule that will mimic the different
+   /// refinement types.
+   IntegrationRule* SetIntRulesFromMesh(Mesh &meshsplit);
+public:
+   TMOPRefinerEstimator(Mesh &mesh_, NonlinearForm &nlf_, int order_,
+                        int amrmetric_) :
+      mesh(&mesh_), nlf(&nlf_), order(order_), amrmetric(amrmetric_),
+      TriIntRule(0), QuadIntRule(0), TetIntRule(0), HexIntRule(0),
+      current_sequence(-1), error_estimates(), aniso_flags(),
+      energy_scaling_factor(1.), spat_gf(NULL), spat_gf_critical(0.)
+   {
+      if (mesh->Dimension() == 2)
+      {
+         SetQuadIntRules();
+         SetTriIntRules();
+      }
+      else
+      {
+         SetHexIntRules();
+         SetTetIntRules();
+      }
+   }
+
+   ~TMOPRefinerEstimator()
+   {
+      for (int i = 0; i < QuadIntRule.Size(); i++) { delete QuadIntRule[i]; }
+      for (int i = 0; i < TriIntRule.Size();  i++) { delete TriIntRule[i]; }
+      for (int i = 0; i < HexIntRule.Size();  i++) { delete HexIntRule[i]; }
+      for (int i = 0; i < TetIntRule.Size();  i++) { delete TetIntRule[i]; }
+   }
+
+   /// Get TMOP-based errors for each element in the mesh computed based on the
+   /// refinement types being considered.
+   virtual const Vector &GetLocalErrors()
+   {
+      if (MeshIsModified()) { ComputeEstimates(); }
+      return error_estimates;
+   }
+   /// For anisotropic refinements, get the refinement type (e.g., x or y)
+   virtual const Array<int> &GetAnisotropicFlags()
+   {
+      if (MeshIsModified()) { ComputeEstimates(); }
+      return aniso_flags;
+   }
+
+   /// Scaling factor for the TMOP refinement energy. An element is refined if
+   /// [mean TMOPEnergy(children)]*energy_scaling_factor < TMOPEnergy(parent)
+   void SetEnergyScalingFactor(double scale) { energy_scaling_factor = scale; }
+
+   /// Spatial indicator function (eta) that can be used to prevent elements
+   /// from being refined even if the energy criterion is met. Using this,
+   /// an element E is not refined if mean(@a spat_gf(E)) < @a spat_gf_critical.
+   void SetSpatialIndicator(GridFunction &spat_gf_,
+                            double spat_gf_critical_ = 0.5)
+   { spat_gf = &spat_gf_; spat_gf_critical = spat_gf_critical_; }
+   void SetSpatialIndicatorCritical(double val_) { spat_gf_critical = val_; }
+
+   /// Reset the error estimator.
+   virtual void Reset() { current_sequence = -1; }
+};
+
+class TMOPDeRefinerEstimator : public ErrorEstimator
+{
+protected:
+   Mesh *mesh;
+   NonlinearForm *nlf;
+#ifdef MFEM_USE_MPI
+   ParMesh *pmesh;
+   ParNonlinearForm *pnlf;
+#endif
+   int order;
+   int amrmetric;
+   long current_sequence;
+   Vector error_estimates;
+   bool serial;
+
+   /// Check if the mesh of the solution was modified.
+   bool MeshIsModified()
+   {
+      long mesh_sequence = mesh->GetSequence();
+      MFEM_ASSERT(mesh_sequence >= current_sequence, "");
+      return (mesh_sequence > current_sequence);
+   }
+
+   /// Compute the element error estimates. For a given element E in the mesh,
+   /// error(E) = TMOPEnergy(parent_of_E)-TMOPEnergy(E). Children element of an
+   /// element are derefined if the mean TMOP energy of children is greated than
+   /// the TMOP energy associated with their parent.
+   void ComputeEstimates();
+
+   void GetTMOPDerefinementEnergy(Mesh &cmesh,
+                                  TMOP_Integrator &tmopi,
+                                  Vector &el_energy_vec);
+
+   bool GetDerefineEnergyForIntegrator(TMOP_Integrator &tmopi,
+                                       Vector &fine_energy);
+public:
+   TMOPDeRefinerEstimator(Mesh &mesh_, NonlinearForm &nlf_) :
+      mesh(&mesh_), nlf(&nlf_),
+      current_sequence(-1), error_estimates(), serial(true)   { }
+#ifdef MFEM_USE_MPI
+   TMOPDeRefinerEstimator(ParMesh &pmesh_, ParNonlinearForm &pnlf_) :
+      mesh(&pmesh_), nlf(&pnlf_), pmesh(&pmesh_), pnlf(&pnlf_),
+      current_sequence(-1), error_estimates(), serial(false) { }
+#endif
+
+   ~TMOPDeRefinerEstimator() { }
+
+   virtual const Vector &GetLocalErrors()
+   {
+      if (MeshIsModified()) { ComputeEstimates(); }
+      return error_estimates;
+   }
+
+   /// Reset the error estimator.
+   virtual void Reset() { current_sequence = -1; }
+};
+
+// hr-adaptivity using TMOP.
+// If hr-adaptivity is disabled, r-adaptivity is done once using the
+// TMOPNewtonSolver.
+// Otherwise, "hr_iter" iterations of r-adaptivity are done followed by
+// "h_per_r_iter" iterations of h-adaptivity after each r-adaptivity iteration.
+// The solver terminates early if an h-adaptivity iteration does not
+// refine/derefine any element in the mesh.
+class TMOPHRSolver
+{
+protected:
+   Mesh *mesh;
+   NonlinearForm *nlf;
+   TMOPNewtonSolver *tmopns;
+   GridFunction *x;
+   Array<GridFunction *> gridfuncarr;
+   Array<FiniteElementSpace *> fespacearr;
+   bool move_bnd, hradaptivity;
+   const int mesh_poly_deg, amr_metric_id;
+#ifdef MFEM_USE_MPI
+   ParMesh *pmesh;
+   ParNonlinearForm *pnlf;
+   Array<ParGridFunction *> pgridfuncarr;
+   Array<ParFiniteElementSpace *> pfespacearr;
+#endif
+   bool serial;
+
+   // All are owned.
+   TMOPRefinerEstimator *tmop_r_est;
+   ThresholdRefiner *tmop_r;
+   TMOPDeRefinerEstimator *tmop_dr_est;
+   ThresholdDerefiner *tmop_dr;
+
+   int hr_iter, h_per_r_iter;
+
+   void Update();
+#ifdef MFEM_USE_MPI
+   void ParUpdate();
+#endif
+   void UpdateNonlinearFormAndBC(Mesh *mesh, NonlinearForm *nlf);
+
+#ifdef MFEM_USE_MPI
+   // Rebalance ParMesh such that all the children elements are moved to the same
+   // MPI rank where the parent will be if the mesh were to be derefined.
+   void RebalanceParNCMesh();
+#endif
+
+public:
+   TMOPHRSolver(Mesh &mesh_, NonlinearForm &nlf_,
+                TMOPNewtonSolver &tmopns_, GridFunction &x_,
+                bool move_bnd_, bool hradaptivity_,
+                int mesh_poly_deg_, int amr_metric_id_,
+                int hr_iter_ = 5, int h_per_r_iter_ = 1);
+#ifdef MFEM_USE_MPI
+   TMOPHRSolver(ParMesh &pmesh_, ParNonlinearForm &pnlf_,
+                TMOPNewtonSolver &tmopns_, ParGridFunction &x_,
+                bool move_bnd_, bool hradaptivity_,
+                int mesh_poly_deg_, int amr_metric_id_,
+                int hr_iter_ = 5, int h_per_r_iter_ = 1);
+#endif
+
+   void Mult();
+
+   /// These are used to update spaces and functions that are not owned by the
+   /// TMOPIntegrator or DiscreteAdaptTC. The owned ones are updated in the
+   /// functions UpdateAfterMeshTopologyChange() of both classes.
+   void AddGridFunctionForUpdate(GridFunction *gf) { gridfuncarr.Append(gf); }
+   void AddFESpaceForUpdate(FiniteElementSpace *fes) { fespacearr.Append(fes); }
+
+#ifdef MFEM_USE_MPI
+   void AddGridFunctionForUpdate(ParGridFunction *pgf_)
+   {
+      pgridfuncarr.Append(pgf_);
+   }
+   void AddFESpaceForUpdate(ParFiniteElementSpace *pfes_)
+   {
+      pfespacearr.Append(pfes_);
+   }
+#endif
+
+   ~TMOPHRSolver()
+   {
+      if (!hradaptivity) { return; }
+      delete tmop_dr;
+      delete tmop_dr_est;
+      delete tmop_r;
+      delete tmop_r_est;
+   }
+
+   /// Total number of hr-adaptivity iterations. At each iteration, we do an
+   /// r-adaptivity iteration followed by a number of h-adaptivity iterations.
+   void SetHRAdaptivityIterations(int iter) { hr_iter = iter; }
+
+   /// Total number of h-adaptivity iterations per r-adaptivity iteration.
+   void SetHAdaptivityIterations(int iter) { h_per_r_iter = iter; }
+};
+
+}
+#endif

fem/tmop_tools.cpp

@@ -407,6 +407,8 @@ double TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
    {
       // Needed for the line search below. The untangling metrics see this
       // reference to detect deteriorations.
+      MFEM_VERIFY(min_det_ptr != NULL, " Initial mesh was valid, but"
+                  " intermediate mesh is invalid. Contact TMOP Developers.");
       *min_det_ptr = untangle_factor * min_detT_in;
    }
 
@@ -576,7 +578,7 @@ void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
          ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
          if (ti)
          {
-            ti->UpdateAfterMeshChange(x_loc);
+            ti->UpdateAfterMeshPositionChange(x_loc);
             ti->ComputeFDh(x_loc, *pfesc);
             UpdateDiscreteTC(*ti, x_loc);
          }
@@ -586,7 +588,7 @@ void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
             Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
             for (int j = 0; j < ati.Size(); j++)
             {
-               ati[j]->UpdateAfterMeshChange(x_loc);
+               ati[j]->UpdateAfterMeshPositionChange(x_loc);
                ati[j]->ComputeFDh(x_loc, *pfesc);
                UpdateDiscreteTC(*ati[j], x_loc);
             }
@@ -613,7 +615,7 @@ void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
          ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
          if (ti)
          {
-            ti->UpdateAfterMeshChange(x_loc);
+            ti->UpdateAfterMeshPositionChange(x_loc);
             ti->ComputeFDh(x_loc, *fesc);
             UpdateDiscreteTC(*ti, x_loc);
          }
@@ -623,7 +625,7 @@ void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
             Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
             for (int j = 0; j < ati.Size(); j++)
             {
-               ati[j]->UpdateAfterMeshChange(x_loc);
+               ati[j]->UpdateAfterMeshPositionChange(x_loc);
                ati[j]->ComputeFDh(x_loc, *fesc);
                UpdateDiscreteTC(*ati[j], x_loc);
             }

general/mem_manager.hpp

@@ -494,8 +494,7 @@ public:
    /// Copy @a size entries from @a *this to @a dest.
    /** The given @a size should not exceed the Capacity() of @a *this and the
        destination, @a dest. */
-   inline void CopyTo(Memory &dest, int size) const
-   { dest.CopyFrom(*this, size); }
+   inline void CopyTo(Memory &dest, int size) const;
 
    /// Copy @a size entries from @a *this to the host pointer @a dest.
    /** The given @a size should not exceed the Capacity() of @a *this. */
@@ -923,6 +922,11 @@ inline void Memory<T>::Wrap(T *ptr, T *d_ptr, int size, MemoryType mt, bool own)
 template <typename T>
 inline void Memory<T>::MakeAlias(const Memory &base, int offset, int size)
 {
+   MFEM_ASSERT(0 <= offset, "invalid offset = " << offset);
+   MFEM_ASSERT(0 <= size, "invalid size = " << size);
+   MFEM_ASSERT(offset + size <= base.capacity,
+               "invalid offset + size = " << offset + size
+               << " > base capacity = " << base.capacity);
    capacity = size;
    h_mt = base.h_mt;
    h_ptr = base.h_ptr + offset;
@@ -1136,6 +1140,7 @@ inline bool Memory<T>::DeviceIsValid() const
 template <typename T>
 inline void Memory<T>::CopyFrom(const Memory &src, int size)
 {
+   MFEM_VERIFY(src.capacity>=size && capacity>=size, "Incorrect size");
    if (!(flags & REGISTERED) && !(src.flags & REGISTERED))
    {
       if (h_ptr != src.h_ptr && size != 0)
@@ -1155,6 +1160,7 @@ inline void Memory<T>::CopyFrom(const Memory &src, int size)
 template <typename T>
 inline void Memory<T>::CopyFromHost(const T *src, int size)
 {
+   MFEM_VERIFY(capacity>=size, "Incorrect size");
    if (!(flags & REGISTERED))
    {
       if (h_ptr != src && size != 0)
@@ -1171,9 +1177,17 @@ inline void Memory<T>::CopyFromHost(const T *src, int size)
    }
 }
 
+template <typename T>
+inline void Memory<T>::CopyTo(Memory &dest, int size) const
+{
+   MFEM_VERIFY(capacity>=size, "Incorrect size");
+   dest.CopyFrom(*this, size);
+}
+
 template <typename T>
 inline void Memory<T>::CopyToHost(T *dest, int size) const
 {
+   MFEM_VERIFY(capacity>=size, "Incorrect size");
    if (!(flags & REGISTERED))
    {
       if (h_ptr != dest && size != 0)

general/version.cpp

@@ -82,9 +82,6 @@ const char *GetConfigStr()
 #ifdef MFEM_USE_CUDA
       "MFEM_USE_CUDA\n"
 #endif
-#ifdef MFEM_USE_EPIC
-      "MFEM_USE_EPIC\n"
-#endif
 #ifdef MFEM_USE_EXCEPTIONS
       "MFEM_USE_EXCEPTIONS\n"
 #endif

linalg/CMakeLists.txt

@@ -83,11 +83,6 @@ if (MFEM_USE_SUNDIALS)
   list(APPEND HDRS sundials.hpp)
 endif()
 
-if (MFEM_USE_EPIC)
-  list(APPEND SRCS epic.cpp)
-  list(APPEND HDRS epic.hpp)
-endif()
-
 if (MFEM_USE_SUPERLU)
   list(APPEND SRCS superlu.cpp)
   # If this list (HDRS -> HEADERS) is used for install, we probably want the

linalg/epic.cpp

@@ -1,171 +0,0 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "epic.hpp"
-
-#ifdef MFEM_USE_EPIC
-
-namespace mfem
-{
-
-EPICSolver::EPICSolver(bool exactJacobian_, EPICNumJacDelta delta)
-{
-   // Allocate an empty serial N_Vector
-   temp = N_VNewEmpty_Serial(0);
-   m[0] = 10;
-   m[1] = 10;
-   MFEM_VERIFY(temp, "error in N_VNewEmpty_Serial()");
-   exactJacobian = exactJacobian_;
-   Jtv = NULL;
-   Delta = delta;
-}
-
-#ifdef MFEM_USE_MPI
-EPICSolver::EPICSolver(MPI_Comm comm)
-{
-   m[0] = 10;
-   m[1] = 10;
-
-   // Allocate an empty vector
-   if (comm == MPI_COMM_NULL)
-   {
-      // Allocate an empty serial N_Vector
-      temp = N_VNewEmpty_Serial(0);
-      MFEM_VERIFY(temp, "error in N_VNewEmpty_Serial()");
-   }
-   else
-   {
-      // Allocate an empty parallel N_Vector
-      temp = N_VNewEmpty_Parallel(comm, 0, 0);  // calls MPI_Allreduce()
-      MFEM_VERIFY(temp, "error in N_VNewEmpty_Parallel()");
-   }
-}
-#endif
-
-int EPICSolver::RHS(realtype t, const N_Vector y, N_Vector ydot, void *user_data)
-{
-   // Get data from N_Vectors
-   const Vector mfem_y(y);
-   Vector mfem_ydot(ydot);
-   EPICSolver *self = static_cast<EPICSolver*>(user_data);
-
-   // Compute y' = f(t, y)
-   self->f->SetTime(t);
-   self->f->Mult(mfem_y, mfem_ydot);
-
-   // Return success
-   return 0;
-}
-
-
-int EPICSolver::Jacobian(N_Vector v, N_Vector Jv, realtype t, N_Vector y, N_Vector fy, void *user_data, N_Vector tmp)
-{
-    // Get data from N_Vectors
-   const Vector mfem_v(v);
-   Vector mfem_Jv(Jv);
-   EPICSolver *self = static_cast<EPICSolver*>(user_data);
-
-   // Compute J(t, y) v
-   self->Jtv->Mult(mfem_v, mfem_Jv);
-
-   return 0;
-}
-
-void EPICSolver::Init(TimeDependentOperator &f)
-{
-    ODESolver::Init(f);
-	    
-    long local_size = f.Height();
-    long global_size = 0;
-#ifdef MFEM_USE_MPI
-    if (Parallel())
-    {
-        MPI_Allreduce(&local_size, &global_size, 1, MPI_LONG, MPI_SUM,
-                      NV_COMM_P(temp));
-    }
-#endif
-
-    Vector mfem_temp(local_size);
-    mfem_temp.ToNVector(temp, global_size);
-}
-
-EPI2::EPI2(bool exactJacobian, EPICNumJacDelta delta) : EPICSolver(exactJacobian, delta) {}
-
-void EPI2::Init(TimeDependentOperator &f)
-{
-    EPICSolver::Init(f);
-    long local_size = f.Height();
-    if (exactJacobian) {
-       integrator = new Epi2_KIOPS(EPICSolver::RHS, EPICSolver::Jacobian, this, 100, temp ,local_size);
-    } else {
-       integrator = new Epi2_KIOPS(EPICSolver::RHS, Delta, this, 100, temp ,local_size);
-    }
-
-}
-
-EPIRK4::EPIRK4(bool exactJacobian, EPICNumJacDelta delta) : EPICSolver(exactJacobian, delta) {}
-
-void EPIRK4::Init(TimeDependentOperator &f)
-{
-    EPICSolver::Init(f);
-    long local_size = f.Height();
-    if (exactJacobian) {
-       integrator = new EpiRK4SC_KIOPS(EPICSolver::RHS, EPICSolver::Jacobian, this, 100, temp ,local_size);
-    } else {
-       integrator = new EpiRK4SC_KIOPS(EPICSolver::RHS, Delta, this, 100, temp ,local_size);
-    }
-}
-
-void EPICSolver::Step(Vector &x, double &t, double &dt)
-{
-   if (!Parallel())
-   {
-      NV_DATA_S(temp) = x.GetData();
-      MFEM_VERIFY(NV_LENGTH_S(temp) == x.Size(), "");
-   }
-   else
-   {
-#ifdef MFEM_USE_MPI
-      NV_DATA_P(temp) = x.GetData();
-      MFEM_VERIFY(NV_LOCLENGTH_P(temp) == x.Size(), "");
-#endif
-   }
-
-   Jtv = &(this->f->GetGradient(x));
-}
-
-void EPI2::Step(Vector &x, double &t, double &dt)
-{
-    EPICSolver::Step(x, t, dt);
-    integrator->Integrate(dt, t, t+dt, 0, temp, 1e-10, m);
-    t += dt;
-}
-
-void EPIRK4::Step(Vector &x, double &t, double &dt)
-{
-    EPICSolver::Step(x, t, dt);
-    integrator->Integrate(dt, t, t+dt, 0, temp, 1e-10, m);
-    t += dt;
-}
-
-EPI2::~EPI2()
-{
-    delete integrator;
-}
-
-EPIRK4::~EPIRK4()
-{
-    delete integrator;
-}
-
-}
-
-#endif

linalg/epic.hpp

@@ -1,98 +0,0 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MFEM_EPIC
-#define MFEM_EPIC
-
-#include "../config/config.hpp"
-
-#ifdef MFEM_USE_EPIC
-
-// SUNDIALS vectors
-#include <nvector/nvector_serial.h>
-#ifdef MFEM_USE_MPI
-#include <mpi.h>
-#include <nvector/nvector_parallel.h>
-#endif
-
-#include "ode.hpp"
-#include "solvers.hpp"
-#include <Epic.h>
-
-namespace mfem
-{
-
-typedef void (*JacobianFun)(const realtype t, const Vector &y, const Vector& v, Vector& Jv, void* user_data);
-// ---------------------------------------------------------------------------
-// Interface to the EPIC library -- exponential methods
-// ---------------------------------------------------------------------------
-
-class EPICSolver : public ODESolver
-{
-protected:
-    EPICNumJacDelta Delta;
-    Operator* Jtv;
-    N_Vector temp;
-    int m[2];
-
-    bool exactJacobian;
-
-    #ifdef MFEM_USE_MPI
-    bool Parallel() const
-    {
-        return (N_VGetVectorID(temp) != SUNDIALS_NVEC_SERIAL);
-    }
-    #else
-    bool Parallel() const { return false; }
-    #endif
-
-public:
-    EPICSolver(bool exactJacobian, EPICNumJacDelta delta=&DefaultDelta);
-    EPICSolver(MPI_Comm comm);
-
-    static int RHS(realtype t, const N_Vector y, N_Vector ydot, void *user_data);
-    static int Jacobian(N_Vector v, N_Vector Jv, realtype t,
-                         N_Vector y, N_Vector fy, void *user_data, N_Vector tmp);
-    virtual void Init(TimeDependentOperator &f);
-    virtual void Step(Vector &x, double &t, double &dt);
-
-    virtual ~EPICSolver() {}
-};
-
-class EPI2 : public EPICSolver
-{
-protected:
-    Epi2_KIOPS* integrator;
-public:
-    EPI2(bool exactJacobian=true, EPICNumJacDelta delta=&DefaultDelta);
-    virtual void Init(TimeDependentOperator &f);
-    virtual void Step(Vector &x, double &t, double &dt);
-
-    virtual ~EPI2();
-};
-
-class EPIRK4 : public EPICSolver
-{
-protected:
-    EpiRK4SC_KIOPS* integrator;
-public:
-    EPIRK4(bool exactJacobian=true, EPICNumJacDelta delta=&DefaultDelta);
-    virtual void Init(TimeDependentOperator &f);
-    virtual void Step(Vector &x, double &t, double &dt);
-
-    virtual ~EPIRK4();
-};
-
-} // namespace mfem
-
-#endif // MFEM_USE_EPIC
-
-#endif // MFEM_EPIC

linalg/imex.hpp

@@ -0,0 +1,244 @@
+#ifndef MFEM_IMEX
+#define MFEM_IMEX
+
+#include "../config/config.hpp"
+#include "operator.hpp"
+#include "ode.hpp"
+#include <vector>
+
+namespace mfem
+{
+
+/** Class for spatial discretizations of a PDE resulting in the time-dependent, 
+    nonlinear set of ODEs with implicit-explicit additive partition
+        M*du/dt = N_E(u,t) + N_I(u,t).
+    MFEM typically treats time integration as 
+        du/dt = F^{-1} G(u),
+    Here F represents what MFEM calls the “implicit” part, and G represents the
+    “explicit” part; in simpler terms, F is typically just a mass matrix. 
+
+    For BDF schemes, the ImplicitSolve function is a bit different, and it is
+    more natural to apply M and M^{-1} separate from the Mult functions, so we
+    include MassMult and MassInv as functions to be provided, and do not include
+    such actions in the Mult functions. */
+class IMEXTimeDependentOperator : public TimeDependentOperator
+{
+protected:
+    mutable Vector temp;    // Auxillary vector
+    
+public:
+    // Sets linearly implicit to false by default
+    IMEXTimeDependentOperator(int n, double t=0.0, Type type=EXPLICIT)
+        : TimeDependentOperator(n, t, type) { };
+    ~IMEXTimeDependentOperator() { };
+
+    /** Apply action of implicit part of operator y <- N_I(x,y). For fully
+        implicit schemes, this just corresponds to applying the time-dependent
+        (nonlinear) operator.
+        PREVIOUSLY CALLED ExplicitMult */
+    virtual void ImplicitMult(const Vector &x, Vector &y) const = 0;
+
+    /** Apply action of explicit part of operator y <- N_E(x,y) */
+    virtual void ExplicitMult(const Vector &x, Vector &y) const { y = 0.0; };
+
+    /** Solve k = f(x+dt*k) for stage k, where f() is the implicit part of
+        the operator. Used in Runge-Kutta methods. */
+    virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
+    { mfem::mfem_error("IMEXTimeDependentOperator::ImplicitSolve() is not overridden!"); };
+
+    /** Solve M*x - dtf(x, t) = b for solution x, where f() is the implicit
+        part of the operator. Used in BDF methods. */
+    virtual void ImplicitSolve2(const double dt, const Vector &b, Vector &x)
+    { mfem::mfem_error("IMEXTimeDependentOperator::ImplicitSolve2() is not overridden!"); };
+
+    /** Apply action mass matrix, y = M*x. 
+        If not re-implemented, this method simply generates an error.
+        PREVIOUSLY CALLED ImplictMult */
+    virtual void MassMult(const Vector &x, Vector &y) const = 0;
+
+    /** Apply action of inverse of mass matrix, y = M^{-1}*x. 
+        If not re-implemented, this method simply generates an error.
+        NOTE : only necessary for PolyIMEX methods. */
+    virtual void MassInv(const Vector &x, Vector &y) const = 0;
+};
+
+/** Class holding RK Butcher tableau, and associated data required by
+    implicit and explicit splitting. */
+class IMEXRKData 
+{
+public:
+    // Implicit Runge Kutta type. Enumeration (s, \sigma, p):
+    // - s = number of implicit stages
+    // - \sigma = number of explicit stages
+    // - p = order
+    // In this notation, when s = \sigma, we satisfy (2.3)/(2.4) in
+    // Ascher et al., and do not need to compute the final explicit
+    // stage. This is represented in the stiffly_accurate boolean.
+    enum Type { 
+        IMEX111 = 111,
+        IMEX121 = 121,
+        IMEX122 = 122,
+        IMEX222 = 222,
+        IMEX232 = 232,
+        IMEX233 = 233,
+        IMEX443 = 443,
+    // ARK ESDIRK-ERK schemes: enumeration (s,p), for total number of
+    // stages s.
+        ARK43 = -43
+    };
+
+    IMEXRKData() : s(-1) { };
+    IMEXRKData(Type ID_) : ID(ID_) { SetData(); };
+    ~IMEXRKData() { };
+    
+    /// Set explicit RK data
+    void SetExplicitData(DenseMatrix Ae_, Vector be_, Vector ce_);
+    /// Set implicit RK data
+    void SetImplicitData(DenseMatrix Ai_, Vector bi_, Vector ci_, bool esdirk_=false);
+    void SetID(Type ID_) { ID=ID_; SetData(); };
+
+    bool esdirk;
+    bool stiffly_accurate;
+    bool use_final_exp_stage;
+    int s;
+
+    DenseMatrix Ai;     // Implicit Butcher matrix
+    Vector bi;          // Implicit Butcher tableau weights
+    DenseMatrix Ae;     // Explicit Butcher matrix
+    Vector be;          // Explicit Butcher tableau weights
+    Vector c0;          // Butcher tableau nodes (same for implicit and explicit!)
+
+private:    
+    Type ID;
+    void SetData();
+    void InitData();
+};
+
+/** Class for two-part additive IMEX RK method, where explicit and implicit
+    stage vectors are stored. Assume same abscissae, {c}, for both schemes.
+    Butcher Data must be provided either in a custom IMEXRKData object, or
+    using the IMEXRKData::Type for predefined tableaux. */
+class IMEXRK : public ODESolver
+{
+protected:
+    IMEXRKData tableaux;
+    std::vector< Vector *> exp_stages;
+    std::vector< Vector *> imp_stages;
+    IMEXTimeDependentOperator *imex;    // Spatial discretization. 
+
+public:
+    IMEXRK(IMEXRKData tableaux_) : ODESolver(), tableaux(tableaux_) { };
+    IMEXRK(IMEXRKData::Type type_) : ODESolver(), tableaux(type_) { };
+    ~IMEXRK();
+    void Init(IMEXTimeDependentOperator &_imex);
+    void Step(Vector &x, double &t, double &dt) override;
+};
+
+/** Class holding BDF integrator data. Setting alpha < 0 (the default
+    constructor) defines alpha = 2/(q-1), corresponding to classical BDF
+    of order q. */
+class BDFData 
+{
+public:
+
+    enum Type {
+        BDF1 = 01, BDF2 = 02, BDF3 = 03, BDF4 = 04,
+        IMEX_BDF1 = 11, IMEX_BDF2 = 12, IMEX_BDF3 = 13,
+        IMEX_BDF4 = 14
+    };
+
+    BDFData() { };
+    BDFData(Type ID_, double alpha_=-1) : ID(ID_), alpha(alpha_) { SetData(); };
+    ~BDFData() { };
+
+    int GetID() { return static_cast<int>(ID); };
+    void SetID(Type ID_, double alpha_=-1) 
+    {
+        ID=ID_;
+        alpha = alpha_;
+        SetData();
+    };
+    void Print()
+    {
+        std::cout << "q     = " << q << "\n";
+        std::cout << "alpha = " << alpha << "\n";
+        std::cout << "A:\n";
+        A.PrintMatlab();
+        std::cout << "Be:\n";
+        Be.PrintMatlab();
+        std::cout << "Bi:\n";
+        Bi.Print();
+        std::cout << "z:\n";
+        z0.Print();
+    };
+
+    double alpha;
+    int q;              // Number of previous values stored
+    bool shifted_nodes; // false = clssical BDF, true = Polynomial BDF w/ shifted nodes
+    DenseMatrix A;      // Previous solution coefficients
+    Vector Bi;          // Implicit coefficients
+    DenseMatrix Be;     // Explicit coefficients
+    Vector z0;
+
+
+private:    
+    Type ID;
+    void SetData();
+    void InitData();
+};
+
+/** Class for IMEX-BDF methods, including classical IMEX-BDF and IMEX-
+    Polynomial-BDF (IMEX-PBDF). IMEX-PBDF methods have an additional 
+    alpha parameter, where larger alpha leads to smaller stability
+    regions and a smaller leading accuracy constant, while smaller
+    alpha leads to larger stabiltiy regions and a larger accuracy
+    constant. For classical methods, there are two implementations:
+        - ClassicalStep() stores previous solutions and the explicit
+        part of the operator evaluated on the solution, and
+        - ClassicalStepNoStore() does not store the explicit 
+        component, but must re-evaluate q times during each time
+        step.
+    This option can be set via the recompute_exp input. The type of
+    scheme must be set through the BDFData structure or BDFData::Type.
+        There is also an option to use pointwise Lagrange interpolating 
+    polynomials to provide an initial guess for the ImplicitSolve. This
+    is set via InterpolateGuess(). This option is only implemented for
+    PBDF. */
+class IMEXBDF : public ODESolver
+{
+private:
+    BDFData data;
+    bool recompute_exp;
+    bool interpolate;
+    int initialized;
+    double dt_prev;
+    std::vector< Vector*> sols;
+    std::vector< Vector*> exp_sols;
+    IMEXTimeDependentOperator *imex;    // Spatial discretization
+    IMEXRK *RKsolver;
+    std::vector<double> exp_nodes;
+
+    void AlphaStep(Vector &x, double &t, double &dt);
+    void ClassicalStep(Vector &x, double &t, double &dt);
+    void ClassicalStepNoStore(Vector &x, double &t, double &dt);
+
+public:
+    IMEXBDF(BDFData data_, bool recompute_exp_=false) :
+      ODESolver(), data(data_), recompute_exp(recompute_exp_),
+      interpolate(false) { };
+    IMEXBDF(BDFData::Type scheme, bool recompute_exp_=false) :
+      ODESolver(), recompute_exp(recompute_exp_), interpolate(false)
+      { data.SetID(scheme); };
+    IMEXBDF(BDFData::Type scheme, double alpha) :
+      ODESolver(), interpolate(false), recompute_exp(false)
+      { data.SetID(scheme, alpha); };
+    ~IMEXBDF();
+
+    void Init(IMEXTimeDependentOperator &_imex);
+    void Step(Vector &x, double &t, double &dt);
+    void InterpolateGuess() {interpolate = true; };
+};
+
+}
+
+#endif

linalg/linalg.hpp

@@ -31,6 +31,7 @@
 #include "invariants.hpp"
 #include "constraints.hpp"
 #include "auxiliary.hpp"
+#include "imex.hpp"
 
 #ifdef MFEM_USE_AMGX
 #include "amgxsolver.hpp"
@@ -40,10 +41,6 @@
 #include "sundials.hpp"
 #endif
 
-#ifdef MFEM_USE_EPIC
-#include "epic.hpp"
-#endif
-
 #ifdef MFEM_USE_HIOP
 #include "hiop.hpp"
 #endif

linalg/sparsemat.cpp

@@ -33,7 +33,12 @@ int SparseMatrix::SparseMatrixCount = 0;
 cusparseHandle_t SparseMatrix::handle = nullptr;
 size_t SparseMatrix::bufferSize = 0;
 void * SparseMatrix::dBuffer = nullptr;
-#endif
+#  if CUSPARSE_VERSION >=  11400
+#    define MFEM_CUSPARSE_ALG CUSPARSE_SPMV_CSR_ALG1
+#  else
+#    define MFEM_CUSPARSE_ALG CUSPARSE_CSRMV_ALG1
+#  endif // CUSPARSE_VERSION >= 11400
+#endif // MFEM_USE_CUDA
 
 void SparseMatrix::InitCuSparse()
 {
@@ -679,25 +684,16 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const double a) const
          cusparseCreateMatDescr(&matA_descr);
          cusparseSetMatIndexBase(matA_descr, CUSPARSE_INDEX_BASE_ZERO);
          cusparseSetMatType(matA_descr, CUSPARSE_MATRIX_TYPE_GENERAL);
-
 #endif
-
          initBuffers = true;
       }
       // Allocate kernel space. Buffer is shared between different sparsemats
       size_t newBufferSize = 0;
 
-#if CUDA_VERSION >= 11020
       cusparseSpMV_bufferSize(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha,
                               matA_descr,
                               vecX_descr, &beta, vecY_descr, CUDA_R_64F,
-                              CUSPARSE_SPMV_CSR_ALG1, &newBufferSize);
-#elif CUDA_VERSION >= 10010
-      cusparseSpMV_bufferSize(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha,
-                              matA_descr,
-                              vecX_descr, &beta, vecY_descr, CUDA_R_64F,
-                              CUSPARSE_CSRMV_ALG1, &newBufferSize);
-#endif
+                              MFEM_CUSPARSE_ALG, &newBufferSize);
 
       // Check if we need to resize
       if (newBufferSize > bufferSize)
@@ -707,30 +703,22 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const double a) const
          CuMemAlloc(&dBuffer, bufferSize);
       }
 
-#if CUDA_VERSION >= 11020
+#if CUDA_VERSION >= 10010
       // Update input/output vectors
       cusparseDnVecSetValues(vecX_descr, const_cast<double *>(d_x));
       cusparseDnVecSetValues(vecY_descr, d_y);
 
       // Y = alpha A * X + beta * Y
       cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA_descr,
-                   vecX_descr, &beta, vecY_descr, CUDA_R_64F, CUSPARSE_SPMV_CSR_ALG1, dBuffer);
-#elif CUDA_VERSION >= 10010
-      // Update input/output vectors
-      cusparseDnVecSetValues(vecX_descr, const_cast<double *>(d_x));
-      cusparseDnVecSetValues(vecY_descr, d_y);
-
-      // Y = alpha A * X + beta * Y
-      cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA_descr,
-                   vecX_descr, &beta, vecY_descr, CUDA_R_64F, CUSPARSE_CSRMV_ALG1, dBuffer);
+                   vecX_descr, &beta, vecY_descr, CUDA_R_64F, MFEM_CUSPARSE_ALG, dBuffer);
 #else
       cusparseDcsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
                      Height(), Width(), J.Capacity(),
                      &alpha, matA_descr,
                      const_cast<double *>(d_A), const_cast<int *>(d_I), const_cast<int *>(d_J),
                      const_cast<double *>(d_x), &beta, d_y);
-#endif
-#endif
+#endif // CUDA_VERSION >= 10010
+#endif // MFEM_USE_CUDA
    }
    else
    {

linalg/sundials.cpp

@@ -425,26 +425,6 @@ static SUNMatrix_ID MatGetID(SUNMatrix)
    return (SUNMATRIX_CUSTOM);
 }
 
-//BEGIN WORKAROUND CODE
-static SUNMatrix MatClone(SUNMatrix A)
-{
-   SUNMatrix B = SUNMatNewEmpty();
-   MFEM_VERIFY(B, "error in MatClone()");
-
-   B->content = A->content;
-
-   int flag = SUNMatCopyOps(A, B);
-   MFEM_VERIFY(flag == SUNMAT_SUCCESS, "error in MatClone()");
-
-   return B;
-}
-
-static int MatCopy(SUNMatrix, SUNMatrix)
-{
-   return 0;
-}
-// END WORKAROUN CODE
-
 static void MatDestroy(SUNMatrix A)
 {
    if (A->content) { A->content = NULL; }
@@ -1513,12 +1493,7 @@ void ARKStepSolver::UseMFEMMassLinearSolver(int tdep)
    MFEM_VERIFY(M, "error in SUNMatNewEmpty()");
 
    M->content      = this;
-// BEGIN WORKAROUND CODE
-   // M->ops->getid   = SUNMatGetID;
-   M->ops->getid = MatGetID;
-   M->ops->clone = MatClone;
-   M->ops->copy = MatCopy;
-// END WORKAROUND CODE
+   M->ops->getid   = SUNMatGetID;
    M->ops->matvec  = ARKStepSolver::MassMult1;
    M->ops->destroy = MatDestroy;
 
@@ -1562,12 +1537,6 @@ void ARKStepSolver::SetSStolerances(double reltol, double abstol)
    MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSStolerances()");
 }
 
-void ARKStepSolver::SetMaxIter(int iterations)
-{
-  flag = ARKStepSetMaxNonlinIters(sundials_mem, iterations);
-  MFEM_VERIFY(flag == ARK_SUCCESS, "error in ARKStepSetMaxNonlinIters()");
-}
-
 void ARKStepSolver::SetMaxStep(double dt_max)
 {
    flag = ARKStepSetMaxStep(sundials_mem, dt_max);

linalg/sundials.hpp

@@ -666,9 +666,6 @@ public:
    /// Set the scalar relative and scalar absolute tolerances.
    void SetSStolerances(double reltol, double abstol);
 
-   /// Set the maximum number of nonlinear iterations per RK stage
-   void SetMaxIter(int iterations);
-
    /// Set the maximum time step.
    void SetMaxStep(double dt_max);
 

linalg/vector.hpp

@@ -82,6 +82,11 @@ public:
    Vector(double *data_, int size_)
    { data.Wrap(data_, size_, false); size = size_; }
 
+   /** @brief Create a Vector referencing a sub-vector of the Vector @a base
+       starting at the given offset, @a base_offset, and size @a size_. */
+   Vector(Vector &base, int base_offset, int size_)
+      : data(base.data, base_offset, size_), size(size_) { }
+
    /// Create a Vector of size @a size_ using MemoryType @a mt.
    Vector(int size_, MemoryType mt)
       : data(size_, mt), size(size_) { }

makefile

@@ -1,4 +1,4 @@
-# Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 # at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 # LICENSE and NOTICE for details. LLNL-CODE-806117.
 #
@@ -119,7 +119,7 @@ $(if $(word 2,$(SRC)),$(error Spaces in SRC = "$(SRC)" are not supported))
 MFEM_GIT_STRING = $(shell [ -d $(MFEM_DIR)/.git ] && git -C $(MFEM_DIR) \
    describe --all --long --abbrev=40 --dirty --always 2> /dev/null)
 
-EXAMPLE_SUBDIRS = amgx ginkgo hiop petsc pumi sundials epic superlu
+EXAMPLE_SUBDIRS = amgx ginkgo hiop petsc pumi sundials superlu
 EXAMPLE_DIRS := examples $(addprefix examples/,$(EXAMPLE_SUBDIRS))
 EXAMPLE_TEST_DIRS := examples
 
@@ -272,7 +272,7 @@ ifeq ($(MFEM_USE_LEGACY_OPENMP),YES)
 endif
 
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
-MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS EPIC MESQUITE\
+MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
  SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP GSLIB\
  OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX CALIPER
 
@@ -333,7 +333,7 @@ endif
 MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_METIS MFEM_USE_METIS_5 MFEM_DEBUG MFEM_USE_EXCEPTIONS MFEM_USE_ZLIB\
  MFEM_USE_LIBUNWIND MFEM_USE_LAPACK MFEM_THREAD_SAFE MFEM_USE_OPENMP\
- MFEM_USE_LEGACY_OPENMP MFEM_USE_MEMALLOC MFEM_TIMER_TYPE MFEM_USE_SUNDIALS MFEM_USE_EPIC\
+ MFEM_USE_LEGACY_OPENMP MFEM_USE_MEMALLOC MFEM_TIMER_TYPE MFEM_USE_SUNDIALS\
  MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE MFEM_USE_GINKGO MFEM_USE_SUPERLU\
  MFEM_USE_STRUMPACK MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC\
  MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT\
@@ -650,7 +650,6 @@ status info:
 	$(info MFEM_USE_MEMALLOC      = $(MFEM_USE_MEMALLOC))
 	$(info MFEM_TIMER_TYPE        = $(MFEM_TIMER_TYPE))
 	$(info MFEM_USE_SUNDIALS      = $(MFEM_USE_SUNDIALS))
-	$(info MFEM_USE_EPIC          = $(MFEM_USE_EPIC))
 	$(info MFEM_USE_MESQUITE      = $(MFEM_USE_MESQUITE))
 	$(info MFEM_USE_SUITESPARSE   = $(MFEM_USE_SUITESPARSE))
 	$(info MFEM_USE_SUPERLU       = $(MFEM_USE_SUPERLU))

mesh/mesh.cpp

@@ -75,7 +75,7 @@ void Mesh::GetElementCenter(int i, Vector &center)
 
 double Mesh::GetElementSize(ElementTransformation *T, int type)
 {
-   DenseMatrix J(Dim);
+   DenseMatrix J(spaceDim,Dim);
 
    Geometry::Type geom = T->GetGeometryType();
    T->SetIntPoint(&Geometries.GetCenter(geom));
@@ -83,7 +83,7 @@ double Mesh::GetElementSize(ElementTransformation *T, int type)
 
    if (type == 0)
    {
-      return pow(fabs(J.Det()), 1./Dim);
+      return pow(fabs(J.Weight()), 1./Dim);
    }
    else if (type == 1)
    {
@@ -102,7 +102,7 @@ double Mesh::GetElementSize(int i, int type)
 
 double Mesh::GetElementSize(int i, const Vector &dir)
 {
-   DenseMatrix J(Dim);
+   DenseMatrix J(spaceDim,Dim);
    Vector d_hat(Dim);
    GetElementJacobian(i, J);
    J.MultTranspose(dir, d_hat);

mesh/mesh.hpp

@@ -503,7 +503,7 @@ public:
    Mesh& operator=(Mesh &&mesh);
 
    /// Explicitly delete the copy assignment operator.
-   Mesh& operator=(Mesh &mesh) = delete;
+   Mesh& operator=(const Mesh &mesh) = delete;
 
    /** @name Named mesh constructors.
 

mesh/ncmesh.cpp

@@ -2283,7 +2283,6 @@ void NCMesh::GetMeshComponents(Mesh &mesh) const
    // left uninitialized here; they will be initialized later by the Mesh from
    // Nodes -- here we just make sure mesh.vertices has the correct size.
 
-   mesh.elements.SetSize(NElements);
    mesh.elements.SetSize(0);
 
    mesh.boundary.SetSize(0);

mesh/ncmesh.hpp

@@ -131,6 +131,7 @@ public:
    int GetNVertices() const { return NVertices; }
    int GetNEdges() const { return NEdges; }
    int GetNFaces() const { return NFaces; }
+   virtual int GetNGhostElements() const { return 0; }
 
    /** Perform the given batch of refinements. Please note that in the presence
        of anisotropic splits additional refinements may be necessary to keep
@@ -160,7 +161,6 @@ public:
        derefinements may have to be skipped to preserve mesh consistency. */
    virtual void Derefine(const Array<int> &derefs);
 
-
    // master/slave lists
 
    /// Identifies a vertex/edge/face in both Mesh and NCMesh.

mesh/pncmesh.cpp

@@ -1339,6 +1339,37 @@ void ParNCMesh::LimitNCLevel(int max_nc_level)
    }
 }
 
+void ParNCMesh::GetFineToCoarsePartitioning(const Array<int> &derefs,
+                                            Array<int> &new_ranks) const
+{
+   new_ranks.SetSize(leaf_elements.Size()-GetNGhostElements());
+   for (int i = 0; i < leaf_elements.Size()-GetNGhostElements(); i++)
+   {
+      new_ranks[i] = elements[leaf_elements[i]].rank;
+   }
+
+   for (int i = 0; i < derefs.Size(); i++)
+   {
+      int row = derefs[i];
+      MFEM_VERIFY(row >= 0 && row < derefinements.Size(),
+                  "invalid derefinement number.");
+
+      const int* fine = derefinements.GetRow(row);
+      int size = derefinements.RowSize(row);
+
+      int coarse_rank = INT_MAX;
+      for (int j = 0; j < size; j++)
+      {
+         int fine_rank = elements[leaf_elements[fine[j]]].rank;
+         coarse_rank = std::min(coarse_rank, fine_rank);
+      }
+      for (int j = 0; j < size; j++)
+      {
+         new_ranks[fine[j]] = coarse_rank;
+      }
+   }
+}
+
 void ParNCMesh::Derefine(const Array<int> &derefs)
 {
    MFEM_VERIFY(Dim < 3 || Iso,

mesh/pncmesh.hpp

@@ -95,6 +95,11 @@ public:
        in sync. The interface is identical. */
    virtual void Derefine(const Array<int> &derefs);
 
+   /** Gets partitioning for the coarse mesh if the current fine mesh were to
+       be derefined. */
+   virtual void GetFineToCoarsePartitioning(const Array<int> &derefs,
+                                            Array<int> &new_ranks) const;
+
    /** Migrate leaf elements of the global refinement hierarchy (including ghost
        elements) so that each processor owns the same number of leaves (+-1).
        The default partitioning strategy is based on equal splitting of the
@@ -243,7 +248,6 @@ public:
        The debug mesh will have element attributes set to element rank + 1. */
    void GetDebugMesh(Mesh &debug_mesh) const;
 
-
 protected: // interface for ParMesh
 
    friend class ParMesh;

miniapps/gslib/pfindpts.cpp

@@ -35,7 +35,7 @@
 //    mpirun -np 2 pfindpts -m ../../data/inline-wedge.mesh -o 3
 //    mpirun -np 2 pfindpts -m ../../data/amr-quad.mesh -o 2
 //    mpirun -np 2 pfindpts -m ../../data/rt-2d-q3.mesh -o 3 -mo 4 -ft 2
-
+//    mpirun -np 2 pfindpts -m ../../data/inline-quad.mesh -ft 1 -no-vis -sr0
 
 #include "mfem.hpp"
 
@@ -74,6 +74,7 @@ int main (int argc, char *argv[])
    bool visualization    = true;
    int fieldtype         = 0;
    int ncomp             = 1;
+   bool search_on_rank_0 = false;
 
    // Parse command-line options.
    OptionsParser args(argc, argv);
@@ -94,6 +95,9 @@ int main (int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+   args.AddOption(&search_on_rank_0, "-sr0", "--search-on-r0", "-no-sr0",
+                  "--no-search-on-r0",
+                  "Enable search only on rank 0 (disable to search points on all tasks).");
    args.Parse();
    if (!args.Good())
    {
@@ -211,7 +215,7 @@ int main (int argc, char *argv[])
    // Note that some points might be outside, if the mesh is not a box. Note
    // also that all tasks search the same points (not mandatory).
    const int pts_cnt_1D = 10;
-   const int pts_cnt = pow(pts_cnt_1D, dim);
+   int pts_cnt = pow(pts_cnt_1D, dim);
    Vector vxyz(pts_cnt * dim);
    if (dim == 2)
    {
@@ -237,6 +241,12 @@ int main (int argc, char *argv[])
       }
    }
 
+   if ( (myid != 0) && (search_on_rank_0) )
+   {
+      pts_cnt = 0;
+      vxyz.Destroy();
+   }
+
    // Find and Interpolate FE function values on the desired points.
    Vector interp_vals(pts_cnt*vec_dim);
    FindPointsGSLIB finder(MPI_COMM_WORLD);
@@ -246,36 +256,37 @@ int main (int argc, char *argv[])
    Array<unsigned int> task_id_out = finder.GetProc();
    Vector dist_p_out = finder.GetDist();
 
-   int face_pts = 0, not_found = 0, found_loc = 0, found_away = 0;
-   double max_err = 0.0, max_dist = 0.0;
-   Vector pos(dim);
-   int npt = 0;
-   for (int j = 0; j < vec_dim; j++)
+   // Print the results for task 0 since either 1) all tasks have the
+   // same set of points or 2) only task 0 has any points.
+   if (myid == 0 )
    {
-      for (int i = 0; i < pts_cnt; i++)
+      int face_pts = 0, not_found = 0, found_loc = 0, found_away = 0;
+      double max_err = 0.0, max_dist = 0.0;
+      Vector pos(dim);
+      int npt = 0;
+      for (int j = 0; j < vec_dim; j++)
       {
-         if (j == 0)
+         for (int i = 0; i < pts_cnt; i++)
          {
-            (task_id_out[i] == (unsigned)myid) ? found_loc++ : found_away++;
-         }
+            if (j == 0)
+            {
+               (task_id_out[i] == (unsigned)myid) ? found_loc++ : found_away++;
+            }
 
-         if (code_out[i] < 2)
-         {
-            for (int d = 0; d < dim; d++) { pos(d) = vxyz(d * pts_cnt + i); }
-            Vector exact_val(vec_dim);
-            F_exact(pos, exact_val);
-            max_err  = std::max(max_err, fabs(exact_val(j) - interp_vals(npt)));
-            max_dist = std::max(max_dist, dist_p_out(i));
-            if (code_out[i] == 1 && j == 0) { face_pts++; }
+            if (code_out[i] < 2)
+            {
+               for (int d = 0; d < dim; d++) { pos(d) = vxyz(d * pts_cnt + i); }
+               Vector exact_val(vec_dim);
+               F_exact(pos, exact_val);
+               max_err  = std::max(max_err, fabs(exact_val(j) - interp_vals(npt)));
+               max_dist = std::max(max_dist, dist_p_out(i));
+               if (code_out[i] == 1 && j == 0) { face_pts++; }
+            }
+            else { if (j == 0) { not_found++; } }
+            npt++;
          }
-         else { if (j == 0) { not_found++; } }
-         npt++;
       }
-   }
 
-   // Print the results for task 0 since all tasks have the same set of points.
-   if (myid == 0)
-   {
       cout << setprecision(16)
            << "Searched unique points: " << pts_cnt
            << "\nFound on local mesh:  " << found_loc

miniapps/meshing/mesh-optimizer.cpp

@@ -39,6 +39,11 @@
 //   Adapted analytic shape+orientation:
 //     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
 //
+//   Adapted analytic shape and/or size with hr-adaptivity:
+//     mesh-optimizer -m square01.mesh -o 2 -tid 9  -ni 50 -li 20 -hmid 55 -mid 7 -hr
+//     mesh-optimizer -m square01.mesh -o 2 -tid 10 -ni 50 -li 20 -hmid 55 -mid 7 -hr
+//     mesh-optimizer -m square01.mesh -o 2 -tid 11 -ni 50 -li 20 -hmid 58 -mid 7 -hr
+//
 //   Adapted discrete size:
 //     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 80 -tid 5 -ni 50 -qo 4 -nor
 //   Adapted discrete size 3D with PA:
@@ -123,6 +128,8 @@ int main(int argc, char *argv[])
    int max_lin_iter      = 100;
    bool move_bnd         = true;
    int combomet          = 0;
+   bool hradaptivity     = false;
+   int h_metric_id       = -1;
    bool normalization    = false;
    bool visualization    = true;
    int verbosity_level   = 0;
@@ -131,6 +138,8 @@ int main(int argc, char *argv[])
    bool exactaction      = false;
    const char *devopt    = "cpu";
    bool pa               = false;
+   int n_hr_iter         = 5;
+   int n_h_iter          = 1;
 
    // 1. Parse command-line options.
    OptionsParser args(argc, argv);
@@ -221,6 +230,12 @@ int main(int argc, char *argv[])
                   "0: Use single metric\n\t"
                   "1: Shape + space-dependent size given analytically\n\t"
                   "2: Shape + adapted size given discretely; shared target");
+   args.AddOption(&hradaptivity, "-hr", "--hr-adaptivity", "-no-hr",
+                  "--no-hr-adaptivity",
+                  "Enable hr-adaptivity.");
+   args.AddOption(&h_metric_id, "-hmid", "--h-metric",
+                  "Same options as metric_id. Used to determine refinement"
+                  " type for each element if h-adaptivity is enabled.");
    args.AddOption(&normalization, "-nor", "--normalization", "-no-nor",
                   "--no-normalization",
                   "Make all terms in the optimization functional unitless.");
@@ -241,6 +256,11 @@ int main(int argc, char *argv[])
                   "Device configuration string, see Device::Configure().");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
                   "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&n_hr_iter, "-nhr", "--n_hr_iter",
+                  "Number of hr-adaptivity iterations.");
+   args.AddOption(&n_h_iter, "-nh", "--n_h_iter",
+                  "Number of h-adaptivity iterations per r-adaptivity"
+                  "iteration.");
    args.Parse();
    if (!args.Good())
    {
@@ -249,6 +269,13 @@ int main(int argc, char *argv[])
    }
    args.PrintOptions(cout);
 
+   if (h_metric_id < 0) { h_metric_id = metric_id; }
+
+   if (hradaptivity)
+   {
+      MFEM_VERIFY(strcmp(devopt,"cpu")==0, "HR-adaptivity is currently only"
+                  " supported on cpus.");
+   }
    Device device(devopt);
    device.Print();
 
@@ -261,7 +288,9 @@ int main(int argc, char *argv[])
    else { cout << "(NONE)"; }
    cout << endl;
 
-   // 3. Define a finite element space on the mesh. Here we use vector finite
+   if (hradaptivity) { mesh->EnsureNCMesh(); }
+
+   // 3. Define a finite element space on the mesh-> Here we use vector finite
    //    elements which are tensor products of quadratic finite elements. The
    //    number of components in the vector finite element space is specified by
    //    the last parameter of the FiniteElementSpace constructor.
@@ -392,9 +421,41 @@ int main(int argc, char *argv[])
          cout << "Unknown metric_id: " << metric_id << endl;
          return 3;
    }
+   TMOP_QualityMetric *h_metric = NULL;
+   if (hradaptivity)
+   {
+      switch (h_metric_id)
+      {
+         case 1: h_metric = new TMOP_Metric_001; break;
+         case 2: h_metric = new TMOP_Metric_002; break;
+         case 7: h_metric = new TMOP_Metric_007; break;
+         case 9: h_metric = new TMOP_Metric_009; break;
+         case 55: h_metric = new TMOP_Metric_055; break;
+         case 56: h_metric = new TMOP_Metric_056; break;
+         case 58: h_metric = new TMOP_Metric_058; break;
+         case 77: h_metric = new TMOP_Metric_077; break;
+         case 315: h_metric = new TMOP_Metric_315; break;
+         case 316: h_metric = new TMOP_Metric_316; break;
+         case 321: h_metric = new TMOP_Metric_321; break;
+         default: cout << "Metric_id not supported for h-adaptivity: " << h_metric_id <<
+                          endl;
+            return 3;
+      }
+   }
+
+   if (metric_id < 300 || h_metric_id < 300)
+   {
+      MFEM_VERIFY(dim == 2, "Incompatible metric for 3D meshes");
+   }
+   if (metric_id >= 300 || h_metric_id >= 300)
+   {
+      MFEM_VERIFY(dim == 3, "Incompatible metric for 2D meshes");
+   }
+
    TargetConstructor::TargetType target_t;
    TargetConstructor *target_c = NULL;
    HessianCoefficient *adapt_coeff = NULL;
+   HRHessianCoefficient *hr_adapt_coeff = NULL;
    H1_FECollection ind_fec(mesh_poly_deg, dim);
    FiniteElementSpace ind_fes(mesh, &ind_fec);
    FiniteElementSpace ind_fesv(mesh, &ind_fec, dim);
@@ -604,6 +665,18 @@ int main(int argc, char *argv[])
          target_c = tc;
          break;
       }
+      // Targets used for hr-adaptivity tests.
+      case 9:  // size target in an annular region.
+      case 10: // size+aspect-ratio in an annular region.
+      case 11: // size+aspect-ratio target for a rotate sine wave
+      {
+         target_t = TargetConstructor::GIVEN_FULL;
+         AnalyticAdaptTC *tc = new AnalyticAdaptTC(target_t);
+         hr_adapt_coeff = new HRHessianCoefficient(dim, target_id - 9);
+         tc->SetAnalyticTargetSpec(NULL, NULL, hr_adapt_coeff);
+         target_c = tc;
+         break;
+      }
       default: cout << "Unknown target_id: " << target_id << endl; return 3;
    }
    if (target_c == NULL)
@@ -611,7 +684,8 @@ int main(int argc, char *argv[])
       target_c = new TargetConstructor(target_t);
    }
    target_c->SetNodes(x0);
-   TMOP_Integrator *he_nlf_integ = new TMOP_Integrator(metric, target_c);
+   TMOP_Integrator *he_nlf_integ = new TMOP_Integrator(metric, target_c,
+                                                       h_metric);
 
    // Finite differences for computations of derivatives.
    if (fdscheme)
@@ -721,10 +795,10 @@ int main(int argc, char *argv[])
             TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE);
          target_c2->SetVolumeScale(0.01);
          target_c2->SetNodes(x0);
-         he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2);
+         he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2, h_metric);
          he_nlf_integ2->SetCoefficient(coeff2);
       }
-      else { he_nlf_integ2 = new TMOP_Integrator(metric2, target_c); }
+      else { he_nlf_integ2 = new TMOP_Integrator(metric2, target_c, h_metric); }
       he_nlf_integ2->SetIntegrationRules(*irules, quad_order);
       if (fdscheme) { he_nlf_integ2->EnableFiniteDifferences(x); }
       he_nlf_integ2->SetExactActionFlag(exactaction);
@@ -778,7 +852,9 @@ int main(int argc, char *argv[])
       tauval -= 0.01 * h0.Min();
    }
 
-   const double init_energy = a.GetGridFunctionEnergy(x);
+   // For HR tests, the energy is normalized by the number of elements.
+   const double init_energy = a.GetGridFunctionEnergy(x) /
+                              (hradaptivity ? mesh->GetNE() : 1);
 
    // Visualize the starting mesh and metric values.
    // Note that for combinations of metrics, this only shows the first metric.
@@ -906,13 +982,25 @@ int main(int argc, char *argv[])
       solver.SetAdaptiveLinRtol(solver_art_type, 0.5, 0.9);
    }
    solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
-   solver.SetOperator(a);
-   solver.Mult(b, x.GetTrueVector());
-   x.SetFromTrueVector();
-   if (solver.GetConverged() == false)
+
+   // hr-adaptivity solver.
+   // If hr-adaptivity is disabled, r-adaptivity is done once using the
+   // TMOPNewtonSolver.
+   // Otherwise, "hr_iter" iterations of r-adaptivity are done followed by
+   // "h_per_r_iter" iterations of h-adaptivity after each r-adaptivity.
+   // The solver terminates if an h-adaptivity iteration does not modify
+   // any element in the mesh.
+   TMOPHRSolver hr_solver(*mesh, a, solver,
+                          x, move_bnd, hradaptivity,
+                          mesh_poly_deg, h_metric_id,
+                          n_hr_iter, n_h_iter);
+   hr_solver.AddGridFunctionForUpdate(&x0);
+   if (adapt_lim_const > 0.)
    {
-      cout << "Nonlinear solver: rtol = " << solver_rtol << " not achieved.\n";
+      hr_solver.AddGridFunctionForUpdate(&zeta_0);
+      hr_solver.AddFESpaceForUpdate(&ind_fes);
    }
+   hr_solver.Mult();
 
    // 15. Save the optimized mesh to a file. This output can be viewed later
    //     using GLVis: "glvis -m optimized.mesh".
@@ -922,14 +1010,15 @@ int main(int argc, char *argv[])
       mesh->Print(mesh_ofs);
    }
 
-   // 16. Compute the amount of energy decrease.
-   const double fin_energy = a.GetGridFunctionEnergy(x);
+   const double fin_energy = a.GetGridFunctionEnergy(x) /
+                             (hradaptivity ? mesh->GetNE() : 1);
    double metric_part = fin_energy;
    if (lim_const > 0.0 || adapt_lim_const > 0.0)
    {
       lim_coeff.constant = 0.0;
       coef_zeta.constant = 0.0;
-      metric_part = a.GetGridFunctionEnergy(x);
+      metric_part = a.GetGridFunctionEnergy(x) /
+                    (hradaptivity ? mesh->GetNE() : 1);
       lim_coeff.constant = lim_const;
       coef_zeta.constant = adapt_lim_const;
    }
@@ -942,7 +1031,7 @@ int main(int argc, char *argv[])
    cout << "The strain energy decreased by: " << setprecision(12)
         << (init_energy - fin_energy) * 100.0 / init_energy << " %." << endl;
 
-   // 17. Visualize the final mesh and metric values.
+   // 16. Visualize the final mesh and metric values.
    if (visualization)
    {
       char title[] = "Final metric values";
@@ -956,13 +1045,13 @@ int main(int argc, char *argv[])
                              600, 600, 300, 300);
    }
 
-   // 18. Visualize the mesh displacement.
+   // 17. Visualize the mesh displacement.
    if (visualization)
    {
-      x0 -= x;
       osockstream sock(19916, "localhost");
       sock << "solution\n";
       mesh->Print(sock);
+      x0 -= x;
       x0.Save(sock);
       sock.send();
       sock << "window_title 'Displacements'\n"
@@ -971,7 +1060,6 @@ int main(int argc, char *argv[])
            << "keys jRmclA" << endl;
    }
 
-   // 19. Free the used memory.
    delete S;
    delete S_prec;
    delete target_c2;
@@ -979,7 +1067,9 @@ int main(int argc, char *argv[])
    delete coeff1;
    delete adapt_evaluator;
    delete target_c;
+   delete hr_adapt_coeff;
    delete adapt_coeff;
+   delete h_metric;
    delete metric;
    delete fespace;
    delete fec;

miniapps/meshing/mesh-optimizer.hpp

@@ -90,8 +90,6 @@ double discrete_aspr_2d(const Vector &x)
    double th = 22.5*M_PI/180.;
    double xn =  cos(th)*xc + sin(th)*yc;
    double yn = -sin(th)*xc + cos(th)*yc;
-   // double th2 = (th > 45.*M_PI/180) ? M_PI/2 - th : th;
-   // double stretch = 1/cos(th2);
    xc = xn; yc = yn;
 
    double tfac = 20;
@@ -223,6 +221,134 @@ public:
    }
 };
 
+class HRHessianCoefficient : public TMOPMatrixCoefficient
+{
+private:
+   int dim;
+   // 0 - size target in an annular region,
+   // 1 - size+aspect-ratio in an annular region,
+   // 2 - size+aspect-ratio target for a rotate sine wave.
+   int hr_target_type;
+
+public:
+   HRHessianCoefficient(int dim_, int hr_target_type_ = 0)
+      : TMOPMatrixCoefficient(dim_), dim(dim_),
+        hr_target_type(hr_target_type_) { }
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      Vector pos(3);
+      T.Transform(ip, pos);
+      if (hr_target_type == 0) // size only circle
+      {
+         double small = 0.001, big = 0.01;
+         if (dim == 3) { small = 0.005, big = 0.1; }
+         const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
+         double zc;
+         if (dim == 3) { zc = pos(2) - 0.5; }
+         double r = sqrt(xc*xc + yc*yc);
+         if (dim == 3) { r = sqrt(xc*xc + yc*yc + zc*zc); }
+         double r1 = 0.15; double r2 = 0.35; double sf=30.0;
+
+         const double tan1 = std::tanh(sf*(r-r1)),
+                      tan2 = std::tanh(sf*(r-r2));
+
+         double ind = (tan1 - tan2);
+         if (ind > 1.0) {ind = 1.;}
+         if (ind < 0.0) {ind = 0.;}
+         double val = ind * small + (1.0 - ind) * big;
+         K = 0.0;
+         K(0, 0) = 1.0;
+         K(0, 1) = 0.0;
+         K(1, 0) = 0.0;
+         K(1, 1) = 1.0;
+         K(0, 0) *= pow(val,0.5);
+         K(1, 1) *= pow(val,0.5);
+         if (dim == 3) { K(2, 2) = pow(val,0.5); }
+      }
+      else if (hr_target_type == 1) // circle with size and AR
+      {
+         const double small = 0.001, big = 0.01;
+         const double xc = pos(0)-0.5, yc = pos(1)-0.5;
+         const double rv = xc*xc + yc*yc;
+         double r = 0;
+         if (rv>0.) {r = sqrt(rv);}
+
+         double r1 = 0.2; double r2 = 0.3; double sf=30.0;
+         const double szfac = 1;
+         const double asfac = 4;
+         const double eps2 = szfac/asfac;
+         const double eps1 = szfac;
+
+         double tan1 = std::tanh(sf*(r-r1)+1),
+                tan2 = std::tanh(sf*(r-r2)-1);
+         double wgt = 0.5*(tan1-tan2);
+
+         tan1 = std::tanh(sf*(r-r1)),
+         tan2 = std::tanh(sf*(r-r2));
+
+         double ind = (tan1 - tan2);
+         if (ind > 1.0) {ind = 1.;}
+         if (ind < 0.0) {ind = 0.;}
+         double szval = ind * small + (1.0 - ind) * big;
+
+         double th = std::atan2(yc,xc)*180./M_PI;
+         if (wgt > 1) { wgt = 1; }
+         if (wgt < 0) { wgt = 0; }
+
+         double maxval = eps2 + eps1*(1-wgt)*(1-wgt);
+         double minval = eps1;
+         double avgval = 0.5*(maxval+minval);
+         double ampval = 0.5*(maxval-minval);
+         double val1 = avgval + ampval*sin(2.*th*M_PI/180.+90*M_PI/180.);
+         double val2 = avgval + ampval*sin(2.*th*M_PI/180.-90*M_PI/180.);
+
+         K(0,1) = 0.0;
+         K(1,0) = 0.0;
+         K(0,0) = val1;
+         K(1,1) = val2;
+
+         K(0,0) *= pow(szval,0.5);
+         K(1,1) *= pow(szval,0.5);
+      }
+      else if (hr_target_type == 2) // sharp rotated sine wave
+      {
+         double xc = pos(0)-0.5, yc = pos(1)-0.5;
+         double th = 15.5*M_PI/180.;
+         double xn =  cos(th)*xc + sin(th)*yc;
+         double yn = -sin(th)*xc + cos(th)*yc;
+         double th2 = (th > 45.*M_PI/180) ? M_PI/2 - th : th;
+         double stretch = 1/cos(th2);
+         xc = xn/stretch;
+         yc = yn;
+         double tfac = 20;
+         double s1 = 3;
+         double s2 = 2;
+         double yl1 = -0.025;
+         double yl2 =  0.025;
+         double wgt = std::tanh((tfac*(yc-yl1) + s2*std::sin(s1*M_PI*xc)) + 1) -
+                      std::tanh((tfac*(yc-yl2) + s2*std::sin(s1*M_PI*xc)) - 1);
+         if (wgt > 1) { wgt = 1; }
+         if (wgt < 0) { wgt = 0; }
+
+         const double eps2 = 25;
+         const double eps1 = 1;
+         K(1,1) = eps1/eps2 + eps1*(1-wgt)*(1-wgt);
+         K(0,0) = eps1;
+         K(0,1) = 0.0;
+         K(1,0) = 0.0;
+      }
+      else { MFEM_ABORT("Unsupported option / wrong input."); }
+   }
+
+   virtual void EvalGrad(DenseMatrix &K, ElementTransformation &T,
+                         const IntegrationPoint &ip, int comp)
+   {
+      K = 0.;
+   }
+};
+
 // Additional IntegrationRules that can be used with the --quad-type option.
 IntegrationRules IntRulesLo(0, Quadrature1D::GaussLobatto);
 IntegrationRules IntRulesCU(0, Quadrature1D::ClosedUniform);
@@ -252,13 +378,13 @@ double adapt_lim_fun(const Vector &x)
 
 void DiffuseField(GridFunction &field, int smooth_steps)
 {
-   //Setup the Laplacian operator
+   // Setup the Laplacian operator
    BilinearForm *Lap = new BilinearForm(field.FESpace());
    Lap->AddDomainIntegrator(new DiffusionIntegrator());
    Lap->Assemble();
    Lap->Finalize();
 
-   //Setup the smoothing operator
+   // Setup the smoothing operator
    DSmoother *S = new DSmoother(0,1.0,smooth_steps);
    S->iterative_mode = true;
    S->SetOperator(Lap->SpMat());
@@ -274,7 +400,7 @@ void DiffuseField(GridFunction &field, int smooth_steps)
 #ifdef MFEM_USE_MPI
 void DiffuseField(ParGridFunction &field, int smooth_steps)
 {
-   //Setup the Laplacian operator
+   // Setup the Laplacian operator
    ParBilinearForm *Lap = new ParBilinearForm(field.ParFESpace());
    Lap->AddDomainIntegrator(new DiffusionIntegrator());
    Lap->Assemble();

miniapps/meshing/pmesh-optimizer.cpp

@@ -39,6 +39,11 @@
 //   Adapted analytic shape+orientation:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
 //
+//   Adapted analytic shape and/or size with hr-adaptivity:
+//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -tid 9  -ni 50 -li 20 -hmid 55 -mid 7 -hr
+//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -tid 10 -ni 50 -li 20 -hmid 55 -mid 7 -hr
+//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -tid 11 -ni 50 -li 20 -hmid 58 -mid 7 -hr
+//
 //   Adapted discrete size:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 80 -tid 5 -ni 50 -qo 4 -nor
 //   Adapted discrete size 3D with PA:
@@ -130,6 +135,8 @@ int main (int argc, char *argv[])
    int max_lin_iter      = 100;
    bool move_bnd         = true;
    int combomet          = 0;
+   bool hradaptivity     = false;
+   int h_metric_id       = -1;
    bool normalization    = false;
    bool visualization    = true;
    int verbosity_level   = 0;
@@ -138,6 +145,8 @@ int main (int argc, char *argv[])
    bool exactaction      = false;
    const char *devopt    = "cpu";
    bool pa               = false;
+   int n_hr_iter         = 5;
+   int n_h_iter          = 1;
 
    // 2. Parse command-line options.
    OptionsParser args(argc, argv);
@@ -230,6 +239,12 @@ int main (int argc, char *argv[])
                   "0: Use single metric\n\t"
                   "1: Shape + space-dependent size given analytically\n\t"
                   "2: Shape + adapted size given discretely; shared target");
+   args.AddOption(&hradaptivity, "-hr", "--hr-adaptivity", "-no-hr",
+                  "--no-hr-adaptivity",
+                  "Enable hr-adaptivity.");
+   args.AddOption(&h_metric_id, "-hmid", "--h-metric",
+                  "Same options as metric_id. Used to determine refinement"
+                  " type for each element if h-adaptivity is enabled.");
    args.AddOption(&normalization, "-nor", "--normalization", "-no-nor",
                   "--no-normalization",
                   "Make all terms in the optimization functional unitless.");
@@ -250,6 +265,11 @@ int main (int argc, char *argv[])
                   "Device configuration string, see Device::Configure().");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
                   "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&n_hr_iter, "-nhr", "--n_hr_iter",
+                  "Number of hr-adaptivity iterations.");
+   args.AddOption(&n_h_iter, "-nh", "--n_h_iter",
+                  "Number of h-adaptivity iterations per r-adaptivity"
+                  "iteration.");
    args.Parse();
    if (!args.Good())
    {
@@ -257,7 +277,13 @@ int main (int argc, char *argv[])
       return 1;
    }
    if (myid == 0) { args.PrintOptions(cout); }
+   if (h_metric_id < 0) { h_metric_id = metric_id; }
 
+   if (hradaptivity)
+   {
+      MFEM_VERIFY(strcmp(devopt,"cpu")==0, "HR-adaptivity is currently only"
+                  " supported on cpus.");
+   }
    Device device(devopt);
    if (myid == 0) { device.Print();}
 
@@ -276,6 +302,7 @@ int main (int argc, char *argv[])
       cout << endl;
    }
 
+   if (hradaptivity) { mesh->EnsureNCMesh(); }
    ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
 
    delete mesh;
@@ -422,9 +449,41 @@ int main (int argc, char *argv[])
          if (myid == 0) { cout << "Unknown metric_id: " << metric_id << endl; }
          return 3;
    }
+   TMOP_QualityMetric *h_metric = NULL;
+   if (hradaptivity)
+   {
+      switch (h_metric_id)
+      {
+         case 1: h_metric = new TMOP_Metric_001; break;
+         case 2: h_metric = new TMOP_Metric_002; break;
+         case 7: h_metric = new TMOP_Metric_007; break;
+         case 9: h_metric = new TMOP_Metric_009; break;
+         case 55: h_metric = new TMOP_Metric_055; break;
+         case 56: h_metric = new TMOP_Metric_056; break;
+         case 58: h_metric = new TMOP_Metric_058; break;
+         case 77: h_metric = new TMOP_Metric_077; break;
+         case 315: h_metric = new TMOP_Metric_315; break;
+         case 316: h_metric = new TMOP_Metric_316; break;
+         case 321: h_metric = new TMOP_Metric_321; break;
+         default: cout << "Metric_id not supported for h-adaptivity: " << h_metric_id <<
+                          endl;
+            return 3;
+      }
+   }
+
+   if (metric_id < 300 || h_metric_id < 300)
+   {
+      MFEM_VERIFY(dim == 2, "Incompatible metric for 3D meshes");
+   }
+   if (metric_id >= 300 || h_metric_id >= 300)
+   {
+      MFEM_VERIFY(dim == 3, "Incompatible metric for 2D meshes");
+   }
+
    TargetConstructor::TargetType target_t;
    TargetConstructor *target_c = NULL;
    HessianCoefficient *adapt_coeff = NULL;
+   HRHessianCoefficient *hr_adapt_coeff = NULL;
    H1_FECollection ind_fec(mesh_poly_deg, dim);
    ParFiniteElementSpace ind_fes(pmesh, &ind_fec);
    ParFiniteElementSpace ind_fesv(pmesh, &ind_fec, dim);
@@ -639,6 +698,18 @@ int main (int argc, char *argv[])
          target_c = tc;
          break;
       }
+      // Targets used for hr-adaptivity tests.
+      case 9:  // size target in an annular region.
+      case 10: // size+aspect-ratio in an annular region.
+      case 11: // size+aspect-ratio target for a rotate sine wave
+      {
+         target_t = TargetConstructor::GIVEN_FULL;
+         AnalyticAdaptTC *tc = new AnalyticAdaptTC(target_t);
+         hr_adapt_coeff = new HRHessianCoefficient(dim, target_id - 9);
+         tc->SetAnalyticTargetSpec(NULL, NULL, hr_adapt_coeff);
+         target_c = tc;
+         break;
+      }
       default:
          if (myid == 0) { cout << "Unknown target_id: " << target_id << endl; }
          return 3;
@@ -649,13 +720,13 @@ int main (int argc, char *argv[])
       target_c = new TargetConstructor(target_t, MPI_COMM_WORLD);
    }
    target_c->SetNodes(x0);
-   TMOP_Integrator *he_nlf_integ= new TMOP_Integrator(metric, target_c);
+   TMOP_Integrator *he_nlf_integ = new TMOP_Integrator(metric, target_c,
+                                                       h_metric);
 
    // Finite differences for computations of derivatives.
    if (fdscheme)
    {
       MFEM_VERIFY(pa == false, "PA for finite differences is not implemented.");
-
       he_nlf_integ->EnableFiniteDifferences(x);
    }
    he_nlf_integ->SetExactActionFlag(exactaction);
@@ -753,7 +824,8 @@ int main (int argc, char *argv[])
       he_nlf_integ->SetCoefficient(*coeff1);
 
       // Second metric.
-      metric2 = new TMOP_Metric_077;
+      if (dim == 2) { metric2 = new TMOP_Metric_077; }
+      else          { metric2 = new TMOP_Metric_315; }
       TMOP_Integrator *he_nlf_integ2 = NULL;
       if (combomet == 1)
       {
@@ -761,10 +833,10 @@ int main (int argc, char *argv[])
             TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE, MPI_COMM_WORLD);
          target_c2->SetVolumeScale(0.01);
          target_c2->SetNodes(x0);
-         he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2);
+         he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2, h_metric);
          he_nlf_integ2->SetCoefficient(coeff2);
       }
-      else { he_nlf_integ2 = new TMOP_Integrator(metric2, target_c); }
+      else { he_nlf_integ2 = new TMOP_Integrator(metric2, target_c, h_metric); }
       he_nlf_integ2->SetIntegrationRules(*irules, quad_order);
       if (fdscheme) { he_nlf_integ2->EnableFiniteDifferences(x); }
       he_nlf_integ2->SetExactActionFlag(exactaction);
@@ -777,7 +849,10 @@ int main (int argc, char *argv[])
 
       a.AddDomainIntegrator(combo);
    }
-   else { a.AddDomainIntegrator(he_nlf_integ); }
+   else
+   {
+      a.AddDomainIntegrator(he_nlf_integ);
+   }
 
    if (pa) { a.Setup(); }
 
@@ -821,7 +896,9 @@ int main (int argc, char *argv[])
       tauval -= 0.01 * h0min_all;
    }
 
-   const double init_energy = a.GetParGridFunctionEnergy(x);
+   // For HR tests, the energy is normalized by the number of elements.
+   const double init_energy = a.GetParGridFunctionEnergy(x) /
+                              (hradaptivity ? pmesh->GetGlobalNE() : 1);
 
    // Visualize the starting mesh and metric values.
    // Note that for combinations of metrics, this only shows the first metric.
@@ -952,13 +1029,25 @@ int main (int argc, char *argv[])
       solver.SetAdaptiveLinRtol(solver_art_type, 0.5, 0.9);
    }
    solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
-   solver.SetOperator(a);
-   solver.Mult(b, x.GetTrueVector());
-   x.SetFromTrueVector();
-   if (myid == 0 && solver.GetConverged() == false)
+
+   // hr-adaptivity solver.
+   // If hr-adaptivity is disabled, r-adaptivity is done once using the
+   // TMOPNewtonSolver.
+   // Otherwise, "hr_iter" iterations of r-adaptivity are done followed by
+   // "h_per_r_iter" iterations of h-adaptivity after each r-adaptivity.
+   // The solver terminates if an h-adaptivity iteration does not modify
+   // any element in the mesh.
+   TMOPHRSolver hr_solver(*pmesh, a, solver,
+                          x, move_bnd, hradaptivity,
+                          mesh_poly_deg, h_metric_id,
+                          n_hr_iter, n_h_iter);
+   hr_solver.AddGridFunctionForUpdate(&x0);
+   if (adapt_lim_const > 0.)
    {
-      cout << "Nonlinear solver: rtol = " << solver_rtol << " not achieved.\n";
+      hr_solver.AddGridFunctionForUpdate(&zeta_0);
+      hr_solver.AddFESpaceForUpdate(&ind_fes);
    }
+   hr_solver.Mult();
 
    // 16. Save the optimized mesh to a file. This output can be viewed later
    //     using GLVis: "glvis -m optimized -np num_mpi_tasks".
@@ -971,13 +1060,15 @@ int main (int argc, char *argv[])
    }
 
    // 17. Compute the amount of energy decrease.
-   const double fin_energy = a.GetParGridFunctionEnergy(x);
+   const double fin_energy = a.GetParGridFunctionEnergy(x) /
+                             (hradaptivity ? pmesh->GetGlobalNE() : 1);
    double metric_part = fin_energy;
    if (lim_const > 0.0 || adapt_lim_const > 0.0)
    {
       lim_coeff.constant = 0.0;
       coef_zeta.constant = 0.0;
-      metric_part = a.GetParGridFunctionEnergy(x);
+      metric_part = a.GetParGridFunctionEnergy(x) /
+                    (hradaptivity ? pmesh->GetGlobalNE() : 1);
       lim_coeff.constant = lim_const;
       coef_zeta.constant = adapt_lim_const;
    }
@@ -1036,7 +1127,9 @@ int main (int argc, char *argv[])
    delete coeff1;
    delete adapt_evaluator;
    delete target_c;
+   delete hr_adapt_coeff;
    delete adapt_coeff;
+   delete h_metric;
    delete metric;
    delete pfespace;
    delete fec;

tests/gitlab/build_and_test

@@ -56,8 +56,10 @@ then
         prefix_opt="--prefix=${prefix}"
     fi
 
+    echo "Fetching uberenv."
     tests/gitlab/get_mfem_uberenv || ( echo "Error fetching Uberenv" && exit 1 );
 
+    echo "Installing dependencies."
     python tests/uberenv/uberenv.py --spec="${spec}" "${upstream_opt}" "${prefix_opt}"
 
 fi

tests/gitlab/get_mfem_uberenv

@@ -18,7 +18,7 @@ set -o errexit
 set -o nounset
 
 uberenv_url="https://github.com/mfem/mfem-uberenv.git"
-uberenv_ref="dbe35eb4dbeaabc7184ad110fcb98377092baeb0"
+uberenv_ref="222e3743377d8c903e044f80ae44624e1668feb4"
 
 [[ ! -d tests/uberenv ]] && git clone ${uberenv_url} tests/uberenv
 cd tests/uberenv