changes in the configuration files

New descriptions for ex71 and ex71p
Serial example
2020-07-18 13:29:18 -07:00 · 2020-07-08 19:13:23 -07:00 · 2020-07-08 18:32:30 -07:00 · 2020-07-08 16:07:25 -07:00 · 2020-07-08 16:04:37 -07:00 · 2020-06-30 22:37:07 -07:00
22 changed files with 4950 additions and 1 deletions
@@ -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 "")
@@ -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
@@ -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.")
+
@@ -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.")
+
@@ -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.")
+
@@ -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
@@ -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.")
@@ -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()

@@ -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;
+}
@@ -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;
+}
+
+
@@ -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
@@ -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;
+}
+
+
@@ -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
  )

@@ -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
+
@@ -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    
+
@@ -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"
@@ -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
@@ -28,6 +28,7 @@
 #include "solvers.hpp"
 #include "handle.hpp"
 #include "invariants.hpp"
+// #include "fdual.hpp"

 #ifdef MFEM_USE_SUNDIALS
 #include "sundials.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
@@ -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
+
@@ -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
@@ -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());
+
+   }
+
+
+
+
+}
+
Author	SHA1	Message	Date
lazarov	28a6c1f96c	changes in the configuration files	2020-07-18 13:29:18 -07:00
lazarov	cda243493a	New descriptions for ex71 and ex71p	2020-07-08 19:13:23 -07:00
lazarov	85e140bfcf	Serial example	2020-07-08 18:32:30 -07:00
lazarov	53ff1a2bf8	Merge branch 'master' into fad	2020-07-08 16:07:25 -07:00
lazarov	4b47d0eb63	Added support for FADBAD++	2020-07-08 16:04:37 -07:00
lazarov	72aeb54227	The name of ADQIntegratorJ/H class is changed to ADQFunctionJ/H	2020-06-30 22:37:07 -07:00
lazarov	22c33cbdf6	Added: AD integrator for pLaplacian Select between AD integrator and hond coded integrator	2020-06-26 10:13:24 -07:00
lazarov	5287c9f509	Added configuration for CODIPACK	2020-06-18 18:13:09 -07:00
lazarov	fa2db9abf2	Intermediate updates	2020-06-18 18:12:22 -07:00
lazarov	a8a7bc4e40	Native implementation before adding adept	2020-06-18 16:20:59 -07:00
lazarov	1e04cf7798	Merge branch 'master' into fad	2020-06-12 19:30:04 -07:00
bslazarov	fa718bab9a	modified: ../../examples/CMakeLists.txt new file: ../../examples/ex23.cpp modified: ../../fem/CMakeLists.txt new file: ../../fem/adnonlininteg.cpp new file: ../../fem/adnonlininteg.hpp modified: ../../fem/fem.hpp modified: ../../linalg/fdual.hpp new file: ../../linalg/taddensemat.hpp new file: ../../linalg/tadvector.hpp	2020-02-25 20:27:28 -08:00
bslazarov	84209babd2	modified: fdual.hpp modified: ../tests/unit/linalg/test_fdual.cpp	2020-02-16 23:34:17 -08:00
bslazarov	e940331e39	new file: ../../linalg/fdual.hpp modified: ../../linalg/linalg.hpp modified: ../../tests/unit/CMakeLists.txt new file: ../../tests/unit/linalg/test_fdual.cpp	2020-02-14 17:48:25 -08:00