changes in the configuration files

Fixing GetValue methods in parallel [bugfix/gf-getvalue-dev]

CHANGELOG

@@ -109,6 +109,10 @@ New and updated examples and miniapps
   for applying Dirichlet, Neumann (both homogeneous and inhomogeneous), Robin,
   and periodic boundary conditions with either H1 or DG discretizations.
 
+- Added a new miniapp, Navier, that solves the time-dependent Navier-Stokes
+  equations of incompressible fluid dynamics. See the miniapps/navier directory
+  for more details.
+
 - Added a simple meshing miniapp, Twist, which demonstrates MFEM's strategy of
   stitching together opposite surfaces of a mesh to create a topologically
   periodic mesh.

CMakeLists.txt

@@ -292,6 +292,25 @@ if (MFEM_USE_HIOP)
   # find_package updates HIOP_FOUND, HIOP_INCLUDE_DIRS, HIOP_LIBRARIES
 endif()
 
+# ADEPT package
+if (MFEM_USE_ADEPT)
+        find_package(ADEPT REQUIRED)
+        # find_package updates ADEPT_FOUND, ADEPT_INCLUDE_DIRS, ADEPT_LIBRARIES
+endif()
+
+# CODIPACK package
+if (MFEM_USE_CODIPACK)
+	find_package(CODIPACK REQUIRED)
+	# find_package updates CODIPACK_FOUND, CODIPACK_INCLUDE_DIRS, CODIPACK_LIBRARIES
+endif()
+
+# FADBAD++ package
+if (MFEM_USE_FADBADPP)
+	find_package(FADBADPP REQUIRED)
+        # find_package updates FADBADPP_FOUND, FADBADPP_INCLUDE_DIRS, FADBADPP_LIBRARIES
+endif()
+
+
 # CUDA
 if (MFEM_USE_CUDA)
    set(CMAKE_CUDA_STANDARD 11)
@@ -353,7 +372,7 @@ endif()
 #    be before SuiteSparse.
 set(MFEM_TPLS MPI_CXX OPENMP BLAS LAPACK METIS HYPRE SuiteSparse SUNDIALS PETSC
     MESQUITE SuperLUDist STRUMPACK AXOM CONDUIT Ginkgo GNUTLS GSLIB NETCDF
-    MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE ADIOS2)
+    MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE ADIOS2 ADEPT CODIPACK FADBADPP) 
 # Add all *_FOUND libraries in the variable TPL_LIBRARIES.
 set(TPL_LIBRARIES "")
 set(TPL_INCLUDE_DIRS "")

config/cmake/config.hpp.in

@@ -153,4 +153,11 @@
 // library.
 #cmakedefine MFEM_USE_SIMMETRIX
 
+#cmakedefine MFEM_USE_ADEPT
+
+#cmakedefine MFEM_USE_CODIPACK
+
+#cmakedefine MFEM_USE_FADBADPP
+
+
 #endif // MFEM_CONFIG_HEADER

config/cmake/modules/FindADEPT.cmake

@@ -0,0 +1,23 @@
+# Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the
+# Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability see http://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the GNU Lesser General Public License (as published by the Free
+# Software Foundation) version 2.1 dated February 1999.
+
+# Sets the following variables:
+#   - ADEPT_FOUND
+#   - ADEPT_INCLUDE_DIRS
+#   - ADEPT_LIBRARIES
+
+include(MfemCmakeUtilities)
+mfem_find_package(ADEPT ADEPT ADEPT_DIR
+  "include" "adept.hpp"
+  "lib" "libadept.so"
+  "Paths to headers required by ADEPT."
+  "Libraries required by ADEPT.")
+

config/cmake/modules/FindCODIPACK.cmake

@@ -0,0 +1,23 @@
+# Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the
+# Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability see http://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the GNU Lesser General Public License (as published by the Free
+# Software Foundation) version 2.1 dated February 1999.
+
+# Sets the following variables:
+#   - CODIPACK_FOUND
+#   - CODIPACK_INCLUDE_DIRS
+#   - CODIPACK_LIBRARIES
+
+include(MfemCmakeUtilities)
+mfem_find_package(CODIPACK CODIPACK CODIPACK_DIR
+  "include" "codi.hpp"
+  "lib" ""
+  "Paths to headers required by CODIPACK."
+  "Libraries required by CODIPACK.")
+

config/cmake/modules/FindFADBADPP.cmake

@@ -0,0 +1,23 @@
+# Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the
+# Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability see http://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the GNU Lesser General Public License (as published by the Free
+# Software Foundation) version 2.1 dated February 1999.
+
+# Sets the following variables:
+#   - FADBADPP_FOUND
+#   - FADBADPP_INCLUDE_DIRS
+#   - FADBADPP_LIBRARIES
+
+include(MfemCmakeUtilities)
+mfem_find_package(FADBADPP FADBADPP FADBADPP_DIR
+  "include" "fadiff.h"
+  "lib" ""
+  "Paths to headers required by FADBADPP."
+  "Libraries required by FADBADPP.")
+

config/defaults.cmake

@@ -51,6 +51,9 @@ option(MFEM_USE_CEED "Enable CEED" OFF)
 option(MFEM_USE_UMPIRE "Enable Umpire" OFF)
 option(MFEM_USE_SIMD "Enable use of SIMD intrinsics" ON)
 option(MFEM_USE_ADIOS2 "Enable ADIOS2" OFF)
+option(MFEM_USE_ADEPT "Enable AD using ADEPT"   OFF)
+option(MFEM_USE_CODIPACK "Enable AD using CoDiPack" OFF)
+option(MFEM_USE_FADBADPP "Enable AD using FADBAD++" OFF)
 
 set(MFEM_MPI_NP 4 CACHE STRING "Number of processes used for MPI tests")
 
@@ -183,6 +186,16 @@ set(BLAS_LIBRARIES "" CACHE STRING "The BLAS library.")
 set(LAPACK_INCLUDE_DIRS "" CACHE STRING "Path to LAPACK headers.")
 set(LAPACK_LIBRARIES "" CACHE STRING "The LAPACK library.")
 
+set(ADEPT_INCLUDE_DIRS "${MFEM_DIR}/../adept-1.1/include" CACHE STRING "Path to ADEPT headers.")
+set(ADEPT_LIBRARIES    "-L${MFEM_DIR}/../adept-1.1/lib -ladept" CACHE STRING "The ADEPT library.")
+
+set(CODIPACK_INCLUDE_DIRS "${MFEM_DIR}/../CoDiPack/include" CACHE STRING "Path to CoDiPack headers.")
+set(CODIPACK_LIBRARIES    "")
+
+set(FADBADPP_INCLUDE_DIRS "${MFEM_DIR}/../FADBAD++" CACHE STRING "Path to FADBAD++ headers.")
+set(FADBADPP_LIBRARIES    "")
+
+
 # Some useful variables:
 set(CMAKE_SKIP_PREPROCESSED_SOURCE_RULES ON) # Skip *.i rules
 set(CMAKE_SKIP_ASSEMBLY_SOURCE_RULES  ON)    # Skip *.s rules

config/user.cmake

@@ -0,0 +1,204 @@
+# 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.
+
+# See the file INSTALL for description of the configuration options.
+
+# Default options. To replace these, copy this file to user.cmake and modify it.
+
+if (NOT CMAKE_BUILD_TYPE)
+set(CMAKE_BUILD_TYPE "Debug" CACHE STRING
+      "Build type: Debug, Release, RelWithDebInfo, or MinSizeRel." FORCE)
+endif()
+
+# MFEM options. Set to mimic the default "defaults.mk" file.
+option(MFEM_USE_MPI "Enable MPI parallel build" ON)
+option(MFEM_USE_METIS "Enable METIS usage" ${MFEM_USE_MPI})
+option(MFEM_USE_EXCEPTIONS "Enable the use of exceptions" OFF)
+option(MFEM_USE_ZLIB "Enable zlib for compressed data streams." OFF)
+option(MFEM_USE_LIBUNWIND "Enable backtrace for errors." ON)
+option(MFEM_USE_LAPACK "Enable LAPACK usage" ON)
+option(MFEM_THREAD_SAFE "Enable thread safety" OFF)
+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_MESQUITE "Enable MESQUITE usage" OFF)
+option(MFEM_USE_SUITESPARSE "Enable SuiteSparse usage" ON)
+option(MFEM_USE_SUPERLU "Enable SuperLU_DIST usage" OFF)
+option(MFEM_USE_STRUMPACK "Enable STRUMPACK usage" OFF)
+option(MFEM_USE_GINKGO "Enable Ginkgo usage" OFF)
+option(MFEM_USE_GNUTLS "Enable GNUTLS usage" OFF)
+option(MFEM_USE_GSLIB "Enable GSLIB usage" OFF)
+option(MFEM_USE_NETCDF "Enable NETCDF usage" OFF)
+option(MFEM_USE_PETSC "Enable PETSc support." ON)
+option(MFEM_USE_MPFR "Enable MPFR usage." OFF)
+option(MFEM_USE_SIDRE "Enable Axom/Sidre usage" OFF)
+option(MFEM_USE_CONDUIT "Enable Conduit usage" OFF)
+option(MFEM_USE_PUMI "Enable PUMI" OFF)
+option(MFEM_USE_HIOP "Enable HiOp" OFF)
+option(MFEM_USE_CUDA "Enable CUDA" OFF)
+option(MFEM_USE_OCCA "Enable OCCA" OFF)
+option(MFEM_USE_RAJA "Enable RAJA" OFF)
+option(MFEM_USE_CEED "Enable CEED" OFF)
+option(MFEM_USE_UMPIRE "Enable Umpire" OFF)
+option(MFEM_USE_SIMD "Enable use of SIMD intrinsics" ON)
+option(MFEM_USE_ADIOS2 "Enable ADIOS2" OFF)
+option(MFEM_USE_ADEPT "Enable AD using ADEPT"   OFF)
+option(MFEM_USE_CODIPACK "Enable AD using CoDiPack" ON)
+option(MFEM_USE_FADBADPP "Enable AD using FADBAD++" OFF)
+
+
+
+set(MFEM_MPI_NP 4 CACHE STRING "Number of processes used for MPI tests")
+
+# Allow a user to disable testing, examples, and/or miniapps at CONFIGURE TIME
+# if they don't want/need them (e.g. if MFEM is "just a dependency" and all they
+# need is the library, building all that stuff adds unnecessary overhead). Note
+# that the examples or miniapps can always be built using the targets 'examples'
+# or 'miniapps', respectively.
+option(MFEM_ENABLE_TESTING "Enable the ctest framework for testing" ON)
+option(MFEM_ENABLE_EXAMPLES "Build all of the examples" OFF)
+option(MFEM_ENABLE_MINIAPPS "Build all of the miniapps" OFF)
+
+# Setting CXX/MPICXX on the command line or in user.cmake will overwrite the
+# autodetected C++ compiler.
+# set(CXX g++)
+# set(MPICXX mpicxx)
+
+# Set the target CUDA architecture
+set(CUDA_ARCH "sm_60" CACHE STRING "Target CUDA architecture.")
+
+set(MFEM_DIR ${CMAKE_CURRENT_SOURCE_DIR})
+
+# The *_DIR paths below will be the first place searched for the corresponding
+# headers and library. If these fail, then standard cmake search is performed.
+# Note: if the variables are already in the cache, they are not overwritten.
+
+set(HYPRE_DIR "/home/blaz/develop/common/dbg/petsc_3.12.5/" CACHE PATH
+    "Path to the hypre library.")
+# If hypre was compiled to depend on BLAS and LAPACK:
+# set(HYPRE_REQUIRED_PACKAGES "BLAS" "LAPACK" CACHE STRING
+#     "Packages that HYPRE depends on.")
+
+set(METIS_DIR "/home/blaz/develop/common/dbg/petsc_3.12.5/" CACHE PATH "Path to the METIS library.")
+
+set(LIBUNWIND_DIR "" CACHE PATH "Path to Libunwind.")
+
+set(SUNDIALS_DIR "/home/blaz/develop/common/dbg/SUNDIALS_5.2.0/" CACHE PATH
+    "Path to the SUNDIALS library.")
+# The following may be necessary, if SUNDIALS was built with KLU:
+# set(SUNDIALS_REQUIRED_PACKAGES "SuiteSparse/KLU/AMD/BTF/COLAMD/config"
+#     CACHE STRING "Additional packages required by SUNDIALS.")
+
+set(MESQUITE_DIR "${MFEM_DIR}/../mesquite-2.99" CACHE PATH
+    "Path to the Mesquite library.")
+
+set(SuiteSparse_DIR "/home/blaz/develop/common/dbg/petsc_3.12.5/" CACHE PATH
+    "Path to the SuiteSparse library.")
+set(SuiteSparse_REQUIRED_PACKAGES "BLAS" "METIS"
+    CACHE STRING "Additional packages required by SuiteSparse.")
+
+set(ParMETIS_DIR "/home/blaz/develop/common/dbg/petsc_3.12.5/" CACHE PATH
+    "Path to the ParMETIS library.")
+set(ParMETIS_REQUIRED_PACKAGES "METIS" CACHE STRING
+    "Additional packages required by ParMETIS.")
+
+set(SuperLUDist_DIR "${MFEM_DIR}/../SuperLU_DIST_5.1.0" CACHE PATH
+    "Path to the SuperLU_DIST library.")
+# SuperLU_DIST may also depend on "OpenMP", depending on how it was compiled.
+set(SuperLUDist_REQUIRED_PACKAGES "MPI" "BLAS" "ParMETIS" CACHE STRING
+    "Additional packages required by SuperLU_DIST.")
+
+set(STRUMPACK_DIR "/home/blaz/develop/common/dbg/petsc_3.12.5/" CACHE PATH
+    "Path to the STRUMPACK library.")
+# STRUMPACK may also depend on "OpenMP", depending on how it was compiled.
+# Starting with v2.2.0 of STRUMPACK, ParMETIS and Scotch are optional.
+set(STRUMPACK_REQUIRED_PACKAGES "MPI" "MPI_Fortran" "ParMETIS" "METIS"
+    "ScaLAPACK" "Scotch/ptscotch/ptscotcherr/scotch/scotcherr" CACHE STRING
+    "Additional packages required by STRUMPACK.")
+# If the MPI package does not find all required Fortran libraries:
+# set(STRUMPACK_REQUIRED_LIBRARIES "gfortran" "mpi_mpifh" CACHE STRING
+#     "Additional libraries required by STRUMPACK.")
+
+# The Scotch library, required by STRUMPACK <= v2.1.0, optional in STRUMPACK >=
+# v2.2.0.
+set(Scotch_DIR "${MFEM_DIR}/../scotch_6.0.4" CACHE PATH
+    "Path to the Scotch and PT-Scotch libraries.")
+set(Scotch_REQUIRED_PACKAGES "Threads" CACHE STRING
+    "Additional packages required by Scotch.")
+# Tell the "Threads" package/module to prefer pthreads.
+set(CMAKE_THREAD_PREFER_PTHREAD TRUE)
+set(Threads_LIB_VARS CMAKE_THREAD_LIBS_INIT)
+
+# The ScaLAPACK library, required by STRUMPACK
+set(ScaLAPACK_DIR "/home/blaz/develop/common/dbg/petsc_3.12.5/"
+    CACHE PATH "Path to the configuration file scalapack-config.cmake")
+set(ScaLAPACK_TARGET_NAMES scalapack)
+# set(ScaLAPACK_TARGET_FORCE)
+# set(ScaLAPACK_IMPORT_CONFIG DEBUG)
+
+set(Ginkgo_DIR "${MFEM_DIR}/../ginkgo" CACHE PATH "Path to the Ginkgo library.")
+
+set(GNUTLS_DIR "" CACHE PATH "Path to the GnuTLS library.")
+
+set(GSLIB_DIR "" CACHE PATH "Path to the GSLIB library.")
+
+set(NETCDF_DIR "" CACHE PATH "Path to the NetCDF library.")
+# May need to add "HDF5" as requirement.
+set(NetCDF_REQUIRED_PACKAGES "" CACHE STRING
+    "Additional packages required by NetCDF.")
+
+set(PETSC_DIR "/home/blaz/develop/common/dbg/petsc_3.12.5/" CACHE PATH
+    "Path to the PETSc main directory.")
+set(PETSC_ARCH "" CACHE STRING "PETSc build architecture.")
+
+set(MPFR_DIR "" CACHE PATH "Path to the MPFR library.")
+
+set(CONDUIT_DIR "${MFEM_DIR}/../conduit" CACHE PATH
+    "Path to the Conduit library.")
+
+set(AXOM_DIR "${MFEM_DIR}/../axom" CACHE PATH "Path to the Axom library.")
+# May need to add "Boost" as requirement.
+set(Axom_REQUIRED_PACKAGES "Conduit/relay/blueprint" CACHE STRING
+    "Additional packages required by Axom.")
+
+set(PUMI_DIR "${MFEM_DIR}/../pumi-2.1.0" CACHE STRING
+    "Directory where PUMI is installed")
+
+set(HIOP_DIR "${MFEM_DIR}/../hiop/install" CACHE STRING
+    "Directory where HiOp is installed")
+set(HIOP_REQUIRED_PACKAGES "BLAS" "LAPACK" CACHE STRING
+    "Packages that HiOp depends on.")
+
+set(OCCA_DIR "${MFEM_DIR}/../occa" CACHE PATH "Path to OCCA")
+set(RAJA_DIR "${MFEM_DIR}/../raja" CACHE PATH "Path to RAJA")
+set(CEED_DIR "${MFEM_DIR}/../libCEED" CACHE PATH "Path to libCEED")
+set(UMPIRE_DIR "${MFEM_DIR}/../umpire" CACHE PATH "Path to Umpire")
+
+set(BLAS_INCLUDE_DIRS "" CACHE STRING "Path to BLAS headers.")
+set(BLAS_LIBRARIES "-L/home/blaz/develop/common/lib -lblas" CACHE STRING "The BLAS library.")
+set(LAPACK_INCLUDE_DIRS "" CACHE STRING "Path to LAPACK headers.")
+set(LAPACK_LIBRARIES "-L/home/blaz/develop/common/lib -llapack" CACHE STRING "The LAPACK library.")
+
+
+set(ADEPT_INCLUDE_DIRS "/home/blaz/develop/common/dbg/adept-1.1/include" CACHE STRING "Path to ADEPT headers.")
+set(ADEPT_LIBRARIES    "/home/blaz/develop/common/dbg/adept-1.1/lib/libadept.so" CACHE STRING "The ADEPT library.")
+
+set(CODIPACK_INCLUDE_DIRS "/home/blaz/develop/common/CoDiPack/include" CACHE STRING "Path to CoDiPack headers.")
+set(CODIPACK_LIBRARIES    "")
+
+set(FADBADPP_INCLUDE_DIRS "/home/blaz/develop/common/FADBAD++" CACHE STRING "Path to FADBAD++ headers.")
+set(FADBADPP_LIBRARIES    "")
+
+# Some useful variables:
+set(CMAKE_SKIP_PREPROCESSED_SOURCE_RULES ON) # Skip *.i rules
+set(CMAKE_SKIP_ASSEMBLY_SOURCE_RULES  ON)    # Skip *.s rules
+# set(CMAKE_VERBOSE_MAKEFILE ON CACHE BOOL "Verbose makefiles.")

examples/CMakeLists.txt

@@ -34,6 +34,8 @@ list(APPEND ALL_EXE_SRCS
   ex25.cpp
   ex26.cpp
   ex27.cpp
+  ex51.cpp
+  ex71.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -64,6 +66,7 @@ if (MFEM_USE_MPI)
     ex25p.cpp
     ex26p.cpp
     ex27p.cpp
+    ex71p.cpp
     )
 endif()
 

examples/ex18.hpp

@@ -88,8 +88,6 @@ private:
    Vector funval2;
    Vector nor;
    Vector fluxN;
-   IntegrationPoint eip1;
-   IntegrationPoint eip2;
 
 public:
    FaceIntegrator(RiemannSolver &rsolver_, const int dim);
@@ -424,19 +422,16 @@ void FaceIntegrator::AssembleFaceVector(const FiniteElement &el1,
    {
       const IntegrationPoint &ip = ir->IntPoint(i);
 
-      Tr.Loc1.Transform(ip, eip1);
-      Tr.Loc2.Transform(ip, eip2);
+      Tr.SetAllIntPoints(&ip); // set face and element int. points
 
       // Calculate basis functions on both elements at the face
-      el1.CalcShape(eip1, shape1);
-      el2.CalcShape(eip2, shape2);
+      el1.CalcShape(Tr.GetElement1IntPoint(), shape1);
+      el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
 
       // Interpolate elfun at the point
       elfun1_mat.MultTranspose(shape1, funval1);
       elfun2_mat.MultTranspose(shape2, funval2);
 
-      Tr.SetIntPoint(&ip);
-
       // Get the normal vector and the flux on the face
       CalcOrtho(Tr.Jacobian(), nor);
       const double mcs = rsolver.Eval(funval1, funval2, nor, fluxN);

examples/ex51.cpp

@@ -0,0 +1,629 @@
+//                                MFEM Example 1
+//
+// Compile with: make ex1
+//
+// Sample runs:  ex1 -m ../data/square-disc.mesh
+//               ex1 -m ../data/star.mesh
+//               ex1 -m ../data/star-mixed.mesh
+//               ex1 -m ../data/escher.mesh
+//               ex1 -m ../data/fichera.mesh
+//               ex1 -m ../data/fichera-mixed.mesh
+//               ex1 -m ../data/toroid-wedge.mesh
+//               ex1 -m ../data/square-disc-p2.vtk -o 2
+//               ex1 -m ../data/square-disc-p3.mesh -o 3
+//               ex1 -m ../data/square-disc-nurbs.mesh -o -1
+//               ex1 -m ../data/star-mixed-p2.mesh -o 2
+//               ex1 -m ../data/disc-nurbs.mesh -o -1
+//               ex1 -m ../data/pipe-nurbs.mesh -o -1
+//               ex1 -m ../data/fichera-mixed-p2.mesh -o 2
+//               ex1 -m ../data/star-surf.mesh
+//               ex1 -m ../data/square-disc-surf.mesh
+//               ex1 -m ../data/inline-segment.mesh
+//               ex1 -m ../data/amr-quad.mesh
+//               ex1 -m ../data/amr-hex.mesh
+//               ex1 -m ../data/fichera-amr.mesh
+//               ex1 -m ../data/mobius-strip.mesh
+//               ex1 -m ../data/mobius-strip.mesh -o -1 -sc
+//
+// Device sample runs:
+//               ex1 -pa -d cuda
+//               ex1 -pa -d raja-cuda
+//               ex1 -pa -d occa-cuda
+//               ex1 -pa -d raja-omp
+//               ex1 -pa -d occa-omp
+//               ex1 -pa -d ceed-cpu
+//               ex1 -pa -d ceed-cuda
+//               ex1 -m ../data/beam-hex.mesh -pa -d cuda
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order, or if order < 1 using an isoparametric/isogeometric
+//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of mesh refinement, finite
+//               element grid functions, as well as linear and bilinear forms
+//               corresponding to the left-hand side and right-hand side of the
+//               discrete linear system. We also cover the explicit elimination
+//               of essential boundary conditions, static condensation, and the
+//               optional connection to the GLVis tool for visualization.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+
+#include "../fem/adnonlininteg.hpp"
+
+
+using namespace std;
+
+namespace mfem{
+    
+class VolNonlinearForm: public NonlinearFormIntegrator
+{
+protected:
+    double eta;
+    double beta;
+public:
+    VolNonlinearForm(double eta_, double beta_){
+        eta=eta_;
+        beta=beta_;}
+    virtual ~VolNonlinearForm(){ }
+    
+    double Project(double inp)
+    {
+        // tanh projection - Wang&Lazarov&Sigmund2011 
+        double a=std::tanh(eta*beta);
+        double b=std::tanh(beta*(1.0-eta));
+        double c=std::tanh(beta*(inp-eta));
+        double rez=(a+c)/(a+b);
+        return rez;
+    }
+
+    double ProjGrad(double inp)
+    {
+        double c=std::tanh(beta*(inp-eta));
+        double a=std::tanh(eta*beta);
+        double b=std::tanh(beta*(1.0-eta));
+        double rez=beta*(1.0-c*c)/(a+b);
+        return rez;
+    }
+
+    
+    double ProjSec(double inp)
+    {
+        double c=std::tanh(beta*(inp-eta));
+        double a=std::tanh(eta*beta);
+        double b=std::tanh(beta*(1.0-eta));
+        double rez=-2.0*beta*beta*c*(1.0-c*c)/(a+b);
+        return rez;
+    }
+    
+    
+    virtual double GetElementEnergy(const mfem::FiniteElement & el,
+                                           mfem::ElementTransformation & trans,
+                                           const mfem::Vector & elfun) override
+    {
+        double energy=0.0;
+        int ndof = el.GetDof();
+        int ndim = el.GetDim();
+        
+        const mfem::IntegrationRule *ir = NULL;
+        int order = 2 * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+        
+        mfem::Vector shapef(ndof);
+        
+        double w;
+        for (int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            el.CalcShape(ip,shapef);
+            w= Project(shapef*elfun);
+            w= ip.weight * trans.Weight() * w;
+            energy = energy + w;
+        }
+        return energy;
+    }
+    
+    virtual void AssembleElementVector(const mfem::FiniteElement & el,
+                                       mfem::ElementTransformation & trans,
+                                       const mfem::Vector & elfun,
+                                       mfem::Vector & elvect) override
+    {
+        
+        
+        int ndof = el.GetDof();
+    
+        const mfem::IntegrationRule *ir = NULL;
+        int order = 2 * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+        
+        elvect.SetSize(ndof);
+        elvect=0.0;
+        
+        mfem::Vector shapef(ndof);
+        double w;
+        for (int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            el.CalcShape(ip,shapef);
+            w= ProjGrad(shapef*elfun);
+            w= ip.weight * trans.Weight() * w;
+            elvect.Add(w,shapef);
+        }
+        
+    }
+    
+    
+    virtual void AssembleElementGrad(const mfem::FiniteElement & el,
+                                     mfem::ElementTransformation & trans,
+                                     const mfem::Vector & elfun,
+                                     mfem::DenseMatrix & elmat) override
+    {
+        int ndof = el.GetDof();
+    
+        const mfem::IntegrationRule *ir = NULL;
+        int order = 2 * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+        
+        elmat.SetSize(ndof);
+        elmat=0.0;
+        
+        mfem::Vector shapef(ndof);
+        double w;
+        for (int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            el.CalcShape(ip,shapef);
+            w= ProjSec(shapef*elfun);
+            w= ip.weight * trans.Weight() * w;
+            AddMult_a_VVt(w, shapef, elmat);
+        }
+        
+   }
+    
+};
+
+
+class VolNonlinearFormADH:public ADNonlinearFormIntegratorH
+{
+private:
+    double eta;
+    double beta;
+    
+    template<typename DType>
+    DType Project(DType inp)
+    {
+        // tanh projection - Wang&Lazarov&Sigmund2011 
+        double a=std::tanh(eta*beta);
+        double b=std::tanh(beta*(1.0-eta));
+        DType c=tanh(beta*(inp-eta));
+        DType rez=(a+c)/(a+b);
+        return rez;
+    }
+    
+public:
+    
+    VolNonlinearFormADH(double eta_, double beta_){
+        eta=eta_;
+        beta=beta_;
+    }
+    
+    
+    virtual ADFType ElementEnergy(const mfem::FiniteElement & el,
+                                       mfem::ElementTransformation & trans,
+                                       const ADFVector & elfun) override
+    {
+        ADFType rez=MyElementEnergy<ADFType,ADFVector>(el,trans,elfun);
+        return rez;
+    }
+    
+    virtual ADSType ElementEnergy(const mfem::FiniteElement & el,
+                                       mfem::ElementTransformation & trans,
+                                       const ADSVector & elfun) override
+    {
+        return MyElementEnergy<ADSType,ADSVector>(el,trans,elfun);
+    }
+    
+    template<typename MDType, typename MVType>
+    MDType MyElementEnergy(const mfem::FiniteElement & el,
+                            mfem::ElementTransformation & trans,
+                            const MVType & elfun)
+    {
+        MDType energy=MDType();
+        int ndof = el.GetDof();
+        
+        const mfem::IntegrationRule *ir = NULL;
+        int order = 2 * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+        
+        mfem::Vector shapef(ndof);
+        
+        MDType w;
+        for (int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            el.CalcShape(ip,shapef);
+            w= Project(elfun*shapef);
+            w= ip.weight * trans.Weight() * w;
+            energy = energy + w;
+        }
+        return energy;
+    }
+    
+    
+    virtual double ElementEnergy(const mfem::FiniteElement & el, 
+                                 mfem::ElementTransformation & Tr, 
+                                 const mfem::Vector & elfun) override
+    {
+        return GetElementEnergy(el,Tr,elfun);
+    }
+    
+    virtual double GetElementEnergy(const mfem::FiniteElement & el,
+                                           mfem::ElementTransformation & trans,
+                                           const mfem::Vector & elfun) override
+    {
+        double rez;
+        rez=MyElementEnergy<double,mfem::Vector>(el,trans,elfun);
+        return rez;
+    
+    }
+};
+
+
+class VolQIntegratorJ: public ADQIntegratorJ
+{
+  private:
+    template<typename DType>
+    DType Project(double eta, double beta, DType inp)
+    {
+        // tanh projection - Wang&Lazarov&Sigmund2011
+        double a=std::tanh(eta*beta);
+        double b=std::tanh(beta*(1.0-eta));
+        DType c=tanh(beta*(inp-eta));
+        DType rez=(a+c)/(a+b);
+        return rez;
+    }
+
+    template<typename DType>
+    DType ProjGrad(double eta, double beta, DType inp)
+    {
+        DType c=tanh(beta*(inp-eta));
+        DType a=tanh(eta*beta);
+        DType b=tanh(beta*(1.0-eta));
+        DType rez=beta*(1.0-c*c)/(a+b);
+        return rez;
+    }
+
+
+  public:
+    VolQIntegratorJ(){}
+    virtual ~VolQIntegratorJ(){}
+
+
+    template<typename MVType>
+    void MyQIntegratorDU(const mfem::Vector& vparam, MVType& uu, MVType& rr)
+    {
+        //implement all evaluations executed at integration point
+        double eta=vparam[0];
+        double beta=vparam[1];
+        rr.SetSize(1); //return the derivative of the projected value
+        rr[0]=ProjGrad(eta,beta,uu[0]);
+        return;
+    }
+
+    virtual void QIntegratorDU(const mfem::Vector& vparam, mfem::Vector& uu, mfem::Vector& rr) override
+    {
+        MyQIntegratorDU<mfem::Vector>(vparam,uu,rr);
+    }
+
+    virtual void QIntegratorDU(const mfem::Vector& vparam, ADFVector& uu, ADFVector& rr) override
+    {
+        MyQIntegratorDU<ADFVector>(vparam,uu,rr);
+    }
+
+
+    virtual double QIntegrator(const Vector &vparam, const Vector &uu) override
+    {
+        //implement all evaluations executed at integration point
+        double eta=vparam[0];
+        double beta=vparam[1];
+        double rez=Project(eta,beta,uu[0]);
+        return rez;
+    }
+
+};
+
+
+class VolNonlinearFormQJ: public NonlinearFormIntegrator
+{
+protected:
+    double eta;
+    double beta;
+    mfem::Vector vparam;
+    VolQIntegratorJ qint;
+public:
+    VolNonlinearFormQJ(double eta_, double beta_){
+        eta=eta_;
+        beta=beta_;
+        vparam.SetSize(2);
+        vparam[0]=eta;
+        vparam[1]=beta;
+    }
+
+    virtual ~VolNonlinearFormQJ(){ }
+
+    virtual double GetElementEnergy(const mfem::FiniteElement & el,
+                                           mfem::ElementTransformation & trans,
+                                           const mfem::Vector & elfun) override
+    {
+        double energy=0.0;
+        int ndof = el.GetDof();
+
+        const mfem::IntegrationRule *ir = NULL;
+        int order = 2 * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+        mfem::Vector shapef(ndof);
+        mfem::Vector uu(1);
+
+        double w;
+        for (int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            el.CalcShape(ip,shapef);
+            uu[0]= shapef*elfun;
+            w= qint.QIntegrator(vparam,uu);
+            w= ip.weight * trans.Weight() * w;
+            energy = energy + w;
+        }
+        return energy;
+    }
+
+    virtual void AssembleElementVector(const mfem::FiniteElement & el,
+                                       mfem::ElementTransformation & trans,
+                                       const mfem::Vector & elfun,
+                                       mfem::Vector & elvect) override
+    {
+        int ndof = el.GetDof();
+        const mfem::IntegrationRule *ir = NULL;
+        int order = 2 * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+        elvect.SetSize(ndof);
+        elvect=0.0;
+
+        mfem::Vector shapef(ndof);
+        mfem::Vector uu(1);
+        mfem::Vector rr(1);
+
+        double w;
+
+        for (int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            el.CalcShape(ip,shapef);
+            uu[0]=shapef*elfun;
+            qint.QIntegratorDU(vparam,uu,rr);
+            w= ip.weight * trans.Weight() * rr[0];
+            elvect.Add(w,shapef);
+        }
+    }
+
+    virtual void AssembleElementGrad(const mfem::FiniteElement & el,
+                                     mfem::ElementTransformation & trans,
+                                     const mfem::Vector & elfun,
+                                     mfem::DenseMatrix & elmat) override
+    {
+        int ndof = el.GetDof();
+
+        const mfem::IntegrationRule *ir = NULL;
+        int order = 2 * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+        elmat.SetSize(ndof);
+        elmat=0.0;
+
+        mfem::DenseMatrix jac(1,1);
+
+        mfem::Vector shapef(ndof);
+        mfem::Vector uu(1);
+
+
+        double w;
+        for (int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            el.CalcShape(ip,shapef);
+            uu[0]=shapef*elfun;
+            qint.QIntegratorDD(vparam,uu,jac);
+            w= ip.weight * trans.Weight() * jac(0,0);
+            AddMult_a_VVt(w, shapef, elmat);
+        }
+
+   }
+
+};
+
+
+}
+
+
+
+double TFunc(const mfem::Vector& a){
+    double sca=4.0;
+    double rez=(std::sin(sca*a[0])*std::sin(sca*a[1])*std::sin(sca*a[2]))*0.5+0.5;
+    return rez;
+}
+
+
+int main(int argc, char *argv[])
+{
+   
+    // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = true;
+
+   mfem::OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   mfem::Device device(device_config);
+   device.Print();
+
+   // 3. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+   
+   
+   // 4. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
+   //    largest number that gives a final mesh with no more than 50,000
+   //    elements.
+   
+   {
+      int ref_levels =
+         (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
+      ref_levels=1;
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+   
+
+   // 5. Define a finite element space on the mesh. Here we use continuous
+   //    Lagrange finite elements of the specified order. If order < 1, we
+   //    instead use an isoparametric/isogeometric space.
+   mfem::FiniteElementCollection *fec;
+   if (order > 0)
+   {
+      fec = new mfem::H1_FECollection(order, dim);
+   }
+   else if (mesh->GetNodes())
+   {
+      fec = mesh->GetNodes()->OwnFEC();
+      cout << "Using isoparametric FEs: " << fec->Name() << endl;
+   }
+   else
+   {
+      fec = new mfem::H1_FECollection(order = 1, dim);
+   }
+   mfem::FiniteElementSpace *fespace = new mfem::FiniteElementSpace(mesh, fec);
+   cout << "Number of finite element unknowns: "
+        << fespace->GetTrueVSize() << endl;
+   
+   mfem::NonlinearForm* nf0=new mfem::NonlinearForm(fespace);
+   mfem::NonlinearForm* nf1=new mfem::NonlinearForm(fespace);
+
+   mfem::FunctionCoefficient ifun(TFunc);
+   //create an input for the NonlinearForm
+   mfem::GridFunction* igf = new mfem::GridFunction(fespace);
+   igf->ProjectCoefficient(ifun);
+   
+   std::cout << "Size of the grid function igf:"<<igf->Size()<<std::endl;
+   
+   
+   mfem::Vector* resv0=new mfem::Vector(fespace->GetTrueVSize());
+   mfem::Vector* resv1=new mfem::Vector(fespace->GetTrueVSize());
+   mfem::Vector* stat=new mfem::Vector(fespace->GetTrueVSize());
+   
+   igf->GetTrueDofs(*stat);
+   
+   
+   //compute the energy - the total volume above 0.5
+   nf0->AddDomainIntegrator(new mfem::VolNonlinearForm(0.5,8.0));
+   //nf1->AddDomainIntegrator(new mfem::VolNonlinearFormADH(0.5,8.0));
+   nf1->AddDomainIntegrator(new mfem::VolNonlinearFormQJ(0.5,8.0));
+   
+   double vol0=nf0->GetEnergy(*stat);
+   double vol1=nf1->GetEnergy(*stat);
+   std::cout<<"The total volume is:("<<vol0<<","<<vol1<<")"<<std::endl;
+   nf0->Mult(*stat,*resv0);
+   nf1->Mult(*stat,*resv1);
+   //project back the gradients to a grid function
+   mfem::GridFunction* ggf0=new mfem::GridFunction(fespace);
+   ggf0->SetFromTrueDofs(*resv0);
+   mfem::GridFunction* ggf1=new mfem::GridFunction(fespace);
+   ggf1->SetFromTrueDofs(*resv1);
+   
+   
+   resv0->Add(-1.0,*resv1);
+   std::cout<<"Norm|v_1-v_0|="<<resv0->Norml2()<<std::endl;
+   
+   mfem::Operator& grad0(nf0->GetGradient(*stat));
+   mfem::SparseMatrix* spmat0=dynamic_cast<mfem::SparseMatrix*>(&grad0);
+   mfem::Operator& grad1(nf1->GetGradient(*stat));
+   mfem::SparseMatrix* spmat1=dynamic_cast<mfem::SparseMatrix*>(&grad1);
+   std::cout<<"Norm mat1="<<spmat0->MaxNorm()<<" mat2="<<spmat1->MaxNorm()<<std::endl;
+   spmat0->Add(-1.0,*spmat1);
+   std::cout<<"Norm diff"<<spmat0->MaxNorm()<<std::endl;
+   {
+       std::fstream mstr; 
+       mstr.open("mat.dat",std::ios::out);
+       spmat0->PrintMatlab(mstr);
+       mstr.close();
+   }
+   
+   mfem::ParaViewDataCollection *dacol=new mfem::ParaViewDataCollection("IGF_OUT",mesh);
+   dacol->SetLevelsOfDetail(2);
+   dacol->SetCycle(1);
+   dacol->SetTime(0.0); // set the time
+   dacol->RegisterField("density",igf);
+   dacol->RegisterField("grads0",ggf0);
+   dacol->RegisterField("grads1",ggf1);
+   dacol->Save();
+   delete dacol;
+   
+   delete ggf0;
+   delete ggf1;
+   delete stat;
+   delete resv0;
+   delete resv1;
+   delete igf;     
+   delete nf0;
+   delete nf1;
+   delete fespace;
+   delete fec;
+   delete mesh;
+   
+   return 0;
+}

examples/ex71.cpp

@@ -0,0 +1,349 @@
+//                       MFEM Example 71 - Serial Version
+//
+// Compile with: make ex71
+//
+// Sample runs:
+//    ex71 -m ../data/beam-quad.mesh
+//    ex71 -m ../data/beam-tri.mesh
+//    ex71 -m ../data/beam-hex.mesh
+//    ex71 -m ../data/beam-tet.mesh
+//    ex71 -m ../data/beam-wedge.mesh
+//
+// Description:  This examples solves a quasi-static nonlinear
+//               pLaplacian problem with zero Dirichlet boundary
+//				 conditions applied on all defined boundaries
+//
+//               The example demonstrates the use of nonlinear operators
+//				 combined with automatic differentiation (AD). The definitions
+//				 of the integrators are written in the ex71.hpp.
+//				 Selecting integrator=0 will use handcoded integrator.
+//				 Selecting integrator=1 will utilize AD integrator.
+//				 The AD integrator can be modifief to use ADQFunctionJ
+//				 or ADQFunctionH by overwritting the class type of qint,
+//				 i.e., pLapIntegrandJ or pLapIntegrandH.
+//
+//				 qint (the integrand) is a function which is evaluated
+//				 at every integration point. For implementations utilizing
+//				 ADQFunctionJ, the user has to implement the function and the
+//				 residual evaluation - all virtual methods. The Jacobian of
+//				 the residual is evaluated using AD
+//
+//				 For implementations utilizing ADQFunctionH, the user has
+//				 to implement only the function evaluation (preferebaly as
+//				 a template) and the first derivative (the residual) and the
+//				 second derivatives (the Hessian) are evaluated using AD.
+//
+//               We recommend viewing examples 1 and 19, before viewing this
+//               example.
+
+#include "ex71.hpp"
+
+#undef MFEM_USE_SUITESPARSE
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options
+   const char *mesh_file = "../data/beam-tet.mesh";
+   int ser_ref_levels = 3;
+   int order = 1;
+   bool visualization = true;
+   double newton_rel_tol = 1e-4;
+   double newton_abs_tol = 1e-6;
+   int newton_iter = 500;
+   int print_level = 0;
+   double pp = 2.0;
+   int integrator=0;
+   mfem::StopWatch* timer=new mfem::StopWatch();
+
+   mfem::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(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&newton_rel_tol, "-rel", "--relative-tolerance",
+                  "Relative tolerance for the Newton solve.");
+   args.AddOption(&newton_abs_tol, "-abs", "--absolute-tolerance",
+                  "Absolute tolerance for the Newton solve.");
+   args.AddOption(&newton_iter, "-it", "--newton-iterations",
+                  "Maximum iterations for the Newton solve.");
+   args.AddOption(&pp, "-pp", "--power-parameter",
+                  "Power parameter (>=2.0) for the p-Laplacian.");
+   args.AddOption((&print_level),"-prt","--print-level",
+                  "Print level.");
+   args.AddOption(&integrator, "-int","--integrator",
+                  "Integrator 0: standard; 1: AD uaing energy; 2: AD using gradients");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(std::cout);
+      return 1;
+   }
+   args.PrintOptions(std::cout);
+
+
+   // 2. Read the (serial) mesh from the given mesh file.
+   mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 3. 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();
+   }
+
+   // 4. Define the power parameter for the p-Laplacian and all other
+   //    coefficients
+   mfem::ConstantCoefficient c_pp(pp);
+   mfem::ConstantCoefficient load(1.000000000);
+   mfem::ConstantCoefficient c_ee(0.000000001);
+
+   // 5. Define the finite element spaces for the solution
+   mfem::H1_FECollection fec(order,dim);
+   mfem::FiniteElementSpace fespace(mesh,&fec,1,mfem::Ordering::byVDIM);
+   int glob_size=fespace.GetTrueVSize();
+   std::cout << "Number of finite element unknowns: " << glob_size << std::endl;
+
+   // 6. Define the Dirichlet conditions
+   mfem::Array<int> ess_bdr(mesh->bdr_attributes.Max());
+   ess_bdr = 1;
+
+   // 7. Define the nonlinear form
+   mfem::NonlinearForm* nf=new mfem::NonlinearForm(&fespace);
+
+   // 8. Define the solution vector x
+   mfem::GridFunction x(&fespace);
+   x = 0.0;
+   mfem::Vector tv(fespace.GetTrueVSize());
+   mfem::Vector sv(fespace.GetTrueVSize());
+   tv=0.0;
+   sv=0.0;
+
+   // 9. Define ParaView DataCollection
+   mfem::ParaViewDataCollection *dacol=new mfem::ParaViewDataCollection("pLap",mesh);
+   dacol->SetLevelsOfDetail(order);
+   dacol->RegisterField("sol",&x);
+
+
+   // 11. Set domain integrators - start with linear diffusion
+   {
+       // the default power coefficient is 2.0
+       mfem::ConstantCoefficient lpp(2.0);
+       if(integrator==0)
+       {
+            nf->AddDomainIntegrator(new mfem::pLaplace(lpp,c_ee,load));
+       }else
+       if(integrator==1)
+       {
+           nf->AddDomainIntegrator(new mfem::pLaplaceAD(lpp,c_ee,load));
+       }
+       nf->SetEssentialBC(ess_bdr);
+       // compute the energy
+       double energy=nf->GetEnergy(tv);
+       std::cout<<"[2] The total energy of the system is E="<<energy<<std::endl;
+       // time the assembly
+       timer->Clear();
+       timer->Start();
+       mfem::Operator &op=nf->GetGradient(sv);
+       timer->Stop();
+       std::cout<<"[2] The assembly time is: "<<timer->RealTime()<<std::endl;
+       mfem::Solver *prec;
+#ifdef MFEM_USE_SUITESPARSE
+       prec=new mfem::UMFPackSolver();
+#else
+       prec=new mfem::GSSmoother();
+#endif
+       mfem::CGSolver *j_pcg = new mfem::CGSolver();
+       j_pcg->SetRelTol(1e-7);
+       j_pcg->SetAbsTol(1e-15);
+       j_pcg->SetMaxIter(500);
+       j_pcg->SetPrintLevel(print_level);
+       j_pcg->SetPreconditioner(*prec);
+
+       mfem::NewtonSolver* ns;
+       ns=new mfem::NewtonSolver();
+       ns->iterative_mode = true;
+       ns->SetSolver(*j_pcg);
+       ns->SetOperator(*nf);
+       ns->SetPrintLevel(print_level);
+       ns->SetRelTol(1e-6);
+       ns->SetAbsTol(1e-12);
+       ns->SetMaxIter(10);
+       //solve the problem
+       timer->Clear();
+       timer->Start();
+       ns->Mult(tv, sv);
+       timer->Stop();
+       std::cout<<"Time for the NewtonSolver: "<<timer->RealTime()<<std::endl;
+
+       energy=nf->GetEnergy(sv);
+       std::cout<<"[pp=2] The total energy of the system is E="<<energy<<std::endl;
+
+       delete ns;
+       delete j_pcg;
+       delete prec;
+
+       x.SetFromTrueDofs(sv);
+       dacol->SetTime(2.0);
+       dacol->SetCycle(2);
+       dacol->Save();
+   }
+
+   // 12. Continue with powers higher than 2
+   for(int i=3;i<pp;i++)
+   {
+       delete nf;
+       nf=new mfem::NonlinearForm(&fespace);
+       mfem::ConstantCoefficient lpp((double)i);
+       if(integrator==0)
+       {
+            nf->AddDomainIntegrator(new mfem::pLaplace(lpp,c_ee,load));
+       }else
+       if(integrator==1)
+       {
+           nf->AddDomainIntegrator(new mfem::pLaplaceAD(lpp,c_ee,load));
+       }
+       nf->SetEssentialBC(ess_bdr);
+       // compute the energy
+       double energy=nf->GetEnergy(sv);
+       std::cout<<"[pp="<<i<<"] The total energy of the system is E="<<energy<<std::endl;
+       // time the assembly
+       timer->Clear();
+       timer->Start();
+       mfem::Operator &op=nf->GetGradient(sv);
+       timer->Stop();
+       std::cout<<"[pp="<<i<<"] The assembly time is: "<<timer->RealTime()<<std::endl;
+       mfem::Solver *prec;
+#ifdef MFEM_USE_SUITESPARSE
+       prec=new mfem::UMFPackSolver();
+#else
+       prec=new mfem::GSSmoother();
+#endif
+       mfem::CGSolver *j_pcg = new mfem::CGSolver();
+       j_pcg->SetRelTol(1e-7);
+       j_pcg->SetAbsTol(1e-15);
+       j_pcg->SetMaxIter(500);
+       j_pcg->SetPrintLevel(print_level);
+       j_pcg->SetPreconditioner(*prec);
+
+       mfem::NewtonSolver* ns;
+       ns=new mfem::NewtonSolver();
+       ns->iterative_mode = true;
+       ns->SetSolver(*j_pcg);
+       ns->SetOperator(*nf);
+       ns->SetPrintLevel(print_level);
+       ns->SetRelTol(1e-6);
+       ns->SetAbsTol(1e-12);
+       ns->SetMaxIter(10);
+       //solve the problem
+       timer->Clear();
+       timer->Start();
+       ns->Mult(tv, sv);
+       timer->Stop();
+       std::cout<<"Time for the NewtonSolver: "<<timer->RealTime()<<std::endl;
+
+       energy=nf->GetEnergy(sv);
+       std::cout<<"[pp="<<i<<"] The total energy of the system is E="<<energy<<std::endl;
+
+       delete ns;
+       delete j_pcg;
+       delete prec;
+
+       x.SetFromTrueDofs(sv);
+       dacol->SetTime(i);
+       dacol->SetCycle(i);
+       dacol->Save();
+   }
+
+   // 13. Continue with the final power
+   if( std::abs(pp-2.0) > std::numeric_limits<double>::epsilon())
+   {
+       delete nf;
+       nf=new mfem::NonlinearForm(&fespace);
+       if(integrator==0)
+       {
+            nf->AddDomainIntegrator(new mfem::pLaplace(c_pp,c_ee,load));
+       }else
+       if(integrator==1)
+       {
+           nf->AddDomainIntegrator(new mfem::pLaplaceAD(c_pp,c_ee,load));
+       }
+       nf->SetEssentialBC(ess_bdr);
+       // compute the energy
+       double energy=nf->GetEnergy(sv);
+       std::cout<<"[pp="<<pp<<"] The total energy of the system is E="<<energy<<std::endl;
+       // time the assembly
+       timer->Clear();
+       timer->Start();
+       mfem::Operator &op=nf->GetGradient(sv);
+       timer->Stop();
+       std::cout<<"[pp="<<pp<<"] The assembly time is: "<<timer->RealTime()<<std::endl;
+       mfem::Solver *prec;
+#ifdef MFEM_USE_SUITESPARSE
+       prec=new mfem::UMFPackSolver();
+#else
+       prec=new mfem::GSSmoother();
+#endif
+       mfem::CGSolver *j_pcg = new mfem::CGSolver();
+       j_pcg->SetRelTol(1e-7);
+       j_pcg->SetAbsTol(1e-15);
+       j_pcg->SetMaxIter(500);
+       j_pcg->SetPrintLevel(print_level);
+       j_pcg->SetPreconditioner(*prec);
+
+       mfem::NewtonSolver* ns;
+       ns=new mfem::NewtonSolver();
+       ns->iterative_mode = true;
+       ns->SetSolver(*j_pcg);
+       ns->SetOperator(*nf);
+       ns->SetPrintLevel(print_level);
+       ns->SetRelTol(1e-6);
+       ns->SetAbsTol(1e-12);
+       ns->SetMaxIter(10);
+       //solve the problem
+       timer->Clear();
+       timer->Start();
+       ns->Mult(tv, sv);
+       timer->Stop();
+       std::cout<<"Time for the NewtonSolver: "<<timer->RealTime()<<std::endl;
+
+       energy=nf->GetEnergy(sv);
+       std::cout<<"[pp="<<pp<<"] The total energy of the system is E="<<energy<<std::endl;
+
+       delete ns;
+       delete j_pcg;
+       delete prec;
+
+       x.SetFromTrueDofs(sv);
+       dacol->SetTime(pp);
+       if(pp<2.0)
+       {
+           dacol->SetCycle(std::floor(pp));
+       }
+       else
+       {
+           dacol->SetCycle(std::ceil(pp));
+       }
+       dacol->Save();
+   }
+
+
+
+   // 19. Free the used memory
+   delete dacol;
+   delete nf;
+   delete mesh;
+   delete timer;
+
+
+   return 0;
+}
+
+

examples/ex71.hpp

@@ -0,0 +1,619 @@
+// shared implementation ex71p/ex71 for the AD integrands and
+// the handconded integrators
+
+
+#ifndef EXAMPLE71_H
+#define EXAMPLE71_H
+
+#include "mfem.hpp"
+#include <memory>
+#include <iostream>
+#include <fstream>
+
+
+namespace mfem {
+
+class pLapIntegrandJ: public ADQFunctionJ
+{
+private:
+
+    template<typename DType, typename MVType>
+    void MyQFunctionDU(const mfem::Vector& vparam, MVType& uu, MVType& rr)
+    {
+        double pp=vparam[0];
+        double ee=vparam[1];
+        double ff=vparam[2];
+
+        DType norm2=uu[0]*uu[0]+uu[1]*uu[1]+uu[2]*uu[2];
+        DType tvar=pow(ee*ee+norm2,(pp-2.0)/2.0);
+
+        rr[0]=tvar*uu[0];
+        rr[1]=tvar*uu[1];
+        rr[2]=tvar*uu[2];
+        rr[3]=-ff;
+
+    }
+
+public:
+    pLapIntegrandJ():ADQFunctionJ(4){} //the residual vector rr has size of 4 elements
+
+    ~pLapIntegrandJ(){}
+
+    virtual double QFunction(const mfem::Vector &vparam,const mfem::Vector &uu) override
+    {
+        double pp=vparam[0];
+        double ee=vparam[1];
+        double ff=vparam[2];
+
+        double  u=uu[3];
+        double  norm2=uu[0]*uu[0]+uu[1]*uu[1]+uu[2]*uu[2];
+
+        double rez= pow(ee*ee+norm2,pp/2.0)/pp-ff*u;
+        return rez;
+    }
+
+
+    virtual void QFunctionDU(const mfem::Vector& vparam, mfem::Vector& uu, mfem::Vector& rr) override
+    {
+        MyQFunctionDU<double,mfem::Vector>(vparam,uu,rr);
+    }
+
+    virtual void QFunctionDU(const mfem::Vector &vparam, ADFVector &uu, ADFVector &rr) override
+    {
+        MyQFunctionDU<ADFType,ADFVector>(vparam,uu,rr);
+    }
+};
+
+
+class pLapIntegrandH: public ADQFunctionH
+{
+private:
+
+    //MVType - vector type taking one of the following
+    // mfem::Vector - scalar double
+    // ADFVector    - scalar ADFType
+    // ADSVector    - scalar ADSType
+    template<typename DType, typename MVType>
+    DType MyQFunction(const mfem::Vector& vparam, MVType& uu)
+    {
+        double pp=vparam[0];
+        double ee=vparam[1];
+        double ff=vparam[2];
+
+        DType  u=uu[3];
+        DType  norm2=uu[0]*uu[0]+uu[1]*uu[1]+uu[2]*uu[2];
+
+        DType rez= pow(ee*ee+norm2,pp/2.0)/pp-ff*u;
+        return rez;
+    }
+
+    public:
+    pLapIntegrandH(){}
+
+    virtual ~pLapIntegrandH(){}
+
+    virtual double QFunction(const mfem::Vector &vparam,const mfem::Vector &uu) override
+    {
+        double rez=MyQFunction<double,const mfem::Vector>(vparam,uu);
+        return rez;
+    }
+
+    virtual ADFType QFunction(const mfem::Vector &vparam, ADFVector& uu) override
+    {
+        ADFType rez=MyQFunction<ADFType,ADFVector>(vparam,uu);
+        return rez;
+    }
+
+    virtual ADSType QFunction(const mfem::Vector &vparam, ADSVector& uu) override
+    {
+        ADSType rez=MyQFunction<ADSType,ADSVector>(vparam,uu);
+        return rez;
+    }
+
+};
+
+
+class pLaplaceAD: public mfem::NonlinearFormIntegrator
+{
+protected:
+    mfem::Coefficient* pp;
+    mfem::Coefficient* coeff;
+    mfem::Coefficient* load;
+
+    pLapIntegrandJ qint;
+public:
+    pLaplaceAD()
+    {
+        coeff=nullptr;
+        pp=nullptr;
+    }
+
+    pLaplaceAD(mfem::Coefficient& pp_):pp(&pp_), coeff(nullptr), load(nullptr)
+    {
+
+    }
+
+    pLaplaceAD(mfem::Coefficient &pp_,mfem::Coefficient& q, mfem::Coefficient& ld_): pp(&pp_), coeff(&q), load(&ld_)
+    {
+
+    }
+
+    virtual ~pLaplaceAD()
+    {
+
+    }
+
+    virtual double GetElementEnergy(const mfem::FiniteElement &el, mfem::ElementTransformation &trans, const mfem::Vector &elfun) override
+    {
+        double energy=0.0;
+        int ndof = el.GetDof();
+        int ndim = el.GetDim();
+        int spaceDim = trans.GetSpaceDim();
+        bool square = (ndim == spaceDim);
+        const mfem::IntegrationRule *ir = NULL;
+        int order = 2 * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+        mfem::Vector shapef(ndof);
+        mfem::DenseMatrix dshape_iso(ndof,ndim);
+        mfem::DenseMatrix dshape_xyz(ndof,spaceDim);
+        mfem::Vector grad(spaceDim);
+
+        mfem::Vector vparam(3);//[power, epsilon, load]
+        mfem::Vector uu(4);//[diff_x,diff_y,diff_z,u]
+
+        uu=0.0;
+        vparam[0]=2.0;  //default power
+        vparam[1]=1e-8; //default epsilon
+        vparam[2]=1.0;  //default load
+
+        double w;
+        double detJ;
+
+        for(int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            w = trans.Weight();
+            detJ = (square ? w : w*w);
+            w = ip.weight *w;
+
+            el.CalcDShape(ip,dshape_iso);
+            el.CalcShape(ip,shapef);
+            // AdjugateJacobian = / adj(J),         if J is square
+            //                    \ adj(J^t.J).J^t, otherwise
+            mfem::Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+            // dshape_xyz should be devided by detJ for obtaining the real value
+            // calculate the gradient
+            dshape_xyz.MultTranspose(elfun,grad);
+
+            //set the power
+            if(pp!=nullptr)
+            {
+                vparam[0]=pp->Eval(trans,ip);
+            }
+
+            //set the coefficient ensuring possitiveness of the tangent matrix
+            if(coeff!=nullptr)
+            {
+                vparam[1]=coeff->Eval(trans,ip);
+            }
+            //add the contribution from the load
+            if(load!=nullptr)
+            {
+                vparam[2]=load->Eval(trans,ip);
+            }
+            //fill the values of vector uu
+            for(int jj=0;jj<spaceDim;jj++)
+            {
+                uu[jj]=grad[jj]/detJ;
+            }
+            uu[3]=shapef*elfun;
+
+            energy = energy + w * (qint.QFunction(vparam,uu));
+
+        }
+        return energy;
+    }
+
+    virtual void AssembleElementVector(const mfem::FiniteElement & el,
+                                       mfem::ElementTransformation & trans,
+                                       const mfem::Vector & elfun,
+                                       mfem::Vector & elvect) override
+    {
+            int ndof = el.GetDof();
+            int ndim = el.GetDim();
+            int spaceDim = trans.GetSpaceDim();
+            bool square = (ndim == spaceDim);
+            const mfem::IntegrationRule *ir = NULL;
+            int order = 2 * trans.OrderGrad(&el) - 1; // correct order?
+            ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+            mfem::Vector shapef(ndof);
+            mfem::DenseMatrix dshape_iso(ndof,ndim);
+            mfem::DenseMatrix dshape_xyz(ndof,spaceDim);
+            mfem::Vector lvec(ndof);
+            elvect.SetSize(ndof);
+            elvect=0.0;
+
+            mfem::DenseMatrix B(ndof,4); //[diff_x,diff_y,diff_z, shape]
+            mfem::Vector vparam(3);//[power, epsilon, load]
+            mfem::Vector uu(4);//[diff_x,diff_y,diff_z,u]
+            mfem::Vector du(4);
+            B=0.0;
+            uu=0.0;
+            //initialize the parameters - keep the same order
+            //utilized in the pLapIntegrator definition
+            vparam[0]=2.0;  //default power
+            vparam[1]=1e-8; //default epsilon
+            vparam[2]=1.0;  //default load
+
+            double w;
+            double detJ;
+
+            for (int i = 0; i < ir -> GetNPoints(); i++)
+            {
+                const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+                trans.SetIntPoint(&ip);
+                w = trans.Weight();
+                detJ = (square ? w : w*w);
+                w = ip.weight * w;
+
+                el.CalcDShape(ip,dshape_iso);
+                el.CalcShape(ip,shapef);
+                mfem::Mult(dshape_iso, trans.InverseJacobian(), dshape_xyz);
+
+                //set the matrix B
+                for(int jj=0;jj<spaceDim;jj++)
+                {
+                    B.SetCol(jj,dshape_xyz.GetColumn(jj));
+                }
+                B.SetCol(3,shapef);
+
+
+                //set the power
+                if(pp!=nullptr)
+                {
+                    vparam[0]=pp->Eval(trans,ip);
+                }
+                //set the coefficient ensuring possitiveness of the tangent matrix
+                if(coeff!=nullptr)
+                {
+                    vparam[1]=coeff->Eval(trans,ip);
+                }
+                //add the contribution from the load
+                if(load!=nullptr)
+                {
+                    vparam[2]=load->Eval(trans,ip);
+                }
+
+                //calculate uu
+                B.MultTranspose(elfun,uu);
+                //calculate derivative of the energy with respect to uu
+                qint.QFunctionDU(vparam,uu,du);
+
+                B.Mult(du,lvec);
+                elvect.Add( w, lvec);
+            }// end integration loop
+    }
+
+    virtual void AssembleElementGrad(const mfem::FiniteElement & el,
+                                         mfem::ElementTransformation & trans,
+                                         const mfem::Vector & elfun, mfem::DenseMatrix & elmat) override
+    {
+        int ndof = el.GetDof();
+        int ndim = el.GetDim();
+        int spaceDim = trans.GetSpaceDim();
+        bool square = (ndim == spaceDim);
+        const mfem::IntegrationRule *ir = NULL;
+        int order = 2 * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+        mfem::Vector shapef(ndof);
+        mfem::DenseMatrix dshape_iso(ndof,ndim);
+        mfem::DenseMatrix dshape_xyz(ndof,spaceDim);
+        elmat.SetSize(ndof,ndof);
+        elmat=0.0;
+
+        mfem::DenseMatrix B(ndof,4); //[diff_x,diff_y,diff_z, shape]
+        mfem::DenseMatrix A(ndof,4);
+        mfem::Vector vparam(3);//[power, epsilon, load]
+        mfem::Vector uu(4);//[diff_x,diff_y,diff_z,u]
+        mfem::DenseMatrix duu(4,4);
+        B=0.0;
+        uu=0.0;
+        //initialize the parameters - keep the same order
+        //utilized in the pLapIntegrator definition
+        vparam[0]=2.0;  //default power
+        vparam[1]=1e-8; //default epsilon
+        vparam[2]=1.0;  //default load
+
+        double w;
+        double detJ;
+
+        for (int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            w = trans.Weight();
+            detJ = (square ? w : w*w);
+            w = ip.weight * w;
+
+            el.CalcDShape(ip,dshape_iso);
+            el.CalcShape(ip,shapef);
+            mfem::Mult(dshape_iso, trans.InverseJacobian(), dshape_xyz);
+
+            //set the matrix B
+            for(int jj=0;jj<spaceDim;jj++)
+            {
+                B.SetCol(jj,dshape_xyz.GetColumn(jj));
+            }
+            B.SetCol(3,shapef);
+
+
+            //set the power
+            if(pp!=nullptr)
+            {
+                vparam[0]=pp->Eval(trans,ip);
+            }
+            //set the coefficient ensuring possitiveness of the tangent matrix
+            if(coeff!=nullptr)
+            {
+                vparam[1]=coeff->Eval(trans,ip);
+            }
+            //add the contribution from the load
+            if(load!=nullptr)
+            {
+                vparam[2]=load->Eval(trans,ip);
+            }
+
+            //calculate uu
+            B.MultTranspose(elfun,uu);
+            //calculate derivative of the energy with respect to uu
+            qint.QFunctionDD(vparam,uu,duu);
+
+            mfem::Mult(B,duu,A);
+            mfem::AddMult_a_ABt(w,A,B,elmat);
+
+        }//end integration loop
+    }
+
+
+};
+
+
+class pLaplace: public mfem::NonlinearFormIntegrator
+{
+protected:
+    mfem::Coefficient* pp;
+    mfem::Coefficient* coeff;
+    mfem::Coefficient* load;
+public:
+    pLaplace()
+    {
+        coeff=nullptr;
+        pp=nullptr;
+    }
+
+    pLaplace(mfem::Coefficient& pp_):pp(&pp_), coeff(nullptr), load(nullptr)
+    {
+
+    }
+
+    pLaplace(mfem::Coefficient &pp_,mfem::Coefficient& q, mfem::Coefficient& ld_): pp(&pp_), coeff(&q), load(&ld_)
+    {
+
+    }
+
+    virtual ~pLaplace()
+    {
+
+    }
+
+    virtual double GetElementEnergy(const mfem::FiniteElement &el, mfem::ElementTransformation &trans, const mfem::Vector &elfun) override
+    {
+        double energy=0.0;
+        int ndof = el.GetDof();
+        int ndim = el.GetDim();
+        int spaceDim = trans.GetSpaceDim();
+        bool square = (ndim == spaceDim);
+        const mfem::IntegrationRule *ir = NULL;
+        int order = 2 * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+        mfem::Vector shapef(ndof);
+        mfem::DenseMatrix dshape_iso(ndof,ndim);
+        mfem::DenseMatrix dshape_xyz(ndof,spaceDim);
+        mfem::Vector grad(spaceDim);
+
+        double w;
+        double detJ;
+        double nrgrad2;
+        double ppp=2.0;
+        double eee=0.0;
+
+        for(int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            w = trans.Weight();
+            detJ = (square ? w : w*w);
+            w = ip.weight *w;
+
+            el.CalcDShape(ip,dshape_iso);
+            el.CalcShape(ip,shapef);
+            // AdjugateJacobian = / adj(J),         if J is square
+            //                    \ adj(J^t.J).J^t, otherwise
+            mfem::Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+            // dshape_xyz should be devided by detJ for obtaining the real value
+            // calculate the gradient
+            dshape_xyz.MultTranspose(elfun,grad);
+            nrgrad2=grad*grad/(detJ*detJ);
+
+            //set the power
+            if(pp!=nullptr)
+            {
+                ppp=pp->Eval(trans,ip);
+            }
+
+            //set the coefficient ensuring possitiveness of the tangent matrix
+            if(coeff!=nullptr)
+            {
+                eee=coeff->Eval(trans,ip);
+            }
+
+            energy = energy + w * std::pow( nrgrad2 + eee * eee , ppp / 2.0 ) / ppp;
+
+            //add the contribution from the load
+            if(load!=nullptr)
+            {
+                energy = energy - w *  (shapef*elfun) * load->Eval(trans,ip);
+            }
+        }
+        return energy;
+    }
+
+    virtual void  AssembleElementVector(const mfem::FiniteElement & el,
+                                                mfem::ElementTransformation & trans,
+                                                const mfem::Vector & elfun,
+                                                mfem::Vector & elvect) override
+    {
+            int ndof = el.GetDof();
+            int ndim = el.GetDim();
+            int spaceDim = trans.GetSpaceDim();
+            bool square = (ndim == spaceDim);
+            const mfem::IntegrationRule *ir = NULL;
+            int order = 2 * trans.OrderGrad(&el) - 1; // correct order?
+            ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+            mfem::Vector shapef(ndof);
+            mfem::DenseMatrix dshape_iso(ndof,ndim);
+            mfem::DenseMatrix dshape_xyz(ndof,spaceDim);
+            mfem::Vector grad(spaceDim);
+            mfem::Vector lvec(ndof);
+            elvect.SetSize(ndof);
+            elvect=0.0;
+
+            double w;
+            double detJ;
+            double nrgrad;
+            double aa;
+            double ppp=2.0;
+            double eee=0.0;
+
+            for (int i = 0; i < ir -> GetNPoints(); i++)
+            {
+                const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+                trans.SetIntPoint(&ip);
+                w = trans.Weight();
+                detJ = (square ? w : w*w);
+                w = ip.weight * w;//w;
+
+                el.CalcDShape(ip,dshape_iso);
+                el.CalcShape(ip,shapef);
+                // AdjugateJacobian = / adj(J),         if J is square
+                //                    \ adj(J^t.J).J^t, otherwise
+                mfem::Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+                // dshape_xyz should be devided by detJ for obtaining the real value
+
+                //calculate the gradient
+                dshape_xyz.MultTranspose(elfun,grad);
+                nrgrad=grad.Norml2()/detJ;
+                //grad is not scaled so far, i.e., grad=grad/detJ
+
+                //set the power
+                if(pp!=nullptr)
+                {
+                    ppp=pp->Eval(trans,ip);
+                }
+
+                //set the coefficient ensuring possitiveness of the tangent matrix
+                if(coeff!=nullptr)
+                {
+                    eee=coeff->Eval(trans,ip);
+                }
+
+                aa = nrgrad * nrgrad + eee * eee;
+                aa=std::pow( aa , ( ppp - 2.0 ) / 2.0 );
+                dshape_xyz.Mult(grad,lvec);
+                elvect.Add( w * aa / ( detJ * detJ ), lvec);
+
+
+                //add loading
+                if(load!=nullptr)
+                {
+                    elvect.Add(-w*load->Eval(trans,ip),shapef);
+                }
+            }// end integration loop
+    }
+
+    virtual void AssembleElementGrad(const mfem::FiniteElement & el,
+                                         mfem::ElementTransformation & trans,
+                                         const mfem::Vector & elfun, mfem::DenseMatrix & elmat) override
+    {
+        int ndof = el.GetDof();
+        int ndim = el.GetDim();
+        int spaceDim = trans.GetSpaceDim();
+        bool square = (ndim == spaceDim);
+        const mfem::IntegrationRule *ir = NULL;
+        int order = 2 * trans.OrderGrad(&el) - 1; // correct order?
+        ir = &mfem::IntRules.Get(el.GetGeomType(), order);
+
+        mfem::DenseMatrix dshape_iso(ndof,ndim);
+        mfem::DenseMatrix dshape_xyz(ndof,spaceDim);
+        mfem::Vector grad(spaceDim);
+        mfem::Vector lvec(ndof);
+        elmat.SetSize(ndof,ndof);
+        elmat=0.0;
+
+        double w;
+        double detJ;
+        double nrgrad;
+        double aa0;
+        double aa1;
+        double ppp=2.0;
+        double eee=0.0;
+
+        for (int i = 0; i < ir -> GetNPoints(); i++)
+        {
+            const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+            trans.SetIntPoint(&ip);
+            w = trans.Weight();
+            detJ = (square ? w : w*w);
+            w = ip.weight * w;
+
+            el.CalcDShape(ip,dshape_iso);
+            // AdjugateJacobian = / adj(J),         if J is square
+            //                    \ adj(J^t.J).J^t, otherwise
+            mfem::Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+            // dshape_xyz should be devided by detJ for obtaining the real value
+            // grad is not scaled so far,i.e., grad=grad/detJ
+
+            //set the power
+            if(pp!=nullptr)
+            {
+                ppp=pp->Eval(trans,ip);
+            }
+            //set the coefficient ensuring possitiveness of the tangent matrix
+            if(coeff!=nullptr)
+            {
+                eee=coeff->Eval(trans,ip);
+            }
+
+            //calculate the gradient
+            dshape_xyz.MultTranspose(elfun,grad);
+            nrgrad = grad.Norml2() / detJ;
+            aa0 = nrgrad * nrgrad + eee * eee;
+            aa1 = std::pow( aa0 , ( ppp - 2.0 ) / 2.0 );
+            aa0 = ( ppp - 2.0 ) * std::pow(aa0, ( ppp - 4.0 ) / 2.0 );
+            dshape_xyz.Mult(grad,lvec);
+            w = w / ( detJ * detJ );
+            mfem::AddMult_a_VVt( w * aa0 / ( detJ * detJ ), lvec, elmat);
+            mfem::AddMult_a_AAt( w * aa1 , dshape_xyz, elmat);
+
+        }//end integration loop
+    }
+};
+
+}
+#endif

examples/ex71p.cpp

@@ -0,0 +1,378 @@
+//                       MFEM Example 71 - Parallel Version
+//
+// Compile with: make ex71p
+//
+// Sample runs:
+//   mpirun -np 2 ex71p -m ../data/beam-quad.mesh
+//   mpirun -np 2 ex71p -m ../data/beam-tri.mesh
+//   mpirun -np 2 ex71p -m ../data/beam-hex.mesh
+//   mpirun -np 2 ex71p -m ../data/beam-tet.mesh
+//   mpirun -np 2 ex71p -m ../data/beam-wedge.mesh
+//
+// Description:  This examples solves a quasi-static nonlinear
+//               pLaplacian problem with zero Dirichlet boundary
+//				 conditions applied on all defined boundaries
+//
+//               The example demonstrates the use of nonlinear operators
+//				 combined with automatic differentiation (AD). The definitions
+//				 of the integrators are written in the ex71.hpp.
+//				 Selecting integrator=0 will use handcoded integrator.
+//				 Selecting integrator=1 will utilize AD integrator.
+//				 The AD integrator can be modifief to use ADQFunctionJ
+//				 or ADQFunctionH by overwritting the class type of qint,
+//				 i.e., pLapIntegrandJ or pLapIntegrandH.
+//
+//				 qint (the integrand) is a function which is evaluated
+//				 at every integration point. For implementations utilizing
+//				 ADQFunctionJ, the user has to implement the function and the
+//				 residual evaluation - all virtual methods. The Jacobian of
+//				 the residual is evaluated using AD
+//
+//				 For implementations utilizing ADQFunctionH, the user has
+//				 to implement only the function evaluation (preferebaly as
+//				 a template) and the first derivative (the residual) and the
+//				 second derivatives (the Hessian) are evaluated using AD.
+//
+//               We recommend viewing examples 1 and 19, before viewing this
+//               example.
+
+
+#include "ex71.hpp"
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI
+   int num_procs, myrank;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+
+   // 2. Parse command-line options
+   const char *mesh_file = "../data/beam-tet.mesh";
+   int ser_ref_levels = 0;
+   int par_ref_levels = 0;
+   int order = 2;
+   bool visualization = true;
+   double newton_rel_tol = 1e-4;
+   double newton_abs_tol = 1e-6;
+   int newton_iter = 500;
+   int print_level = 0;
+   double pp = 2.0;
+   int integrator=0;
+   mfem::StopWatch* timer=new mfem::StopWatch();
+
+   mfem::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(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&newton_rel_tol, "-rel", "--relative-tolerance",
+                  "Relative tolerance for the Newton solve.");
+   args.AddOption(&newton_abs_tol, "-abs", "--absolute-tolerance",
+                  "Absolute tolerance for the Newton solve.");
+   args.AddOption(&newton_iter, "-it", "--newton-iterations",
+                  "Maximum iterations for the Newton solve.");
+   args.AddOption(&pp, "-pp", "--power-parameter",
+                  "Power parameter (>=2.0) for the p-Laplacian.");
+   args.AddOption((&print_level),"-prt","--print-level",
+                  "Print level.");
+   args.AddOption(&integrator, "-int","--integrator",
+                  "Integrator 0: standard; 1: AD uaing energy; 2: AD using gradients");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myrank == 0)
+      {
+         args.PrintUsage(std::cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myrank == 0)
+   {
+      args.PrintOptions(std::cout);
+   }
+
+   // 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.
+   mfem::Mesh *mesh = new mfem::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.
+   mfem::ParMesh *pmesh = new mfem::ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // 6. Define the power parameter for the p-Laplacian and all other
+   //    coefficients
+   mfem::ConstantCoefficient c_pp(pp);
+   mfem::ConstantCoefficient load(1.000000000);
+   mfem::ConstantCoefficient c_ee(0.000000001);
+
+   // 7. Define the finite element spaces for the solution
+   mfem::H1_FECollection fec(order,dim);
+   mfem::ParFiniteElementSpace fespace(pmesh,&fec,1,mfem::Ordering::byVDIM);
+   HYPRE_Int glob_size=fespace.GlobalTrueVSize();
+   if (myrank == 0)
+   {
+        std::cout << "Number of finite element unknowns: " << glob_size << std::endl;
+   }
+
+   // 8. Define the Dirichlet conditions
+   mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+   ess_bdr = 1;
+
+   // 9. Define the nonlinear form
+   mfem::ParNonlinearForm* nf=new mfem::ParNonlinearForm(&fespace);
+
+   // 10. Define the solution vector x as a parallel finite element grid function
+   //     corresponding to fespace. Initialize x with initial guess of zero,
+   //     which satisfies the boundary conditions.
+   mfem::ParGridFunction x(&fespace);
+   x = 0.0;
+   mfem::HypreParVector* tv=x.GetTrueDofs();
+   mfem::HypreParVector* sv=x.GetTrueDofs();
+
+   // 11. Define ParaView DataCollection
+   mfem::ParaViewDataCollection *dacol=new mfem::ParaViewDataCollection("pLap",pmesh);
+   dacol->SetLevelsOfDetail(order);
+   dacol->RegisterField("sol",&x);
+
+
+   // 11. Set domain integrators - start with linear diffusion
+   {
+       // the default power coefficient is 2.0
+       mfem::ConstantCoefficient lpp(2.0);
+       if(integrator==0)
+       {
+            nf->AddDomainIntegrator(new mfem::pLaplace(lpp,c_ee,load));
+       }else
+       if(integrator==1)
+       {
+           nf->AddDomainIntegrator(new mfem::pLaplaceAD(lpp,c_ee,load));
+       }
+       nf->SetEssentialBC(ess_bdr);
+       // compute the energy
+       double energy=nf->GetEnergy(*tv);
+       if(myrank==0){
+           std::cout<<"[2] The total energy of the system is E="<<energy<<std::endl;}
+       // time the assembly
+       timer->Clear();
+       timer->Start();
+       mfem::Operator &op=nf->GetGradient(*sv);
+       timer->Stop();
+       if(myrank==0){
+           std::cout<<"[2] The assembly time is: "<<timer->RealTime()<<std::endl;}
+       mfem::Solver *prec=new mfem::HypreBoomerAMG();
+       mfem::GMRESSolver *j_gmres = new mfem::GMRESSolver(MPI_COMM_WORLD);
+       j_gmres->SetRelTol(1e-7);
+       j_gmres->SetAbsTol(1e-15);
+       j_gmres->SetMaxIter(300);
+       j_gmres->SetPrintLevel(print_level);
+       j_gmres->SetPreconditioner(*prec);
+
+       mfem::NewtonSolver* ns;
+       ns=new mfem::NewtonSolver(MPI_COMM_WORLD);
+       ns->iterative_mode = true;
+       ns->SetSolver(*j_gmres);
+       ns->SetOperator(*nf);
+       ns->SetPrintLevel(print_level);
+       ns->SetRelTol(1e-6);
+       ns->SetAbsTol(1e-12);
+       ns->SetMaxIter(3);
+       //solve the problem
+       timer->Clear();
+       timer->Start();
+       ns->Mult(*tv, *sv);
+       timer->Stop();
+       if(myrank==0){
+           std::cout<<"Time for the NewtonSolver: "<<timer->RealTime()<<std::endl;}
+
+       energy=nf->GetEnergy(*sv);
+       if(myrank==0){
+                 std::cout<<"[pp=2] The total energy of the system is E="<<energy<<std::endl;}
+
+       delete ns;
+       delete j_gmres;
+       delete prec;
+
+       x.SetFromTrueDofs(*sv);
+       dacol->SetTime(2.0);
+       dacol->SetCycle(2);
+       dacol->Save();
+   }
+
+   // 12. Continue with powers higher than 2
+   for(int i=3;i<pp;i++)
+   {
+       delete nf;
+       nf=new mfem::ParNonlinearForm(&fespace);
+       mfem::ConstantCoefficient lpp((double)i);
+       if(integrator==0)
+       {
+            nf->AddDomainIntegrator(new mfem::pLaplace(lpp,c_ee,load));
+       }else
+       if(integrator==1)
+       {
+           nf->AddDomainIntegrator(new mfem::pLaplaceAD(lpp,c_ee,load));
+       }
+       nf->SetEssentialBC(ess_bdr);
+       // compute the energy
+       double energy=nf->GetEnergy(*sv);
+       if(myrank==0){
+           std::cout<<"[pp="<<i<<"] The total energy of the system is E="<<energy<<std::endl;}
+       // time the assembly
+       timer->Clear();
+       timer->Start();
+       mfem::Operator &op=nf->GetGradient(*sv);
+       timer->Stop();
+       if(myrank==0){
+           std::cout<<"[pp="<<i<<"] The assembly time is: "<<timer->RealTime()<<std::endl;}
+       mfem::Solver *prec=new mfem::HypreBoomerAMG();
+       mfem::GMRESSolver *j_gmres = new mfem::GMRESSolver(MPI_COMM_WORLD);
+       j_gmres->SetRelTol(1e-7);
+       j_gmres->SetAbsTol(1e-15);
+       j_gmres->SetMaxIter(300);
+       j_gmres->SetPrintLevel(print_level);
+       j_gmres->SetPreconditioner(*prec);
+
+       mfem::NewtonSolver* ns;
+       ns=new mfem::NewtonSolver(MPI_COMM_WORLD);
+       ns->iterative_mode = true;
+       ns->SetSolver(*j_gmres);
+       ns->SetOperator(*nf);
+       ns->SetPrintLevel(print_level);
+       ns->SetRelTol(1e-6);
+       ns->SetAbsTol(1e-12);
+       ns->SetMaxIter(3);
+       //solve the problem
+       timer->Clear();
+       timer->Start();
+       ns->Mult(*tv, *sv);
+       timer->Stop();
+       if(myrank==0){
+           std::cout<<"Time for the NewtonSolver: "<<timer->RealTime()<<std::endl;}
+
+       energy=nf->GetEnergy(*sv);
+       if(myrank==0){
+                 std::cout<<"[pp="<<i<<"] The total energy of the system is E="<<energy<<std::endl;}
+
+       delete ns;
+       delete j_gmres;
+       delete prec;
+
+       x.SetFromTrueDofs(*sv);
+       dacol->SetTime(i);
+       dacol->SetCycle(i);
+       dacol->Save();
+   }
+
+   // 13. Continue with the final power
+   if( std::abs(pp-2.0) > std::numeric_limits<double>::epsilon())
+   {
+       delete nf;
+       nf=new mfem::ParNonlinearForm(&fespace);
+       if(integrator==0)
+       {
+            nf->AddDomainIntegrator(new mfem::pLaplace(c_pp,c_ee,load));
+       }else
+       if(integrator==1)
+       {
+           nf->AddDomainIntegrator(new mfem::pLaplaceAD(c_pp,c_ee,load));
+       }
+       nf->SetEssentialBC(ess_bdr);
+       // compute the energy
+       double energy=nf->GetEnergy(*sv);
+       if(myrank==0){
+           std::cout<<"[pp="<<pp<<"] The total energy of the system is E="<<energy<<std::endl;}
+       // time the assembly
+       timer->Clear();
+       timer->Start();
+       mfem::Operator &op=nf->GetGradient(*sv);
+       timer->Stop();
+       if(myrank==0){
+           std::cout<<"[pp="<<pp<<"] The assembly time is: "<<timer->RealTime()<<std::endl;}
+       mfem::Solver *prec=new mfem::HypreBoomerAMG();
+       mfem::GMRESSolver *j_gmres = new mfem::GMRESSolver(MPI_COMM_WORLD);
+       j_gmres->SetRelTol(1e-8);
+       j_gmres->SetAbsTol(1e-15);
+       j_gmres->SetMaxIter(300);
+       j_gmres->SetPrintLevel(print_level);
+       j_gmres->SetPreconditioner(*prec);
+
+       mfem::NewtonSolver* ns;
+       ns=new mfem::NewtonSolver(MPI_COMM_WORLD);
+       ns->iterative_mode = true;
+       ns->SetSolver(*j_gmres);
+       ns->SetOperator(*nf);
+       ns->SetPrintLevel(print_level);
+       ns->SetRelTol(1e-6);
+       ns->SetAbsTol(1e-12);
+       ns->SetMaxIter(3);
+       //solve the problem
+       timer->Clear();
+       timer->Start();
+       ns->Mult(*tv, *sv);
+       timer->Stop();
+       if(myrank==0){
+           std::cout<<"Time for the NewtonSolver: "<<timer->RealTime()<<std::endl;}
+
+       energy=nf->GetEnergy(*sv);
+       if(myrank==0){
+                 std::cout<<"[pp="<<pp<<"] The total energy of the system is E="<<energy<<std::endl;}
+
+       delete ns;
+       delete j_gmres;
+       delete prec;
+
+       x.SetFromTrueDofs(*sv);
+       dacol->SetTime(pp);
+       if(pp<2.0)
+       {
+           dacol->SetCycle(std::floor(pp));
+       }
+       else
+       {
+           dacol->SetCycle(std::ceil(pp));
+       }
+       dacol->Save();
+   }
+
+
+
+   // 19. Free the used memory
+   delete dacol;
+   delete sv;
+   delete tv;
+   delete nf;
+   delete pmesh;
+   delete timer;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+

fem/CMakeLists.txt

@@ -56,6 +56,7 @@ set(SRCS
   tmop.cpp
   tmop_tools.cpp
   gslib.cpp
+  adnonlininteg.cpp
   transfer.cpp
   )
 
@@ -98,6 +99,7 @@ set(HDRS
   tmop.hpp
   tmop_tools.hpp
   gslib.hpp
+  adnonlininteg.hpp
   transfer.hpp
   )
 

fem/adnonlininteg.cpp

@@ -0,0 +1,408 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "fem.hpp"
+#include "../general/forall.hpp"
+#include "adnonlininteg.hpp"
+
+namespace mfem
+{
+
+
+void ADQFunctionJ::QFunctionDD(const Vector &vparam, const Vector &uu, DenseMatrix &jac)
+{
+#if defined MFEM_USE_ADEPT
+    //use ADEPT package
+    adept::Stack* p_stack=adept::active_stack();
+    p_stack->deactivate();
+
+    int n=uu.Size();
+    jac.SetSize(m,n);
+    jac=0.0;
+    m_stack.activate();
+    {
+        ADFVector aduu(uu);
+        ADFVector rr(m); //residual vector
+        m_stack.new_recording();
+        this->QFunctionDU(vparam,aduu,rr);
+        m_stack.independent(aduu.GetData(), n);//independent variables
+        m_stack.dependent(rr.GetData(), m);//dependent variables
+        m_stack.jacobian(jac.Data());
+    }
+    m_stack.deactivate();
+#elif defined MFEM_USE_CODIPACK
+#if defined MFEM_USE_ADFORWARD
+    //use CoDipack
+    int n=uu.Size();
+    jac.SetSize(m,n);
+    jac=0.0;
+    {
+        ADFVector aduu(n);
+        ADFVector rr(m);
+        for(int i=0;i<n;i++)
+        {
+            aduu[i]=uu[i];
+            aduu[i].setGradient(0.0);
+        }
+
+        for(int ii=0;ii<n;ii++){
+            aduu[ii].setGradient(1.0);
+            this->QFunctionDU(vparam,aduu,rr);
+            for(int jj=0;jj<m;jj++)
+            {
+                jac(jj,ii)=rr[jj].getGradient();
+            }
+            aduu[ii].setGradient(0.0);
+        }
+
+    }
+#else
+    //use CoDiPack in reverse mode
+    int n=uu.Size();
+    jac.SetSize(m,n);
+    jac=0.0;
+    {
+        ADFVector aduu(n);
+        ADFVector rr(m);
+        for(int i=0;i<n;i++)
+        {
+            aduu[i]=uu[i];
+        }
+
+        ADFType::TapeType& tape= ADFType::getGlobalTape();
+        typename ADFType::TapeType::Position pos=tape.getPosition();
+
+        tape.setActive();
+        for(int ii=0;ii<n;ii++){ tape.registerInput(aduu[ii]); }
+        this->QFunctionDU(vparam,aduu,rr);
+        for(int ii=0;ii<m;ii++){ tape.registerOutput(rr[ii]); }
+        tape.setPassive();
+
+        for(int jj=0;jj<m;jj++){
+            rr[jj].setGradient(1.0);
+            tape.evaluate();
+            for(int ii=0;ii<n;ii++){
+                jac(jj,ii)=aduu[ii].getGradient();
+            }
+            rr[jj].setGradient(0.0);
+        }
+        tape.reset(pos);
+    }
+#endif
+#elif defined MFEM_USE_FADBADPP
+    //use FADBAD++
+    #ifdef MFEM_USE_ADFORWARD
+        int n=uu.Size();
+        jac.SetSize(m,n);
+        jac=0.0;
+        {
+            ADFVector aduu(uu);
+            ADFVector rr(m);
+
+            for(int ii=0;ii<n;ii++){
+                aduu[ii].diff(ii,n);
+            }
+            this->QFunctionDU(vparam,aduu,rr);
+            for(int ii=0;ii<n;ii++){
+                for(int jj=0;jj<m;jj++)
+                {
+                    jac(jj,ii)=rr[jj].d(ii);
+                }
+            }
+        }
+    #else
+        int n=uu.Size();
+        jac.SetSize(m,n);
+        jac=0.0;
+        {
+            ADFVector aduu(uu);
+            ADFVector rr(m);
+            this->QFunctionDU(vparam,aduu,rr);
+            for(int ii=0;ii<m;ii++)
+            {
+                rr[ii].diff(ii,m);
+            }
+            for(int ii=0;ii<n;ii++){
+                for(int jj=0;jj<m;jj++)
+                {
+                    jac(jj,ii)=aduu[ii].d(jj);
+                }
+            }
+
+        }
+    #endif
+#else
+    //use native AD package
+   int n=uu.Size();
+   jac.SetSize(m,n);
+   jac=0.0;
+   {
+       ADFVector aduu(uu); //all dual numbers are initialized to zero
+       ADFVector rr(m);
+
+       for(int ii=0;ii<n;ii++){
+           aduu[ii].dual(1.0);
+           this->QFunctionDU(vparam,aduu,rr);
+           for(int jj=0;jj<m;jj++)
+           {
+               jac(jj,ii)=rr[jj].dual();
+           }
+           aduu[ii].dual(0.0);
+       }
+   }
+#endif
+}
+
+void ADQFunctionH::QFunctionDU(const Vector &vparam, Vector &uu, Vector &rr)
+{
+
+#if defined MFEM_USE_CODIPACK
+    int n=uu.Size();
+    rr.SetSize(n);
+    ADFVector aduu(n);
+    ADFType   rez;
+    for(int ii=0;ii<n;ii++)
+    {
+        aduu[ii].setValue(uu[ii]);
+        aduu[ii].setGradient(0.0);
+    }
+    for(int ii=0;ii<n;ii++)
+    {
+        aduu[ii].setGradient(1.0);
+        rez=this->QFunction(vparam,aduu);
+        rr[ii]=rez.getGradient();
+        aduu[ii].setGradient(0.0);
+    }
+#elif defined MFEM_USE_FADBADPP
+    int n=uu.Size();
+    rr.SetSize(n);
+    ADFVector aduu(uu);
+    ADFType   rez;
+    rez=this->QFunction(vparam,aduu);
+    rez.diff(0,1);
+    for(int ii=0;ii<n;ii++)
+    {
+        rr[ii]=aduu[ii].d(0);
+    }
+#else
+    int n=uu.Size();
+    rr.SetSize(n);
+    ADFVector aduu(uu);
+    ADFType   rez;
+    for(int ii=0;ii<n;ii++)
+    {
+        aduu[ii].dual(1.0);
+        rez=this->QFunction(vparam,aduu);
+        rr[ii]=rez.dual();
+        aduu[ii].dual(0.0);
+    }
+#endif
+}
+
+void ADQFunctionH::QFunctionDD(const Vector &vparam, const Vector &uu, DenseMatrix &jac)
+{
+
+#if defined MFEM_USE_CODIPACK
+#if defined MFEM_USE_ADFORWARD
+    int n=uu.Size();
+    jac.SetSize(n);
+    jac=0.0;
+    {
+        ADSVector aduu(n);
+        for(int ii = 0;  ii < n ; ii++)
+        {
+            aduu[ii].value().value()=uu[ii];
+            aduu[ii].value().gradient()=0.0;
+            aduu[ii].gradient().value()=0.0;
+            aduu[ii].gradient().gradient()=0.0;
+        }
+
+        for(int ii = 0; ii < n ; ii++)
+        {
+            aduu[ii].value().gradient()=1.0;
+            for(int jj=0; jj<(ii+1); jj++)
+            {
+                aduu[ii].gradient().value()=1.0;
+                ADSType rez= this->QFunction(vparam,aduu);
+                jac(ii,jj)=rez.gradient().gradient();
+                jac(jj,ii)=jac(ii,jj);
+                aduu[jj].gradient().value()=0.0;
+            }
+            aduu[ii].value().gradient()=0.0;
+        }
+    }
+#else
+    int n=uu.Size();
+    jac.SetSize(n);
+    jac=0.0;
+    {
+        ADSVector aduu(n);
+        for(int ii=0;ii < n ; ii++)
+        {
+            aduu[ii].value().value()=uu[ii];
+        }
+
+        ADSType rez;
+        ADSType::TapeType& tape = ADSType::getGlobalTape();
+        typename ADSType::TapeType::Position pos;
+
+        for(int ii = 0; ii < n ; ii++)
+        {
+            pos=tape.getPosition();
+            tape.setActive();
+
+            for(int jj=0;jj < n; jj++) {
+                if(jj==ii) {aduu[jj].value().gradient()=1.0;}
+                else {aduu[jj].value().gradient()=0.0;}
+                tape.registerInput(aduu[jj]);
+            }
+
+            rez=this->QFunction(vparam,aduu);
+            tape.registerOutput(rez);
+            tape.setPassive();
+
+            rez.gradient().value()=1.0;
+            tape.evaluate();
+
+            for(int jj=0; jj<(ii+1); jj++)
+            {
+                jac(ii,jj)=aduu[jj].gradient().gradient();
+                jac(jj,ii)=jac(ii,jj);
+            }
+            tape.reset(pos);
+        }
+    }
+#endif
+#elif defined MFEM_USE_FADBADPP
+    int n=uu.Size();
+    jac.SetSize(n);
+    jac=0.0;
+    {
+        ADSVector aduu(n);
+        for(int ii = 0;  ii < n ; ii++)
+        {
+            aduu[ii]=uu[ii];
+            aduu[ii].x().diff(ii,n);
+        }
+        ADSType rez= this->QFunction(vparam,aduu);
+        rez.diff(0,1);
+        for(int ii = 0; ii < n ; ii++)
+        {
+            for(int jj=0; jj<ii; jj++)
+            {
+                jac(ii,jj)=aduu[ii].d(0).d(jj);
+                jac(jj,ii)=aduu[jj].d(0).d(ii);
+            }
+            jac(ii,ii)=aduu[ii].d(0).d(ii);
+        }
+    }
+#else
+    int n=uu.Size();
+    jac.SetSize(n);
+    jac=0.0;
+    {
+        ADSVector aduu(n);
+        for(int ii = 0;  ii < n ; ii++)
+        {
+            aduu[ii].real(ADFType(uu[ii],0.0));
+            aduu[ii].dual(ADFType(0.0,0.0));
+        }
+
+        for(int ii = 0; ii < n ; ii++)
+        {
+            aduu[ii].real(ADFType(uu[ii],1.0));
+            for(int jj=0; jj<(ii+1); jj++)
+            {
+                aduu[jj].dual(ADFType(1.0,0.0));
+                ADSType rez= this->QFunction(vparam,aduu);
+                jac(ii,jj)=rez.dual().dual();
+                jac(jj,ii)=rez.dual().dual();
+                aduu[jj].dual(ADFType(0.0,0.0));
+            }
+            aduu[ii].real(ADFType(uu[ii],0.0));
+        }
+    }
+#endif
+}
+
+
+
+
+
+double ADNonlinearFormIntegratorH::GetElementEnergy(const mfem::FiniteElement & el,
+                                    mfem::ElementTransformation & Tr, 
+                                    const mfem::Vector & elfun)
+{
+    return this->ElementEnergy(el,Tr,elfun);
+}
+
+void ADNonlinearFormIntegratorH::AssembleElementVector(const mfem::FiniteElement & el,
+                                            mfem::ElementTransformation & Tr, 
+                                            const mfem::Vector & elfun, mfem::Vector & elvect)
+{
+    
+    int ndof = el.GetDof();
+    elvect.SetSize(ndof);
+    
+    {
+        ADFVector adelfun(elfun);
+        //all dual numbers in adelfun are initialized to 0.0
+        for(int ii = 0; ii < adelfun.Size(); ii++)
+        {
+            //set the dual for the ii^th element to 1.0
+            adelfun[ii].dual(1.0);
+            ADFType rez= this->ElementEnergy(el,Tr, adelfun);
+            elvect[ii]=rez.dual();
+            //return it back to zero
+            adelfun[ii].dual(0.0);
+        }
+    }
+    
+}
+
+void ADNonlinearFormIntegratorH::AssembleElementGrad(const mfem::FiniteElement & el,
+                                          mfem::ElementTransformation & Tr, 
+                                          const mfem::Vector & elfun, 
+                                          mfem::DenseMatrix & elmat)
+{
+    
+    int ndof = el.GetDof();
+    elmat.SetSize(ndof);
+    elmat=0.0;
+    {
+        ADSVector adelfun(ndof);
+        for(int ii = 0; ii < ndof; ii++)
+        {
+            adelfun[ii].real(ADFType(elfun[ii],0.0));
+            adelfun[ii].dual(ADFType(0.0,0.0));
+        }
+        
+        for(int ii = 0; ii < adelfun.Size(); ii++)
+        {
+            adelfun[ii].real(ADFType(elfun[ii],1.0));
+            for(int jj = 0; jj < (ii+1); jj++)
+            {
+                adelfun[jj].dual(ADFType(1.0,0.0));
+                ADSType rez= this->ElementEnergy(el,Tr, adelfun);
+                elmat(ii,jj)=rez.dual().dual();
+                elmat(jj,ii)=rez.dual().dual();
+                adelfun[jj].dual(ADFType(0.0,0.0));
+            }
+            adelfun[ii].real(ADFType(elfun[ii],0.0));
+        }
+        
+    }
+    
+}
+
+    
+} //end namespace mfem
+

fem/adnonlininteg.hpp

@@ -0,0 +1,204 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+
+
+#ifndef MFEM_ADNONLININTEG
+#define MFEM_ADNONLININTEG
+
+#include "../config/config.hpp"
+#include "fe.hpp"
+#include "coefficient.hpp"
+#include "fespace.hpp"
+#include "nonlininteg.hpp"
+#include "../linalg/tadvector.hpp"
+#include "../linalg/taddensemat.hpp"
+#include "../linalg/fdual.hpp"
+
+#if defined MFEM_USE_ADEPT
+#include <adept.h>
+#elif defined  MFEM_USE_CODIPACK
+#include <codi.hpp>
+#elif defined  MFEM_USE_FADBADPP
+#include <fadiff.h>
+#include <badiff.h>
+#endif
+
+//define Forward AD mode
+//#define MFEM_USE_ADFORWARD
+
+namespace mfem
+{
+
+class ADQFunctionJ
+{
+  private:
+    int m;  //dimension of the residual vector
+            //the Jacobian will have dimensions [m,length(uu)]
+  protected:protected:
+#ifdef MFEM_USE_ADEPT
+    adept::Stack  m_stack;
+#endif
+
+  public:
+#if defined MFEM_USE_ADEPT
+    typedef adept::adouble 			   ADFType;
+    typedef TADVector<ADFType>         ADFVector;
+    typedef TADDenseMatrix<ADFType>    ADFDenseMatrix;
+#elif defined MFEM_USE_CODIPACK
+    #if defined MFEM_USE_ADFORWARD
+        typedef codi::RealForward 		   ADFType;
+        typedef TADVector<ADFType>         ADFVector;
+        typedef TADDenseMatrix<ADFType>    ADFDenseMatrix;
+    #else
+        typedef codi::RealRevers 		   ADFType;
+        typedef TADVector<ADFType>         ADFVector;
+        typedef TADDenseMatrix<ADFType>    ADFDenseMatrix;
+    #endif
+#elif defined MFEM_USE_FADBADPP
+    #ifdef MFEM_USE_ADFORWARD
+        typedef fadbad::F<double>		   ADFType;
+        typedef TADVector<ADFType>         ADFVector;
+        typedef TADDenseMatrix<ADFType>    ADFDenseMatrix;
+    #else
+        typedef fadbad::B<double>		   ADFType;
+        typedef TADVector<ADFType>         ADFVector;
+        typedef TADDenseMatrix<ADFType>    ADFDenseMatrix;
+    #endif
+#else
+    typedef mfem::ad::FDual<double>    ADFType;
+    typedef TADVector<ADFType>         ADFVector;
+    typedef TADDenseMatrix<ADFType>    ADFDenseMatrix;
+#endif
+
+#ifdef MFEM_USE_ADEPT
+    ADQFunctionJ(int m_=1):m_stack(false)
+    {
+        m=m_;
+    }
+#else
+    ADQFunctionJ(int m_=1){ m=m_;}
+#endif
+
+    virtual ~ADQFunctionJ(){}
+
+    virtual double QFunction(const mfem::Vector& vparam, const mfem::Vector& uu)=0;
+    virtual void QFunctionDU(const mfem::Vector& vparam, ADFVector& uu, ADFVector& rr)=0;
+    virtual void QFunctionDU(const mfem::Vector& vparam, mfem::Vector& uu, mfem::Vector& rr)=0;
+
+    void QFunctionDD(const mfem::Vector& vparam, const mfem::Vector& uu, mfem::DenseMatrix& jj);
+
+};
+
+
+class ADQFunctionH
+{
+  public:
+#if defined MFEM_USE_CODIPACK
+#if defined MFEM_USE_ADFORWARD
+    //use forward mode for both the first and the second derivatives
+    typedef codi::RealForwardGen<double> ADFType;
+    typedef TADVector<ADFType>        	 ADFVector;
+    typedef TADDenseMatrix<ADFType>    	 ADFDenseMatrix;
+
+    typedef codi::RealForwardGen<ADFType> 	ADSType;
+    typedef TADVector<ADSType>         		ADSVector;
+    typedef TADDenseMatrix<ADSType>    		ADSDenseMatrix;
+#else
+    //use mixed forward and reverse mode
+    typedef codi::RealForwardGen<double> ADFType;
+    typedef TADVector<ADFType>        	 ADFVector;
+    typedef TADDenseMatrix<ADFType>    	 ADFDenseMatrix;
+
+    typedef codi::RealReverseGen<ADFType>   ADSType;
+    typedef TADVector<ADSType>         		ADSVector;
+    typedef TADDenseMatrix<ADSType>    		ADSDenseMatrix;
+#endif
+#elif defined MFEM_USE_FADBADPP
+        typedef fadbad::B<double> 			 ADFType;
+        typedef TADVector<ADFType>        	 ADFVector;
+        typedef TADDenseMatrix<ADFType>    	 ADFDenseMatrix;
+
+        typedef fadbad::B<fadbad::F<double>>    ADSType;
+        typedef TADVector<ADSType>         		ADSVector;
+        typedef TADDenseMatrix<ADSType>    		ADSDenseMatrix;
+#else
+    typedef mfem::ad::FDual<double>    ADFType;
+    typedef TADVector<ADFType>         ADFVector;
+    typedef TADDenseMatrix<ADFType>    ADFDenseMatrix;
+
+    typedef mfem::ad::FDual<ADFType>   ADSType;
+    typedef TADVector<ADSType>         ADSVector;
+    typedef TADDenseMatrix<ADSType>    ADSDenseMatrix;
+#endif
+
+    ADQFunctionH(){}
+
+    virtual ~ADQFunctionH(){}
+
+    virtual double QFunction(const mfem::Vector& vparam, const mfem::Vector& uu)=0;
+    virtual ADFType QFunction(const mfem::Vector& vparam, ADFVector& uu)=0;
+    virtual ADSType QFunction(const mfem::Vector &vparam, ADSVector& uu)=0;
+
+    virtual void QFunctionDU(const mfem::Vector& vparam, mfem::Vector& uu, mfem::Vector& rr);
+    void QFunctionDD(const mfem::Vector& vparam, const mfem::Vector& uu, mfem::DenseMatrix& jj);
+};
+
+
+
+class ADNonlinearFormIntegratorH: public NonlinearFormIntegrator
+{
+public:
+    
+    typedef mfem::ad::FDual<double>    ADFType;
+    typedef TADVector<ADFType>         ADFVector;
+    typedef TADDenseMatrix<ADFType>    ADFDenseMatrix;
+    
+    typedef mfem::ad::FDual<ADFType>   ADSType;
+    typedef TADVector<ADSType>         ADSVector;
+    typedef TADDenseMatrix<ADSType>    ADSDenseMatrix;
+    
+    ADNonlinearFormIntegratorH(){}
+    
+    virtual ~ADNonlinearFormIntegratorH(){}
+    
+    virtual ADSType ElementEnergy(const mfem::FiniteElement & el,
+                                 mfem::ElementTransformation & Tr,
+                                 const ADSVector & elfun)=0;
+                                 
+    virtual ADFType ElementEnergy(const mfem::FiniteElement & el,
+                                 mfem::ElementTransformation & Tr,
+                                 const ADFVector & elfun)=0;
+                                
+    virtual double ElementEnergy(const mfem::FiniteElement & el, 
+                                 mfem::ElementTransformation & Tr,
+                                 const mfem::Vector & elfun)=0;
+                                    
+                                    
+    virtual double GetElementEnergy(const mfem::FiniteElement & el, 
+                                    mfem::ElementTransformation & Tr, 
+                                    const mfem::Vector & elfun) override;
+    
+    virtual void AssembleElementVector(const mfem::FiniteElement & el,
+                                            mfem::ElementTransformation & Tr, 
+                                            const mfem::Vector & elfun, mfem::Vector & elvect) override;
+    
+    virtual void AssembleElementGrad(const mfem::FiniteElement & el, 
+                                          mfem::ElementTransformation & Tr, 
+                                          const mfem::Vector & elfun, 
+                                          mfem::DenseMatrix & elmat) override;
+};
+
+
+}
+
+#endif    
+

fem/bilininteg.cpp

@@ -926,11 +926,14 @@ void BoundaryMassIntegrator::AssembleFaceMatrix(
    for (int i = 0; i < ir->GetNPoints(); i++)
    {
       const IntegrationPoint &ip = ir->IntPoint(i);
-      IntegrationPoint eip;
-      Trans.Loc1.Transform(ip, eip);
+
+      // Set the integration point in the face and the neighboring element
+      Trans.SetAllIntPoints(&ip);
+
+      // Access the neighboring element's integration point
+      const IntegrationPoint &eip = Trans.GetElement1IntPoint();
       el1.CalcShape(eip, shape);
 
-      Trans.SetIntPoint(&ip);
       w = Trans.Weight() * ip.weight;
       if (Q)
       {
@@ -2571,15 +2574,16 @@ void DGTraceIntegrator::AssembleFaceMatrix(const FiniteElement &el1,
    for (int p = 0; p < ir->GetNPoints(); p++)
    {
       const IntegrationPoint &ip = ir->IntPoint(p);
-      IntegrationPoint eip1, eip2;
-      Trans.Loc1.Transform(ip, eip1);
-      if (ndof2)
-      {
-         Trans.Loc2.Transform(ip, eip2);
-      }
-      el1.CalcShape(eip1, shape1);
 
-      Trans.SetIntPoint(&ip);
+      // Set the integration point in the face and the neighboring elements
+      Trans.SetAllIntPoints(&ip);
+
+      // Access the neighboring elements' integration points
+      // Note: eip2 will only contain valid data if Elem2 exists
+      const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
+      const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
+
+      el1.CalcShape(eip1, shape1);
 
       u->Eval(vu, *Trans.Elem1, eip1);
 
@@ -2727,10 +2731,15 @@ void DGDiffusionIntegrator::AssembleFaceMatrix(
    for (int p = 0; p < ir->GetNPoints(); p++)
    {
       const IntegrationPoint &ip = ir->IntPoint(p);
-      IntegrationPoint eip1, eip2;
 
-      Trans.Loc1.Transform(ip, eip1);
-      Trans.SetIntPoint(&ip);
+      // Set the integration point in the face and the neighboring elements
+      Trans.SetAllIntPoints(&ip);
+
+      // Access the neighboring elements' integration points
+      // Note: eip2 will only contain valid data if Elem2 exists
+      const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
+      const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
+
       if (dim == 1)
       {
          nor(0) = 2*eip1.x - 1.0;
@@ -2787,7 +2796,6 @@ void DGDiffusionIntegrator::AssembleFaceMatrix(
 
       if (ndof2)
       {
-         Trans.Loc2.Transform(ip, eip2);
          el2.CalcShape(eip2, shape2);
          el2.CalcDShape(eip2, dshape2);
          w = ip.weight/2/Trans.Elem2->Weight();
@@ -3005,9 +3013,14 @@ void DGElasticityIntegrator::AssembleFaceMatrix(
    for (int pind = 0; pind < ir->GetNPoints(); ++pind)
    {
       const IntegrationPoint &ip = ir->IntPoint(pind);
-      IntegrationPoint eip1, eip2; // integration point in the reference space
-      Trans.Loc1.Transform(ip, eip1);
-      Trans.SetIntPoint(&ip);
+
+      // Set the integration point in the face and the neighboring elements
+      Trans.SetAllIntPoints(&ip);
+
+      // Access the neighboring elements' integration points
+      // Note: eip2 will only contain valid data if Elem2 exists
+      const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
+      const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
 
       el1.CalcShape(eip1, shape1);
       el1.CalcDShape(eip1, dshape1);
@@ -3027,7 +3040,6 @@ void DGElasticityIntegrator::AssembleFaceMatrix(
       double w, wLM;
       if (ndofs2)
       {
-         Trans.Loc2.Transform(ip, eip2);
          el2.CalcShape(eip2, shape2);
          el2.CalcDShape(eip2, dshape2);
          CalcAdjugate(Trans.Elem2->Jacobian(), adjJ);
@@ -3165,17 +3177,22 @@ void TraceJumpIntegrator::AssembleFaceMatrix(
    for (int p = 0; p < ir->GetNPoints(); p++)
    {
       const IntegrationPoint &ip = ir->IntPoint(p);
-      IntegrationPoint eip1, eip2;
+
+      // Set the integration point in the face and the neighboring elements
+      Trans.SetAllIntPoints(&ip);
+
+      // Access the neighboring elements' integration points
+      // Note: eip2 will only contain valid data if Elem2 exists
+      const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
+      const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
+
       // Trace finite element shape function
-      Trans.SetIntPoint(&ip);
       trial_face_fe.CalcShape(ip, face_shape);
       // Side 1 finite element shape function
-      Trans.Loc1.Transform(ip, eip1);
       test_fe1.CalcShape(eip1, shape1);
       if (ndof2)
       {
          // Side 2 finite element shape function
-         Trans.Loc2.Transform(ip, eip2);
          test_fe2.CalcShape(eip2, shape2);
       }
       w = ip.weight;

fem/coefficient.cpp

@@ -209,18 +209,24 @@ void GradientGridFunctionCoefficient::Eval(
    GridFunc->GetGradients(T, ir, M);
 }
 
-CurlGridFunctionCoefficient::CurlGridFunctionCoefficient (
+CurlGridFunctionCoefficient::CurlGridFunctionCoefficient(
    const GridFunction *gf)
-   : VectorCoefficient ((gf) ?
-                        gf -> FESpace() -> GetMesh() -> SpaceDimension() : 0)
+   : VectorCoefficient(0)
 {
-   GridFunc = gf;
+   SetGridFunction(gf);
 }
 
 void CurlGridFunctionCoefficient::SetGridFunction(const GridFunction *gf)
 {
-   GridFunc = gf; vdim = (gf) ?
-                         gf -> FESpace() -> GetMesh() -> SpaceDimension() : 0;
+   if (gf)
+   {
+      int sdim = gf -> FESpace() -> GetMesh() -> SpaceDimension();
+      MFEM_VERIFY(sdim == 2 || sdim == 3,
+                  "CurlGridFunctionCoefficient "
+                  "only defind for spaces of dimension 2 or 3.");
+   }
+   GridFunc = gf;
+   vdim = (gf) ? (2 * gf -> FESpace() -> GetMesh() -> SpaceDimension() - 3) : 0;
 }
 
 void CurlGridFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,

fem/eltrans.cpp

@@ -552,26 +552,32 @@ void IntegrationPointTransformation::Transform (const IntegrationRule &ir1,
    }
 }
 
-void FaceElementTransformations::SetIntPoint(const IntegrationPoint *ip)
+void FaceElementTransformations::SetIntPoint(const IntegrationPoint *face_ip)
 {
-   IsoparametricTransformation::SetIntPoint(ip);
+   IsoparametricTransformation::SetIntPoint(face_ip);
 
-   if (Elem1)
+   if (mask & 4)
    {
-      Loc1.Transform(*ip, eip1);
-      Elem1->SetIntPoint(&eip1);
+      Loc1.Transform(*face_ip, eip1);
+      if (Elem1)
+      {
+         Elem1->SetIntPoint(&eip1);
+      }
    }
-   if (Elem2)
+   if (mask & 8)
    {
-      Loc2.Transform(*ip, eip2);
-      Elem2->SetIntPoint(&eip2);
+      Loc2.Transform(*face_ip, eip2);
+      if (Elem2)
+      {
+         Elem2->SetIntPoint(&eip2);
+      }
    }
 }
 
 ElementTransformation &
 FaceElementTransformations::GetElement1Transformation()
 {
-   MFEM_VERIFY(mask & 1 && Elem1 != NULL, "The ElementTransformation "
+   MFEM_VERIFY(mask & HAVE_ELEM1 && Elem1 != NULL, "The ElementTransformation "
                "for the element has not been configured for side 1.");
    return *Elem1;
 }
@@ -579,7 +585,7 @@ FaceElementTransformations::GetElement1Transformation()
 ElementTransformation &
 FaceElementTransformations::GetElement2Transformation()
 {
-   MFEM_VERIFY(mask & 2 && Elem2 != NULL, "The ElementTransformation "
+   MFEM_VERIFY(mask & HAVE_ELEM2 && Elem2 != NULL, "The ElementTransformation "
                "for the element has not been configured for side 2.");
    return *Elem2;
 }
@@ -587,7 +593,7 @@ FaceElementTransformations::GetElement2Transformation()
 IntegrationPointTransformation &
 FaceElementTransformations::GetIntPoint1Transformation()
 {
-   MFEM_VERIFY(mask & 4, "The IntegrationPointTransformation "
+   MFEM_VERIFY(mask & HAVE_LOC1, "The IntegrationPointTransformation "
                "for the element has not been configured for side 1.");
    return Loc1;
 }
@@ -595,7 +601,7 @@ FaceElementTransformations::GetIntPoint1Transformation()
 IntegrationPointTransformation &
 FaceElementTransformations::GetIntPoint2Transformation()
 {
-   MFEM_VERIFY(mask & 8, "The IntegrationPointTransformation "
+   MFEM_VERIFY(mask & HAVE_LOC2, "The IntegrationPointTransformation "
                "for the element has not been configured for side 2.");
    return Loc2;
 }
@@ -603,7 +609,7 @@ FaceElementTransformations::GetIntPoint2Transformation()
 void FaceElementTransformations::Transform(const IntegrationPoint &ip,
                                            Vector &trans)
 {
-   MFEM_VERIFY(mask & 16, "The ElementTransformation "
+   MFEM_VERIFY(mask & HAVE_FACE, "The ElementTransformation "
                "for the face has not been configured.");
    IsoparametricTransformation::Transform(ip, trans);
 }
@@ -611,7 +617,7 @@ void FaceElementTransformations::Transform(const IntegrationPoint &ip,
 void FaceElementTransformations::Transform(const IntegrationRule &ir,
                                            DenseMatrix &tr)
 {
-   MFEM_VERIFY(mask & 16, "The ElementTransformation "
+   MFEM_VERIFY(mask & HAVE_FACE, "The ElementTransformation "
                "for the face has not been configured.");
    IsoparametricTransformation::Transform(ir, tr);
 }
@@ -619,7 +625,7 @@ void FaceElementTransformations::Transform(const IntegrationRule &ir,
 void FaceElementTransformations::Transform(const DenseMatrix &matrix,
                                            DenseMatrix &result)
 {
-   MFEM_VERIFY(mask & 16, "The ElementTransformation "
+   MFEM_VERIFY(mask & HAVE_FACE, "The ElementTransformation "
                "for the face has not been configured.");
    IsoparametricTransformation::Transform(matrix, result);
 }

fem/eltrans.hpp

@@ -439,15 +439,57 @@ public:
    void Transform (const IntegrationRule  &, IntegrationRule  &);
 };
 
+/** @brief A specialized ElementTransformation class representing a face and
+    its two neighboring elements.
 
+    This class can be used as a container for the element transformation data
+    needed for integrating discontinuous fields on element interfaces in a
+    Discontinuous Galerkin (DG) context.
+
+    The secondary purpose of this class is to enable the
+    GridFunction::GetValue function, and various related functions, to properly
+    evaluate fields with limited continuity on boundary elements.
+*/
 class FaceElementTransformations : public IsoparametricTransformation
 {
 private:
+
+   // Bitwise OR of ConfigMasks
    int mask;
 
    IntegrationPoint eip1, eip2;
 
+protected: // interface for Mesh to be able to configure this object.
+
+   friend class Mesh;
+#ifdef MFEM_USE_MPI
+   friend class ParMesh;
+#endif
+
+   /// Set the mask indicating which portions of the object have been setup
+   /** The argument @a m is a bitmask used in
+       Mesh::GetFaceElementTransformations to indicate which portions of the
+       FaceElementTransformations object have been configured.
+
+       mask &  1: Elem1 is configured
+       mask &  2: Elem2 is configured
+       mask &  4: Loc1 is configured
+       mask &  8: Loc2 is configured
+       mask & 16: The Face transformation itself is configured
+   */
+   void SetConfigurationMask(int m) { mask = m; }
+
 public:
+
+   enum ConfigMasks
+   {
+      HAVE_ELEM1 =  1, ///< Element on side 1 is configured
+      HAVE_ELEM2 =  2, ///< Element on side 2 is configured
+      HAVE_LOC1  =  4, ///< Point transformation for side 1 is configured
+      HAVE_LOC2  =  8, ///< Point transformation for side 2 is configured
+      HAVE_FACE  = 16  ///< Face transformation is configured
+   };
+
    int Elem1No, Elem2No;
    Geometry::Type &FaceGeom; ///< @deprecated Use GetGeometryType instead
    ElementTransformation *Elem1, *Elem2;
@@ -466,10 +508,10 @@ public:
    */
    void SetGeometryType(Geometry::Type g) { geom = g; }
 
-   /// Set the mask indicating which portions of the object have been setup
-   /** The argument @a m is a bitmask used in
-       Mesh::GetFaceElementTransformations to indicate which portions of the
-       FaceElement Transformations object have been configured.
+   /** @brief Return the mask defining the configuration state.
+
+       The mask value indicates which portions of FaceElementTransformations
+       object have been configured.
 
        mask &  1: Elem1 is configured
        mask &  2: Elem2 is configured
@@ -477,12 +519,45 @@ public:
        mask &  8: Loc2 is configured
        mask & 16: The Face transformation itself is configured
    */
-   void SetConfigurationMask(int m) { mask = m; }
-   int  GetConfigurationMask() const { return mask; }
+   int GetConfigurationMask() const { return mask; }
 
    /** @brief Set the integration point in the Face and the two neighboring
-       elements, if present. */
-   void SetIntPoint(const IntegrationPoint *ip);
+       elements, if present.
+
+       The point @a face_ip must be in the reference coordinate system of the
+       face.
+   */
+   void SetIntPoint(const IntegrationPoint *face_ip);
+
+   /** @brief Set the integration point in the Face and the two neighboring
+       elements, if present.
+
+       This is a more expressive member function name than SetIntPoint, which
+       in this special case, does the same thing. This function can be used for
+       greater code clarity.
+   */
+   inline void SetAllIntPoints(const IntegrationPoint *face_ip)
+   { FaceElementTransformations::SetIntPoint(face_ip); }
+
+   /** @brief Get a const reference to the integration point in neighboring
+       element 1 corresponding to the currently set integration point on the
+       face.
+
+       This IntegrationPoint object will only contain up-to-date data if
+       SetIntPoint or SetAllIntPoints has been called with the latest
+       integration point for the face and the appropriate point transformation
+       has been configured. */
+   const IntegrationPoint &GetElement1IntPoint() { return eip1; }
+
+   /** @brief Get a const reference to the integration point in neighboring
+       element 2 corresponding to the currently set integration point on the
+       face.
+
+       This IntegrationPoint object will only contain up-to-date data if
+       SetIntPoint or SetAllIntPoints has been called with the latest
+       integration point for the face and the appropriate point transformation
+       has been configured. */
+   const IntegrationPoint &GetElement2IntPoint() { return eip2; }
 
    virtual void Transform(const IntegrationPoint &, Vector &);
    virtual void Transform(const IntegrationRule &, DenseMatrix &);

fem/fem.hpp

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

fem/gridfunc.cpp

@@ -752,19 +752,21 @@ double GridFunction::GetValue(ElementTransformation &T,
             IntegrationPoint fip;
             be_to_bfe(FET->GetGeometryType(), o, ip, fip);
 
-            FET->SetIntPoint(&fip);
+            // Compute and set the point in element 1 from fip
+            FET->SetAllIntPoints(&fip);
             ElementTransformation & T1 = FET->GetElement1Transformation();
             return GetValue(T1, T1.GetIntPoint(), comp);
          }
-         break;
       }
+      break;
       case ElementTransformation::BDR_FACE:
       {
          FaceElementTransformations * FET =
             dynamic_cast<FaceElementTransformations *>(&T);
 
          // Evaluate in neighboring element for both continuous and
-         // discontinuous fields.
+         // discontinuous fields (the integration point in T1 should have
+         // already been set).
          ElementTransformation & T1 = FET->GetElement1Transformation();
          return GetValue(T1, T1.GetIntPoint(), comp);
       }
@@ -888,19 +890,21 @@ void GridFunction::GetVectorValue(ElementTransformation &T,
             IntegrationPoint fip;
             be_to_bfe(FET->GetGeometryType(), o, ip, fip);
 
-            FET->SetIntPoint(&fip);
+            // Compute and set the point in element 1 from fip
+            FET->SetAllIntPoints(&fip);
             ElementTransformation & T1 = FET->GetElement1Transformation();
             return GetVectorValue(T1, T1.GetIntPoint(), val);
          }
-         break;
       }
+      break;
       case ElementTransformation::BDR_FACE:
       {
          FaceElementTransformations * FET =
             dynamic_cast<FaceElementTransformations *>(&T);
 
          // Evaluate in neighboring element for both continuous and
-         // discontinuous fields.
+         // discontinuous fields (the integration point in T1 should have
+         // already been set).
          ElementTransformation & T1 = FET->GetElement1Transformation();
          return GetVectorValue(T1, T1.GetIntPoint(), val);
       }
@@ -1338,107 +1342,262 @@ void GridFunction::GetVectorGradientHat(
    MultAtB(loc_data_mat, dshape, gh);
 }
 
-double GridFunction::GetDivergence(ElementTransformation &tr) const
+double GridFunction::GetDivergence(ElementTransformation &T) const
 {
-   double div_v;
-   int elNo = tr.ElementNo;
-   const FiniteElement *FElem = fes->GetFE(elNo);
-   if (FElem->GetRangeType() == FiniteElement::SCALAR)
+   switch (T.ElementType)
    {
-      MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
-                  "invalid FE map type");
-      DenseMatrix grad_hat;
-      GetVectorGradientHat(tr, grad_hat);
-      const DenseMatrix &Jinv = tr.InverseJacobian();
-      div_v = 0.0;
-      for (int i = 0; i < Jinv.Width(); i++)
+      case ElementTransformation::ELEMENT:
       {
-         for (int j = 0; j < Jinv.Height(); j++)
+         int elNo = T.ElementNo;
+         const FiniteElement *fe = fes->GetFE(elNo);
+         if (fe->GetRangeType() == FiniteElement::SCALAR)
          {
-            div_v += grad_hat(i, j) * Jinv(j, i);
+            MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE,
+                        "invalid FE map type");
+            DenseMatrix grad_hat;
+            GetVectorGradientHat(T, grad_hat);
+            const DenseMatrix &Jinv = T.InverseJacobian();
+            double div_v = 0.0;
+            for (int i = 0; i < Jinv.Width(); i++)
+            {
+               for (int j = 0; j < Jinv.Height(); j++)
+               {
+                  div_v += grad_hat(i, j) * Jinv(j, i);
+               }
+            }
+            return div_v;
+         }
+         else
+         {
+            // Assuming RT-type space
+            Array<int> dofs;
+            fes->GetElementDofs(elNo, dofs);
+            Vector loc_data, divshape(fe->GetDof());
+            GetSubVector(dofs, loc_data);
+            fe->CalcDivShape(T.GetIntPoint(), divshape);
+            return (loc_data * divshape) / T.Weight();
          }
       }
+      break;
+      case ElementTransformation::BDR_ELEMENT:
+      {
+         // In order to properly capture the derivative of the normal component
+         // of the field (as well as the transverse divergence of the
+         // tangential compoents) we must evaluate it in the neighboring
+         // element.
+         FaceElementTransformations * FET =
+            fes->GetMesh()->GetBdrFaceTransformations(T.ElementNo);
+
+         // Boundary elements and Boundary Faces may have different
+         // orientations so adjust the integration point if necessary.
+         int o = 0;
+         if (fes->GetMesh()->Dimension() == 3)
+         {
+            int f;
+            fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+         }
+
+         IntegrationPoint fip;
+         be_to_bfe(FET->GetGeometryType(), o, T.GetIntPoint(), fip);
+
+         // Compute and set the point in element 1 from fip
+         FET->SetAllIntPoints(&fip);
+         ElementTransformation & T1 = FET->GetElement1Transformation();
+
+         return GetDivergence(T1);
+      }
+      break;
+      case ElementTransformation::BDR_FACE:
+      {
+         // This must be a DG context so this dynamic cast must succeed.
+         FaceElementTransformations * FET =
+            dynamic_cast<FaceElementTransformations *>(&T);
+
+         // Evaluate in neighboring element (the integration point in T1 should
+         // have already been set).
+         ElementTransformation & T1 = FET->GetElement1Transformation();
+         return GetDivergence(T1);
+      }
+      break;
+      default:
+      {
+         MFEM_ABORT("GridFunction::GetDivergence: Unsupported element type \""
+                    << T.ElementType << "\"");
+      }
    }
-   else
-   {
-      // Assuming RT-type space
-      Array<int> dofs;
-      fes->GetElementDofs(elNo, dofs);
-      Vector loc_data, divshape(FElem->GetDof());
-      GetSubVector(dofs, loc_data);
-      FElem->CalcDivShape(tr.GetIntPoint(), divshape);
-      div_v = (loc_data * divshape) / tr.Weight();
-   }
-   return div_v;
+   return 0.0; // never reached
 }
 
-void GridFunction::GetCurl(ElementTransformation &tr, Vector &curl) const
+void GridFunction::GetCurl(ElementTransformation &T, Vector &curl) const
 {
-   int elNo = tr.ElementNo;
-   const FiniteElement *FElem = fes->GetFE(elNo);
-   if (FElem->GetRangeType() == FiniteElement::SCALAR)
+   switch (T.ElementType)
    {
-      MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
-                  "invalid FE map type");
-      DenseMatrix grad_hat;
-      GetVectorGradientHat(tr, grad_hat);
-      const DenseMatrix &Jinv = tr.InverseJacobian();
-      DenseMatrix grad(grad_hat.Height(), Jinv.Width()); // vdim x FElem->Dim
-      Mult(grad_hat, Jinv, grad);
-      MFEM_ASSERT(grad.Height() == grad.Width(), "");
-      if (grad.Height() == 3)
+      case ElementTransformation::ELEMENT:
       {
-         curl.SetSize(3);
-         curl(0) = grad(2,1) - grad(1,2);
-         curl(1) = grad(0,2) - grad(2,0);
-         curl(2) = grad(1,0) - grad(0,1);
+         int elNo = T.ElementNo;
+         const FiniteElement *fe = fes->GetFE(elNo);
+         if (fe->GetRangeType() == FiniteElement::SCALAR)
+         {
+            MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE,
+                        "invalid FE map type");
+            DenseMatrix grad_hat;
+            GetVectorGradientHat(T, grad_hat);
+            const DenseMatrix &Jinv = T.InverseJacobian();
+            // Dimensions of grad are vdim x FElem->Dim
+            DenseMatrix grad(grad_hat.Height(), Jinv.Width());
+            Mult(grad_hat, Jinv, grad);
+            MFEM_ASSERT(grad.Height() == grad.Width(), "");
+            if (grad.Height() == 3)
+            {
+               curl.SetSize(3);
+               curl(0) = grad(2,1) - grad(1,2);
+               curl(1) = grad(0,2) - grad(2,0);
+               curl(2) = grad(1,0) - grad(0,1);
+            }
+            else if (grad.Height() == 2)
+            {
+               curl.SetSize(1);
+               curl(0) = grad(1,0) - grad(0,1);
+            }
+         }
+         else
+         {
+            // Assuming ND-type space
+            Array<int> dofs;
+            fes->GetElementDofs(elNo, dofs);
+            Vector loc_data;
+            GetSubVector(dofs, loc_data);
+            DenseMatrix curl_shape(fe->GetDof(), fe->GetDim() == 3 ? 3 : 1);
+            fe->CalcCurlShape(T.GetIntPoint(), curl_shape);
+            curl.SetSize(curl_shape.Width());
+            if (curl_shape.Width() == 3)
+            {
+               double curl_hat[3];
+               curl_shape.MultTranspose(loc_data, curl_hat);
+               T.Jacobian().Mult(curl_hat, curl);
+            }
+            else
+            {
+               curl_shape.MultTranspose(loc_data, curl);
+            }
+            curl /= T.Weight();
+         }
       }
-      else if (grad.Height() == 2)
+      break;
+      case ElementTransformation::BDR_ELEMENT:
       {
-         curl.SetSize(1);
-         curl(0) = grad(1,0) - grad(0,1);
+         // In order to capture the tangential components of the curl we
+         // must evaluate it in the neighboring element.
+         FaceElementTransformations * FET =
+            fes->GetMesh()->GetBdrFaceTransformations(T.ElementNo);
+
+         // Boundary elements and Boundary Faces may have different
+         // orientations so adjust the integration point if necessary.
+         int o = 0;
+         if (fes->GetMesh()->Dimension() == 3)
+         {
+            int f;
+            fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+         }
+
+         IntegrationPoint fip;
+         be_to_bfe(FET->GetGeometryType(), o, T.GetIntPoint(), fip);
+
+         // Compute and set the point in element 1 from fip
+         FET->SetAllIntPoints(&fip);
+         ElementTransformation & T1 = FET->GetElement1Transformation();
+
+         GetCurl(T1, curl);
       }
-   }
-   else
-   {
-      // Assuming ND-type space
-      Array<int> dofs;
-      fes->GetElementDofs(elNo, dofs);
-      Vector loc_data;
-      GetSubVector(dofs, loc_data);
-      DenseMatrix curl_shape(FElem->GetDof(), FElem->GetDim() == 3 ? 3 : 1);
-      FElem->CalcCurlShape(tr.GetIntPoint(), curl_shape);
-      curl.SetSize(curl_shape.Width());
-      if (curl_shape.Width() == 3)
+      break;
+      case ElementTransformation::BDR_FACE:
       {
-         double curl_hat[3];
-         curl_shape.MultTranspose(loc_data, curl_hat);
-         tr.Jacobian().Mult(curl_hat, curl);
+         // This must be a DG context so this dynamic cast must succeed.
+         FaceElementTransformations * FET =
+            dynamic_cast<FaceElementTransformations *>(&T);
+
+         // Evaluate in neighboring element (the integration point in T1 should
+         // have already been set).
+         ElementTransformation & T1 = FET->GetElement1Transformation();
+         GetCurl(T1, curl);
       }
-      else
+      break;
+      default:
       {
-         curl_shape.MultTranspose(loc_data, curl);
+         MFEM_ABORT("GridFunction::GetCurl: Unsupported element type \""
+                    << T.ElementType << "\"");
       }
-      curl /= tr.Weight();
    }
 }
 
-void GridFunction::GetGradient(ElementTransformation &tr, Vector &grad) const
+void GridFunction::GetGradient(ElementTransformation &T, Vector &grad) const
 {
-   int elNo = tr.ElementNo;
-   const FiniteElement *fe = fes->GetFE(elNo);
-   MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE, "invalid FE map type");
-   int dim = fe->GetDim(), dof = fe->GetDof();
-   DenseMatrix dshape(dof, dim);
-   Vector lval, gh(dim);
-   Array<int> dofs;
+   switch (T.ElementType)
+   {
+      case ElementTransformation::ELEMENT:
+      {
+         const FiniteElement * fe = fes->GetFE(T.ElementNo);
+         MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE,
+                     "invalid FE map type");
+         int spaceDim = fes->GetMesh()->SpaceDimension();
+         int dim = fe->GetDim(), dof = fe->GetDof();
+         DenseMatrix dshape(dof, dim);
+         Vector lval, gh(dim);
+         Array<int> dofs;
 
-   grad.SetSize(dim);
-   fes->GetElementDofs(elNo, dofs);
-   GetSubVector(dofs, lval);
-   fe->CalcDShape(tr.GetIntPoint(), dshape);
-   dshape.MultTranspose(lval, gh);
-   tr.InverseJacobian().MultTranspose(gh, grad);
+         grad.SetSize(spaceDim);
+         fes->GetElementDofs(T.ElementNo, dofs);
+         GetSubVector(dofs, lval);
+         fe->CalcDShape(T.GetIntPoint(), dshape);
+         dshape.MultTranspose(lval, gh);
+         T.InverseJacobian().MultTranspose(gh, grad);
+      }
+      break;
+      case ElementTransformation::BDR_ELEMENT:
+      {
+         // In order to properly capture the normal component of the gradient
+         // as well as its tangential components we must evaluate it in the
+         // neighboring element.
+         FaceElementTransformations * FET =
+            fes->GetMesh()->GetBdrFaceTransformations(T.ElementNo);
+
+         // Boundary elements and Boundary Faces may have different
+         // orientations so adjust the integration point if necessary.
+         int o = 0;
+         if (fes->GetMesh()->Dimension() == 3)
+         {
+            int f;
+            fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+         }
+
+         IntegrationPoint fip;
+         be_to_bfe(FET->GetGeometryType(), o, T.GetIntPoint(), fip);
+
+         // Compute and set the point in element 1 from fip
+         FET->SetAllIntPoints(&fip);
+         ElementTransformation & T1 = FET->GetElement1Transformation();
+
+         GetGradient(T1, grad);
+      }
+      break;
+      case ElementTransformation::BDR_FACE:
+      {
+         // This must be a DG context so this dynamic cast must succeed.
+         FaceElementTransformations * FET =
+            dynamic_cast<FaceElementTransformations *>(&T);
+
+         // Evaluate in neighboring element (the integration point in T1 should
+         // have already been set).
+         ElementTransformation & T1 = FET->GetElement1Transformation();
+         GetGradient(T1, grad);
+      }
+      break;
+      default:
+      {
+         MFEM_ABORT("GridFunction::GetGradient: Unsupported element type \""
+                    << T.ElementType << "\"");
+      }
+   }
 }
 
 void GridFunction::GetGradients(ElementTransformation &tr,
@@ -1467,15 +1626,65 @@ void GridFunction::GetGradients(ElementTransformation &tr,
 }
 
 void GridFunction::GetVectorGradient(
-   ElementTransformation &tr, DenseMatrix &grad) const
+   ElementTransformation &T, DenseMatrix &grad) const
 {
-   MFEM_ASSERT(fes->GetFE(tr.ElementNo)->GetMapType() == FiniteElement::VALUE,
-               "invalid FE map type");
-   DenseMatrix grad_hat;
-   GetVectorGradientHat(tr, grad_hat);
-   const DenseMatrix &Jinv = tr.InverseJacobian();
-   grad.SetSize(grad_hat.Height(), Jinv.Width());
-   Mult(grad_hat, Jinv, grad);
+   switch (T.ElementType)
+   {
+      case ElementTransformation::ELEMENT:
+      {
+         MFEM_ASSERT(fes->GetFE(T.ElementNo)->GetMapType() ==
+                     FiniteElement::VALUE, "invalid FE map type");
+         DenseMatrix grad_hat;
+         GetVectorGradientHat(T, grad_hat);
+         const DenseMatrix &Jinv = T.InverseJacobian();
+         grad.SetSize(grad_hat.Height(), Jinv.Width());
+         Mult(grad_hat, Jinv, grad);
+      }
+      break;
+      case ElementTransformation::BDR_ELEMENT:
+      {
+         // In order to capture the normal component of the gradient we
+         // must evaluate it in the neighboring element.
+         FaceElementTransformations * FET =
+            fes->GetMesh()->GetBdrFaceTransformations(T.ElementNo);
+
+         // Boundary elements and Boundary Faces may have different
+         // orientations so adjust the integration point if necessary.
+         int o = 0;
+         if (fes->GetMesh()->Dimension() == 3)
+         {
+            int f;
+            fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+         }
+
+         IntegrationPoint fip;
+         be_to_bfe(FET->GetGeometryType(), o, T.GetIntPoint(), fip);
+
+         // Compute and set the point in element 1 from fip
+         FET->SetAllIntPoints(&fip);
+         ElementTransformation & T1 = FET->GetElement1Transformation();
+
+         GetVectorGradient(T1, grad);
+      }
+      break;
+      case ElementTransformation::BDR_FACE:
+      {
+         // This must be a DG context so this dynamic cast must succeed.
+         FaceElementTransformations * FET =
+            dynamic_cast<FaceElementTransformations *>(&T);
+
+         // Evaluate in neighboring element (the integration point in T1 should
+         // have already been set).
+         ElementTransformation & T1 = FET->GetElement1Transformation();
+         GetVectorGradient(T1, grad);
+      }
+      break;
+      default:
+      {
+         MFEM_ABORT("GridFunction::GetVectorGradient: "
+                    "Unsupported element type \"" << T.ElementType << "\"");
+      }
+   }
 }
 
 void GridFunction::GetElementAverages(GridFunction &avgs) const

fem/gridfunc.hpp

@@ -208,13 +208,14 @@ public:
    ///@{
    /** Return a scalar value from within the element indicated by the
        ElementTransformation Object. */
-   double GetValue(ElementTransformation &T, const IntegrationPoint &ip,
-                   int comp = 0, Vector *tr = NULL) const;
+   virtual double GetValue(ElementTransformation &T, const IntegrationPoint &ip,
+                           int comp = 0, Vector *tr = NULL) const;
 
    /** Return a vector value from within the element indicated by the
        ElementTransformation Object. */
-   void GetVectorValue(ElementTransformation &T, const IntegrationPoint &ip,
-                       Vector &val, Vector *tr = NULL) const;
+   virtual void GetVectorValue(ElementTransformation &T,
+                               const IntegrationPoint &ip,
+                               Vector &val, Vector *tr = NULL) const;
    ///@}
 
    /** @name ElementTransformation Get Values Methods

fem/gslib.cpp

@@ -192,6 +192,7 @@ void FindPointsGSLIB::Interpolate(Array<unsigned int> &codes,
    const int ncomp      = field_in.FESpace()->GetVDim(),
              points_fld = field_in.Size() / ncomp,
              points_cnt = codes.Size();
+   field_out.SetSize(points_cnt*ncomp);
 
    for (int i = 0; i < ncomp; i++)
    {

fem/lininteg.cpp

@@ -159,10 +159,13 @@ void BoundaryLFIntegrator::AssembleRHSElementVect(
    for (int i = 0; i < ir->GetNPoints(); i++)
    {
       const IntegrationPoint &ip = ir->IntPoint(i);
-      IntegrationPoint eip;
-      Tr.Loc1.Transform(ip, eip);
 
-      Tr.Face->SetIntPoint (&ip);
+      // Set the integration point in the face and the neighboring element
+      Tr.SetAllIntPoints(&ip);
+
+      // Access the neighboring element's integration point
+      const IntegrationPoint &eip = Tr.GetElement1IntPoint();
+
       double val = Tr.Face->Weight() * ip.weight * Q.Eval(*Tr.Face, ip);
 
       el.CalcShape(eip, shape);
@@ -359,10 +362,12 @@ void VectorBoundaryLFIntegrator::AssembleRHSElementVect(
    for (int i = 0; i < ir->GetNPoints(); i++)
    {
       const IntegrationPoint &ip = ir->IntPoint(i);
-      IntegrationPoint eip;
-      Tr.Loc1.Transform(ip, eip);
 
-      Tr.SetIntPoint(&ip);
+      // Set the integration point in the face and the neighboring element
+      Tr.SetAllIntPoints(&ip);
+
+      // Access the neighboring element's integration point
+      const IntegrationPoint &eip = Tr.GetElement1IntPoint();
 
       // Use Tr transformation in case Q depends on boundary attribute
       Q.Eval(vec, Tr, ip);
@@ -683,11 +688,13 @@ void BoundaryFlowIntegrator::AssembleRHSElementVect(
    for (int p = 0; p < ir->GetNPoints(); p++)
    {
       const IntegrationPoint &ip = ir->IntPoint(p);
-      IntegrationPoint eip;
-      Tr.Loc1.Transform(ip, eip);
-      el.CalcShape(eip, shape);
 
-      Tr.SetIntPoint(&ip);
+      // Set the integration point in the face and the neighboring element
+      Tr.SetAllIntPoints(&ip);
+
+      // Access the neighboring element's integration point
+      const IntegrationPoint &eip = Tr.GetElement1IntPoint();
+      el.CalcShape(eip, shape);
 
       // Use Tr.Elem1 transformation for u so that it matches the coefficient
       // used with the ConvectionIntegrator and/or the DGTraceIntegrator.
@@ -752,10 +759,13 @@ void DGDirichletLFIntegrator::AssembleRHSElementVect(
    for (int p = 0; p < ir->GetNPoints(); p++)
    {
       const IntegrationPoint &ip = ir->IntPoint(p);
-      IntegrationPoint eip;
 
-      Tr.Loc1.Transform(ip, eip);
-      Tr.SetIntPoint(&ip);
+      // Set the integration point in the face and the neighboring element
+      Tr.SetAllIntPoints(&ip);
+
+      // Access the neighboring element's integration point
+      const IntegrationPoint &eip = Tr.GetElement1IntPoint();
+
       if (dim == 1)
       {
          nor(0) = 2*eip.x - 1.0;
@@ -774,14 +784,14 @@ void DGDirichletLFIntegrator::AssembleRHSElementVect(
       {
          if (Q)
          {
-            w *= Q->Eval(Tr, ip);
+            w *= Q->Eval(*Tr.Elem1, eip);
          }
          ni.Set(w, nor);
       }
       else
       {
          nh.Set(w, nor);
-         MQ->Eval(mq, Tr, ip);
+         MQ->Eval(mq, *Tr.Elem1, eip);
          mq.MultTranspose(nh, ni);
       }
       CalcAdjugate(Tr.Elem1->Jacobian(), adjJ);
@@ -845,9 +855,12 @@ void DGElasticityDirichletLFIntegrator::AssembleRHSElementVect(
    for (int pi = 0; pi < ir->GetNPoints(); ++pi)
    {
       const IntegrationPoint &ip = ir->IntPoint(pi);
-      IntegrationPoint eip;
-      Tr.Loc1.Transform(ip, eip);
-      Tr.SetIntPoint(&ip);
+
+      // Set the integration point in the face and the neighboring element
+      Tr.SetAllIntPoints(&ip);
+
+      // Access the neighboring element's integration point
+      const IntegrationPoint &eip = Tr.GetElement1IntPoint();
 
       // Evaluate the Dirichlet b.c. using the face transformation.
       uD.Eval(u_dir, Tr, ip);

fem/pbilinearform.cpp

@@ -198,8 +198,9 @@ void ParBilinearForm::AssembleSharedFaces(int skip_zeros)
    for (int i = 0; i < nfaces; i++)
    {
       T = pmesh->GetSharedFaceTransformations(i);
+      int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
       pfes->GetElementVDofs(T->Elem1No, vdofs1);
-      pfes->GetFaceNbrElementVDofs(T->Elem2No, vdofs2);
+      pfes->GetFaceNbrElementVDofs(Elem2NbrNo, vdofs2);
       vdofs1.Copy(vdofs_all);
       for (int j = 0; j < vdofs2.Size(); j++)
       {
@@ -216,7 +217,7 @@ void ParBilinearForm::AssembleSharedFaces(int skip_zeros)
       for (int k = 0; k < fbfi.Size(); k++)
       {
          fbfi[k]->AssembleFaceMatrix(*pfes->GetFE(T->Elem1No),
-                                     *pfes->GetFaceNbrFE(T->Elem2No),
+                                     *pfes->GetFaceNbrFE(Elem2NbrNo),
                                      *T, elemmat);
          if (keep_nbr_block)
          {

fem/pgridfunc.cpp

@@ -214,7 +214,7 @@ void ParGridFunction::ExchangeFaceNbrData()
    ParMesh *pmesh = pfes->GetParMesh();
 
    face_nbr_data.SetSize(pfes->GetFaceNbrVSize());
-   Vector send_data(pfes->send_face_nbr_ldof.Size_of_connections());
+   send_data.SetSize(pfes->send_face_nbr_ldof.Size_of_connections());
 
    int *send_offset = pfes->send_face_nbr_ldof.GetI();
    const int *d_send_ldof = mfem::Read(pfes->send_face_nbr_ldof.GetJMemory(),
@@ -299,6 +299,108 @@ const
    return (DofVal * LocVec);
 }
 
+double ParGridFunction::GetValue(ElementTransformation &T,
+                                 const IntegrationPoint &ip,
+                                 int comp, Vector *tr) const
+{
+   // We can assume faces and edges are local
+   if (T.ElementType != ElementTransformation::ELEMENT)
+   {
+      return GridFunction::GetValue(T, ip, comp, tr);
+   }
+
+   // Check for evaluation in a local element
+   int nbr_el_no = T.ElementNo - pfes->GetParMesh()->GetNE();
+   if (nbr_el_no < 0)
+   {
+      return GridFunction::GetValue(T, ip, comp, tr);
+   }
+
+   // Evaluate using DoFs from a neighboring element
+   if (tr)
+   {
+      T.SetIntPoint(&ip);
+      T.Transform(ip, *tr);
+   }
+
+   Array<int> dofs;
+   const FiniteElement * fe = pfes->GetFaceNbrFE(nbr_el_no);
+   pfes->GetFaceNbrElementVDofs(nbr_el_no, dofs);
+
+   pfes->DofsToVDofs(comp-1, dofs);
+   Vector DofVal(dofs.Size()), LocVec;
+   if (fe->GetMapType() == FiniteElement::VALUE)
+   {
+      fe->CalcShape(ip, DofVal);
+   }
+   else
+   {
+      fe->CalcPhysShape(T, DofVal);
+   }
+   face_nbr_data.GetSubVector(dofs, LocVec);
+
+   return (DofVal * LocVec);
+}
+
+void ParGridFunction::GetVectorValue(ElementTransformation &T,
+                                     const IntegrationPoint &ip,
+                                     Vector &val, Vector *tr) const
+{
+   // We can assume faces and edges are local
+   if (T.ElementType != ElementTransformation::ELEMENT)
+   {
+      return GridFunction::GetVectorValue(T, ip, val, tr);
+   }
+
+   // Check for evaluation in a local element
+   int nbr_el_no = T.ElementNo - pfes->GetParMesh()->GetNE();
+   if (nbr_el_no < 0)
+   {
+      return GridFunction::GetVectorValue(T, ip, val, tr);
+   }
+
+   // Evaluate using DoFs from a neighboring element
+   if (tr)
+   {
+      T.SetIntPoint(&ip);
+      T.Transform(ip, *tr);
+   }
+
+   Array<int> vdofs;
+   pfes->GetFaceNbrElementVDofs(nbr_el_no, vdofs);
+   const FiniteElement *fe = pfes->GetFaceNbrFE(nbr_el_no);
+
+   int dof = fe->GetDof();
+   Vector loc_data;
+   face_nbr_data.GetSubVector(vdofs, loc_data);
+   if (fe->GetRangeType() == FiniteElement::SCALAR)
+   {
+      Vector shape(dof);
+      if (fe->GetMapType() == FiniteElement::VALUE)
+      {
+         fe->CalcShape(ip, shape);
+      }
+      else
+      {
+         fe->CalcPhysShape(T, shape);
+      }
+      int vdim = pfes->GetVDim();
+      val.SetSize(vdim);
+      for (int k = 0; k < vdim; k++)
+      {
+         val(k) = shape * ((const double *)loc_data + dof * k);
+      }
+   }
+   else
+   {
+      int spaceDim = pfes->GetMesh()->SpaceDimension();
+      DenseMatrix vshape(dof, spaceDim);
+      fe->CalcVShape(T, vshape);
+      val.SetSize(spaceDim);
+      vshape.MultTranspose(loc_data, val);
+   }
+}
+
 void ParGridFunction::ProjectCoefficient(Coefficient &coeff)
 {
    DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);

fem/pgridfunc.hpp

@@ -38,6 +38,11 @@ protected:
        initialized by ExchangeFaceNbrData(). */
    Vector face_nbr_data;
 
+   /** @brief Vector used as an MPI buffer to send face-neighbor data
+       in ExchangeFaceNbrData() to neighboring processors. */
+   //TODO: Use temporary memory to avoid CUDA malloc allocation cost.
+   Vector send_data;
+
    void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff,
                               Array<int> &attr);
 
@@ -204,6 +209,15 @@ public:
    double GetValue(ElementTransformation &T)
    { return GetValue(T.ElementNo, T.GetIntPoint()); }
 
+   // Redefine to handle the case when T describes a face-neighbor element
+   virtual double GetValue(ElementTransformation &T, const IntegrationPoint &ip,
+                           int comp = 0, Vector *tr = NULL) const;
+
+   // Redefine to handle the case when T describes a face-neighbor element
+   virtual void GetVectorValue(ElementTransformation &T,
+                               const IntegrationPoint &ip,
+                               Vector &val, Vector *tr = NULL) const;
+
    using GridFunction::ProjectCoefficient;
    virtual void ProjectCoefficient(Coefficient &coeff);
 

fem/pnonlinearform.cpp

@@ -64,12 +64,13 @@ void ParNonlinearForm::Mult(const Vector &x, Vector &y) const
       for (int i = 0; i < n_shared_faces; i++)
       {
          tr = pmesh->GetSharedFaceTransformations(i, true);
+         int Elem2NbrNo = tr->Elem2No - pmesh->GetNE();
 
          fe1 = pfes->GetFE(tr->Elem1No);
-         fe2 = pfes->GetFaceNbrFE(tr->Elem2No);
+         fe2 = pfes->GetFaceNbrFE(Elem2NbrNo);
 
          pfes->GetElementVDofs(tr->Elem1No, vdofs1);
-         pfes->GetFaceNbrElementVDofs(tr->Elem2No, vdofs2);
+         pfes->GetFaceNbrElementVDofs(Elem2NbrNo, vdofs2);
 
          el_x.SetSize(vdofs1.Size() + vdofs2.Size());
          X.GetSubVector(vdofs1, el_x.GetData());

fem/tmop.cpp

@@ -177,24 +177,6 @@ double TMOP_Metric_SSA2D::EvalW(const DenseMatrix &Jpt) const
    return Mat.FNorm2();
 }
 
-// mu_85 = |T-T'|^2, where T'= |T|*I/sqrt(2)
-double TMOP_Metric_SS2D::EvalW(const DenseMatrix &Jpt) const
-{
-   MFEM_VERIFY(Jtr != NULL,
-               "Requires a target Jacobian, use SetTargetJacobian().");
-
-   DenseMatrix Id(2,2);
-   DenseMatrix Mat(2,2);
-   Mat = Jpt;
-
-   Id(0,0) = 1; Id(0,1) = 0;
-   Id(1,0) = 0; Id(1,1) = 1;
-   Id *= Mat.FNorm()/pow(2,0.5);
-
-   Mat.Add(-1.,Id);
-   return Mat.FNorm2();
-}
-
 double TMOP_Metric_002::EvalW(const DenseMatrix &Jpt) const
 {
    ie.SetJacobian(Jpt.GetData());
@@ -484,6 +466,24 @@ void TMOP_Metric_077::AssembleH(const DenseMatrix &Jpt,
    ie.Assemble_TProd(weight * I2inv_sq / I2, ie.Get_dI2(), A.GetData());
 }
 
+// mu_85 = |T-T'|^2, where T'= |T|*I/sqrt(2)
+double TMOP_Metric_085::EvalW(const DenseMatrix &Jpt) const
+{
+   MFEM_VERIFY(Jtr != NULL,
+               "Requires a target Jacobian, use SetTargetJacobian().");
+
+   DenseMatrix Id(2,2);
+   DenseMatrix Mat(2,2);
+   Mat = Jpt;
+
+   Id(0,0) = 1; Id(0,1) = 0;
+   Id(1,0) = 0; Id(1,1) = 1;
+   Id *= Mat.FNorm()/pow(2,0.5);
+
+   Mat.Add(-1.,Id);
+   return Mat.FNorm2();
+}
+
 double TMOP_Metric_211::EvalW(const DenseMatrix &Jpt) const
 {
    // mu_211 = (det(J) - 1)^2 - det(J) + (det(J)^2 + eps)^{1/2}

fem/tmop.hpp

@@ -162,21 +162,6 @@ public:
    { MFEM_ABORT("Not implemented"); }
 };
 
-/// Shape+Size metric, 2D.
-class TMOP_Metric_SS2D : public TMOP_QualityMetric
-{
-public:
-   // W = 0.5 (1 - cos(theta_Jpr - theta_Jtr)).
-   virtual double EvalW(const DenseMatrix &Jpt) const;
-
-   virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
-   { MFEM_ABORT("Not implemented"); }
-
-   virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
-                          const double weight, DenseMatrix &A) const
-   { MFEM_ABORT("Not implemented"); }
-};
-
 /// Shape, ideal barrier metric, 2D
 class TMOP_Metric_002 : public TMOP_QualityMetric
 {
@@ -331,6 +316,21 @@ public:
 
 };
 
+/// Shape & orientation metric, 2D.
+class TMOP_Metric_085 : public TMOP_QualityMetric
+{
+public:
+   // W = |T-T'|^2, where T'= |T|*I/sqrt(2).
+   virtual double EvalW(const DenseMatrix &Jpt) const;
+
+   virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
+   { MFEM_ABORT("Not implemented"); }
+
+   virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
+                          const double weight, DenseMatrix &A) const
+   { MFEM_ABORT("Not implemented"); }
+};
+
 /// Untangling metric, 2D
 class TMOP_Metric_211 : public TMOP_QualityMetric
 {

fem/tmop_tools.cpp

@@ -353,20 +353,59 @@ double TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
       energy_in = nlf->GetEnergy(x);
    }
 
-   const bool have_b = (b.Size() == Height());
-
    const int NE = fes->GetMesh()->GetNE(), dim = fes->GetFE(0)->GetDim(),
              dof = fes->GetFE(0)->GetDof(), nsp = ir.GetNPoints();
    Array<int> xdofs(dof * dim);
    DenseMatrix Jpr(dim), dshape(dof, dim), pos(dof, dim);
    Vector posV(pos.Data(), dof * dim);
+   Vector x_out_loc(fes->GetVSize());
 
-   Vector x_out(x.Size()), x_out_loc(fes->GetVSize());
+   if (serial)
+   {
+      const SparseMatrix *cP = fes->GetConformingProlongation();
+      if (!cP) { x_out_loc = x; }
+      else     { cP->Mult(x, x_out_loc); }
+   }
+#ifdef MFEM_USE_MPI
+   else
+   {
+      fes->GetProlongationMatrix()->Mult(x, x_out_loc);
+   }
+#endif
+
+   double min_detJ = infinity();
+   for (int i = 0; i < NE; i++)
+   {
+      fes->GetElementVDofs(i, xdofs);
+      x_out_loc.GetSubVector(xdofs, posV);
+
+      for (int j = 0; j < nsp; j++)
+      {
+         fes->GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
+         MultAtB(pos, dshape, Jpr);
+         min_detJ = std::min(min_detJ, Jpr.Det());
+      }
+   }
+   double min_detJ_all = min_detJ;
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      MPI_Allreduce(&min_detJ, &min_detJ_all, 1, MPI_DOUBLE, MPI_MIN,
+                    p_nlf->ParFESpace()->GetComm());
+   }
+#endif
+   bool untangling = false;
+   if (min_detJ_all <= 0) { untangling = true; }
+
+   const bool have_b = (b.Size() == Height());
+
+   Vector x_out(x.Size());
    bool x_out_ok = false;
    double scale = 1.0, energy_out = 0.0;
    double norm0 = Norm(r);
 
-   // Decreases the scaling of the update until the new mesh is valid.
+   const double detJ_factor = (solver_type == 1) ? 0.25 : 0.5;
+
    for (int i = 0; i < 12; i++)
    {
       add(x, -scale, c, x_out);
@@ -384,35 +423,39 @@ double TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
       }
 #endif
 
-      int jac_ok = 1;
-      for (int i = 0; i < NE; i++)
+      // Check det(Jpr) > 0.
+      if (!untangling)
       {
-         fes->GetElementVDofs(i, xdofs);
-         x_out_loc.GetSubVector(xdofs, posV);
-         for (int j = 0; j < nsp; j++)
+         int jac_ok = 1;
+         for (int i = 0; i < NE; i++)
          {
-            fes->GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
-            MultAtB(pos, dshape, Jpr);
-            if (Jpr.Det() <= 0.0) { jac_ok = 0; goto break2; }
+            fes->GetElementVDofs(i, xdofs);
+            x_out_loc.GetSubVector(xdofs, posV);
+            for (int j = 0; j < nsp; j++)
+            {
+               fes->GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
+               MultAtB(pos, dshape, Jpr);
+               if (Jpr.Det() <= 0.0) { jac_ok = 0; goto break2; }
+            }
          }
-      }
 
-   break2:
-      int jac_ok_all = jac_ok;
+      break2:
+         int jac_ok_all = jac_ok;
 #ifdef MFEM_USE_MPI
-      if (parallel)
-      {
-         MPI_Allreduce(&jac_ok, &jac_ok_all, 1, MPI_INT, MPI_LAND,
-                       p_nlf->ParFESpace()->GetComm());
-      }
+         if (parallel)
+         {
+            MPI_Allreduce(&jac_ok, &jac_ok_all, 1, MPI_INT, MPI_LAND,
+                          p_nlf->ParFESpace()->GetComm());
+         }
 #endif
 
-      if (jac_ok_all == 0)
-      {
-         if (print_level >= 0)
-         { mfem::out << "Scale = " << scale << " Neg det(J) found.\n"; }
-         scale *= 0.5; continue;
-      }
+         if (jac_ok_all == 0)
+         {
+            if (print_level >= 0)
+            { mfem::out << "Scale = " << scale << " Neg det(J) found.\n"; }
+            scale *= detJ_factor; continue;
+         }
+      } // endif(!untangling)
 
       ProcessNewState(x_out);
       if (serial)
@@ -425,25 +468,37 @@ double TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
          energy_out = p_nlf->GetParGridFunctionEnergy(x_out_loc);
       }
 #endif
-      if (energy_out > 1.2*energy_in || std::isnan(energy_out) != 0)
-      {
-         if (print_level >= 0)
-         { mfem::out << "Scale = " << scale << " Increasing energy.\n"; }
-         scale *= 0.5; continue;
-      }
 
-      oper->Mult(x_out, r);
-      if (have_b) { r -= b; }
-      double norm = Norm(r);
-
-      if (norm > 1.2*norm0)
+      if (untangling)
       {
-         if (print_level >= 0)
-         { mfem::out << "Scale = " << scale << " Norm increased.\n"; }
-         scale *= 0.5; continue;
+         if (energy_out > energy_in || std::isnan(energy_out) != 0)
+         {
+            scale *= 0.5;
+         }
+         else { x_out_ok = true; break; }
       }
-      else { x_out_ok = true; break; }
-   }
+      else
+      {
+         if (energy_out > 1.2*energy_in || std::isnan(energy_out) != 0)
+         {
+            if (print_level >= 0)
+            { mfem::out << "Scale = " << scale << " Increasing energy.\n"; }
+            scale *= 0.5; continue;
+         }
+
+         oper->Mult(x_out, r);
+         if (have_b) { r -= b; }
+         double norm = Norm(r);
+
+         if (norm > 1.2*norm0)
+         {
+            if (print_level >= 0)
+            { mfem::out << "Scale = " << scale << " Norm increased.\n"; }
+            scale *= 0.5; continue;
+         }
+         else { x_out_ok = true; break; }
+      } // endif (untangling)
+   } // enddo (i)
 
    if (print_level >= 0)
    {
@@ -570,105 +625,6 @@ void TMOPNewtonSolver::UpdateDiscreteTC(const TMOP_Integrator &ti,
    }
 }
 
-double TMOPDescentNewtonSolver::ComputeScalingFactor(const Vector &x,
-                                                     const Vector &b) const
-{
-   const FiniteElementSpace *fes = NULL;
-   double energy_in = 0.0;
-#ifdef MFEM_USE_MPI
-   const ParNonlinearForm *p_nlf = dynamic_cast<const ParNonlinearForm *>(oper);
-   MFEM_VERIFY(!(parallel && p_nlf == NULL), "Invalid Operator subclass.");
-   if (parallel)
-   {
-      fes = p_nlf->FESpace();
-      energy_in = p_nlf->GetEnergy(x);
-   }
-#endif
-   const bool serial = !parallel;
-   const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
-   MFEM_VERIFY(!(serial && nlf == NULL), "Invalid Operator subclass.");
-   if (serial)
-   {
-      fes = nlf->FESpace();
-      energy_in = nlf->GetEnergy(x);
-   }
-
-   const int NE = fes->GetMesh()->GetNE(), dim = fes->GetFE(0)->GetDim(),
-             dof = fes->GetFE(0)->GetDof(), nsp = ir.GetNPoints();
-   Array<int> xdofs(dof * dim);
-   DenseMatrix Jpr(dim), dshape(dof, dim), pos(dof, dim);
-   Vector posV(pos.Data(), dof * dim);
-   Vector x_loc(fes->GetVSize());
-
-   double min_detJ = infinity();
-   for (int i = 0; i < NE; i++)
-   {
-      fes->GetElementVDofs(i, xdofs);
-      // TODO x_loc doesn't have valid values here!
-      MFEM_ABORT("This function has to be fixed!");
-      x_loc.GetSubVector(xdofs, posV);
-
-      for (int j = 0; j < nsp; j++)
-      {
-         fes->GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
-         MultAtB(pos, dshape, Jpr);
-         min_detJ = std::min(min_detJ, Jpr.Det());
-      }
-   }
-   double min_detJ_all = min_detJ;
-#ifdef MFEM_USE_MPI
-   if (parallel)
-   {
-      MPI_Allreduce(&min_detJ, &min_detJ_all, 1, MPI_DOUBLE, MPI_MIN,
-                    p_nlf->ParFESpace()->GetComm());
-   }
-#endif
-   if (print_level >= 0)
-   {
-      mfem::out << "Minimum det(J) = " << min_detJ_all << '\n';
-   }
-
-   Vector x_out(x.Size());
-   bool x_out_ok = false;
-   double scale = 1.0, energy_out = 0.0;
-
-   for (int i = 0; i < 7; i++)
-   {
-      add(x, -scale, c, x_out);
-      if (serial)
-      {
-         const SparseMatrix *cP = fes->GetConformingProlongation();
-         if (!cP) { x_loc = x_out; }
-         else     { cP->Mult(x_out,x_loc); }
-         energy_out = nlf->GetGridFunctionEnergy(x_loc);
-      }
-#ifdef MFEM_USE_MPI
-      else
-      {
-         fes->GetProlongationMatrix()->Mult(x_out, x_loc);
-         energy_out = p_nlf->GetParGridFunctionEnergy(x_loc);
-      }
-#endif
-
-      if (energy_out > energy_in || std::isnan(energy_out) != 0)
-      {
-         scale *= 0.5;
-      }
-      else { x_out_ok = true; break; }
-   }
-
-   if (print_level >= 0)
-   {
-      mfem::out << "Energy decrease: "
-                << (energy_in - energy_out) / energy_in * 100.0
-                << "% with " << scale << " scaling.\n";
-   }
-
-   if (x_out_ok == false) { return 0.0; }
-
-   return scale;
-}
-
 #ifdef MFEM_USE_MPI
 // Metric values are visualized by creating an L2 finite element functions and
 // computing the metric values at the nodes.

fem/tmop_tools.hpp

@@ -109,9 +109,11 @@ public:
 };
 #endif
 
-class TMOPNewtonSolver : public NewtonSolver
+class TMOPNewtonSolver : public LBFGSSolver
 {
 protected:
+   // 0 - Newton, 1 - LBFGS.
+   int solver_type;
    bool parallel;
 
    // Quadrature points that are checked for negative Jacobians etc.
@@ -121,29 +123,42 @@ protected:
 
 public:
 #ifdef MFEM_USE_MPI
-   TMOPNewtonSolver(MPI_Comm comm, const IntegrationRule &irule)
-      : NewtonSolver(comm), parallel(true), ir(irule) { }
+   TMOPNewtonSolver(MPI_Comm comm, const IntegrationRule &irule, int type = 0)
+      : LBFGSSolver(comm), solver_type(type), parallel(true), ir(irule) { }
 #endif
-   TMOPNewtonSolver(const IntegrationRule &irule)
-      : NewtonSolver(), parallel(false), ir(irule) { }
+   TMOPNewtonSolver(const IntegrationRule &irule, int type = 0)
+      : LBFGSSolver(), solver_type(type), parallel(false), ir(irule) { }
 
    virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
 
    virtual void ProcessNewState(const Vector &x) const;
-};
 
-/// Allows negative Jacobians. Used for untangling.
-class TMOPDescentNewtonSolver : public TMOPNewtonSolver
-{
-public:
-#ifdef MFEM_USE_MPI
-   TMOPDescentNewtonSolver(MPI_Comm comm, const IntegrationRule &irule)
-      : TMOPNewtonSolver(comm, irule) { }
-#endif
-   TMOPDescentNewtonSolver(const IntegrationRule &irule)
-      : TMOPNewtonSolver(irule) { }
+   virtual void Mult(const Vector &b, Vector &x) const
+   {
+      if (solver_type == 0)
+      {
+         NewtonSolver::Mult(b, x);
+      }
+      else if (solver_type == 1)
+      {
+         LBFGSSolver::Mult(b, x);
+      }
+      else { MFEM_ABORT("Invalid type"); }
+   }
 
-   virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
+   virtual void SetSolver(Solver &solver)
+   {
+      if (solver_type == 0)
+      {
+         NewtonSolver::SetSolver(solver);
+      }
+      else if (solver_type == 1)
+      {
+         LBFGSSolver::SetSolver(solver);
+      }
+      else { MFEM_ABORT("Invalid type"); }
+   }
+   virtual void SetPreconditioner(Solver &pr) { SetSolver(pr); }
 };
 
 void vis_tmop_metric_s(int order, TMOP_QualityMetric &qm,

general/mem_manager.hpp

@@ -450,16 +450,16 @@ private:
 
    template <std::size_t align_bytes, bool dummy = true> struct Alloc
    {
+#if __cplusplus < 201703L
       static inline T *New(std::size_t)
       {
-#if __cplusplus < 201703L
          // Generate an error in debug mode
          MFEM_ASSERT(false, "overaligned type cannot use MemoryType::HOST");
          return nullptr;
-#else
-         return new T[size];
-#endif
       }
+#else
+      static inline T *New(std::size_t size) { return new T[size]; }
+#endif
    };
 
 #if __cplusplus < 201703L

linalg/fdual.hpp

@@ -0,0 +1,628 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+
+#ifndef FDUAL_H
+#define FDUAL_H
+
+#include <cmath>
+#include <type_traits>
+
+namespace mfem
+{
+
+namespace ad
+{
+// Forward AD - simple class for automatic differentiation
+template<typename tbase>
+class FDual
+{
+private:
+   tbase pr;
+   tbase du;
+
+public:
+
+   FDual():pr(0),du(0)
+   {
+       
+   }
+
+   
+   template <class fltyp, class = typename 
+   std::enable_if<std::is_arithmetic<fltyp>::value>::type>
+   FDual(fltyp& f):pr(f),du(0) 
+   {
+
+   }
+    
+   template <class fltyp, class = typename 
+   std::enable_if<std::is_arithmetic<fltyp>::value>::type>
+   FDual(const fltyp& f):pr(f),du(0) 
+   {
+
+   }
+   
+   FDual(tbase& pr_,tbase& du_):pr(pr_),du(du_)
+   {
+
+   }
+
+   FDual(const tbase& pr_,const tbase& du_):pr(pr_),du(du_)
+   {
+       
+   }
+   
+   FDual(FDual<tbase>& nm):pr(nm.pr),du(nm.du)
+   {
+       
+   }
+   
+   FDual(const FDual<tbase>& nm):pr(nm.pr),du(nm.du)
+   {
+       
+   }
+   
+   tbase prim() const
+   {
+       return pr;
+   }
+   
+   tbase real() const
+   {
+      return pr;
+   }
+
+   tbase dual() const
+   {
+      return du;
+   }
+
+   void set(const tbase& pr_,const tbase& du_)
+   {
+       pr=pr_;
+       du=du_;
+   }
+   
+   void prim(const tbase& pr_)
+   {
+       pr=pr_;
+   }
+
+   void real(const tbase& pr_)
+   {
+       pr=pr_;
+   }
+   
+   void dual(const tbase& du_)
+   {
+       du=du_;
+   }
+   
+   FDual<tbase> & operator=(tbase sc_)
+   {
+      pr=sc_;
+      du=tbase(0);
+      return *this;
+   }
+
+   FDual<tbase> & operator+=(tbase sc_)
+   {
+      pr=pr+sc_;
+      return *this;
+   }
+
+   FDual<tbase> & operator-=(tbase sc_)
+   {
+      pr=pr-sc_;
+      return *this;
+
+   }
+
+   FDual<tbase> & operator*=(tbase sc_)
+   {
+      pr=pr*sc_;
+      du=du*sc_;
+      return *this;
+   }
+
+   FDual<tbase>& operator/=(tbase sc_)
+   {
+      pr=pr/sc_;
+      du=du/sc_;
+      return *this;
+   }
+
+
+   FDual<tbase>& operator=(const FDual<tbase> & f)
+   {
+      pr = f.real();
+      du = f.dual();
+      return *this;
+   }
+
+   FDual<tbase>& operator+=(const FDual<tbase>& f)
+   {
+      pr += f.real();
+      du += f.dual();
+      return *this;
+   }
+
+   
+   FDual<tbase>& operator-=(const FDual<tbase>& f)
+   {
+      pr -= f.real();
+      du -= f.dual();
+      return *this;
+   }
+
+   FDual<tbase>& operator*=(const FDual<tbase>& f)
+   {
+      du = du * f.real();
+      du = du+ pr * f.dual();
+      pr = pr * f.real();
+      return *this;
+   }
+
+   FDual<tbase>& operator/=(const FDual<tbase>& f_)
+   {
+      pr = pr / f_.real();
+      du = du - pr * f_.dual();
+      du = du / f_.real();
+      return *this;
+   }
+};
+
+
+// non-member functions
+// boolean operations
+template <typename tbase>
+inline
+bool operator==(const FDual<tbase>& a1, const FDual<tbase>& a2)
+{
+   return a1.real() == a2.real();
+}
+
+template <typename tbase>
+inline
+bool operator==(tbase a, const FDual<tbase>& f_)
+{
+   return a == f_.real();
+}
+
+template <typename tbase>
+inline
+bool operator==(const FDual<tbase>& a, tbase b)
+{
+   return a.real() == b;
+}
+
+
+template <typename tbase>
+inline
+bool operator<(const FDual<tbase>& f1, const FDual<tbase>& f2)
+{
+   return f1.real() < f2.real();
+}
+
+template <typename tbase>
+inline
+bool operator<(const FDual<tbase>& f, tbase a)
+{
+   return f.real() < a;
+}
+
+template <typename tbase>
+inline
+bool operator<(tbase a, const FDual<tbase>& f)
+{
+   return a < f.real();
+}
+
+template <typename tbase>
+inline
+bool operator>(const FDual<tbase>& f1, const FDual<tbase>& f2)
+{
+   return f1.real() > f2.real();
+}
+
+template <typename tbase>
+inline
+bool operator>(const FDual<tbase>& f, tbase a)
+{
+   return f.real() > a;
+}
+
+template <typename tbase>
+inline
+bool operator>(tbase a, const FDual<tbase>& f)
+{
+   return (a > f.real());
+}
+
+template <typename tbase>
+inline
+FDual<tbase> operator-(const FDual<tbase>& f)
+{
+   return FDual<tbase>(-f.real(), -f.dual());
+}
+
+template <typename tbase>
+inline
+FDual<tbase> operator-(const FDual<tbase>& f, tbase a)
+{
+   return FDual<tbase>(f.real() - a, f.dual());
+}
+
+template <typename tbase>
+inline
+FDual<FDual<tbase>> operator-(const FDual<FDual<tbase>>& f, tbase a)
+{
+   return FDual<FDual<tbase>>(f.real() - a, f.dual());
+}
+
+
+template <typename tbase>
+inline
+FDual<tbase> operator+(const FDual<tbase>& f, tbase a)
+{
+   return FDual<tbase>(f.real() + a, f.dual());
+}
+
+template <typename tbase>
+inline
+FDual<FDual<tbase>> operator+(const FDual<FDual<tbase>>& f, tbase a)
+{
+   return FDual<FDual<tbase>>(f.real() + a, f.dual());
+}
+
+
+template <typename tbase>
+inline
+FDual<tbase> operator*(const FDual<tbase>& f, tbase a)
+{
+   return FDual<tbase>(f.real() * a, f.dual() * a);
+}
+
+
+template <typename tbase>
+inline
+FDual<tbase> operator/(const FDual<tbase>& f, tbase a)
+{
+   return FDual<tbase>(f.real() / a, f.dual() / a);
+}
+
+template <typename tbase>
+inline
+FDual<FDual<tbase>> operator/(const FDual<FDual<tbase>>& f, tbase a)
+{
+   return FDual<FDual<tbase>>(f.real() / a, f.dual() / a);
+}
+
+template <typename tbase>
+inline
+FDual<tbase> operator+(tbase a, const FDual<tbase>& f)
+{
+   return FDual<tbase>(a + f.real(), f.dual());
+}
+
+template <typename tbase>
+inline
+FDual<FDual<tbase>> operator+(tbase a, const FDual<FDual<tbase>>& f)
+{
+   return FDual<FDual<tbase>>(a + f.real(), f.dual());
+}
+
+
+template <typename tbase>
+inline
+FDual<tbase> operator-(tbase a, const FDual<tbase>& f)
+{
+   return FDual<tbase>(a - f.real(), -f.dual());
+}
+
+template <typename tbase>
+inline
+FDual<FDual<tbase>> operator-(tbase a, const FDual<FDual<tbase>>& f)
+{
+   return FDual<FDual<tbase>>(a - f.real(), -f.dual());
+}
+
+template <typename tbase>
+inline
+FDual<tbase> operator*(tbase a, const FDual<tbase>& f)
+{
+   return FDual<tbase>(f.real() * a, f.dual() *a);
+}
+
+template <typename tbase>
+inline
+FDual<FDual<tbase>> operator*(tbase a, const FDual<FDual<tbase>>& f)
+{
+   return FDual<FDual<tbase>>(f.real() * a, f.dual() *a);
+}
+
+
+template <typename tbase>
+inline
+FDual<tbase> operator/(tbase a, const FDual<tbase>& f)
+{
+   a = a / f.real();
+   return FDual<tbase>(a, -a * f.dual() / f.real());
+}
+
+template <typename tbase>
+inline
+FDual<tbase> operator+(const FDual<tbase>& f1, const FDual<tbase>& f2)
+{
+   return FDual<tbase>(f1.real() + f2.real(), f1.dual() + f2.dual());
+}
+
+template <typename tbase>
+inline
+FDual<tbase> operator-(const FDual<tbase>& f1, const FDual<tbase>& f2)
+{
+   return FDual<tbase>(f1.real() - f2.real(), f1.dual() - f2.dual());
+}
+
+template <typename tbase>
+inline
+FDual<tbase> operator*(const FDual<tbase>& f1, const FDual<tbase>& f2)
+{
+   return FDual<tbase>(f1.real() * f2.real(),
+                       f1.real() * f2.dual() + f1.dual() * f2.real());
+}
+
+template <typename tbase>
+inline
+FDual<tbase> operator/(const FDual<tbase>& f1, const FDual<tbase>& f2)
+{
+   tbase a=tbase(1)/f2.real();
+   tbase b=f1.real()*a;
+   return FDual<tbase>(b, (f1.dual() - f2.dual()*b)*a);
+}
+
+template <typename tbase>
+inline
+FDual<tbase> acos(const FDual<tbase>& f)
+{
+   return FDual<tbase>(acos(f.real()),
+                       -f.dual() / sqrt(tbase(1) - f.real() * f.real()));
+}
+
+template <>
+inline
+FDual<double> acos(const FDual<double>& f)
+{
+   return FDual<double>(std::acos(f.real()),
+                        -f.dual() / std::sqrt(double(1) - f.real() * f.real()));
+}
+
+template <typename tbase>
+inline
+FDual<tbase> asin(const FDual<tbase>& f)
+{
+   return FDual<tbase>(asin(f.real()),
+                       f.dual() / sqrt(tbase(1) - f.real() * f.real()));
+}
+
+template <>
+inline
+FDual<double> asin(const FDual<double>& f)
+{
+   return FDual<double>(std::asin(f.real()),
+                        f.dual() / std::sqrt(double(1) - f.real() * f.real()));
+}
+
+template <typename tbase>
+inline
+FDual<tbase> atan(const FDual<tbase>& f)
+{
+   return FDual<tbase>(atan(f.real()),
+                       f.dual() / (tbase(1) + f.real() * f.real()));
+}
+
+template <>
+inline
+FDual<double> atan(const FDual<double>& f)
+{
+   return FDual<double>(std::atan(f.real()),
+                        f.dual() / (double(1) + f.real() * f.real()));
+}
+
+template <typename tbase>
+inline
+FDual<tbase> cos(const FDual<tbase>& f)
+{
+   return FDual<tbase>(cos(f.real()), -f.dual() * sin(f.real()));
+}
+
+template <>
+inline
+FDual<double> cos(const FDual<double>& f)
+{
+   return FDual<double>(std::cos(f.real()), -f.dual() * std::sin(f.real()));
+}
+
+template <typename tbase>
+inline
+FDual<tbase> cosh(const FDual<tbase>& f)
+{
+   return FDual<tbase>(cosh(f.real()), f.dual() * sinh(f.real()));
+}
+
+template <>
+inline
+FDual<double> cosh(const FDual<double>& f)
+{
+   return FDual<double>(std::cosh(f.real()), f.dual() * std::sinh(f.real()));
+}
+
+template <typename tbase>
+inline
+FDual<tbase> exp(const FDual<tbase>& f)
+{
+   tbase x = exp(f.real());
+   return FDual<tbase>(x, f.dual() * x);
+}
+
+template <>
+inline
+FDual<double> exp(const FDual<double>& f)
+{
+   double x = std::exp(f.real());
+   return FDual<double>(x, f.dual() * x);
+}
+
+template <typename tbase>
+inline
+FDual<tbase> log(const FDual<tbase>& f)
+{
+   return FDual<tbase>(log(f.real()), f.dual() / f.real());
+}
+
+template <>
+inline
+FDual<double> log(const FDual<double>& f)
+{
+   return FDual<double>(std::log(f.real()), f.dual() / f.real());
+}
+
+template <typename tbase>
+inline
+FDual<tbase> log10(const FDual<tbase>& f)
+{
+   return log(f) / log(tbase(10));
+}
+
+template <>
+inline
+FDual<double> log10(const FDual<double>& f)
+{
+   return log(f) / std::log(double(10));
+}
+
+template <typename tbase>
+inline
+FDual<tbase> pow(const FDual<tbase>& a, const FDual<tbase>& b)
+{
+   return exp(log(a) * b);
+}
+
+template <typename tbase, typename tbase1>
+inline
+FDual<tbase> pow(const FDual<tbase>& a, const tbase1& b)
+{
+   return exp(log(a) * tbase(b));
+}
+
+
+template <typename tbase, typename tbase1>
+inline
+FDual<tbase> pow(const tbase1& a, const FDual<tbase>& b)
+{
+   return exp(log(tbase(a)) * b);
+}
+
+template <>
+inline
+FDual<double> pow(const double& a, const FDual<double>& b)
+{
+   return exp(std::log(a) * b);
+}
+
+template <typename tbase>
+inline
+FDual<tbase> sin(const FDual<tbase>& f)
+{
+   return FDual<tbase>(sin(f.real()), f.dual() * cos(f.real()));
+}
+
+template <>
+inline
+FDual<double> sin(const FDual<double>& f)
+{
+   return FDual<double>(std::sin(f.real()), f.dual() * std::cos(f.real()));
+}
+
+template <typename tbase>
+inline
+FDual<tbase> sinh(const FDual<tbase>& f)
+{
+   return FDual<tbase>(sinh(f.real()), f.dual() * cosh(f.real()));
+}
+
+template <>
+inline
+FDual<double> sinh(const FDual<double>& f)
+{
+   return FDual<double>(std::sinh(f.real()), f.dual() * std::cosh(f.real()));
+}
+
+
+template <typename tbase>
+inline
+FDual<tbase> sqrt(const FDual<tbase>& f)
+{
+   tbase a = sqrt(f.real());
+   return FDual<tbase>(a, f.dual() / (tbase(2) * a));
+}
+
+template <>
+inline
+FDual<double> sqrt(const FDual<double>& f)
+{
+   double a = std::sqrt(f.real());
+   return FDual<double>(a, f.dual() / (double(2) * a));
+}
+
+
+template <typename tbase>
+inline
+FDual<tbase> tan(const FDual<tbase>& f)
+{
+   tbase a = tan(f.real());
+   return FDual<tbase>(a,f.dual() * (tbase(1) + a * a));
+}
+
+template <>
+inline
+FDual<double> tan(const FDual<double>& f)
+{
+   double a = std::tan(f.real());
+   return FDual<double>(a,f.dual() * (double(1) + a * a));
+}
+
+
+template <typename tbase>
+inline
+FDual<tbase> tanh(const FDual<tbase>& f)
+{
+   tbase a = tanh(f.real());
+   return FDual<tbase>(a, f.dual() * (tbase(1) - a * a));
+}
+
+template <>
+inline
+FDual<double> tanh(const FDual<double>& f)
+{
+   double a = std::tanh(f.real());
+   return FDual<double>(a, f.dual() * (double(1) - a * a));
+}
+
+}
+
+}
+
+
+
+
+#endif

linalg/linalg.hpp

@@ -28,6 +28,7 @@
 #include "solvers.hpp"
 #include "handle.hpp"
 #include "invariants.hpp"
+// #include "fdual.hpp"
 
 #ifdef MFEM_USE_SUNDIALS
 #include "sundials.hpp"

linalg/solvers.cpp

@@ -1641,6 +1641,138 @@ void NewtonSolver::Mult(const Vector &b, Vector &x) const
    final_norm = norm;
 }
 
+void LBFGSSolver::Mult(const Vector &b, Vector &x) const
+{
+   MFEM_VERIFY(oper != NULL, "the Operator is not set (use SetOperator).");
+
+   // Quadrature points that are checked for negative Jacobians etc.
+   Vector sk, rk, yk, rho, alpha;
+   DenseMatrix skM(width, m), ykM(width, m);
+
+   //r - r_{k+1}, c - descent direction
+   sk.SetSize(width);    //x_{k+1}-x_k
+   rk.SetSize(width);    //nabla(f(x_{k}))
+   yk.SetSize(width);    //r_{k+1}-r_{k}
+   rho.SetSize(m);       //1/(dot(yk,sk)
+   alpha.SetSize(m);    //rhok*sk'*c
+   int last_saved_id = -1;
+
+   int it;
+   double norm0, norm, norm_goal;
+   const bool have_b = (b.Size() == Height());
+
+   if (!iterative_mode)
+   {
+      x = 0.0;
+   }
+
+   // r = F(x)-b
+   oper->Mult(x, r);
+   if (have_b) { r -= b; }
+
+   c = r;           // initial descent direction
+
+   norm0 = norm = Norm(r);
+   norm_goal = std::max(rel_tol*norm, abs_tol);
+   for (it = 0; true; it++)
+   {
+      MFEM_ASSERT(IsFinite(norm), "norm = " << norm);
+      if (print_level >= 0)
+      {
+         mfem::out << "LBFGS iteration " <<  it
+                   << " : ||r|| = " << norm;
+         if (it > 0)
+         {
+            mfem::out << ", ||r||/||r_0|| = " << norm/norm0;
+         }
+         mfem::out << '\n';
+      }
+
+      if (norm <= norm_goal)
+      {
+         converged = 1;
+         break;
+      }
+
+      if (it >= max_iter)
+      {
+         converged = 0;
+         break;
+      }
+
+      rk = r;
+      const double c_scale = ComputeScalingFactor(x, b);
+      if (c_scale == 0.0)
+      {
+         converged = 0;
+         break;
+      }
+      add(x, -c_scale, c, x); //x_{k+1} = x_k - c_scale*c
+
+      ProcessNewState(x);
+
+      oper->Mult(x, r);
+      if (have_b)
+      {
+         r -= b;
+      }
+
+      //    LBFGS - construct descent direction
+      subtract(r, rk, yk);   // yk = r_{k+1} - r_{k}
+      sk = c; sk *= -c_scale; //sk = x_{k+1} - x_{k} = -c_scale*c
+      const double gamma = Dot(sk, yk)/Dot(yk, yk);
+
+      //  Save last m vectors
+      last_saved_id = (last_saved_id == m-1) ? 0 : last_saved_id+1;
+      skM.SetCol(last_saved_id, sk);
+      ykM.SetCol(last_saved_id, yk);
+
+      c = r;
+      for (int i = last_saved_id; i > -1; i--)
+      {
+         skM.GetColumn(i, sk);
+         ykM.GetColumn(i, yk);
+         rho(i) = 1./Dot(sk, yk);
+         alpha(i) = rho(i)*Dot(sk,c);
+         add(c, -alpha(i), yk, c);
+      }
+      if (it > m-1)
+      {
+         for (int i = m-1; i > last_saved_id; i--)
+         {
+            skM.GetColumn(i, sk);
+            ykM.GetColumn(i, yk);
+            rho(i) = 1./Dot(sk, yk);
+            alpha(i) = rho(i)*Dot(sk,c);
+            add(c, -alpha(i), yk, c);
+         }
+      }
+
+      c *= gamma;   // scale search direction
+      if (it > m-1)
+      {
+         for (int i = last_saved_id+1; i < m ; i++)
+         {
+            skM.GetColumn(i,sk);
+            ykM.GetColumn(i,yk);
+            double betai = rho(i)*Dot(yk, c);
+            add(c, alpha(i)-betai, sk, c);
+         }
+      }
+      for (int i = 0; i < last_saved_id+1 ; i++)
+      {
+         skM.GetColumn(i,sk);
+         ykM.GetColumn(i,yk);
+         double betai = rho(i)*Dot(yk, c);
+         add(c, alpha(i)-betai, sk, c);
+      }
+
+      norm = Norm(r);
+   }
+
+   final_iter = it;
+   final_norm = norm;
+}
 
 int aGMRES(const Operator &A, Vector &x, const Vector &b,
            const Operator &M, int &max_iter,

linalg/solvers.hpp

@@ -416,6 +416,32 @@ public:
    virtual void ProcessNewState(const Vector &x) const { }
 };
 
+/** L-BFGS method for solving F(x)=b for a given operator F, by minimizing
+    the norm of F(x) - b. Requires only the action of the operator F. */
+class LBFGSSolver : public NewtonSolver
+{
+protected:
+   int m = 10;
+
+public:
+   LBFGSSolver() : NewtonSolver() { }
+
+#ifdef MFEM_USE_MPI
+   LBFGSSolver(MPI_Comm _comm) : NewtonSolver(_comm) { }
+#endif
+
+   void SetHistorySize(int dim) { m = dim; }
+
+   /// Solve the nonlinear system with right-hand side @a b.
+   /** If `b.Size() != Height()`, then @a b is assumed to be zero. */
+   virtual void Mult(const Vector &b, Vector &x) const;
+
+   virtual void SetPreconditioner(Solver &pr)
+   { MFEM_WARNING("L-BFGS won't use the given preconditioner."); }
+   virtual void SetSolver(Solver &solver)
+   { MFEM_WARNING("L-BFGS won't use the given solver."); }
+};
+
 /** Adaptive restarted GMRES.
     m_max and m_min(=1) are the maximal and minimal restart parameters.
     m_step(=1) is the step to use for going from m_max and m_min.

linalg/taddensemat.hpp

@@ -0,0 +1,532 @@
+// Copyright (c) 2020, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+
+#ifndef TADDENSEMATRIX_H
+#define TADDENSEMATRIX_H
+
+
+#include "../config/config.hpp"
+#include "../general/globals.hpp"
+
+#include "tadvector.hpp"
+#include "densemat.hpp"
+
+
+
+namespace mfem
+{
+
+template<typename dtype>
+class TADDenseMatrix
+{
+private:
+   int height; ///< Dimension of the output / number of rows in the matrix.
+   int width;  ///< Dimension of the input / number of columns in the matrix.
+   dtype *data;
+   int capacity; // zero or negative capacity means we do not own the data.
+
+
+public:
+   /// Get the height (size of output) of the Operator. Synonym with NumRows().
+   inline int Height() const { return height; }
+   /** @brief Get the number of rows (size of output) of the Operator. Synonym
+       with Height(). */
+   inline int NumRows() const { return height; }
+
+   /// Get the width (size of input) of the Operator. Synonym with NumCols().
+   inline int Width() const { return width; }
+   /** @brief Get the number of columns (size of input) of the Operator. Synonym
+       with Width(). */
+   inline int NumCols() const { return width; }
+
+   /** Default constructor for TADDenseMatrix.
+      Sets data = NULL and height = width = 0. */
+   TADDenseMatrix()
+   {
+      data=nullptr;
+      capacity=0;
+      height=0;
+      width=0;
+   }
+
+   /// Copy constructor
+   template<typename idtype>
+   TADDenseMatrix(const TADDenseMatrix<idtype> &m)
+   {
+      height=m.GetHeight();
+      width=m.GetWidth();
+      const int hw = height * width;
+      if (hw > 0)
+      {
+         idtype* mdata=m.Data();
+         MFEM_ASSERT(mdata, "invalid source matrix");
+         data = new dtype[hw];
+         capacity = hw;
+         for (int i=0; i<hw; i++)
+         {
+            data[i]=mdata[i];
+         }
+      }
+      else
+      {
+         data = nullptr;
+         capacity = 0;
+         width=0;
+         height=0;
+      }
+   }
+
+   TADDenseMatrix(const DenseMatrix &m)
+   {
+      height=m.Height();
+      width=m.Width();
+      const int hw = height * width;
+      if (hw > 0)
+      {
+         double* mdata=m.Data();
+         MFEM_ASSERT(mdata, "invalid source matrix");
+         data = new dtype[hw];
+         capacity = hw;
+         for (int i=0; i<hw; i++)
+         {
+            data[i]=mdata[i];
+         }
+      }
+      else
+      {
+         data = nullptr;
+         capacity = 0;
+         width=0;
+         height=0;
+      }
+   }
+
+
+   /// Creates square matrix of size s.
+   explicit TADDenseMatrix(int s)
+   {
+      MFEM_ASSERT(s >= 0, "invalid DenseMatrix size: " << s);
+      height=s;
+      width=s;
+      capacity = s*s;
+      if (capacity > 0)
+      {
+         data = new dtype[capacity](); // init with zeroes
+      }
+      else
+      {
+         data = NULL;
+      }
+   }
+
+   /// Creates rectangular matrix of size m x n.
+   TADDenseMatrix(int m, int n)
+   {
+      MFEM_ASSERT(m >= 0 && n >= 0,
+                  "invalid DenseMatrix size: " << m << " x " << n);
+      height=m;
+      width=n;
+      capacity = m*n;
+      if (capacity > 0)
+      {
+         data = new dtype[capacity](); // init with zeroes
+      }
+      else
+      {
+         data = NULL;
+      }
+   }
+
+   TADDenseMatrix(const TADDenseMatrix<dtype> &mat, char ch)
+   {
+      height=mat.Width();
+      width=mat.Height();
+      capacity = height*width;
+      if (capacity > 0)
+      {
+         data = new dtype[capacity];
+
+         for (int i = 0; i < height; i++)
+         {
+            for (int j = 0; j < width; j++)
+            {
+               (*this)(i,j) = mat(j,i);
+            }
+         }
+      }
+      else
+      {
+         data = NULL;
+      }
+   }
+
+
+   /// Change the size of the DenseMatrix to s x s.
+   void SetSize(int s) { SetSize(s, s); }
+
+   /// Change the size of the DenseMatrix to h x w.
+   void SetSize(int h, int w)
+   {
+      MFEM_ASSERT(h >= 0 && w >= 0,
+                  "invalid DenseMatrix size: " << h << " x " << w);
+      if (Height() == h && Width() == w)
+      {
+         return;
+      }
+      height = h;
+      width = w;
+      const int hw = h*w;
+      if (hw > std::abs(capacity))
+      {
+         if (capacity > 0)
+         {
+            delete [] data;
+         }
+         capacity = hw;
+         data = new dtype[hw](); // init with zeroes
+      }
+   }
+
+   /// Returns the matrix data array.
+   inline dtype *Data() const { return data; }
+   /// Returns the matrix data array.
+   inline dtype *GetData() const { return data; }
+
+   inline bool OwnsData() const { return (capacity > 0); }
+
+   /// Returns reference to a_{ij}.
+   dtype& operator()(int i, int j)
+   {
+      MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
+      return data[i+j*height];
+   }
+
+   const dtype& operator()(int i, int j) const
+   {
+      MFEM_ASSERT(data && i >= 0 && i < height && j >= 0 && j < width, "");
+      return data[i+j*height];
+   }
+
+
+   dtype& Elem(int i, int j)
+   {
+      return (*this)(i,j);
+   }
+
+   const dtype& Elem(int i, int j) const
+   {
+      return (*this)(i,j);
+   }
+
+   void Mult(const dtype *x, dtype *y) const
+   {
+      if (width == 0)
+      {
+         for (int row = 0; row < height; row++)
+         {
+            y[row] = 0.0;
+         }
+         return;
+      }
+      dtype *d_col = data;
+      dtype x_col = x[0];
+      for (int row = 0; row < height; row++)
+      {
+         y[row] = x_col*d_col[row];
+      }
+      d_col += height;
+      for (int col = 1; col < width; col++)
+      {
+         x_col = x[col];
+         for (int row = 0; row < height; row++)
+         {
+            y[row] += x_col*d_col[row];
+         }
+         d_col += height;
+      }
+   }
+
+   void Mult(const TADVector<dtype> &x, TADVector<dtype> &y) const
+   {
+      MFEM_ASSERT(height == y.Size() && width == x.Size(),
+                  "incompatible dimensions");
+
+      Mult((const dtype *)x, (dtype *)y);
+   }
+
+   dtype operator *(const TADDenseMatrix<dtype> &m) const
+   {
+      MFEM_ASSERT(Height() == m.Height() && Width() == m.Width(),
+                  "incompatible dimensions");
+
+      const int hw = height * width;
+      dtype a = 0.0;
+      for (int i = 0; i < hw; i++)
+      {
+         a += data[i] * m.data[i];
+      }
+
+      return a;
+   }
+
+   void MultTranspose(const dtype *x, dtype *y) const
+   {
+      dtype *d_col = data;
+      for (int col = 0; col < width; col++)
+      {
+         double y_col = 0.0;
+         for (int row = 0; row < height; row++)
+         {
+            y_col += x[row]*d_col[row];
+         }
+         y[col] = y_col;
+         d_col += height;
+      }
+   }
+
+   void MultTranspose(const  TADVector<dtype> &x,  TADVector<dtype> &y) const
+   {
+      MFEM_ASSERT(height == x.Size() && width == y.Size(),
+                  "incompatible dimensions");
+
+      MultTranspose((const dtype *)x, (dtype *)y);
+   }
+
+
+   void Randomize(int seed)
+   {
+      // static unsigned int seed = time(0);
+      const double max = (double)(RAND_MAX) + 1.;
+
+      if (seed == 0)
+      {
+         seed = (int)time(0);
+      }
+
+      // srand(seed++);
+      srand((unsigned)seed);
+
+      for (int i = 0; i < capacity; i++)
+      {
+         data[i] = (dtype)(std::abs(rand()/max));
+      }
+   }
+
+   void RandomizeDiag(int seed)
+   {
+      // static unsigned int seed = time(0);
+      const double max = (double)(RAND_MAX) + 1.;
+
+      if (seed == 0)
+      {
+         seed = (int)time(0);
+      }
+
+      // srand(seed++);
+      srand((unsigned)seed);
+
+      for (int i = 0; i < std::min(height,width); i++)
+      {
+         Elem(i,i) = (dtype)(std::abs(rand()/max));
+      }
+   }
+
+
+   /// Creates n x n diagonal matrix with diagonal elements c
+   void Diag(dtype c, int n)
+   {
+      SetSize(n);
+
+      const int N = n*n;
+      for (int i = 0; i < N; i++)
+      {
+         data[i] = (dtype)0.0;
+      }
+      for (int i = 0; i < n; i++)
+      {
+         data[i*(n+1)] = c;
+      }
+   }
+   /// Creates n x n diagonal matrix with diagonal given by diag
+   template<typename itype>
+   void Diag(itype *diag, int n)
+   {
+      SetSize(n);
+
+      int i, N = n*n;
+      for (i = 0; i < N; i++)
+      {
+         data[i] = 0.0;
+      }
+      for (i = 0; i < n; i++)
+      {
+         data[i*(n+1)] = (dtype) diag[i];
+      }
+   }
+
+   /// (*this) = (*this)^t
+   void Transpose()
+   {
+      int i, j;
+      dtype t;
+
+      if (Width() == Height())
+      {
+         for (i = 0; i < Height(); i++)
+            for (j = i+1; j < Width(); j++)
+            {
+               t = (*this)(i,j);
+               (*this)(i,j) = (*this)(j,i);
+               (*this)(j,i) = t;
+            }
+      }
+      else
+      {
+         TADDenseMatrix<dtype> T(*this,'t');
+         (*this) = T;
+      }
+   }
+   /// (*this) = A^t
+   template<typename itype>
+   void Transpose(const TADDenseMatrix<itype> &A)
+   {
+      SetSize(A.Width(),A.Height());
+
+      for (int i = 0; i < Height(); i++)
+         for (int j = 0; j < Width(); j++)
+         {
+            (*this)(i,j) = (dtype) A(j,i);
+         }
+   }
+
+
+   /// (*this) = 1/2 ((*this) + (*this)^t)
+   void Symmetrize()
+   {
+#ifdef MFEM_DEBUG
+      if (Width() != Height())
+      {
+         mfem_error("DenseMatrix::Symmetrize() : not a square matrix!");
+      }
+#endif
+
+      for (int i = 0; i < Height(); i++)
+         for (int j = 0; j < i; j++)
+         {
+            dtype a = 0.5 * ((*this)(i,j) + (*this)(j,i));
+            (*this)(j,i) = (*this)(i,j) = a;
+         }
+   }
+
+   void Lump()
+   {
+      for (int i = 0; i < Height(); i++)
+      {
+         dtype L = 0.0;
+         for (int j = 0; j < Width(); j++)
+         {
+            L += (*this)(i, j);
+            (*this)(i, j) = (dtype) 0.0;
+         }
+         (*this)(i, i) = L;
+      }
+   }
+
+
+
+
+};
+
+
+
+template<typename dtype>
+void CalcAdjugate(const TADDenseMatrix<dtype> &a, TADDenseMatrix<dtype> &adja)
+{
+#ifdef MFEM_DEBUG
+   if (a.Width() > a.Height() || a.Width() < 1 || a.Height() > 3)
+   {
+      mfem_error("CalcAdjugate(...)");
+   }
+   if (a.Width() != adja.Height() || a.Height() != adja.Width())
+   {
+      mfem_error("CalcAdjugate(...)");
+   }
+#endif
+
+   if (a.Width() < a.Height())
+   {
+      const dtype *d = a.Data();
+      dtype *ad = adja.Data();
+      if (a.Width() == 1)
+      {
+         // N x 1, N = 2,3
+         ad[0] = d[0];
+         ad[1] = d[1];
+         if (a.Height() == 3)
+         {
+            ad[2] = d[2];
+         }
+      }
+      else
+      {
+         // 3 x 2
+         double e, g, f;
+         e = d[0]*d[0] + d[1]*d[1] + d[2]*d[2];
+         g = d[3]*d[3] + d[4]*d[4] + d[5]*d[5];
+         f = d[0]*d[3] + d[1]*d[4] + d[2]*d[5];
+
+         ad[0] = d[0]*g - d[3]*f;
+         ad[1] = d[3]*e - d[0]*f;
+         ad[2] = d[1]*g - d[4]*f;
+         ad[3] = d[4]*e - d[1]*f;
+         ad[4] = d[2]*g - d[5]*f;
+         ad[5] = d[5]*e - d[2]*f;
+      }
+      return;
+   }
+
+   if (a.Width() == 1)
+   {
+      adja(0,0) = (dtype)1.0;
+   }
+   else if (a.Width() == 2)
+   {
+      adja(0,0) =  a(1,1);
+      adja(0,1) = -a(0,1);
+      adja(1,0) = -a(1,0);
+      adja(1,1) =  a(0,0);
+   }
+   else
+   {
+      adja(0,0) = a(1,1)*a(2,2)-a(1,2)*a(2,1);
+      adja(0,1) = a(0,2)*a(2,1)-a(0,1)*a(2,2);
+      adja(0,2) = a(0,1)*a(1,2)-a(0,2)*a(1,1);
+
+      adja(1,0) = a(1,2)*a(2,0)-a(1,0)*a(2,2);
+      adja(1,1) = a(0,0)*a(2,2)-a(0,2)*a(2,0);
+      adja(1,2) = a(0,2)*a(1,0)-a(0,0)*a(1,2);
+
+      adja(2,0) = a(1,0)*a(2,1)-a(1,1)*a(2,0);
+      adja(2,1) = a(0,1)*a(2,0)-a(0,0)*a(2,1);
+      adja(2,2) = a(0,0)*a(1,1)-a(0,1)*a(1,0);
+   }
+}
+
+
+
+
+
+
+
+}
+
+#endif

linalg/tadvector.hpp

@@ -0,0 +1,687 @@
+// Copyright (c) 2020, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_TADVECTOR
+#define MFEM_TADVECTOR
+
+#include "../general/mem_manager.hpp"
+#include "vector.hpp"
+
+#include <cmath>
+#include <iostream>
+#include <limits>
+#if defined(_MSC_VER) && (_MSC_VER < 1800)
+#include <float.h>
+#define isfinite _finite
+#endif
+
+
+namespace mfem
+{
+
+
+/// Vector data type.
+template<typename dtype>
+class TADVector
+{
+protected:
+
+   Memory<dtype> data;
+   int size;
+
+public:
+
+   /// Default constructor for Vector. Sets size = 0 and data = NULL.
+   TADVector() { data.Reset(); size = 0; }
+
+   /// Copy constructor. Allocates a new data array and copies the data.
+   TADVector(const TADVector<dtype> &v)
+   {
+      const int s = v.Size();
+      if (s > 0)
+      {
+         size = s;
+         data.New(s);
+         for (int i=0; i<s; i++)
+         {
+            data[i]=v[i];
+         }
+      }
+      else
+      {
+         size = 0;
+         data.Reset();
+      }
+   }
+
+   TADVector(const Vector &v)
+   {
+      const int s = v.Size();
+      if (s > 0)
+      {
+         size = s;
+         data.New(s);
+         for (int i=0; i<s; i++)
+         {
+            data[i]=v[i];
+         }
+      }
+      else
+      {
+         size = 0;
+         data.Reset();
+      }
+   }
+
+   /// @brief Creates vector of size s.
+   /// @warning Entries are not initialized to zero!
+   explicit TADVector(int s)
+   {
+      if (s > 0)
+      {
+         size = s;
+         data.New(size);
+      }
+      else
+      {
+         size = 0;
+         data.Reset();
+      }
+   }
+
+   /// Creates a vector referencing an array of doubles, owned by someone else.
+   /** The pointer @a _data can be NULL. The data array can be replaced later
+       with SetData(). */
+   TADVector(dtype *_data, int _size)
+   { data.Wrap(_data, _size, false); size = _size; }
+
+   /// Create a Vector of size @a size_ using MemoryType @a mt.
+   TADVector(int size_, MemoryType mt)
+      : data(size_, mt), size(size_) { }
+
+   /// Enable execution of Vector operations using the mfem::Device.
+   /** The default is to use Backend::CPU (serial execution on each MPI rank),
+       regardless of the mfem::Device configuration.
+
+       When appropriate, MFEM functions and class methods will enable the use
+       of the mfem::Device for their Vector parameters.
+
+       Some derived classes, e.g. GridFunction, enable the use of the
+       mfem::Device by default. */
+   void UseDevice(bool use_dev) const { data.UseDevice(use_dev); }
+
+   /// Return the device flag of the Memory object used by the Vector
+   bool UseDevice() const { return data.UseDevice(); }
+
+   /// Reads a vector from multiple files
+   void Load(std::istream ** in, int np, int * dim)
+   {
+      int i, j, s;
+
+      s = 0;
+      for (i = 0; i < np; i++)
+      {
+         s += dim[i];
+      }
+      SetSize(s);
+
+      int p = 0;
+      double tmpd;
+      for (i = 0; i < np; i++)
+      {
+         for (j = 0; j < dim[i]; j++)
+         {
+            *in[i] >> tmpd;
+            data[p++]=dtype(tmpd);
+         }
+      }
+   }
+
+   /// Load a vector from an input stream.
+   void Load(std::istream &in, int Size)
+   {
+      SetSize(Size);
+      double tmpd;
+      for (int i = 0; i < size; i++)
+      {
+         in >> tmpd;
+         data[i]=dtype(tmpd);
+      }
+   }
+
+   /// Load a vector from an input stream, reading the size from the stream.
+   void Load(std::istream &in) { int s; in >> s; Load(in, s); }
+
+   /// @brief Resize the vector to size @a s.
+   /** If the new size is less than or equal to Capacity() then the internal
+       data array remains the same. Otherwise, the old array is deleted, if
+       owned, and a new array of size @a s is allocated without copying the
+       previous content of the Vector.
+       @warning In the second case above (new size greater than current one),
+       the vector will allocate new data array, even if it did not own the
+       original data! Also, new entries are not initialized! */
+   void SetSize(int s)
+   {
+      if (s == size)
+      {
+         return;
+      }
+      if (s <= data.Capacity())
+      {
+         size = s;
+         return;
+      }
+      // preserve a valid MemoryType and device flag
+      const MemoryType mt = data.GetMemoryType();
+      const bool use_dev = data.UseDevice();
+      data.Delete();
+      size = s;
+      data.New(s, mt);
+      data.UseDevice(use_dev);
+   }
+
+   /// Resize the vector to size @a s using MemoryType @a mt.
+   void SetSize(int s, MemoryType mt)
+   {
+      if (mt == data.GetMemoryType())
+      {
+         if (s == size)
+         {
+            return;
+         }
+         if (s <= data.Capacity())
+         {
+            size = s;
+            return;
+         }
+      }
+      const bool use_dev = data.UseDevice();
+      data.Delete();
+      if (s > 0)
+      {
+         data.New(s, mt);
+         size = s;
+      }
+      else
+      {
+         data.Reset();
+         size = 0;
+      }
+      data.UseDevice(use_dev);
+   }
+
+   /// Set the Vector data.
+   /// @warning This method should be called only when OwnsData() is false.
+   void SetData(dtype *d) { data.Wrap(d, data.Capacity(), false); }
+
+   /// Set the Vector data and size.
+   /** The Vector does not assume ownership of the new data. The new size is
+       also used as the new Capacity().
+       @warning This method should be called only when OwnsData() is false.
+       @sa NewDataAndSize(). */
+   void SetDataAndSize(dtype *d, int s)
+   { data.Wrap(d, s, false); size = s; }
+
+   /// Set the Vector data and size, deleting the old data, if owned.
+   /** The Vector does not assume ownership of the new data. The new size is
+       also used as the new Capacity().
+       @sa SetDataAndSize(). */
+   void NewDataAndSize(dtype *d, int s)
+   {
+      data.Delete();
+      SetDataAndSize(d, s);
+   }
+
+   /// Reset the Vector to use the given external Memory @a mem and size @a s.
+   /** If @a own_mem is false, the Vector will not own any of the pointers of
+       @a mem.
+       @sa NewDataAndSize(). */
+   void NewMemoryAndSize(const Memory<dtype> &mem, int s, bool own_mem)
+   {
+      data.Delete();
+      size = s;
+      data = mem;
+      if (!own_mem) { data.ClearOwnerFlags(); }
+
+   }
+
+   /// Reset the Vector to be a reference to a sub-vector of @a base.
+   inline void MakeRef(TADVector<dtype> &base, int offset, int size_)
+   {
+      data.Delete();
+      size = size_;
+      data.MakeAlias(base.GetMemory(), offset, size_);
+   }
+
+   /** @brief Reset the Vector to be a reference to a sub-vector of @a base
+       without changing its current size. */
+   inline void MakeRef(TADVector<dtype> &base, int offset)
+   {
+      data.Delete();
+      data.MakeAlias(base.GetMemory(), offset, size);
+   }
+
+   /// Set the Vector data (host pointer) ownership flag.
+   inline void MakeDataOwner() const { data.SetHostPtrOwner(true); }
+
+   /// Destroy a vector
+   void Destroy()
+   {
+      data.Delete();
+      size = 0;
+      data.Reset();
+   }
+
+   /// Returns the size of the vector.
+   inline int Size() const { return size; }
+
+   /// Return the size of the currently allocated data array.
+   /** It is always true that Capacity() >= Size(). */
+   inline int Capacity() const { return data.Capacity(); }
+
+   /// Return a pointer to the beginning of the Vector data.
+   /** @warning This method should be used with caution as it gives write access
+       to the data of const-qualified Vector%s. */
+   inline dtype *GetData() const
+   { return const_cast<dtype*>((const dtype*)data); }
+
+   /// Conversion to `double *`.
+   /** @note This conversion function makes it possible to use [] for indexing
+       in addition to the overloaded operator()(int). */
+   inline operator dtype *() { return data; }
+
+   /// Conversion to `const double *`.
+   /** @note This conversion function makes it possible to use [] for indexing
+       in addition to the overloaded operator()(int). */
+   inline operator const dtype *() const { return data; }
+
+   /// Return a reference to the Memory object used by the Vector.
+   Memory<dtype> &GetMemory() { return data; }
+
+   /** @brief Return a reference to the Memory object used by the Vector, const
+       version. */
+   const Memory<dtype> &GetMemory() const { return data; }
+
+   /// Update the memory location of the vector to match @a v.
+   void SyncMemory(const TADVector<dtype> &v) { GetMemory().Sync(v.GetMemory()); }
+
+   /// Update the alias memory location of the vector to match @a v.
+   void SyncAliasMemory(const TADVector<dtype> &v)
+   { GetMemory().SyncAlias(v.GetMemory(),Size()); }
+
+   /// Read the Vector data (host pointer) ownership flag.
+   inline bool OwnsData() const { return data.OwnsHostPtr(); }
+
+   /// Changes the ownership of the data; after the call the Vector is empty
+   inline void StealData(dtype **p)
+   { *p = data; data.Reset(); size = 0; }
+
+   /// Changes the ownership of the data; after the call the Vector is empty
+   inline dtype *StealData() { dtype *p; StealData(&p); return p; }
+
+   /// Access Vector entries. Index i = 0 .. size-1.
+   dtype &Elem(int i)
+   {
+      return operator()(i);
+   }
+   /// Read only access to Vector entries. Index i = 0 .. size-1.
+   const double &Elem(int i) const
+   {
+      return operator()(i);
+   }
+
+   /// Access Vector entries using () for 0-based indexing.
+   /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
+   inline double &operator()(int i)
+   {
+      MFEM_ASSERT(data && i >= 0 && i < size,
+                  "index [" << i << "] is out of range [0," << size << ")");
+
+      return data[i];
+   }
+
+   /// Read only access to Vector entries using () for 0-based indexing.
+   /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
+   inline const double &operator()(int i) const
+   {
+      MFEM_ASSERT(data && i >= 0 && i < size,
+                  "index [" << i << "] is out of range [0," << size << ")");
+
+      return data[i];
+   }
+
+   /// Dot product with a `dtype *` array.
+   dtype operator*(const dtype *v) const
+   {
+      dtype dot = 0.0;
+#ifdef MFEM_USE_LEGACY_OPENMP
+      #pragma omp parallel for reduction(+:dot)
+#endif
+      for (int i = 0; i < size; i++)
+      {
+         dot += data[i] * v[i];
+      }
+      return dot;
+   }
+
+
+   /// Return the inner-product.
+   dtype operator*(const TADVector<dtype> &v) const
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      dtype dot = 0.0;
+      for (int i = 0; i < size; i++)
+      {
+         dot += data[i] * v[i];
+      }
+      return dot;
+   }
+
+   dtype operator*(const Vector &v) const
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      dtype dot = 0.0;
+      for (int i = 0; i < size; i++)
+      {
+         dot += data[i] * v[i];
+      }
+      return dot;
+   }
+
+   /// Copy Size() entries from @a v.
+   TADVector<dtype> &operator=(const dtype *v)
+   {
+      for (int i=0; i<size; i++)
+      {
+         data[i]=v[i];
+      }
+      return *this;
+   }
+
+   /// Copy assignment.
+   /** @note Defining this method overwrites the implicitly defined copy
+       assignemnt operator. */
+   TADVector<dtype> &operator=(const TADVector<dtype> &v)
+   {
+      SetSize(v.Size());
+      for (int i=0; i<size; i++)
+      {
+         data[i]=v[i];
+      }
+      return *this;
+   }
+
+   TADVector<dtype> &operator=(const Vector &v)
+   {
+      SetSize(v.Size());
+      for (int i=0; i<size; i++)
+      {
+         data[i]=v[i];
+      }
+      return *this;
+   }
+
+   /// Redefine '=' for vector = constant.
+   template<typename ivtype>
+   TADVector &operator=(ivtype value)
+   {
+      for (int i=0; i<size; i++)
+      {
+         data[i]=value;
+      }
+      return *this;
+   }
+
+   template<typename ivtype>
+   TADVector &operator*=(ivtype c)
+   {
+      for (int i=0; i<size; i++)
+      {
+         data[i]=data[i]*c;
+      }
+      return *this;
+   }
+
+   template<typename ivtype>
+   TADVector &operator/=(ivtype c)
+   {
+      for (int i=0; i<size; i++)
+      {
+         data[i]=data[i]/c;
+      }
+      return *this;
+   }
+
+   template<typename ivtype>
+   TADVector &operator-=(ivtype c)
+   {
+      for (int i=0; i<size; i++)
+      {
+         data[i]=data[i]-c;
+      }
+      return *this;
+   }
+
+   TADVector &operator-=(const TADVector<dtype> &v)
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      for (int i=0; i<size; i++)
+      {
+         data[i]=data[i]-v[i];
+      }
+      return *this;
+   }
+
+   TADVector &operator+=(const TADVector<dtype> &v)
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      for (int i=0; i<size; i++)
+      {
+         data[i]=data[i]+v[i];
+      }
+      return *this;
+   }
+
+   /// (*this) += a * Va
+   template<typename ivtype, typename vtype>
+   TADVector &Add(const ivtype a, const vtype &v)
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      for (int i=0; i<size; i++)
+      {
+         data[i]=data[i]+a*v[i];
+      }
+      return *this;
+   }
+
+   /// (*this) = a * x
+   template<typename ivtype, typename vtype>
+   TADVector &Set(const ivtype a, const vtype &v)
+   {
+      MFEM_ASSERT(size == v.Size(), "incompatible Vectors!");
+      for (int i=0; i<size; i++)
+      {
+         data[i]=a*v[i];
+      }
+      return *this;
+   }
+
+   template<typename vtype>
+   void SetVector(const vtype &v, int offset)
+   {
+      MFEM_ASSERT(v.Size() + offset <= size, "invalid sub-vector");
+      for (int i = 0; i < size; i++)
+      {
+         data[i+offset] = v[i];
+      }
+   }
+
+   /// (*this) = -(*this)
+   void Neg()
+   {
+      for (int i = 0; i < size; i++)
+      {
+         data[i]=-data[i];
+      }
+   }
+
+   /// Swap the contents of two Vectors
+   inline void Swap(TADVector &other)
+   {
+      Swap(data, other.data);
+      Swap(size, other.size);
+   }
+
+   /// Set v = v1 + v2.
+   template<typename vtype1, typename vtype2>
+   friend void add(const vtype1 &v1, const vtype2 &v2, TADVector<dtype> &v)
+   {
+      MFEM_ASSERT(v1.Size() == v.Size(), "incompatible Vectors!");
+      MFEM_ASSERT(v2.Size() == v.Size(), "incompatible Vectors!");
+      for (int i=0; i<v.Size(); i++)
+      {
+         v[i]=v1[i]+v2[i];
+      }
+   }
+
+   /// Set v = v1 + alpha * v2.
+   template<typename vtype1, typename ivtype, typename vtype2>
+   friend void add(const vtype1 &v1, ivtype alpha, const vtype2 &v2,
+                   TADVector<dtype> &v)
+   {
+      MFEM_ASSERT(v1.Size() == v.Size(), "incompatible Vectors!");
+      MFEM_ASSERT(v2.Size() == v.Size(), "incompatible Vectors!");
+      for (int i=0; i<v.Size(); i++)
+      {
+         v[i]=v1[i]+alpha*v2[i];
+      }
+   }
+
+
+   /// Destroys vector.
+   ~TADVector()
+   {
+      data.Delete();
+   }
+
+
+   /// Prints vector to stream out.
+   void Print(std::ostream &out = mfem::out, int width = 8) const
+   {
+      if (!size) { return; }
+      data.Read(MemoryClass::HOST, size);
+      for (int i = 0; 1; )
+      {
+         out << data[i];
+         i++;
+         if (i == size)
+         {
+            break;
+         }
+         if ( i % width == 0 )
+         {
+            out << '\n';
+         }
+         else
+         {
+            out << ' ';
+         }
+      }
+      out << '\n';
+   }
+
+   /// Set random values in the vector.
+   void Randomize(int seed = 0)
+   {
+      // static unsigned int seed = time(0);
+      const double max = (double)(RAND_MAX) + 1.;
+
+      if (seed == 0)
+      {
+         seed = (int)time(0);
+      }
+
+      // srand(seed++);
+      srand((unsigned)seed);
+
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = std::abs(rand()/max);
+      }
+   }
+   /// Returns the l2 norm of the vector.
+   dtype Norml2() const
+   {
+      // Scale entries of Vector on the fly, using algorithms from
+      // std::hypot() and LAPACK's drm2. This scaling ensures that the
+      // argument of each call to std::pow is <= 1 to avoid overflow.
+      if (0 == size)
+      {
+         return 0.0;
+      } // end if 0 == size
+
+      if (1 == size)
+      {
+         return std::abs(data[0]);
+      } // end if 1 == size
+
+      dtype scale = 0.0;
+      dtype sum = 0.0;
+
+      for (int i = 0; i < size; i++)
+      {
+         if (data[i] != 0.0)
+         {
+            const dtype absdata = abs(data[i]);
+            if (scale <= absdata)
+            {
+               const dtype sqr_arg = scale / absdata;
+               sum = 1.0 + sum * (sqr_arg * sqr_arg);
+               scale = absdata;
+               continue;
+            } // end if scale <= absdata
+            const dtype sqr_arg = absdata / scale;
+            sum += (sqr_arg * sqr_arg); // else scale > absdata
+         } // end if data[i] != 0
+      }
+      return scale * sqrt(sum);
+   }
+
+   /// Returns the l_infinity norm of the vector.
+   dtype Normlinf() const
+   {
+      dtype max = 0.0;
+      for (int i = 0; i < size; i++)
+      {
+         max = max(abs(data[i]), max);
+      }
+      return max;
+   }
+   /// Returns the l_1 norm of the vector.
+   dtype Norml1() const
+   {
+      dtype sum = 0.0;
+      for (int i = 0; i < size; i++)
+      {
+         sum += abs(data[i]);
+      }
+      return sum;
+   }
+
+
+};
+
+
+
+} // namespace mfem
+
+#endif
+

mesh/mesh.cpp

@@ -859,35 +859,41 @@ FaceElementTransformations *Mesh::GetFaceElementTransformations(int FaceNo,
 {
    FaceInfo &face_info = faces_info[FaceNo];
 
-   FaceElemTr.SetConfigurationMask(0);
+   int cmask = 0;
+   FaceElemTr.SetConfigurationMask(cmask);
    FaceElemTr.Elem1 = NULL;
    FaceElemTr.Elem2 = NULL;
 
    // setup the transformation for the first element
    FaceElemTr.Elem1No = face_info.Elem1No;
-   if (mask & 1)
+   if (mask & FaceElementTransformations::HAVE_ELEM1)
    {
       GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
       FaceElemTr.Elem1 = &Transformation;
+      cmask |= 1;
    }
 
    //  setup the transformation for the second element
    //     return NULL in the Elem2 field if there's no second element, i.e.
    //     the face is on the "boundary"
    FaceElemTr.Elem2No = face_info.Elem2No;
-   if ((mask & 2) && FaceElemTr.Elem2No >= 0)
+   if ((mask & FaceElementTransformations::HAVE_ELEM2) &&
+       FaceElemTr.Elem2No >= 0)
    {
 #ifdef MFEM_DEBUG
-      if (NURBSext && (mask & 1)) { MFEM_ABORT("NURBS mesh not supported!"); }
+      if (NURBSext && (mask & FaceElementTransformations::HAVE_ELEM1))
+      { MFEM_ABORT("NURBS mesh not supported!"); }
 #endif
       GetElementTransformation(FaceElemTr.Elem2No, &Transformation2);
       FaceElemTr.Elem2 = &Transformation2;
+      cmask |= 2;
    }
 
    // setup the face transformation
-   if (mask & 16)
+   if (mask & FaceElementTransformations::HAVE_FACE)
    {
       GetFaceTransformation(FaceNo, &FaceElemTr);
+      cmask |= 16;
    }
    else
    {
@@ -896,13 +902,15 @@ FaceElementTransformations *Mesh::GetFaceElementTransformations(int FaceNo,
 
    // setup Loc1 & Loc2
    int face_type = GetFaceElementType(FaceNo);
-   if (mask & 4)
+   if (mask & FaceElementTransformations::HAVE_LOC1)
    {
       int elem_type = GetElementType(face_info.Elem1No);
       GetLocalFaceTransformation(face_type, elem_type,
                                  FaceElemTr.Loc1.Transf, face_info.Elem1Inf);
+      cmask |= 4;
    }
-   if ((mask & 8) && FaceElemTr.Elem2No >= 0)
+   if ((mask & FaceElementTransformations::HAVE_LOC2) &&
+       FaceElemTr.Elem2No >= 0)
    {
       int elem_type = GetElementType(face_info.Elem2No);
       GetLocalFaceTransformation(face_type, elem_type,
@@ -913,9 +921,10 @@ FaceElementTransformations *Mesh::GetFaceElementTransformations(int FaceNo,
       {
          ApplyLocalSlaveTransformation(FaceElemTr, face_info, false);
       }
+      cmask |= 8;
    }
 
-   FaceElemTr.SetConfigurationMask(mask);
+   FaceElemTr.SetConfigurationMask(cmask);
 
    // This check can be useful for internal debugging, however it will fail on
    // periodic boundary faces, so we keep it disabled in general.
@@ -1000,7 +1009,7 @@ FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
    {
       return NULL;
    }
-   tr = GetFaceElementTransformations(fn);
+   tr = GetFaceElementTransformations(fn, 21);
    tr->Attribute = boundary[BdrElemNo]->GetAttribute();
    tr->ElementNo = BdrElemNo;
    tr->ElementType = ElementTransformation::BDR_FACE;

mesh/mesh.hpp

@@ -961,7 +961,7 @@ public:
    /// Returns the transformation defining the given face element
    ElementTransformation *GetEdgeTransformation(int EdgeNo);
 
-   /// Returns (a pointer to a structure containing) the following data:
+   /// Returns (a pointer to an object containing) the following data:
    ///
    /// 1) Elem1No - the index of the first element that contains this face this
    ///    is the element that has the same outward unit normal vector as the
@@ -989,6 +989,8 @@ public:
    /// The mask specifies which fields in the structure to return:
    ///    mask & 1 - Elem1, mask & 2 - Elem2
    ///    mask & 4 - Loc1, mask & 8 - Loc2, mask & 16 - Face.
+   /// These mask values are defined in the ConfigMasks enum type as part of the
+   /// FaceElementTransformations class in fem/eltrans.hpp.
    FaceElementTransformations *GetFaceElementTransformations(int FaceNo,
                                                              int mask = 31);
 

mesh/pmesh.cpp

@@ -1690,6 +1690,7 @@ void ParMesh::GetFaceNbrElementTransformation(
 
    ElTr->Attribute = elem->GetAttribute();
    ElTr->ElementNo = NumOfElements + i;
+   ElTr->ElementType = ElementTransformation::ELEMENT;
 
    if (Nodes == NULL)
    {
@@ -2402,6 +2403,11 @@ GetSharedFaceTransformations(int sf, bool fill2)
    bool is_slave = Nonconforming() && IsSlaveFace(face_info);
    bool is_ghost = Nonconforming() && FaceNo >= GetNumFaces();
 
+   int mask = 0;
+   FaceElemTr.SetConfigurationMask(0);
+   FaceElemTr.Elem1 = NULL;
+   FaceElemTr.Elem2 = NULL;
+
    NCFaceInfo* nc_info = NULL;
    if (is_slave) { nc_info = &nc_faces_info[face_info.NCFace]; }
 
@@ -2413,13 +2419,21 @@ GetSharedFaceTransformations(int sf, bool fill2)
    FaceElemTr.Elem1No = face_info.Elem1No;
    GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
    FaceElemTr.Elem1 = &Transformation;
+   mask |= FaceElementTransformations::HAVE_ELEM1;
 
    // setup the transformation for the second (neighbor) element
+   int Elem2NbrNo;
    if (fill2)
    {
-      FaceElemTr.Elem2No = -1 - face_info.Elem2No;
-      GetFaceNbrElementTransformation(FaceElemTr.Elem2No, &Transformation2);
+      Elem2NbrNo = -1 - face_info.Elem2No;
+      // Store the "shifted index" for element 2 in FaceElemTr.Elem2No.
+      // `Elem2NbrNo` is the index of the face neighbor (starting from 0),
+      // and `FaceElemTr.Elem2No` will be offset by the number of (local)
+      // elements in the mesh.
+      FaceElemTr.Elem2No = NumOfElements + Elem2NbrNo;
+      GetFaceNbrElementTransformation(Elem2NbrNo, &Transformation2);
       FaceElemTr.Elem2 = &Transformation2;
+      mask |= FaceElementTransformations::HAVE_ELEM2;
    }
    else
    {
@@ -2431,6 +2445,7 @@ GetSharedFaceTransformations(int sf, bool fill2)
    {
       GetFaceTransformation(FaceNo, &FaceElemTr);
       // NOTE: The above call overwrites FaceElemTr.Loc1
+      mask |= FaceElementTransformations::HAVE_FACE;
    }
    else
    {
@@ -2441,12 +2456,14 @@ GetSharedFaceTransformations(int sf, bool fill2)
    int elem_type = GetElementType(face_info.Elem1No);
    GetLocalFaceTransformation(face_type, elem_type, FaceElemTr.Loc1.Transf,
                               face_info.Elem1Inf);
+   mask |= FaceElementTransformations::HAVE_LOC1;
 
    if (fill2)
    {
-      elem_type = face_nbr_elements[FaceElemTr.Elem2No]->GetType();
+      elem_type = face_nbr_elements[Elem2NbrNo]->GetType();
       GetLocalFaceTransformation(face_type, elem_type, FaceElemTr.Loc2.Transf,
                                  face_info.Elem2Inf);
+      mask |= FaceElementTransformations::HAVE_LOC2;
    }
 
    // adjust Loc1 or Loc2 of the master face if this is a slave face
@@ -2463,9 +2480,10 @@ GetSharedFaceTransformations(int sf, bool fill2)
    if (is_ghost)
    {
       GetGhostFaceTransformation(&FaceElemTr, face_type, face_geom);
+      mask |= FaceElementTransformations::HAVE_FACE;
    }
 
-   FaceElemTr.SetConfigurationMask(fill2 ? 31 : 21);
+   FaceElemTr.SetConfigurationMask(mask);
 
    // This check can be useful for internal debugging, however it will fail on
    // periodic boundary faces, so we keep it disabled in general.

miniapps/electromagnetics/joule_solver.cpp

@@ -906,7 +906,7 @@ double JouleHeatingCoefficient::Eval(ElementTransformation &T,
 {
    Vector E;
    double thisSigma;
-   E_gf.GetVectorValue(T.ElementNo, ip, E);
+   E_gf.GetVectorValue(T, ip, E);
    thisSigma = sigma.Eval(T, ip);
    return thisSigma*(E*E);
 }

miniapps/meshing/mesh-optimizer.cpp

@@ -32,27 +32,29 @@
 // Compile with: make mesh-optimizer
 //
 // Sample runs:
-//   Adapted analytic Hessian:
+//   Adapted analytic shape:
 //     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 4 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
-//   Adapted analytic Hessian with size+orientation:
+//   Adapted analytic size+orientation:
 //     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd
-//   Adapted analytic Hessian with shape+size+orientation
-//     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 87 -tid 4 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd
+//   Adapted analytic shape+orientation:
+//     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd
+//
 //   Adapted discrete size:
 //     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 5 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
 //     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 5 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8 -cmb 2 -nor
-//
-//   Adapted size+aspect ratio to discrete material indicator
+//   Adapted discrete size+aspect_ratio:
 //     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100  -ls 2 -li 100 -bnd -qt 1 -qo 8
-//   Adapted discrete size+orientation (requires GSLIB)
+//   Adapted discrete size+orientation (requires GSLIB):
 //   * mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 8 -ni 100  -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd -ae 1
-//   Adapted discrete aspect-ratio+orientation (requires GSLIB)
-//   * mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 87 -tid 8 -ni 10  -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd -ae 1
-//   Adapted discrete aspect ratio (3D)
+//   Adapted discrete aspect-ratio+orientation (requires GSLIB):
+//   * mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 8 -ni 10  -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd -ae 1
+//   Adapted discrete aspect ratio (3D):
 //     mesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 302 -tid 7 -ni 20  -ls 2 -li 100 -bnd -qt 1 -qo 8
 //
 //   Adaptive limiting:
-//     mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5 -ae 0
+//     mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5
+//   Adaptive limiting through the L-BFGS solver:
+//     mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 400 -qo 5 -nor -vl 1 -alc 0.5 -st 1
 //   Adaptive limiting through FD (requires GSLIB):
 //   * mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5 -fd -ae 1
 //
@@ -100,8 +102,9 @@ int main(int argc, char *argv[])
    double adapt_lim_const = 0.0;
    int quad_type         = 1;
    int quad_order        = 8;
-   int newton_iter       = 10;
-   double newton_rtol    = 1e-10;
+   int solver_type       = 0;
+   int solver_iter       = 10;
+   double solver_rtol    = 1e-10;
    int lin_solver        = 2;
    int max_lin_iter      = 100;
    bool move_bnd         = true;
@@ -161,9 +164,11 @@ int main(int argc, char *argv[])
                   "3: Closed uniform points");
    args.AddOption(&quad_order, "-qo", "--quad_order",
                   "Order of the quadrature rule.");
-   args.AddOption(&newton_iter, "-ni", "--newton-iters",
+   args.AddOption(&solver_type, "-st", "--solver-type",
+                  " Type of solver: (default) 0: Newton, 1: LBFGS");
+   args.AddOption(&solver_iter, "-ni", "--newton-iters",
                   "Maximum number of Newton iterations.");
-   args.AddOption(&newton_rtol, "-rtol", "--newton-rel-tolerance",
+   args.AddOption(&solver_rtol, "-rtol", "--newton-rel-tolerance",
                   "Relative tolerance for the Newton solver.");
    args.AddOption(&lin_solver, "-ls", "--lin-solver",
                   "Linear solver: 0 - l1-Jacobi, 1 - CG, 2 - MINRES.");
@@ -313,7 +318,7 @@ int main(int argc, char *argv[])
       case 56: metric = new TMOP_Metric_056; break;
       case 58: metric = new TMOP_Metric_058; break;
       case 77: metric = new TMOP_Metric_077; break;
-      case 87: metric = new TMOP_Metric_SS2D; break;
+      case 85: metric = new TMOP_Metric_085; break;
       case 211: metric = new TMOP_Metric_211; break;
       case 252: metric = new TMOP_Metric_252(tauval); break;
       case 301: metric = new TMOP_Metric_301; break;
@@ -511,7 +516,7 @@ int main(int argc, char *argv[])
             tc->SetSerialDiscreteTargetSize(size);
          }
 
-         if (metric_id == 87)
+         if (metric_id == 85)
          {
             FunctionCoefficient aspr_coeff(discrete_aspr_2d);
             aspr.ProjectCoefficient(aspr_coeff);
@@ -737,43 +742,26 @@ int main(int argc, char *argv[])
       }
    }
    cout << "Minimum det(J) of the original mesh is " << tauval << endl;
+   tauval -= 0.01 * h0.Min(); // Slightly below minJ0 to avoid div by 0.
 
-   // 19. Finally, perform the nonlinear optimization.
-   NewtonSolver *newton = NULL;
-   if (tauval > 0.0)
+   // Perform the nonlinear optimization.
+   TMOPNewtonSolver solver(*ir, solver_type);
+   if (solver_type == 0)
    {
-      tauval = 0.0;
-      TMOPNewtonSolver *tns = new TMOPNewtonSolver(*ir);
-      newton = tns;
-      cout << "TMOPNewtonSolver is used (as all det(J) > 0).\n";
+      // Specify linear solver when we use a Newton-based solver.
+      solver.SetPreconditioner(*S);
    }
-   else
-   {
-      if ( (dim == 2 && metric_id != 22 && metric_id != 252) ||
-           (dim == 3 && metric_id != 352) )
-      {
-         cout << "The mesh is inverted. Use an untangling metric." << endl;
-         return 3;
-      }
-      tauval -= 0.01 * h0.Min(); // Slightly below minJ0 to avoid div by 0.
-      newton = new TMOPDescentNewtonSolver(*ir);
-      cout << "The TMOPDescentNewtonSolver is used (as some det(J) < 0).\n";
-   }
-   newton->SetPreconditioner(*S);
-   newton->SetMaxIter(newton_iter);
-   newton->SetRelTol(newton_rtol);
-   newton->SetAbsTol(0.0);
-   newton->SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
-   newton->SetOperator(a);
-   newton->Mult(b, x.GetTrueVector());
+   solver.SetMaxIter(solver_iter);
+   solver.SetRelTol(solver_rtol);
+   solver.SetAbsTol(0.0);
+   solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
+   solver.SetOperator(a);
+   solver.Mult(b, x.GetTrueVector());
    x.SetFromTrueVector();
-
-   if (newton->GetConverged() == false)
+   if (solver.GetConverged() == false)
    {
-      cout << "NewtonIteration: rtol = " << newton_rtol << " not achieved."
-           << endl;
+      cout << "Nonlinear solver: rtol = " << solver_rtol << " not achieved.\n";
    }
-   delete newton;
 
    // 20. Save the optimized mesh to a file. This output can be viewed later
    //     using GLVis: "glvis -m optimized.mesh".

miniapps/meshing/mesh-optimizer.hpp

@@ -103,7 +103,7 @@ public:
    {
       Vector pos(3);
       T.Transform(ip, pos);
-      if (metric != 14 && metric != 87)
+      if (metric != 14 && metric != 85)
       {
          const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
          const double r = sqrt(xc*xc + yc*yc);
@@ -131,7 +131,7 @@ public:
 
          K *= alpha_bar;
       }
-      else if (metric == 87) // Shape + Alignment
+      else if (metric == 85) // Shape + Alignment
       {
          Vector x = pos;
          double xc = x(0)-0.5, yc = x(1)-0.5;

miniapps/meshing/pmesh-optimizer.cpp

@@ -32,27 +32,29 @@
 // Compile with: make pmesh-optimizer
 //
 // Sample runs:
-//   Adapted analytic Hessian:
+//   Adapted analytic shape:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 4 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
-//   Adapted analytic Hessian with size+orientation:
+//   Adapted analytic size+orientation:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd
-//   Adapted analytic Hessian with Shape+size+orientation
-//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 87 -tid 4 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd
+//   Adapted analytic shape+orientation:
+//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd
+//
 //   Adapted discrete size:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 5 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 5 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8 -cmb 2 -nor
-//
-//   Adapted size+aspect ratio to discrete material indicator
+//   Adapted discrete size+aspect_ratio:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100  -ls 2 -li 100 -bnd -qt 1 -qo 8
-//   Adapted discrete size+orientation (requires GSLIB)
+//   Adapted discrete size+orientation (requires GSLIB):
 //   * mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 8 -ni 100  -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd -ae 1
-//   Adapted discrete aspect-ratio+orientation (requires GSLIB)
-//   * mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 87 -tid 8 -ni 10  -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd -ae 1
-//   Adapted discrete aspect ratio (3D)
+//   Adapted discrete aspect_ratio+orientation (requires GSLIB):
+//   * mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 8 -ni 10  -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd -ae 1
+//   Adapted discrete aspect ratio (3D):
 //     mpirun -np 4 pmesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 302 -tid 7 -ni 20  -ls 2 -li 100 -bnd -qt 1 -qo 8
 //
 //   Adaptive limiting:
-//     mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5 -ae 0
+//     mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5
+//   Adaptive limiting through the L-BFGS solver:
+//     mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 400 -qo 5 -nor -vl 1 -alc 0.5 -st 1
 //   Adaptive limiting through FD (requires GSLIB):
 //   * mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5 -fd -ae 1
 //
@@ -106,8 +108,9 @@ int main (int argc, char *argv[])
    double adapt_lim_const = 0.0;
    int quad_type         = 1;
    int quad_order        = 8;
-   int newton_iter       = 10;
-   double newton_rtol    = 1e-10;
+   int solver_type       = 0;
+   int solver_iter       = 10;
+   double solver_rtol    = 1e-10;
    int lin_solver        = 2;
    int max_lin_iter      = 100;
    bool move_bnd         = true;
@@ -168,9 +171,11 @@ int main (int argc, char *argv[])
                   "3: Closed uniform points");
    args.AddOption(&quad_order, "-qo", "--quad_order",
                   "Order of the quadrature rule.");
-   args.AddOption(&newton_iter, "-ni", "--newton-iters",
+   args.AddOption(&solver_type, "-st", "--solver-type",
+                  " Type of solver: (default) 0: Newton, 1: LBFGS");
+   args.AddOption(&solver_iter, "-ni", "--newton-iters",
                   "Maximum number of Newton iterations.");
-   args.AddOption(&newton_rtol, "-rtol", "--newton-rel-tolerance",
+   args.AddOption(&solver_rtol, "-rtol", "--newton-rel-tolerance",
                   "Relative tolerance for the Newton solver.");
    args.AddOption(&lin_solver, "-ls", "--lin-solver",
                   "Linear solver: 0 - l1-Jacobi, 1 - CG, 2 - MINRES.");
@@ -341,7 +346,7 @@ int main (int argc, char *argv[])
       case 56: metric = new TMOP_Metric_056; break;
       case 58: metric = new TMOP_Metric_058; break;
       case 77: metric = new TMOP_Metric_077; break;
-      case 87: metric = new TMOP_Metric_SS2D; break;
+      case 85: metric = new TMOP_Metric_085; break;
       case 211: metric = new TMOP_Metric_211; break;
       case 252: metric = new TMOP_Metric_252(tauval); break;
       case 301: metric = new TMOP_Metric_301; break;
@@ -548,7 +553,7 @@ int main (int argc, char *argv[])
             tc->SetParDiscreteTargetSize(size);
          }
 
-         if (metric_id == 87)
+         if (metric_id == 85)
          {
             FunctionCoefficient aspr_coeff(discrete_aspr_2d);
             aspr.ProjectCoefficient(aspr_coeff);
@@ -783,47 +788,28 @@ int main (int argc, char *argv[])
    tauval = minJ0;
    if (myid == 0)
    { cout << "Minimum det(J) of the original mesh is " << tauval << endl; }
+   double h0min = h0.Min(), h0min_all;
+   MPI_Allreduce(&h0min, &h0min_all, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+   tauval -= 0.01 * h0min_all; // Slightly below minJ0 to avoid div by 0.
 
-   // 20. Finally, perform the nonlinear optimization.
-   NewtonSolver *newton = NULL;
-   if (tauval > 0.0)
+   // Perform the nonlinear optimization.
+   TMOPNewtonSolver solver(pfespace->GetComm(), *ir, solver_type);
+   if (solver_type == 0)
    {
-      tauval = 0.0;
-      TMOPNewtonSolver *tns = new TMOPNewtonSolver(pfespace->GetComm(), *ir);
-      newton = tns;
-      if (myid == 0)
-      { cout << "TMOPNewtonSolver is used (as all det(J) > 0)." << endl; }
+      // Specify linear solver when we use a Newton-based solver.
+      solver.SetPreconditioner(*S);
    }
-   else
-   {
-      if ( (dim == 2 && metric_id != 22 && metric_id != 252) ||
-           (dim == 3 && metric_id != 352) )
-      {
-         if (myid == 0)
-         { cout << "The mesh is inverted. Use an untangling metric.\n"; }
-         return 3;
-      }
-      double h0min = h0.Min(), h0min_all;
-      MPI_Allreduce(&h0min, &h0min_all, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
-      tauval -= 0.01 * h0min_all; // Slightly below minJ0 to avoid div by 0.
-      newton = new TMOPDescentNewtonSolver(pfespace->GetComm(), *ir);
-      if (myid == 0)
-      { cout << "TMOPDescentNewtonSolver is used (as some det(J) < 0).\n"; }
-   }
-   newton->SetPreconditioner(*S);
-   newton->SetMaxIter(newton_iter);
-   newton->SetRelTol(newton_rtol);
-   newton->SetAbsTol(0.0);
-   newton->SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
-   newton->SetOperator(a);
-   newton->Mult(b, x.GetTrueVector());
+   solver.SetMaxIter(solver_iter);
+   solver.SetRelTol(solver_rtol);
+   solver.SetAbsTol(0.0);
+   solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
+   solver.SetOperator(a);
+   solver.Mult(b, x.GetTrueVector());
    x.SetFromTrueVector();
-   if (myid == 0 && newton->GetConverged() == false)
+   if (myid == 0 && solver.GetConverged() == false)
    {
-      cout << "NewtonIteration: rtol = " << newton_rtol << " not achieved."
-           << endl;
+      cout << "Nonlinear solver: rtol = " << solver_rtol << " not achieved.\n";
    }
-   delete newton;
 
    // 21. Save the optimized mesh to a file. This output can be viewed later
    //     using GLVis: "glvis -m optimized -np num_mpi_tasks".

miniapps/navier/makefile

@@ -12,7 +12,7 @@
 # Use the MFEM build directory
 MFEM_DIR ?= ../..
 MFEM_BUILD_DIR ?= ../..
-# SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/navier/,)
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/navier/,)
 CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
 # Use the MFEM install directory
 # MFEM_INSTALL_DIR = ../../mfem
@@ -40,10 +40,10 @@ endif
 %: %.cpp
 %.o: %.cpp
 
-%: %.cpp $(NAVIER_COMMON_OBJ)
+%: $(SRC)%.cpp $(NAVIER_COMMON_OBJ)
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -o $@ $(NAVIER_COMMON_OBJ) $(MFEM_LIBS)
 
-%.o: %.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+%.o: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
 	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
 
 all: $(MINIAPPS)

tests/unit/CMakeLists.txt

@@ -28,6 +28,7 @@ set(UNIT_TESTS_SRCS
   linalg/test_matrix_rectangular.cpp
   linalg/test_matrix_square.cpp
   linalg/test_ode.cpp
+  linalg/test_fdual.cpp
   linalg/test_ode2.cpp
   linalg/test_operator.cpp
   linalg/test_cg_indefinite.cpp

tests/unit/linalg/test_fdual.cpp

@@ -0,0 +1,195 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "mfem.hpp"
+#include "catch.hpp"
+
+
+using namespace mfem;
+
+
+template<typename tbase>
+tbase exprp02(tbase x,tbase y)
+{
+   return sin(x)*cos(y)+tan(x*y);
+}
+
+template<typename tbase>
+tbase exprp02x(tbase x,tbase y)
+{
+   return cos(x)*cos(y)+y*(1.0+pow(tan(x*y),2.0));
+}
+
+template<typename tbase>
+tbase exprp02y(tbase x,tbase y)
+{
+   return -sin(x)*sin(y)+x*(1.0+pow(tan(x*y),2.0));
+}
+
+
+TEST_CASE("Simple AD tests", "[Simple_AD_tests]")
+{
+   SECTION("sin")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x,1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::sin(xx);
+      d = std::cos(x);
+      REQUIRE(std::abs(d-lrez.dual())<std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("cos")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x,1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::cos(xx);
+      d = -std::sin(x);
+      REQUIRE(std::abs(d-lrez.dual())<std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("tan")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x,1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::tan(xx);
+      d = 1.0+std::tan(x)*std::tan(x);
+      REQUIRE(std::abs(d-lrez.dual())<std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("exp")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x,1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::exp(xx);
+      d = exp(x);
+      REQUIRE(std::abs(d-lrez.dual())<std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("log")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x,1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::log(xx);
+      d = 1.0/x;
+      REQUIRE(std::abs(d-lrez.dual())<std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("pow")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x,1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::pow(xx,1.5);
+      d = 1.5*std::pow(x,0.5);
+      REQUIRE(std::abs(d-lrez.dual())<std::numeric_limits<double>::epsilon());
+
+   }
+
+   SECTION("atan")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x,1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::atan(xx);
+      d = 1.0/(1.0+x*x);
+      REQUIRE(std::abs(d-lrez.dual())<std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("asin")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x,1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::asin(xx);
+      d = 1.0/std::sqrt(1.0-x*x);
+      REQUIRE(std::abs(d-lrez.dual())<std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("acos")
+   {
+      double x = 0.5;
+      double d;
+      ad::FDual<double> xx(x,1.0);
+      ad::FDual<double> lrez;
+      lrez = ad::acos(xx);
+      d = -1.0/std::sqrt(1.0-x*x);
+      REQUIRE(std::abs(d-lrez.dual())<std::numeric_limits<double>::epsilon());
+   }
+
+   SECTION("general")
+   {
+      double x = 1.0;
+      double y = 1.5;
+
+      double pr = exprp02(x,y);
+      double dx = exprp02x(x,y);
+      double dy = exprp02y(x,y);
+
+      {
+         mfem::ad::FDual<double> xx(x,1.0);
+         mfem::ad::FDual<double> yy(y,0.0);
+         mfem::ad::FDual<double> rr=exprp02(xx,yy);
+         REQUIRE(std::abs(rr.real()-pr)<std::numeric_limits<double>::epsilon());
+         REQUIRE(std::abs(rr.dual()-dx)<std::numeric_limits<double>::epsilon());
+      }
+
+      {
+         mfem::ad::FDual<double> xx(x,0.0);
+         mfem::ad::FDual<double> yy(y,1.0);
+         mfem::ad::FDual<double> rr=exprp02(xx,yy);
+         REQUIRE(std::abs(rr.real()-pr)<std::numeric_limits<double>::epsilon());
+         REQUIRE(std::abs(rr.dual()-dy)<std::numeric_limits<double>::epsilon());
+      }
+   }
+
+   SECTION("second_derivative")
+   {
+
+      double x = 0.5;
+      double d;
+      mfem::ad::FDual<mfem::ad::FDual<double>> xxx(mfem::ad::FDual<double>(x,1.0),
+                                                   mfem::ad::FDual<double>(1.0,0.0));
+      mfem::ad::FDual<mfem::ad::FDual<double>> drez=mfem::ad::exp(xxx);
+      d=exp(x);
+      REQUIRE(std::abs(d-drez.dual().dual())<std::numeric_limits<double>::epsilon());
+
+      drez = mfem::ad::log(xxx);
+      d = -1.0/(x*x);
+      REQUIRE(std::abs(d-drez.dual().dual())<std::numeric_limits<double>::epsilon());
+
+      drez = mfem::ad::sin(xxx);
+      d = -sin(x);
+      REQUIRE(std::abs(d-drez.dual().dual())<std::numeric_limits<double>::epsilon());
+
+      drez = mfem::ad::cos(xxx);
+      d = -cos(x);
+      REQUIRE(std::abs(d-drez.dual().dual())<std::numeric_limits<double>::epsilon());
+
+   }
+
+
+
+
+}
+