wip

CMakeLists.txt

@@ -351,6 +351,30 @@ 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()
+
+# FADBAD++ package
+if (MFEM_USE_FADBADPP)
+        find_package(FADBADPP REQUIRED)
+        # find_package updates FADBADPP_FOUND, FADBADPP_INCLUDE_DIRS, FADBADPP_LIBRARIES
+endif()
+
+# CoDiPack package
+if (MFEM_USE_CODIPACK)
+	find_package(CODIPACK REQUIRED)
+	# find_package updates CODIPACK_FOUND, CODIPACK_INCLUDE_DIRS, CODIPACK_LIBRARIES
+endif()
+
+# Eigen package
+if (MFEM_USE_EIGEN)
+	find_package(EIGEN REQUIRED)
+	# find_package updates EIGEN_FOUND, EIGEN_INCLUDE_DIRS, EIGEN_LIBRARIES
+endif()
+
 # OCCA
 if (MFEM_USE_OCCA)
    find_package(OCCA REQUIRED)
@@ -417,7 +441,7 @@ endif()
 set(MFEM_TPLS MPI_CXX OPENMP HYPRE BLAS LAPACK SuperLUDist METIS SuiteSparse SUNDIALS PETSC
     SLEPC MESQUITE MUMPS STRUMPACK AXOM CONDUIT Ginkgo GNUTLS GSLIB NETCDF
     MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE ADIOS2
-    CUSPARSE MKL_CPARDISO AMGX CALIPER)
+    CUSPARSE MKL_CPARDISO AMGX CALIPER FADBADPP ADEPT CODIPACK EIGEN)
 # Add all *_FOUND libraries in the variable TPL_LIBRARIES.
 set(TPL_LIBRARIES "")
 set(TPL_INCLUDE_DIRS "")

INSTALL

@@ -462,6 +462,27 @@ MFEM_USE_HIOP = YES/NO
    Enable the usage of HiOp (https://github.com/LLNL/hiop) in MFEM. HiOp is an
    HPC solver for nonlinear optimization problems.
 
+MFEM_USE_ADEPT = YES/NO
+  Enable automatic differentiation using the ADEPT library.
+  (http://www.met.reading.ac.uk/clouds/adept)
+  Please, compile the library with flag --disable-openmp.
+
+MFEM_USE_FADBADPP = YES/NO
+  Enable automatic differentiation using the FADBAD++ library.
+  www.fadbad.com/fadbad.html
+
+MFEM_USE_CODIPACK = YES/NO
+  Enable automatic differentiation using the CoDiPack library.
+  www.scicomp.uni-kl.de/codi/
+
+MFEM_USE_EIGEN = YES/NO
+  Enable the Eigen library.
+  https://gitlab.com/libeigen/eigen
+
+MFEM_USE_ADFORWARD = YES/NO
+  Enable forward mode for AD packages. This option is valid
+  only if the AD package supports two modes (backward/forward).
+
 MFEM_USE_CUDA = YES/NO
    Enables support for CUDA devices in MFEM. CUDA is a parallel computing
    platform and programming model for general computing on graphical processing
@@ -684,6 +705,26 @@ The specific libraries and their options are:
   Options: HIOP_OPT, HIOP_LIB.
   Versions: HIOP >= 0.4.
 
+- ADEPT (optional), used with MFEM_USE_ADEPT = YES
+  URL: www.met.reading.ac.uk/clouds/adept/
+  Options: ADEPT_OPT, ADEPT_LIB
+  Versions: 1.1 and 2.0.5
+
+- FADBAD++ (optiobal), used with MFEM_USE_FADBADPP = YES
+  URL: www.fadbad.com/fadbad.html
+  Options: FADBADPP_OPT
+  Versions: 2.1
+
+- CoDiPack (optiobal), used with MFEM_USE_CODIPACK = YES
+  URL: https://www.scicomp.uni-kl.de/codi/
+  Options: CODIPACK_OPT
+  Versions: 1.9.3
+
+- Eigen (optiobal), used with MFEM_USE_EIGEN = YES
+  URL: https://gitlab.com/libeigen/eigen
+  Options: EIGEN_OPT
+  Versions: 3.3.9
+
 - GSLIB (optional), used when MFEM_USE_GSLIB = YES. The gslib library must be
   built prior to the MFEM build, as follows: download gslib-1.0.7, untar it at
   the same level as MFEM and create a symbolic link: "ln -s gslib-1.0.7 gslib".
@@ -875,6 +916,11 @@ MFEM_USE_MPFR
 MFEM_USE_ZLIB
 MFEM_USE_PUMI
 MFEM_USE_HIOP
+MFEM_USE_ADEPT
+MFEM_USE_FADBADPP
+MFEM_USE_CODIPACK
+MFEM_USE_EIGEN
+MFEM_USE_ADFORWARD
 MFEM_USE_CUDA
 MFEM_USE_OCCA
 MFEM_USE_CEED
@@ -931,6 +977,10 @@ The CMake build system adds auto-detection for the following packages/libraries:
  - POSIXCLOCKS
  - PUMI
  - HIOP
+ - ADEPT
+ - FADBAD++
+ - CoDiPack
+ - Eigen
  - OCCA
  - RAJA
  - UMPIRE

config/cmake/MFEMConfig.cmake.in

@@ -53,6 +53,11 @@ set(MFEM_USE_CEED @MFEM_USE_CEED@)
 set(MFEM_USE_UMPIRE @MFEM_USE_UMPIRE@)
 set(MFEM_USE_SIMD @MFEM_USE_SIMD@)
 set(MFEM_USE_ADIOS2 @MFEM_USE_ADIOS2@)
+set(MFEM_USE_ADEPT @MFEM_USE_ADEPT@)
+set(MFEM_USE_FADBADPP @MFEM_USE_FADBADPP@)
+set(MFEM_USE_CODIPACK @MFEM_USE_CODIPACK@)
+set(MFEM_USE_EIGEN @MFEM_USE_EIGEN@)
+set(MFEM_USE_ADFORWARD @MFEM_USE_ADFORWARD@)
 set(MFEM_USE_CALIPER @MFEM_USE_CALIPER@)
 
 set(MFEM_CXX_COMPILER "@CMAKE_CXX_COMPILER@")

config/cmake/config.hpp.in

@@ -169,6 +169,21 @@
 // library.
 #cmakedefine MFEM_USE_SIMMETRIX
 
+// use ADEPT library for AD
+#cmakedefine MFEM_USE_ADEPT
+
+// use FADBAD++ library for AD
+#cmakedefine MFEM_USE_FADBADPP
+
+// use CoDiPack library for AD
+#cmakedefine MFEM_USE_CODIPACK
+
+// use Eigen library
+#cmakedefine MFEM_USE_EIGEN
+
+// use forward mode for automatic differentiation
+#cmakedefine MFEM_USE_ADFORWARD
+
 // Enable interface to the MKL CPardiso library.
 #cmakedefine MFEM_USE_MKL_CPARDISO
 

config/cmake/modules/FindADEPT.cmake

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

config/cmake/modules/FindCODIPACK.cmake

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

config/cmake/modules/FindEIGEN.cmake

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

config/cmake/modules/FindFADBADPP.cmake

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

config/config.hpp.in

@@ -174,6 +174,21 @@
 // library.
 // #define MFEM_USE_SIMMETRIX
 
+// use ADEPT library for AD
+// #define MFEM_USE_ADEPT
+
+// use FADBAD++ library for AD
+// #define MFEM_USE_FADBADPP
+
+// use CoDiPack library for AD
+// #define MFEM_USE_CODIPACK
+
+// use Eigen library 
+// #define MFEM_USE_EIGEN
+
+// use forward mode for automatic differentiation
+// #define MFEM_USE_ADFORWARD
+
 // Enable interface to the MKL CPardiso library.
 // #define MFEM_USE_MKL_CPARDISO
 

config/config.mk.in

@@ -46,6 +46,11 @@ MFEM_USE_SIDRE         = @MFEM_USE_SIDRE@
 MFEM_USE_CONDUIT       = @MFEM_USE_CONDUIT@
 MFEM_USE_PUMI          = @MFEM_USE_PUMI@
 MFEM_USE_HIOP          = @MFEM_USE_HIOP@
+MFEM_USE_ADEPT         = @MFEM_USE_ADEPT@
+MFEM_USE_FADBADPP      = @MFEM_USE_FADBADPP@
+MFEM_USE_CODIPACK      = @MFEM_USE_CODIPACK@
+MFEM_USE_EIGEN         = @MFEM_USE_EIGEN@
+MFEM_USE_ADFORWARD     = @MFEM_USE_ADFORWARD@
 MFEM_USE_GSLIB         = @MFEM_USE_GSLIB@
 MFEM_USE_CUDA          = @MFEM_USE_CUDA@
 MFEM_USE_HIP           = @MFEM_USE_HIP@

config/defaults.cmake

@@ -55,6 +55,11 @@ option(MFEM_USE_CEED "Enable CEED" OFF)
 option(MFEM_USE_UMPIRE "Enable Umpire" OFF)
 option(MFEM_USE_SIMD "Enable use of SIMD intrinsics" OFF)
 option(MFEM_USE_ADIOS2 "Enable ADIOS2" OFF)
+option(MFEM_USE_ADEPT "Enable AD using ADEPT" OFF)
+option(MFEM_USE_FADBADPP "Enable AD using FADBAD++" OFF)
+option(MFEM_USE_CODIPACK "Enable AD using CoDiPack" OFF)
+option(MFEM_USE_EIGEN "Enable Eigen" OFF)
+option(MFEM_USE_ADFORWARD "Enable forward mode for AD" OFF)
 option(MFEM_USE_CALIPER "Enable Caliper support" OFF)
 option(MFEM_USE_MKL_CPARDISO "Enable MKL CPardiso" OFF)
 
@@ -218,6 +223,18 @@ 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(FADBADPP_INCLUDE_DIRS "${MFEM_DIR}/../FADBAD++" CACHE STRING "Path to FADBAD++ headers.")
+set(FADBADPP_LIBRARIES "")
+
+set(CODIPACK_INCLUDE_DIRS "${MFEM_DIR}/../CoDiPack/inlude" CACHE STRING "Path to CoDoPack headers.")
+set(CODIPACK_LIBRARIES "")
+
+set(EIGEN_INCLUDE_DIRS "${MFEM_DIR}/../Eigen/" CACHE STRING "Path to Eigen headers.")
+set(EIGEN_LIBRARIES "")
+
 # Some useful variables:
 set(CMAKE_SKIP_PREPROCESSED_SOURCE_RULES ON) # Skip *.i rules
 set(CMAKE_SKIP_ASSEMBLY_SOURCE_RULES  ON)    # Skip *.s rules

config/defaults.mk

@@ -146,6 +146,11 @@ MFEM_USE_CALIPER       = NO
 MFEM_USE_UMPIRE        = NO
 MFEM_USE_SIMD          = NO
 MFEM_USE_ADIOS2        = NO
+MFEM_USE_ADEPT         = NO
+MFEM_USE_FADBADPP      = NO
+MFEM_USE_CODIAPCK      = NO
+MFEM_USE_EIGEN         = NO
+MFEM_USE_ADFORWARD     = NO
 MFEM_USE_MKL_CPARDISO  = NO
 
 # MPI library compile and link flags
@@ -379,6 +384,26 @@ HIOP_DIR = @MFEM_DIR@/../hiop/install
 HIOP_OPT = -I$(HIOP_DIR)/include
 HIOP_LIB = -L$(HIOP_DIR)/lib -lhiop $(LAPACK_LIB)
 
+# ADEPT
+ADEPT_DIR = @MFEM_DIR@/../adept-1.1
+ADEPT_OPT = -I$(ADEPT_DIR)/include
+ADEPT_LIB = -L$(ADEPT_DIR)/lib -ladept
+
+# FADBAD++
+FADBADPP_DIR = @MFEM_DIR@/../FADBAD++
+FADBADPP_OPT = -I$(FADBADPP_DIR)
+FADBADPP_LIB = -L.
+
+# CoDiPack
+CODIPACK_DIR = @MFEM_DIR@/../CoDiPack
+CODIPACK_OPT = -I$(CODIPACK_DIR)
+CODIPACK_LIB = -L.
+
+# Eigen
+EIGEN_DIR = @MFEM_DIR@/../Eigen
+EIGEN_OPT = -I$(EIGEN_DIR)
+EIGEN_LIB = -L.
+
 # GSLIB library
 GSLIB_DIR = @MFEM_DIR@/../gslib/build
 GSLIB_OPT = -I$(GSLIB_DIR)/include

fem/nonlinearform.cpp

@@ -979,6 +979,8 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
 
             for (int k = 0; k < bfnfi.Size(); ++k)
             {
+               if (bfnfi_marker[k] &&
+                   (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }	    
                bfnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
                for (int l=0; l<fes.Size(); ++l)
                {

makefile

@@ -267,7 +267,7 @@ endif
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
 MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
  SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP GSLIB\
- OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX CALIPER
+ OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX CALIPER ADEPT FADBADPP CODIPACK EIGEN
 
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
 SLEPC_ERROR_MSG = $(if $(SLEPC_FOUND),,. SLEPC config not found: $(SLEPC_VARS))
@@ -333,7 +333,8 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP MFEM_USE_OCCA\
  MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE MFEM_USE_SIMD MFEM_USE_ADIOS2\
  MFEM_USE_MKL_CPARDISO MFEM_USE_AMGX MFEM_USE_MUMPS MFEM_USE_CALIPER\
- MFEM_SOURCE_DIR MFEM_INSTALL_DIR
+ MFEM_USE_ADEPT MFEM_USE_FADBADPP MFEM_USE_ADFORWARD MFEM_USE_CODIPACK\
+ MFEM_SOURCE_DIR MFEM_INSTALL_DIR MFEM_USE_EIGEN
 
 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -658,6 +659,11 @@ status info:
 	$(info MFEM_USE_UMPIRE        = $(MFEM_USE_UMPIRE))
 	$(info MFEM_USE_SIMD          = $(MFEM_USE_SIMD))
 	$(info MFEM_USE_ADIOS2        = $(MFEM_USE_ADIOS2))
+	$(info MFEM_USE_ADEPT         = $(MFEM_USE_ADEPT))
+	$(info MFEM_USE_FADBADPP      = $(MFEM_USE_FADBADPP))
+	$(info MFEM_USE_CODIPACK      = $(MFEM_USE_CODIPACK))
+	$(info MFEM_USE_EIGEN         = $(MFEM_USE_EIGEN))
+	$(info MFEM_USE_ADFORWARD     = $(MFEM_USE_ADFORWARD))
 	$(info MFEM_USE_MKL_CPARDISO  = $(MFEM_USE_MKL_CPARDISO))
 	$(info MFEM_CXX               = $(value MFEM_CXX))
 	$(info MFEM_HOST_CXX          = $(value MFEM_HOST_CXX))

miniapps/CMakeLists.txt

@@ -29,3 +29,4 @@ add_subdirectory(gslib)
 add_subdirectory(solvers)
 add_subdirectory(shifted)
 add_subdirectory(mtop)
+add_subdirectory(autodiff)

miniapps/autodiff/CMakeLists.txt

@@ -0,0 +1,75 @@
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+list(APPEND SEQADIFF_COMMON_SOURCES
+        stokes.cpp
+	advdiff.cpp)
+
+list(APPEND SEQADIFF_COMMON_HEADERS
+        fdual.hpp
+	tadvector.hpp
+	taddensemat.hpp
+	adnonlininteg.hpp
+	admfem.hpp
+	stokes.hpp
+	advdiff.hpp)
+
+
+
+convert_filenames_to_full_paths(SEQADIFF_COMMON_SOURCES)
+convert_filenames_to_full_paths(SEQADIFF_COMMON_HEADERS)
+
+set(SEQADIFF_COMMON_FILES
+	EXTRA_SOURCES ${SEQADIFF_COMMON_SOURCES}
+	EXTRA_HEADERS ${SEQADIFF_COMMON_HEADERS})
+
+add_mfem_miniapp(seqadiff
+	MAIN seq_example.cpp
+	${SEQADIFF_COMMON_FILES}
+	LIBRARIES mfem)
+
+add_mfem_miniapp(seqtest
+	MAIN seq_test.cpp
+	${SEQADIFF_COMMON_FILES}
+        LIBRARIES mfem)
+
+if(MFEM_USE_MPI)
+
+list(APPEND PARADIFF_COMMON_SOURCES
+	)
+list(APPEND PARADIFF_COMMON_HEADERS
+	)
+
+convert_filenames_to_full_paths(PARADIFF_COMMON_SOURCES)
+convert_filenames_to_full_paths(PARADIFF_COMMON_HEADERS)
+
+set(PARADIFF_COMMON_FILES
+	EXTRA_SOURCES ${PARADIFF_COMMON_SOURCES} ${SEQADIFF_COMMON_SOURCES}
+	EXTRA_HEADERS ${PARADIFF_COMMON_HEADERS} ${SEQADIFF_COMMON_HEADERS})
+
+message(STATUS "PARADIFF_COMMON_FILES: ${PARADIFF_COMMON_FILES}")
+message(STATUS "SEQADIFF_COMMON_FILES: ${SEQADIFF_COMMON_FILES}")
+
+add_mfem_miniapp(paradiff
+     MAIN par_example.cpp
+     ${PARADIFF_COMMON_FILES}
+     LIBRARIES mfem)
+
+add_mfem_miniapp(stokes
+     MAIN stokes_test.cpp
+     LIBRARIES mfem)
+
+add_mfem_miniapp(advdiff
+     MAIN advdiff_test.cpp
+     LIBRARIES mfem)
+
+
+endif ()

miniapps/autodiff/admfem.hpp

@@ -0,0 +1,576 @@
+#ifndef ADMFEM_HPP 
+#define ADMFEM_HPP
+
+#include "mfem.hpp"
+#include "fdual.hpp"
+#include "tadvector.hpp"
+#include "taddensemat.hpp"
+
+#ifdef MFEM_USE_CODIPACK
+#include <codi.hpp>
+
+namespace mfem {
+
+namespace ad {
+#ifdef MFEM_USE_ADFORWARD
+typedef codi::RealForward ADFloatType;
+typedef TADVector<ADFloatType> ADVectorType;
+typedef TADDenseMatrix<ADFloatType> ADMatrixType;
+#else
+typedef codi::RealReverse ADFloatType;
+typedef TADVector<ADFloatType> ADVectorType;
+typedef TADDenseMatrix<ADFloatType> ADMatrixType;
+#endif
+}
+
+
+template<int vector_size=1, int state_size=1, int param_size=0>
+class VectorFuncAutoDiff
+{
+public:
+    VectorFuncAutoDiff(std::function<void(mfem::Vector&, ad::ADVectorType&, ad::ADVectorType&)> F_)
+    {
+            F=F_;
+    }
+
+    void QJacobian(mfem::Vector &vparam, mfem::Vector &uu, mfem::DenseMatrix &jac)
+    {
+#ifdef MFEM_USE_ADFORWARD
+        // use forward mode
+        jac.SetSize(vector_size, state_size);
+        jac = 0.0;
+        {
+            ad::ADVectorType aduu(state_size);
+            ad::ADVectorType rr(vector_size);
+            for(int i=0;i<state_size;i++){
+                aduu[i].setValue(uu[i]);
+                aduu[i].setGradient(0.0);
+            }
+            for(int ii=0;ii<state_size;ii++){
+                aduu[ii].setGradient(1.0);
+                F(vparam,aduu,rr);
+                for(int jj=0;jj<vector_size;jj++)
+                {
+                    jac(jj,ii)=rr[jj].getGradient();
+                }
+                aduu[ii].setGradient(0.0);
+            }
+        }
+#else // use reverse mode
+        jac.SetSize(vector_size, state_size);
+        jac = 0.0;
+        {
+            ad::ADVectorType aduu(state_size);
+            ad::ADVectorType rr(vector_size);
+            for(int i=0;i<state_size;i++){
+                aduu[i]=uu[i];
+            }
+
+            ad::ADFloatType::TapeType& tape =ad::ADFloatType::getGlobalTape();
+            typename ad::ADFloatType::TapeType::Position pos=tape.getPosition();
+
+            tape.setActive();
+            for(int ii=0;ii<state_size;ii++){ tape.registerInput(aduu[ii]); }
+            F(vparam,aduu,rr);
+            for(int ii=0;ii<vector_size;ii++){ tape.registerOutput(rr[ii]); }
+            tape.setPassive();
+
+            for(int jj=0;jj<vector_size;jj++){
+                rr[jj].setGradient(1.0);
+                tape.evaluate();
+                for(int ii=0;ii<state_size;ii++){
+                    jac(jj,ii)=aduu[ii].getGradient();
+                }
+                tape.clearAdjoints();
+                rr[jj].setGradient(0.0);
+            }
+            tape.reset(pos);
+        }
+#endif
+    }
+private:
+    std::function<void(mfem::Vector&, ad::ADVectorType&, ad::ADVectorType&)> F;
+}; //VectorFuncAutoDiff
+
+
+template<template<typename, typename, typename, int, int, int> class CTD
+         , int vector_size=1, int state_size=1, int param_size=0>
+class QVectorFuncAutoDiff
+{
+public:
+    /// Evaluates the vector function for given set of parameters and state
+    /// values in vector uu. The result is returned in vector rr.
+    void QVectorFunc(const mfem::Vector &vparam, mfem::Vector &uu, mfem::Vector& rr)
+    {
+        CTD<double,const mfem::Vector, mfem::Vector,
+                        vector_size, state_size, param_size> tf;
+        tf(vparam,uu,rr);
+    }
+
+    /// Returns the gradient of CTD(...) in the dense matrix jac.
+    /// The dimensions of jac are vector_size x state_size, where state_size is the
+    /// length of vector uu.
+    void QJacobian(mfem::Vector &vparam, mfem::Vector &uu, mfem::DenseMatrix &jac)
+    {
+#ifdef MFEM_USE_ADFORWARD
+        // use forward mode
+
+        CTD<ad::ADFloatType, const Vector, ad::ADVectorType,
+                        vector_size, state_size, param_size> tf;
+        jac.SetSize(vector_size, state_size);
+        jac = 0.0;
+        {
+            ad::ADVectorType aduu(state_size);
+            ad::ADVectorType rr(vector_size);
+            for(int i=0;i<state_size;i++){
+                aduu[i].setValue(uu[i]);
+                aduu[i].setGradient(0.0);
+            }
+
+            for(int ii=0;ii<state_size;ii++){
+                aduu[ii].setGradient(1.0);
+                tf(vparam,aduu,rr);
+                for(int jj=0;jj<vector_size;jj++)
+                {
+                    jac(jj,ii)=rr[jj].getGradient();
+                }
+                aduu[ii].setGradient(0.0);
+            }
+        }
+#else //end MFEM_USE_ADFORWARD
+        // use reverse mode
+        CTD<ad::ADFloatType, const Vector, ad::ADVectorType,
+                vector_size, state_size, param_size> tf;
+
+        jac.SetSize(vector_size, state_size);
+        jac = 0.0;
+        {
+            ad::ADVectorType aduu(state_size);
+            ad::ADVectorType rr(vector_size);
+            for(int i=0;i<state_size;i++){
+                aduu[i]=uu[i];
+            }
+
+            ad::ADFloatType::TapeType& tape =ad::ADFloatType::getGlobalTape();
+            typename ad::ADFloatType::TapeType::Position pos=tape.getPosition();
+
+            tape.setActive();
+            for(int ii=0;ii<state_size;ii++){ tape.registerInput(aduu[ii]); }
+            tf(vparam,aduu,rr);
+            for(int ii=0;ii<vector_size;ii++){ tape.registerOutput(rr[ii]); }
+            tape.setPassive();
+            for(int jj=0;jj<vector_size;jj++){
+                rr[jj].setGradient(1.0);
+                tape.evaluate();
+                for(int ii=0;ii<state_size;ii++){
+                    jac(jj,ii)=aduu[ii].getGradient();
+                }
+                tape.clearAdjoints();
+                rr[jj].setGradient(0.0);
+            }
+            tape.reset(pos);
+        }
+#endif
+    }
+};
+
+template<template<typename, typename, typename, int, int> class CTD
+         , int state_size=1, int param_size=0>
+class QFunctionAutoDiff
+{
+public:
+
+    /// Evaluates a function for arguments vparam and uu.
+    /// The evaluatin is based on the operator() in the
+    /// user provided functor CTD.
+    double QEval(const Vector &vparam, Vector &uu)
+    {
+        CTD<double, const Vector, Vector, state_size, param_size> tf;
+        return tf(vparam,uu);
+    }
+
+    /// Provides the same functionality as QGrad.
+    void QVectorFunc(const Vector &vparam, Vector &uu, Vector &rr)
+    {
+        QGrad(vparam,uu,rr);
+    }
+
+    /// Returns the first derivative of CTD(...) with
+    /// respect to the active arguments proved in vector uu.
+    /// The length of rr is the same as for uu.
+    void QGrad(const Vector &vparam, Vector &uu, Vector &rr)
+    {
+
+#ifdef MFEM_USE_ADFORWARD
+        // use forward mode
+        CTD<ad::ADFloatType, const Vector, ad::ADVectorType,
+                        state_size, param_size> tf;
+        rr.SetSize(state_size);
+        {
+            ad::ADVectorType aduu(state_size);
+            for(int i=0;i<state_size;i++){
+                aduu[i].setValue(uu[i]);
+                aduu[i].setGradient(0.0);
+            }
+
+            ad::ADFloatType rez;
+
+            for(int ii=0;ii<state_size;ii++){
+                aduu[ii].setGradient(1.0);
+                rez=tf(vparam,aduu);
+                rr[ii]=rez.getGradient();
+                aduu[ii].setGradient(0.0);
+            }
+        }
+
+#else
+        typedef codi::RealReverse ADFType;
+        typedef TADVector<ADFType> ADFVector;
+
+        CTD<ADFType, const Vector, ADFVector, state_size, param_size> tf;
+        {
+            ADFVector aduu(state_size);
+            ADFType rez;
+            for(int i=0;i<state_size;i++){
+                aduu[i]=uu[i];
+            }
+
+            ADFType::TapeType& tape =ADFType::getGlobalTape();
+            typename ADFType::TapeType::Position pos=tape.getPosition();
+
+            tape.setActive();
+            for(int ii=0;ii<state_size;ii++){ tape.registerInput(aduu[ii]); }
+
+            rez=tf(vparam,aduu);
+            tape.registerOutput(rez);
+            tape.setPassive();
+
+            rez.setGradient(1.0);
+            tape.evaluate();
+            for(int i=0;i<state_size;i++){
+                rr[i]=aduu[i].getGradient();
+            }
+            tape.reset(pos);
+        }
+#endif
+    }
+
+    /// Provides same functionality as QHessian.
+    void QJacobian(mfem::Vector &vparam, mfem::Vector &uu, mfem::DenseMatrix &jac)
+    {
+        QHessian(vparam,uu,jac);
+    }
+
+    /// Returns the Hessian of CTD(...) in the dense matrix hh.
+    /// The dimensions of jac are state_size x state_size, where state_size is the
+    /// length of vector uu.
+    void QHessian(mfem::Vector &vparam, mfem::Vector &uu, mfem::DenseMatrix &jac)
+    {
+#ifdef MFEM_USE_ADFORWARD
+        // use forward-forward mode
+        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;
+
+        CTD<ADSType, const Vector, ADSVector, state_size, param_size> tf;
+
+        jac.SetSize(state_size);
+        jac=0.0;
+        {
+            ADSVector aduu(state_size);
+            for(int ii = 0;  ii < state_size; 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 < state_size; ii++)
+            {
+                aduu[ii].value().gradient()=1.0;
+                for(int jj=0; jj<(ii+1); jj++)
+                {
+                    aduu[jj].gradient().value()=1.0;
+                    ADSType rez=tf(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
+        //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;
+
+        CTD<ADSType, const Vector, ADSVector, state_size, param_size> tf;
+
+        jac.SetSize(state_size);
+        jac=0.0;
+        {
+            ADSVector aduu(state_size);
+            for(int ii=0;ii < state_size ; ii++)
+            {
+                aduu[ii].value().value()=uu[ii];
+            }
+
+            ADSType rez;
+
+            ADSType::TapeType& tape = ADSType::getGlobalTape();
+            typename ADSType::TapeType::Position pos;
+            for(int ii = 0; ii < state_size ; ii++)
+            {
+                pos=tape.getPosition();
+                tape.setActive();
+
+                for(int jj=0;jj < state_size; jj++) {
+                    if(jj==ii) {aduu[jj].value().gradient()=1.0;}
+                    else {aduu[jj].value().gradient()=0.0;}
+                    tape.registerInput(aduu[jj]);
+                }
+
+                rez=tf(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
+    }
+};
+
+}
+#else //end MFEM_USE_CODIPACK
+
+//USE NATIVE IMPLEMENTATION
+namespace mfem {
+
+namespace ad {
+typedef FDual<double> ADFloatType;
+typedef TADVector<ADFloatType> ADVectorType;
+typedef TADDenseMatrix<ADFloatType> ADMatrixType;
+}
+
+template<int vector_size=1, int state_size=1, int param_size=0>
+class VectorFuncAutoDiff
+{
+public:
+    VectorFuncAutoDiff(std::function<void(mfem::Vector&, ad::ADVectorType&, ad::ADVectorType&)> F_)
+    {
+            F=F_;
+    }
+
+public:
+    void QJacobian(mfem::Vector &vparam, mfem::Vector &uu, mfem::DenseMatrix &jac)
+    {
+        jac.SetSize(vector_size, state_size);
+        jac = 0.0;
+        {
+           ad::ADVectorType aduu(uu); //all dual numbers are initialized to zero
+           ad::ADVectorType rr(vector_size);
+
+           for (int ii = 0; ii < state_size; ii++)
+           {
+              aduu[ii].dual(1.0);
+              F(vparam,aduu,rr);
+              for (int jj = 0; jj < vector_size; jj++)
+              {
+                 jac(jj, ii) = rr[jj].dual();
+              }
+              aduu[ii].dual(0.0);
+           }
+        }
+    }
+
+private:
+    std::function<void(mfem::Vector&, ad::ADVectorType&, ad::ADVectorType&)> F;
+
+};
+
+
+template<template<typename, typename, typename, int, int, int> class CTD
+         , int vector_size=1, int state_size=1, int param_size=0>
+class QVectorFuncAutoDiff
+{
+private:
+    /// MFEM native forward AD-type
+    typedef ad::FDual<double> ADFType;
+    typedef TADVector<ADFType> ADFVector;
+    typedef TADDenseMatrix<ADFType> ADFDenseMatrix;
+
+public:
+    /// Returns a vector valued function rr for supplied passive arguments
+    /// vparam and active arguments uu. The evaluation is based on the
+    /// user supplied CTD template class.
+    void QVectorFunc(const Vector &vparam, Vector &uu, Vector &rr)
+    {
+       CTD<double, const Vector, Vector,
+               vector_size, state_size, param_size> func;
+       func(vparam, uu, rr);
+    }
+
+    /// Returns the gradient of CTD(...) residual in the dense matrix jac.
+    /// The dimensions of jac are vector_size x state_size, where state_size is the
+    /// length of vector uu.
+    void QJacobian(mfem::Vector &vparam, mfem::Vector &uu, mfem::DenseMatrix &jac)
+    {
+        //use native AD package
+        jac.SetSize(vector_size, state_size);
+        jac = 0.0;
+        {
+           ADFVector aduu(uu); //all dual numbers are initialized to zero
+           ADFVector rr(vector_size);
+
+           for (int ii = 0; ii < state_size; ii++)
+           {
+              aduu[ii].dual(1.0);
+              QEval(vparam, aduu, rr);
+              for (int jj = 0; jj < vector_size; jj++)
+              {
+                 jac(jj, ii) = rr[jj].dual();
+              }
+              aduu[ii].dual(0.0);
+           }
+        }
+    }
+
+private:
+    /// Evaluates the residual from CTD(...).
+    /// Intended for internal use only.
+    void QEval(const Vector &vparam, ADFVector &uu, ADFVector &rr)
+    {
+       CTD<ADFType, const Vector, ADFVector,
+               vector_size, state_size, param_size> tf;
+       tf(vparam, uu, rr);
+    }
+};
+
+template<template<typename, typename, typename, int, int> class CTD
+         , int state_size=1, int param_size=0>
+class QFunctionAutoDiff
+{
+private:
+    ///MFEM native AD-type for first derivatives
+    typedef ad::FDual<double> ADFType;
+    typedef TADVector<ADFType> ADFVector;
+    typedef TADDenseMatrix<ADFType> ADFDenseMatrix;
+    ///MFEM native AD-type for second derivatives
+    typedef ad::FDual<ADFType> ADSType;
+    typedef TADVector<ADSType> ADSVector;
+    typedef TADDenseMatrix<ADSType> ADSDenseMatrix;
+
+public:
+    /// Evaluates a function for arguments vparam and uu.
+    /// The evaluatin is based on the operator() in the
+    /// user provided functor CTD.
+    double QEval(const Vector &vparam, Vector &uu)
+    {
+        CTD<double, const Vector, Vector, state_size, param_size> tf;
+        return tf(vparam,uu);
+    }
+
+    /// Provides the same functionality as QGrad.
+    void QVectorFunc(const Vector &vparam, Vector &uu, Vector &rr)
+    {
+        QGrad(vparam,uu,rr);
+    }
+
+    /// Returns the first derivative of CTD(...) with
+    /// respect to the active arguments proved in vector uu.
+    /// The length of rr is the same as for uu.
+    void QGrad(const Vector &vparam, Vector &uu, Vector &rr)
+    {
+        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 = QEval(vparam, aduu);
+           rr[ii] = rez.dual();
+           aduu[ii].dual(0.0);
+        }
+    }
+
+    /// Provides same functionality as QHessian.
+    void QJacobian(mfem::Vector &vparam, mfem::Vector &uu, mfem::DenseMatrix &jac)
+    {
+        QHessian(vparam,uu,jac);
+    }
+
+    /// Returns the Hessian of CTD(...) in the dense matrix hh.
+    /// The dimensions of jac are state_size x state_size, where state_size is the
+    /// length of vector uu.
+    void QHessian(mfem::Vector &vparam, mfem::Vector &uu, mfem::DenseMatrix &jac)
+    {
+        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 = QEval(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));
+           }
+        }
+    }
+
+private:
+    /// Evaluates the first derivative of CTD(...).
+    /// Intended for internal use only.
+    ADFType QEval(const Vector &vparam, ADFVector &uu)
+    {
+       CTD<ADFType, const Vector, ADFVector, state_size, param_size> tf;
+       return tf(vparam, uu);
+    }
+
+    /// Evaluates the second derivative of CTD(...).
+    /// Intended for internal use only.
+    ADSType QEval(const Vector &vparam, ADSVector &uu)
+    {
+       CTD<ADSType, const Vector, ADSVector, state_size, param_size> tf;
+       return tf(vparam, uu);
+    }
+
+};
+}//end namespace mfem
+#endif // NATIVE
+#endif // ADMFEM_HPP
+

miniapps/autodiff/adnonlininteg.hpp

@@ -0,0 +1,503 @@
+// 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_ADNONLININTEG
+#define MFEM_ADNONLININTEG
+
+#include "mfem.hpp"
+#include "tadvector.hpp"
+#include "taddensemat.hpp"
+#include "fdual.hpp"
+
+#if defined MFEM_USE_ADEPT
+#include <adept.h>
+#elif defined MFEM_USE_FADBADPP
+#include <fadiff.h>
+#include <badiff.h>
+#endif
+
+
+namespace mfem
+{
+
+///Automatic differentiation class - for user coded functor CTD returning  PDE's
+/// energy and residual at a point, provides the derivative of the residual
+/// with respect to the active arguments of the functor.
+///
+/**
+   The two templated parameters are
+   template<typename, typename> class CTD and int m.
+   The templated class CTD is a user-provided functor which implements
+   the function of interest. The CTD class should have the following signature:
+
+template<typename DType, typename MVType>
+
+class MyQFunctorJ
+
+{
+
+    DType operator()(const Vector &vparam, MVType &uu)
+
+        {...}
+
+    void operator()(const Vector &vparam, MVType &uu, MVType &rr)
+
+    {...}
+
+};
+
+
+DType - represents the scalars in the functor, and  MVType - represents
+the vectors in the functor. In ADQFunctionTJ,  DType will be replaced
+either with double or AD-type.
+The first operator()(const Vector &vparam, MVType &uu) takes as input
+standard MFEM Vector vparam (passive arguments), which holds all parameters
+ supplied to the functor and MVType vector uu, which holds all active
+arguments provided to the functor. The operator returns a scalar
+representing the value of the function(energy) for passive parameters vparam
+and active arguments uu.
+
+
+The second void operator()(const Vector &vparam, MVType &uu, MVType &rr)
+returns a vector-valued function rr with passive input arguments vparam
+and active (AD) arguments supplied in vector uu. Since the size of the
+return vector rr is not known in advance, it should be provided explicitly
+to the ADQFunctionTJ class, i.e., the template parameter m in ADQFunctionTJ
+represents the dimensions of the return vector rr.
+
+
+For PDE discretized problem, the first operator should return the energy
+evaluated a point, and the second operator should return the residual
+evaluated at the point. The derivative of the residual vector rr with
+respect to the active arguments uu is evaluated automatically by
+ADQFunctionTJ and returned by the ADQFunctionTJ:: QFunctionDD(..) method.
+*/
+template<template<typename, typename> class CTD, int m>
+class ADQFunctionTJ
+{
+   // m - dimension of the residual vector
+   // the Jacobian will have dimensions [m,length(uu)]
+protected:
+#ifdef MFEM_USE_ADEPT
+   adept::Stack m_stack;
+#endif
+
+public:
+#if defined MFEM_USE_ADEPT
+   /// AD-type based on the Adept AD library.
+   typedef adept::adouble ADFType;
+   typedef TADVector<ADFType> ADFVector;
+   typedef TADDenseMatrix<ADFType> ADFDenseMatrix;
+#elif defined MFEM_USE_FADBADPP
+#ifdef MFEM_USE_ADFORWARD
+   /// AD-type based on the FADBAD++ library -
+   /// used only for forward differentiation.
+   typedef fadbad::F<double> ADFType;
+   typedef TADVector<ADFType> ADFVector;
+   typedef TADDenseMatrix<ADFType> ADFDenseMatrix;
+#else
+   /// AD-type based on the FADBAD++ library -
+   /// used only for reverse AD mode.
+   typedef fadbad::B<double> ADFType;
+   typedef TADVector<ADFType> ADFVector;
+   typedef TADDenseMatrix<ADFType> ADFDenseMatrix;
+#endif
+#else
+   /// MFEM native forward AD-type
+   typedef ad::FDual<double> ADFType;
+   typedef TADVector<ADFType> ADFVector;
+   typedef TADDenseMatrix<ADFType> ADFDenseMatrix;
+#endif
+
+#ifdef MFEM_USE_ADEPT
+   ADQFunctionTJ() : m_stack(false) {}
+#else
+   ADQFunctionTJ() {}
+#endif
+
+   ~ADQFunctionTJ() {}
+
+   ///Returns the energy for passive arguments vparam and
+   /// active arguments uu. The evaluation is based on the
+   /// first operator in the user-supplied CTD template class.
+   double QFunction(const Vector &vparam, Vector &uu)
+   {
+      CTD<double, Vector> func;
+      return func(vparam, uu);
+   }
+
+   ///Returns vector valued function rr for passive input parametrs vparam
+   /// and active arguments provided in uu.
+   void QFunctionDU(const Vector &vparam, ADFVector &uu, ADFVector &rr)
+   {
+      CTD<ADFType, ADFVector> func;
+      func(vparam, uu, rr);
+   }
+
+   ///Evaluates automatically the first derivative of
+   /// QFunction(const Vector &vparam, Vector &uu) with respect to
+   /// the active vector arguments uu. The dimension or rr is the same
+   /// as the dimension of uu. The method can be used for testing
+   /// the correctness of hand-coded gradients of  QFunction(...).
+   void QFunctionAU(const Vector &vparam, Vector &uu, Vector &rr)
+   {
+      //the result is computed automatically by differentiating
+      //QFunction with respect to uu
+      CTD<ADFType, ADFVector> func;
+      int n = uu.Size();
+      rr.SetSize(n);
+
+#if defined MFEM_USE_ADEPT
+      //use ADEPT package
+      adept::Stack *p_stack = adept::active_stack();
+      p_stack->deactivate();
+
+      m_stack.activate();
+      {
+         ADFVector aduu(uu);
+         ADFType rez;
+         m_stack.new_recording();
+         rez = func(vparam, aduu);
+         m_stack.independent(aduu.GetData(), n); //independent variables
+         m_stack.dependent(&rez, 1);             //dependent variables
+         m_stack.jacobian(rr.GetData());
+      }
+      m_stack.deactivate();
+#elif defined MFEM_USE_FADBADPP
+      //use FADBAD++
+#ifdef MFEM_USE_ADFORWARD
+      {
+         ADFVector aduu(uu);
+         ADFType rez;
+         for (int ii = 0; ii < n; ii++)
+         {
+            aduu[ii].diff(ii, n);
+         }
+         rez = func(vparam, aduu);
+         for (int ii = 0; ii < n; ii++)
+         {
+            rr[ii] = rez.d(ii);
+         }
+      }
+#else
+      {
+         ADFVector aduu(uu);
+         ADFType rez;
+         rez = func(vparam, aduu);
+         rez.diff(0, 1);
+         for (int ii = 0; ii < n; ii++)
+         {
+            rr[ii] = aduu[ii].d(0);
+         }
+      }
+#endif
+#else
+      //use native AD package
+      {
+         ADFVector aduu(uu); //all dual numbers are initialized to zero
+         ADFType rez;
+
+         for (int ii = 0; ii < n; ii++)
+         {
+            aduu[ii].dual(1.0);
+            rez = func(vparam, aduu);
+            rr[ii] = rez.dual();
+            aduu[ii].dual(0.0);
+         }
+      }
+#endif
+   }
+
+   ///Returns vector valued function rr for supplied passive arguments
+   /// vparam and active arguments uu. The evaluation is based on the
+   /// user supplied CTD template class (second operator).
+   void QFunctionDU(const Vector &vparam, Vector &uu, Vector &rr)
+   {
+      CTD<double, Vector> func;
+      func(vparam, uu, rr);
+   }
+
+   ///Evaluates automatically the derivative or the vector function
+   /// QFunctionDU(...), i.e., the retuned vector rr, with respect to
+   /// the active arguments uu. The dimensions of the the dense matix
+   /// jac are [m,n] where m is the size of the vector rr and n is the
+   /// size of the active vector uu. The parameter m should be supplied as
+   /// template parameter int m in ADQFunctionTJ.
+   void QFunctionDD(const Vector &vparam, 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();
+         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_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);
+         }
+         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);
+         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);
+            QFunctionDU(vparam, aduu, rr);
+            for (int jj = 0; jj < m; jj++)
+            {
+               jac(jj, ii) = rr[jj].dual();
+            }
+            aduu[ii].dual(0.0);
+         }
+      }
+#endif
+   }
+};
+
+//template class for differentiation; the function
+//for differentiation is supplied as a functor
+//the operator()(scalar,vector) defines the actual function
+
+///ADQFunctionTH is a templated class evaluating the first and
+/// the second derivatives of a user supplied function encoded
+/// in a functor class CTD. The signature of CTD is as follows:
+/**
+
+template<typename DType, typename MVType> class CTD
+
+{
+
+    DType operator()(const Vector &vparam, MVType &uu)
+    {...}
+
+};
+
+The operator returns a scalar representing the value of the
+function for passive arguments vparam and active arguments uu.
+The first derivative of the operator is evaluated automatically
+by ADQFunctionTH::QFunctionDU(const Vector &vparam, Vector &uu, Vector &rr)
+method. The length of the return vector rr is the same as for
+the input vector uu.  The second derivate (the Hessian) of the function
+is evaluated again automatically by
+QFunctionDD(const Vector &vparam, const Vector &uu, DenseMatrix &jac).
+
+Hessian evaluation is an expensive process. The provided AD-functionality
+is intended to be utilized for development purposes. The performance will
+increase significantly by replacing the AD-provided derivatives with hand-coded.
+*/
+template<template<typename, typename> class CTD>
+class ADQFunctionTH
+{
+public:
+#if defined MFEM_USE_FADBADPP
+   ///AD-type derived from FADBAD++
+   typedef fadbad::B<double> ADFType;
+   ///AD vector type for evaluation of first derivatives
+   typedef TADVector<ADFType> ADFVector;
+   ///AD dense matrix type for evaluation of first derivatives
+   typedef TADDenseMatrix<ADFType> ADFDenseMatrix;
+   ///AD-type for the second derivatives derived from
+   ///FADBAD++
+   typedef fadbad::B<fadbad::F<double>> ADSType;
+   /// AD vector type for evaluation of second derivatives
+   typedef TADVector<ADSType> ADSVector;
+   /// AD dense matrix type for  evaluation of second derivatives
+   typedef TADDenseMatrix<ADSType> ADSDenseMatrix;
+#else
+   ///MFEM native AD-type for first derivatives
+   typedef ad::FDual<double> ADFType;
+   typedef TADVector<ADFType> ADFVector;
+   typedef TADDenseMatrix<ADFType> ADFDenseMatrix;
+   ///MFEM native AD-type for second derivatives
+   typedef ad::FDual<ADFType> ADSType;
+   typedef TADVector<ADSType> ADSVector;
+   typedef TADDenseMatrix<ADSType> ADSDenseMatrix;
+#endif
+
+   ADQFunctionTH() {}
+
+   ~ADQFunctionTH() {}
+
+   ///Evaluates a function for arguments vparam and uu.
+   /// The evaluatin is based on the operator() in the
+   /// user provided functor CTD.
+   double QFunction(const Vector &vparam, Vector &uu)
+   {
+      CTD<double, Vector> tf;
+      return tf(vparam, uu);
+   }
+
+   ///Evaluates the first derivative of QFunction(...).
+   /// Intended for internal use only.
+   ADFType QFunction(const Vector &vparam, ADFVector &uu)
+   {
+      CTD<ADFType, ADFVector> tf;
+      return tf(vparam, uu);
+   }
+
+   ///Evaluates the second derivative of QFunction(...).
+   /// Intended for internal use only.
+   ADSType QFunction(const Vector &vparam, ADSVector &uu)
+   {
+      CTD<ADSType, ADSVector> tf;
+      return tf(vparam, uu);
+   }
+
+   ///Returns the first derivative of QFunction(...) with
+   /// respect to the active arguments proved in vector uu.
+   /// The length of rr is the same as for uu.
+   void QFunctionDU(const Vector &vparam, Vector &uu, Vector &rr)
+   {
+#if defined MFEM_USE_FADBADPP
+      int n = uu.Size();
+      rr.SetSize(n);
+      ADFVector aduu(uu);
+      ADFType rez;
+      rez = 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 = QFunction(vparam, aduu);
+         rr[ii] = rez.dual();
+         aduu[ii].dual(0.0);
+      }
+#endif
+   }
+
+   ///Returns the Hessian of QFunction(...) in the dense matrix jac.
+   /// The dimensions of jac are m x m, where m is the length of vector uu.
+   void QFunctionDD(const Vector &vparam, const Vector &uu, DenseMatrix &jac)
+   {
+#if 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 = 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 = 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
+   }
+
+};
+
+} // namespace mfem
+
+#endif

miniapps/autodiff/advdiff.cpp

@@ -0,0 +1,108 @@
+#include "advdiff.hpp"
+#include "petsc.h"
+
+namespace mfem {
+
+
+void AdvectionDiffusionMXSolver::DirectSolver(mfem::BlockOperator& A)
+{
+    delete psol; psol=nullptr;
+    delete prec; prec=nullptr;
+    delete pmat; pmat=nullptr;
+
+    mfem::HypreParMatrix* A00=static_cast<mfem::HypreParMatrix*>(&(A.GetBlock(0,0)));
+    mfem::HypreParMatrix* A01=static_cast<mfem::HypreParMatrix*>(&(A.GetBlock(0,1)));
+    mfem::HypreParMatrix* A02=static_cast<mfem::HypreParMatrix*>(&(A.GetBlock(0,2)));
+    mfem::HypreParMatrix* A10=static_cast<mfem::HypreParMatrix*>(&(A.GetBlock(1,0)));
+    mfem::HypreParMatrix* A11=static_cast<mfem::HypreParMatrix*>(&(A.GetBlock(1,1)));
+    mfem::HypreParMatrix* A12=static_cast<mfem::HypreParMatrix*>(&(A.GetBlock(1,2)));
+    mfem::HypreParMatrix* A20=static_cast<mfem::HypreParMatrix*>(&(A.GetBlock(2,0)));
+    mfem::HypreParMatrix* A21=static_cast<mfem::HypreParMatrix*>(&(A.GetBlock(2,1)));
+    mfem::HypreParMatrix* A22=static_cast<mfem::HypreParMatrix*>(&(A.GetBlock(2,2)));
+
+    Array2D< HypreParMatrix * > bm(3,3);
+    bm(0,0)=A00; bm(0,1)=A01; bm(0,2)=A02;
+    bm(1,0)=A10; bm(1,1)=A11; bm(1,2)=A12;
+    bm(2,0)=A20; bm(2,1)=A21; bm(2,2)=A22;
+
+    HypreParMatrix* MM=mfem::HypreParMatrixFromBlocks(bm);
+
+    mfem::PetscParMatrix* pmat=new mfem::PetscParMatrix(MM,mfem::Operator::PETSC_MATAIJ);
+    mfem::PetscLinearSolver*   psol=new mfem::PetscLinearSolver(pmesh->GetComm());
+    psol->SetOperator(*pmat);
+    psol->SetAbsTol(abs_tol);
+    psol->SetRelTol(rel_tol);
+    psol->SetMaxIter(max_iter);
+    psol->SetPrintLevel(print_level);
+    psol->Mult(rhs, sol);
+
+    delete MM;
+}
+
+void AdvectionDiffusionMXSolver::PETSCSolver(mfem::BlockOperator& A)
+{
+    delete psol; psol=nullptr;
+    delete prec; prec=nullptr;
+    delete pmat; pmat=nullptr;
+
+    pmat= new mfem::PetscParMatrix(pmesh->GetComm(),&A, mfem::Operator::PETSC_MATAIJ);
+    // construct the preconditioner
+    prec = new mfem::PetscFieldSplitSolver(pmesh->GetComm(),*pmat,"prec_");
+    // construct the linear solver
+    psol = new mfem::PetscLinearSolver(pmesh->GetComm());
+
+    psol->SetOperator(*pmat);
+    psol->SetPreconditioner(*prec);
+    psol->SetAbsTol(abs_tol);
+    psol->SetRelTol(rel_tol);
+    psol->SetMaxIter(max_iter);
+    psol->SetPrintLevel(print_level);
+    psol->Mult(rhs, sol);
+
+}
+
+void AdvectionDiffusionMXSolver::FSolve()
+{
+    sol=0.0;
+
+    if(nfin.size()==0)
+    {
+        //add the domain integrator
+        nfin.push_back(new mfem::AdvectionDiffusionMX(dicoef,vecoef,mucoef,incoef));
+        nf->AddDomainIntegrator(nfin[nfin.size()-1]);
+
+
+        //add BC face integrators
+        for(auto it=bcc.begin();it!=bcc.end();it++)
+        {
+            mfem::AdvectionDiffusionMX* iin=new mfem::AdvectionDiffusionMX(dicoef,vecoef,mucoef,incoef);
+            nfin.push_back(iin);
+            iin->SetDirichletBCCoeficient(it->first);
+            nf->AddBdrFaceIntegrator(iin,it->second);
+        }
+
+
+        nf->SetGradientType(mfem::Operator::Type::Hypre_ParCSR);
+    }
+
+    // set the RHS
+    nf->Mult(sol,rhs);
+    rhs.Neg();
+
+    mfem::BlockOperator& A=nf->GetGradient(sol);
+
+    DirectSolver(A);
+
+
+}
+
+
+
+void AdvectionDiffusionMXSolver::ASolve(BlockVector &rhs)
+{
+
+}
+
+
+
+}

miniapps/autodiff/advdiff.hpp

@@ -0,0 +1,981 @@
+#ifndef ADVDIFF_H
+#define ADVDIFF_H
+
+#include<map>
+#include<vector>
+#include "mfem.hpp"
+
+namespace mfem {
+
+
+
+class AdvectionDiffusionMX:public mfem::BlockNonlinearFormIntegrator
+{
+public:
+    AdvectionDiffusionMX()
+    {
+        diffc=nullptr;
+        veloc=nullptr;
+        inpuc=nullptr;
+        mucoe=nullptr;
+
+        dbc=nullptr;
+        gamma=10;
+
+
+    }
+
+    AdvectionDiffusionMX(mfem::Coefficient* diffusion_, mfem::VectorCoefficient* veloc_,
+                         mfem::Coefficient* mu_, mfem::Coefficient* inp_)
+    {
+        diffc=diffusion_;
+        veloc=veloc_;
+        inpuc=inp_;
+        mucoe=mu_;
+
+        dbc=nullptr;
+        gamma=10;
+    }
+
+
+    void SetDirichletBCCoeficient(mfem::Coefficient* bc)
+    {
+        dbc=bc;
+    }
+
+    void SetVelocity(mfem::VectorCoefficient* veloc_)
+    {
+        veloc=veloc_;
+    }
+
+    void SetDiffusion(mfem::Coefficient* diffusion_)
+    {
+        diffc=diffusion_;
+    }
+
+    void SetReactionCoefficient(mfem::Coefficient* mu_)
+    {
+        mucoe=mu_;
+    }
+
+    void SetVolInput(mfem::Coefficient* imp_)
+    {
+        inpuc=imp_;
+    }
+
+    void SetDirichletBCPenalization(double penal)
+    {
+        gamma=penal;
+    }
+
+    virtual ~AdvectionDiffusionMX(){}
+
+    virtual
+    double GetElementEnergy(const Array<const FiniteElement *> &el,
+                            ElementTransformation &Tr,
+                            const Array<const Vector *> &elfun) override
+    {
+        return 0.0;
+    }
+
+    virtual
+    void AssembleElementVector(const Array<const FiniteElement *> &el,
+                               ElementTransformation &Tr,
+                               const Array<const Vector *> &elfun,
+                               const Array<Vector *> &elvec) override
+    {
+        int dof_u = el[0]->GetDof();
+        int dof_p = el[1]->GetDof();
+        int dof_z = el[2]->GetDof();
+
+        int dim = el[0]->GetDim();
+
+        elvec[0]->SetSize(dof_u);
+        elvec[1]->SetSize(dof_p);
+        elvec[2]->SetSize(dof_z);
+
+        *(elvec[0])=0.0;
+        *(elvec[1])=0.0;
+        *(elvec[2])=0.0;
+
+
+        int spaceDim = Tr.GetDimension();
+        if (dim != spaceDim)
+        {
+           mfem::mfem_error("AdvectionDiffusionMX::AssembleElementVector"
+                            " is not defined on manifold meshes");
+        }
+
+        mfem::DenseMatrix bsu;
+        mfem::DenseMatrix bsp;
+        mfem::DenseMatrix bsz;
+
+        //set B-matrices
+        bsu.SetSize(dof_u,4); // [u, ux,uy,uz]
+        bsp.SetSize(dof_p,4); // [px,py,pz, div(p)]
+        bsz.SetSize(dof_z,1); // [z]
+
+        bsu=0.0;
+        bsp=0.0;
+        bsz=0.0;
+
+        Vector sh;
+        DenseMatrix dh;
+
+        const IntegrationRule *ir = nullptr;
+        int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0]);
+        ir=&IntRules.Get(Tr.GetGeometryType(),order);
+
+        double mmu;
+        DenseMatrix kap(3,3); kap=0.0;
+        Vector vel(3); vel=0.0;
+        double rhc;
+        double w;
+
+        Vector ss(9);
+        Vector rr(9);
+
+        for (int i = 0; i < ir->GetNPoints(); i++)
+        {
+           const IntegrationPoint &ip = ir->IntPoint(i);
+           Tr.SetIntPoint(&ip);
+           w=Tr.Weight();
+           w = ip.weight * w;
+
+           sh.SetDataAndSize(bsu.GetData(),dof_u);
+           el[0]->CalcPhysShape(Tr,sh);
+
+           dh.UseExternalData(bsu.GetData()+dof_u,dof_u,dim);
+           el[0]->CalcPhysDShape(Tr,dh);
+
+           sh.SetDataAndSize(ss.GetData(),4);
+           bsu.MultTranspose(*(elfun[0]),sh);
+
+
+           dh.UseExternalData(bsp.GetData(), dof_p, dim);
+           el[1]->CalcPhysVShape(Tr,dh);
+
+           sh.SetDataAndSize(bsp.GetData()+3*dof_p, dof_p);
+           el[1]->CalcPhysDivShape(Tr,sh);
+
+           sh.SetDataAndSize(ss.GetData()+4,4);
+           bsp.MultTranspose(*(elfun[1]),sh);
+
+           sh.SetDataAndSize(bsz.GetData(),dof_z);
+           el[2]->CalcPhysShape(Tr,sh);
+
+           sh.SetDataAndSize(ss.GetData()+8,1);
+           bsz.MultTranspose(*(elfun[2]),sh);
+
+
+           mmu=0.0;
+           if(mucoe!=nullptr)
+           {
+               mmu=mucoe->Eval(Tr,ip);
+           }
+
+           if(diffc!=nullptr)
+           {
+               kap(0,0)=diffc->Eval(Tr,ip);
+               kap(1,1)=kap(0,0);
+               kap(2,2)=kap(0,0);
+           }
+
+           vel=0.0;
+           if(veloc!=nullptr)
+           {
+               veloc->Eval(vel,Tr,ip);
+           }
+
+           rhc=0.0;
+           if(inpuc!=nullptr)
+           {
+               rhc=inpuc->Eval(Tr,ip);
+           }
+
+           EvalQRes(kap.GetData(), vel.GetData(), mmu, rhc, ss.GetData(), rr.GetData());
+
+           sh.SetDataAndSize(rr.GetData(),4);
+           bsu.AddMult_a(w,sh, *(elvec[0]));
+           sh.SetDataAndSize(rr.GetData()+4,4);
+           bsp.AddMult_a(w,sh, *(elvec[1]));
+           sh.SetDataAndSize(rr.GetData()+8,1);
+           bsz.AddMult_a(w,sh, *(elvec[2]));
+
+
+        }//end integration loop
+
+    }
+
+    virtual
+    void AssembleElementGrad(const Array<const FiniteElement *> &el,
+                             ElementTransformation &Tr,
+                             const Array<const Vector *> &elfun,
+                             const Array2D<DenseMatrix *> &elmats) override
+    {
+        int dof_u = el[0]->GetDof();
+        int dof_p = el[1]->GetDof();
+        int dof_z = el[2]->GetDof();
+
+        int dim = el[0]->GetDim();
+        int spaceDim = Tr.GetDimension();
+        if (dim != spaceDim)
+        {
+           mfem::mfem_error("AdvectionDiffusionMX::AssembleElementVector"
+                            " is not defined on manifold meshes");
+        }
+
+        elmats(0,0)->SetSize(dof_u,dof_u);
+        elmats(0,1)->SetSize(dof_u,dof_p);
+        elmats(0,2)->SetSize(dof_u,dof_z);
+
+        elmats(1,0)->SetSize(dof_p,dof_u);
+        elmats(1,1)->SetSize(dof_p,dof_p);
+        elmats(1,2)->SetSize(dof_p,dof_z);
+
+        elmats(2,0)->SetSize(dof_z,dof_u);
+        elmats(2,1)->SetSize(dof_z,dof_p);
+        elmats(2,2)->SetSize(dof_z,dof_z);
+
+        (*elmats(0,0))=0.0;
+        (*elmats(0,1))=0.0;
+        (*elmats(0,2))=0.0;
+
+        (*elmats(1,0))=0.0;
+        (*elmats(1,1))=0.0;
+        (*elmats(1,2))=0.0;
+
+        (*elmats(2,0))=0.0;
+        (*elmats(2,1))=0.0;
+        (*elmats(2,2))=0.0;
+
+        mfem::DenseMatrix bsu;
+        mfem::DenseMatrix bsp;
+        mfem::DenseMatrix bsz;
+
+        //set B-matrices
+        bsu.SetSize(dof_u,4); // [u, ux,uy,uz]
+        bsp.SetSize(dof_p,4); // [px,py,pz, div(p)]
+        bsz.SetSize(dof_z,1); // [z]
+
+        bsu=0.0;
+        bsp=0.0;
+        bsz=0.0;
+
+        Vector sh;
+        DenseMatrix dh;
+
+        DenseMatrix th;
+        DenseMatrix mh;
+        DenseMatrix rh;
+
+        const IntegrationRule *ir = nullptr;
+        int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0]);
+        ir=&IntRules.Get(Tr.GetGeometryType(),order);
+
+        double mmu;
+        DenseMatrix kap(3,3); kap=0.0;
+        Vector vel(3); vel=0.0;
+        double w;
+
+        DenseMatrix mm; //state matrix
+        mm.SetSize(9,9); //set the size of the state matrix
+
+        for (int i = 0; i < ir->GetNPoints(); i++)
+        {
+           const IntegrationPoint &ip = ir->IntPoint(i);
+           Tr.SetIntPoint(&ip);
+           w=Tr.Weight();
+           w = ip.weight * w;
+
+           sh.SetDataAndSize(bsu.GetData(),dof_u);
+           el[0]->CalcPhysShape(Tr,sh);
+
+           dh.UseExternalData(bsu.GetData()+dof_u,dof_u,dim);
+           el[0]->CalcPhysDShape(Tr,dh);
+
+           dh.UseExternalData(bsp.GetData(), dof_p, dim);
+           el[1]->CalcPhysVShape(Tr,dh);
+
+           sh.SetDataAndSize(bsp.GetData()+3*dof_p, dof_p);
+           el[1]->CalcPhysDivShape(Tr,sh);
+
+           sh.SetDataAndSize(bsz.GetData(),dof_z);
+           el[2]->CalcPhysShape(Tr,sh);
+
+           mmu=0.0;
+           if(mucoe!=nullptr)
+           {
+               mmu=mucoe->Eval(Tr,ip);
+           }
+
+           if(diffc!=nullptr)
+           {
+               kap(0,0)=diffc->Eval(Tr,ip);
+               kap(1,1)=kap(0,0);
+               kap(2,2)=kap(0,0);
+           }
+
+           vel=0.0;
+           if(veloc!=nullptr)
+           {
+               veloc->Eval(vel,Tr,ip);
+           }
+
+           EvalQMat(kap.GetData(),vel.GetData(),mmu,mm.GetData());
+
+           mh.SetSize(4,4);
+           mh.CopyMN(mm,4,4,0,0);
+           mh.Transpose();
+           th.SetSize(dof_u,4);
+           rh.SetSize(dof_u,dof_u);
+           MultABt(bsu,mh,th);
+           MultABt(th,bsu,rh);
+           elmats(0,0)->AddMatrix(w,rh,0,0);
+
+           mh.SetSize(4,4);
+           mh.CopyMN(mm,4,4,0,4);
+           mh.Transpose();
+           th.SetSize(dof_u,4);
+           rh.SetSize(dof_u,dof_p);
+           MultABt(bsu,mh,th);
+           MultABt(th,bsp,rh);
+           elmats(0,1)->AddMatrix(w,rh,0,0);
+
+           mh.SetSize(4,1);
+           mh.CopyMN(mm,4,1,0,8);
+           mh.Transpose();
+           th.SetSize(dof_u,1);
+           rh.SetSize(dof_u,dof_z);
+           MultABt(bsu,mh,th);
+           MultABt(th,bsz,rh);
+           elmats(0,2)->AddMatrix(w,rh,0,0);
+
+
+           mh.SetSize(4,4);
+           mh.CopyMN(mm,4,4,4,0);
+           mh.Transpose();
+           th.SetSize(dof_p,4);
+           rh.SetSize(dof_p,dof_u);
+           MultABt(bsp,mh,th);
+           MultABt(th,bsu,rh);
+           elmats(1,0)->AddMatrix(w,rh,0,0);
+
+           mh.SetSize(4,4);
+           mh.CopyMN(mm,4,4,4,4);
+           mh.Transpose();
+           th.SetSize(dof_p,4);
+           rh.SetSize(dof_p,dof_p);
+           MultABt(bsp,mh,th);
+           MultABt(th,bsp,rh);
+           elmats(1,1)->AddMatrix(w,rh,0,0);
+
+           mh.SetSize(4,1);
+           mh.CopyMN(mm,4,1,4,8);
+           mh.Transpose();
+           th.SetSize(dof_p,1);
+           rh.SetSize(dof_p,dof_z);
+           MultABt(bsp,mh,th);
+           MultABt(th,bsz,rh);
+           elmats(1,2)->AddMatrix(w,rh,0,0);
+
+           mh.SetSize(1,4);
+           mh.CopyMN(mm,1,4,8,0);
+           mh.Transpose();
+           th.SetSize(dof_z,4);
+           rh.SetSize(dof_z,dof_u);
+           MultABt(bsz,mh,th);
+           MultABt(th,bsu,rh);
+           elmats(2,0)->AddMatrix(w,rh,0,0);
+
+           mh.SetSize(1,4);
+           mh.CopyMN(mm,1,4,8,4);
+           mh.Transpose();
+           th.SetSize(dof_z,4);
+           rh.SetSize(dof_z,dof_p);
+           MultABt(bsz,mh,th);
+           MultABt(th,bsp,rh);
+           elmats(2,1)->AddMatrix(w,rh,0,0);
+
+           mh.SetSize(1,1);
+           mh.CopyMN(mm,1,1,8,8);
+           mh.Transpose();
+           th.SetSize(dof_z,1);
+           rh.SetSize(dof_z,dof_z);
+           MultABt(bsz,mh,th);
+           MultABt(th,bsz,rh);
+           elmats(2,2)->AddMatrix(w,rh,0,0);
+        }
+    }
+
+    virtual
+    void AssembleFaceVector(const Array<const FiniteElement *> &el1,
+                            const Array<const FiniteElement *> &el2,
+                            FaceElementTransformations &Tr,
+                            const Array<const Vector *> &elfun,
+                            const Array<Vector *> &elvec)
+    {
+        int dom_id=Tr.Attribute;
+
+        int dof_u = el1[0]->GetDof();
+        int dof_p = el1[1]->GetDof();
+        int dof_z = el1[2]->GetDof();
+
+        int dim = el1[0]->GetDim();
+
+        elvec[0]->SetSize(dof_u);
+        elvec[1]->SetSize(dof_p);
+        elvec[2]->SetSize(dof_z);
+
+        *(elvec[0])=0.0;
+        *(elvec[1])=0.0;
+        *(elvec[2])=0.0;
+
+
+
+        mfem::Vector bsu;
+        bsu.SetSize(dof_u);
+
+        mfem::Vector nor; //normal vector
+        mfem::Vector nir; //unit normal vector
+
+        nor.SetSize(dim);
+        nir.SetSize(dim);
+
+        const IntegrationRule *ir = nullptr;
+        int order= 2 * el1[0]->GetOrder();
+        ir=&IntRules.Get(Tr.GetGeometryType(),order);
+
+        double ih=0.0;  //inverse of the element characteristic length
+        double nr=0.0;  //norm of the normal vector
+        double ek=0.0;  //the smallest eigenvalue of the diffusion tensor
+        double gg=0.0;  //boundary value
+        double bp=0.0;
+        double w;
+
+        mfem::Vector ev(3);
+        Vector vel(3); vel=0.0;
+        DenseMatrix kap(3,3); kap=0.0;
+
+
+        for (int i = 0; i < ir->GetNPoints(); i++)
+        {
+           const IntegrationPoint &ipg = ir->IntPoint(i);
+           Tr.SetAllIntPoints(&ipg);
+           const mfem::IntegrationPoint &ip=Tr.GetElement1IntPoint();
+           mfem::CalcOrtho(Tr.Jacobian(),nor);
+
+           w = Tr.Weight();
+           w = ipg.weight * w;
+
+           nr=nor.Norml2();
+           ih=nr/Tr.Elem1->Weight();
+           nir.Set(1.0/nr,nor);
+
+
+           if(diffc!=nullptr)
+           {
+               kap(0,0)=diffc->Eval(Tr,ip);
+               kap(1,1)=kap(0,0);
+               kap(2,2)=kap(0,0);
+               ek=kap(0,0);
+
+               //general case
+               //kap.Eigenvalues(ev);
+               //ek=std::min(ev(0),ev(1));
+               //ek=std::min(ek,ev(2));
+           }
+
+           if(veloc!=nullptr)
+           {
+               veloc->Eval(vel,Tr,ip);
+           }
+
+           if(dbc!=nullptr)
+           {
+               gg=dbc->Eval(Tr,ip);
+           }
+
+           el1[0]->CalcShape(ip,bsu);
+
+           bp=0.0;
+           for(int ii=0;ii<dim;ii++)
+           {
+               bp=bp+nir(ii)*vel(ii);
+           }
+           bp=std::min(0.0,bp);
+
+           w=w*gg*(bp*bp/ih+gamma*ek*ek*ih);
+           elvec[0]->Add(-w,bsu);
+        }
+
+    }
+
+    virtual
+    void AssembleFaceGrad(const Array<const FiniteElement *> &el1,
+                          const Array<const FiniteElement *> &el2,
+                          FaceElementTransformations &Tr,
+                          const Array<const Vector *> &elfun,
+                          const Array2D<DenseMatrix *> &elmats)
+    {
+        int dom_id=Tr.Attribute;
+
+        int dof_u = el1[0]->GetDof();
+        int dof_p = el1[1]->GetDof();
+        int dof_z = el1[2]->GetDof();
+
+        int dim = el1[0]->GetDim();
+
+        elmats(0,0)->SetSize(dof_u,dof_u);
+        elmats(0,1)->SetSize(dof_u,dof_p);
+        elmats(0,2)->SetSize(dof_u,dof_z);
+
+        elmats(1,0)->SetSize(dof_p,dof_u);
+        elmats(1,1)->SetSize(dof_p,dof_p);
+        elmats(1,2)->SetSize(dof_p,dof_z);
+
+        elmats(2,0)->SetSize(dof_z,dof_u);
+        elmats(2,1)->SetSize(dof_z,dof_p);
+        elmats(2,2)->SetSize(dof_z,dof_z);
+
+        for(int i=0;i<3;i++){
+            for(int j=0;j<3;j++){
+                (*elmats(i,j))=0.0;
+            }
+        }
+
+        mfem::Vector bsu;
+        bsu.SetSize(dof_u);
+
+        mfem::Vector nor; //normal vector
+        mfem::Vector nir; //unit normal vector
+
+        nor.SetSize(dim);
+        nir.SetSize(dim);
+
+        const IntegrationRule *ir = nullptr;
+        int order= 2 * el1[0]->GetOrder();
+        ir=&IntRules.Get(Tr.GetGeometryType(),order);
+
+        double ih=0.0;  //inverse of the element characteristic length
+        double nr=0.0;  //norm of the normal vector
+        double ek=0.0;  //the smallest eigenvalue of the diffusion tensor
+        double bp;
+        double w;
+
+        mfem::Vector ev(3);
+        Vector vel(3); vel=0.0;
+        DenseMatrix kap(3,3); kap=0.0;
+
+        for (int i = 0; i < ir->GetNPoints(); i++)
+        {
+           const IntegrationPoint &ipg = ir->IntPoint(i);
+           Tr.SetAllIntPoints(&ipg);
+           const mfem::IntegrationPoint &ip=Tr.GetElement1IntPoint();
+           mfem::CalcOrtho(Tr.Jacobian(),nor);
+
+           w = Tr.Weight();
+           w = ipg.weight * w;
+
+           nr=nor.Norml2();
+           ih=nr/Tr.Elem1->Weight();
+           nir.Set(1.0/nr,nor);
+
+           if(diffc!=nullptr)
+           {
+               kap(0,0)=diffc->Eval(Tr,ip);
+               kap(1,1)=kap(0,0);
+               kap(2,2)=kap(0,0);
+
+               ek=kap(0,0);
+
+               //general case
+               //kap.Eigenvalues(ev);
+               //ek=std::min(ev(0),ev(1));
+               //ek=std::min(ek,ev(2));
+           }
+
+           if(veloc!=nullptr)
+           {
+               veloc->Eval(vel,Tr,ip);
+           }
+
+           el1[0]->CalcShape(ip,bsu);
+
+           bp=0.0;
+           for(int ii=0;ii<nir.Size();ii++)
+           {
+               bp=bp+nir(ii)*vel(ii);
+           }
+           bp=std::min(0.0,bp);
+
+           w=w*(bp*bp/ih+gamma*ek*ek*ih);
+
+           mfem::AddMult_a_VVt(w,bsu,*elmats(0,0));
+
+        }
+    }
+
+
+private:
+    mfem::Coefficient* diffc;
+    mfem::VectorCoefficient* veloc;
+    mfem::Coefficient* inpuc;
+    mfem::Coefficient* mucoe;
+
+    //boundary faces
+    double gamma; //penalization for the Nitsche method
+    mfem::Coefficient* dbc; //Dirichlet BC
+
+
+    //aa[3,3] - diffusion matrix
+    //bb[3] - velocity
+    //mmu - reaction coeficient
+    void EvalQRes(double* kap, double* bb, double mmu, double inp, double* uu, double* rr)
+    {
+        double t5,t11,t17;
+        t5 = bb[0]*uu[0]-kap[0]*uu[1]-kap[3]*uu[2]-kap[6]*uu[3]-uu[4];
+        t11 = bb[1]*uu[0]-kap[1]*uu[1]-kap[4]*uu[2]-kap[7]*uu[3]-uu[5];
+        t17 = bb[2]*uu[0]-kap[2]*uu[1]-kap[5]*uu[2]-kap[8]*uu[3]-uu[6];
+        rr[0] = mmu*uu[8]+t11*bb[1]+t17*bb[2]+t5*bb[0];
+        rr[1] = -t11*kap[1]-t17*kap[2]-t5*kap[0];
+        rr[2] = -t11*kap[4]-t17*kap[5]-t5*kap[3];
+        rr[3] = -t11*kap[7]-t17*kap[8]-t5*kap[6];
+        rr[4] = -t5;
+        rr[5] = -t11;
+        rr[6] = -t17;
+        rr[7] = uu[8];
+        rr[8] = mmu*uu[0]-inp+uu[7];
+    }
+
+    void EvalQMat(double* kap, double* bb, double mmu, double* kmat)
+    {
+        double t1,t2,t3,t8,t12,t16,t17,t18,t19,t24,t28,t29,t30;
+        double t31,t36,t37,t38,t39;
+        t1 = bb[0]*bb[0];
+        t2 = bb[1]*bb[1];
+        t3 = bb[2]*bb[2];
+        t8 = -bb[0]*kap[0]-bb[1]*kap[1]-bb[2]*kap[2];
+        t12 = -bb[0]*kap[3]-bb[1]*kap[4]-bb[2]*kap[5];
+        t16 = -bb[0]*kap[6]-bb[1]*kap[7]-bb[2]*kap[8];
+        t17 = kap[0]*kap[0];
+        t18 = kap[1]*kap[1];
+        t19 = kap[2]*kap[2];
+        t24 = kap[0]*kap[3]+kap[1]*kap[4]+kap[2]*kap[5];
+        t28 = kap[0]*kap[6]+kap[1]*kap[7]+kap[2]*kap[8];
+        t29 = kap[3]*kap[3];
+        t30 = kap[4]*kap[4];
+        t31 = kap[5]*kap[5];
+        t36 = kap[3]*kap[6]+kap[4]*kap[7]+kap[5]*kap[8];
+        t37 = kap[6]*kap[6];
+        t38 = kap[7]*kap[7];
+        t39 = kap[8]*kap[8];
+        kmat[0] = t1+t2+t3;
+        kmat[1] = t8;
+        kmat[2] = t12;
+        kmat[3] = t16;
+        kmat[4] = -bb[0];
+        kmat[5] = -bb[1];
+        kmat[6] = -bb[2];
+        kmat[7] = 0.0;
+        kmat[8] = mmu;
+        kmat[9] = t8;
+        kmat[10] = t17+t18+t19;
+        kmat[11] = t24;
+        kmat[12] = t28;
+        kmat[13] = kap[0];
+        kmat[14] = kap[1];
+        kmat[15] = kap[2];
+        kmat[16] = 0.0;
+        kmat[17] = 0.0;
+        kmat[18] = t12;
+        kmat[19] = t24;
+        kmat[20] = t29+t30+t31;
+        kmat[21] = t36;
+        kmat[22] = kap[3];
+        kmat[23] = kap[4];
+        kmat[24] = kap[5];
+        kmat[25] = 0.0;
+        kmat[26] = 0.0;
+        kmat[27] = t16;
+        kmat[28] = t28;
+        kmat[29] = t36;
+        kmat[30] = t37+t38+t39;
+        kmat[31] = kap[6];
+        kmat[32] = kap[7];
+        kmat[33] = kap[8];
+        kmat[34] = 0.0;
+        kmat[35] = 0.0;
+        kmat[36] = -bb[0];
+        kmat[37] = kap[0];
+        kmat[38] = kap[3];
+        kmat[39] = kap[6];
+        kmat[40] = 1.0;
+        kmat[41] = 0.0;
+        kmat[42] = 0.0;
+        kmat[43] = 0.0;
+        kmat[44] = 0.0;
+        kmat[45] = -bb[1];
+        kmat[46] = kap[1];
+        kmat[47] = kap[4];
+        kmat[48] = kap[7];
+        kmat[49] = 0.0;
+        kmat[50] = 1.0;
+        kmat[51] = 0.0;
+        kmat[52] = 0.0;
+        kmat[53] = 0.0;
+        kmat[54] = -bb[2];
+        kmat[55] = kap[2];
+        kmat[56] = kap[5];
+        kmat[57] = kap[8];
+        kmat[58] = 0.0;
+        kmat[59] = 0.0;
+        kmat[60] = 1.0;
+        kmat[61] = 0.0;
+        kmat[62] = 0.0;
+        kmat[63] = 0.0;
+        kmat[64] = 0.0;
+        kmat[65] = 0.0;
+        kmat[66] = 0.0;
+        kmat[67] = 0.0;
+        kmat[68] = 0.0;
+        kmat[69] = 0.0;
+        kmat[70] = 0.0;
+        kmat[71] = 1.0;
+        kmat[72] = mmu;
+        kmat[73] = 0.0;
+        kmat[74] = 0.0;
+        kmat[75] = 0.0;
+        kmat[76] = 0.0;
+        kmat[77] = 0.0;
+        kmat[78] = 0.0;
+        kmat[79] = 1.0;
+        kmat[80] = 0.0;
+    }
+};
+
+
+
+class AdvectionDiffusionMXSolver
+{
+public:
+    AdvectionDiffusionMXSolver(mfem::ParMesh* pmesh_, int order_=1)
+    {
+        pmesh=pmesh_;
+        order=order_;
+
+        int dim=pmesh->Dimension();
+
+        ufec=new mfem::H1_FECollection(order,dim);
+        pfec=new mfem::RT_FECollection(order,dim);
+        zfec=new mfem::L2_FECollection(order,dim);
+
+
+        ufes=new mfem::ParFiniteElementSpace(pmesh,ufec);
+        pfes=new mfem::ParFiniteElementSpace(pmesh,pfec);
+        zfes=new mfem::ParFiniteElementSpace(pmesh,zfec);
+
+        sfes.Append(ufes);
+        sfes.Append(pfes);
+        sfes.Append(zfes);
+
+        nf=new mfem::ParBlockNonlinearForm(sfes);
+
+        rhs.Update(nf->GetBlockTrueOffsets()); rhs=0.0;
+        sol.Update(nf->GetBlockTrueOffsets()); sol=0.0;
+        adj.Update(nf->GetBlockTrueOffsets()); adj=0.0;
+
+        fprim.SetSpace(ufes);
+        fflux.SetSpace(pfes);
+        fmult.SetSpace(zfes);
+
+        SetSolver();
+
+        dicoef=nullptr;
+        mucoef=nullptr;
+        vecoef=nullptr;
+        incoef=nullptr;
+
+        pmat=nullptr;
+        prec=nullptr;
+        psol=nullptr;
+
+    }
+
+    ~AdvectionDiffusionMXSolver()
+    {
+
+        delete psol;
+        delete prec;
+        delete pmat;
+
+        delete nf;
+
+        delete ufes;
+        delete pfes;
+        delete zfes;
+
+        delete ufec;
+        delete pfec;
+        delete zfec;
+
+        for(auto it=bc.begin();it!=bc.end();it++)
+        {
+            delete *it;
+        }
+    }
+
+    void SetSolver(double rtol=1e-8, double atol=1e-12,int miter=1000, int prt_level=1)
+    {
+        rel_tol=rtol;
+        abs_tol=atol;
+        max_iter=miter;
+        print_level=prt_level;
+    }
+
+
+    void SetDiffusion(mfem::Coefficient* coef_)
+    {
+        dicoef=coef_;
+    }
+
+    void SetReaction(mfem::Coefficient* coef_)
+    {
+        mucoef=coef_;
+    }
+
+    void SetVelocity(mfem::VectorCoefficient* coef_)
+    {
+        vecoef=coef_;
+    }
+
+    void SetLoad(mfem::Coefficient* coef_)
+    {
+        incoef=coef_;
+    }
+
+    /// Solves the forward problem.
+    void FSolve();
+
+    /// Solves the adjoint with the provided rhs.
+    void ASolve(mfem::BlockVector& rhs);
+
+
+    mfem::ParGridFunction& GetPrimField()
+    {
+        fprim.SetFromTrueDofs(sol.GetBlock(0));
+        return fprim;
+    }
+
+    mfem::ParGridFunction& GetFluxField()
+    {
+        fflux.SetFromTrueDofs(sol.GetBlock(1));
+        return fflux;
+    }
+
+    mfem::ParGridFunction& GetMultField()
+    {
+        fmult.SetFromTrueDofs(sol.GetBlock(2));
+        return fmult;
+    }
+
+    void AddDirichletBC(int mark, double val)
+    {
+        int ni=bc.size();
+        bc.push_back(new mfem::ConstantCoefficient(val));
+        mfem::Array<int> markers(pmesh->bdr_attributes.Max());
+        markers=0;
+        markers[mark-1]=1;
+        bcc[bc[ni]]=markers;
+    }
+
+    void AddDirichletBC(int mark, mfem::Coefficient* cc)
+    {
+        //check if cc is already in
+        auto it=bcc.find(cc);
+        if(it!=bcc.end())
+        {
+            (it->second)[mark-1]=1;
+        }
+        else{
+            mfem::Array<int> markers(pmesh->bdr_attributes.Max());
+            markers=0;
+            markers[mark-1]=1;
+            bcc[cc]=markers;
+        }
+    }
+
+private:
+
+    mfem::Coefficient* dicoef; //diffusion
+    mfem::Coefficient* mucoef; //reaction coefficient
+    mfem::VectorCoefficient* vecoef; //velocity
+    mfem::Coefficient* incoef; //input
+
+    mfem::ParMesh* pmesh;
+    int order;
+
+    std::vector<mfem::ConstantCoefficient*> bc;
+
+    std::map<mfem::Coefficient*, mfem::Array<int>> bcc;
+
+
+    mfem::ParFiniteElementSpace* ufes;
+    mfem::ParFiniteElementSpace* pfes;
+    mfem::ParFiniteElementSpace* zfes;
+
+    mfem::Array<mfem::ParFiniteElementSpace*> sfes;
+
+    mfem::FiniteElementCollection* ufec;
+    mfem::FiniteElementCollection* pfec;
+    mfem::FiniteElementCollection* zfec;
+
+    std::vector<mfem::AdvectionDiffusionMX*> nfin;
+    mfem::ParBlockNonlinearForm* nf;
+
+
+    mfem::BlockVector rhs;
+    mfem::BlockVector sol;
+    mfem::BlockVector adj;
+
+    // forward fields
+    mfem::ParGridFunction fprim;
+    mfem::ParGridFunction fflux;
+    mfem::ParGridFunction fmult;
+
+    // adjoint fields
+    mfem::ParGridFunction aprim;
+    mfem::ParGridFunction agrad;
+    mfem::ParGridFunction amult;
+
+    /// The PETSc objects are allocated once the problem is
+    /// assembled. They are utilized in computing the adjoint
+    /// solutions.
+    mfem::PetscParMatrix* pmat;
+    mfem::PetscPreconditioner* prec;
+    mfem::PetscLinearSolver*   psol;
+
+
+    double abs_tol;
+    double rel_tol;
+    int print_level;
+    int max_iter;
+
+    void DirectSolver(mfem::BlockOperator& A);
+    void PETSCSolver(mfem::BlockOperator& A);
+};
+
+
+
+
+
+
+
+
+}
+
+
+
+
+
+#endif

miniapps/autodiff/advdiff_test.cpp

@@ -0,0 +1,162 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+#include "advdiff.hpp"
+
+int main(int argc, char* argv[])
+{
+
+    // Initialize MPI.
+    int nprocs, myrank;
+    MPI_Init(&argc, &argv);
+    MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+
+
+
+    // Parse command-line options.
+    const char *mesh_file = "../../data/star.mesh";
+    int order = 1;
+    bool static_cond = false;
+    int ser_ref_levels = 1;
+    int par_ref_levels = 1;
+    double newton_rel_tol = 1e-7;
+    double newton_abs_tol = 1e-12;
+    int newton_iter = 10;
+    int print_level = 1;
+    bool visualization = false;
+
+    const char *petscrc_file = "advdiff_fieldsplit";
+
+    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",
+                   "Finite element order (polynomial degree) or -1 for"
+                   " isoparametric space.");
+    args.AddOption(&visualization,
+                   "-vis",
+                   "--visualization",
+                   "-no-vis",
+                   "--no-visualization",
+                   "Enable or disable GLVis visualization.");
+    args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                   "--no-static-condensation", "Enable static condensation.");
+    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(&petscrc_file, "-petscopts", "--petscopts",
+                                  "PetscOptions file to use.");
+    args.Parse();
+    if (!args.Good())
+    {
+       if (myrank == 0)
+       {
+          args.PrintUsage(std::cout);
+       }
+       MPI_Finalize();
+       return 1;
+    }
+
+    if (myrank == 0)
+    {
+       args.PrintOptions(std::cout);
+    }
+
+    mfem::MFEMInitializePetsc(NULL,NULL,petscrc_file,NULL);
+
+    // Read the (serial) mesh from the given mesh file on all processors.  We
+    // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+    // and volume meshes with the same code.
+    mfem::Mesh mesh(mesh_file, 1, 1);
+    int dim = mesh.Dimension();
+
+    // Refine the serial mesh on all processors 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 10,000 elements.
+    {
+       int ref_levels =
+          (int)floor(log(100./mesh.GetNE())/log(2.)/dim);
+       for (int l = 0; l < ref_levels; l++)
+       {
+          mesh.UniformRefinement();
+       }
+    }
+
+    // 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(MPI_COMM_WORLD, mesh);
+    mesh.Clear();
+    /*
+    {
+       for (int l = 0; l < par_ref_levels; l++)
+       {
+          pmesh.UniformRefinement();
+       }
+    }*/
+
+    mfem::AdvectionDiffusionMXSolver* solver=new mfem::AdvectionDiffusionMXSolver(&pmesh,order);
+
+    mfem::ConstantCoefficient dicoef(1.0);
+    mfem::ConstantCoefficient mucoef(0.0);
+    mfem::ConstantCoefficient incoef(1.0);
+
+    mfem::Vector veloc(3); veloc=10.0; veloc(2)=0.0;
+    mfem::VectorConstantCoefficient vecoef(veloc);
+
+    solver->AddDirichletBC(1,0.0);
+    solver->AddDirichletBC(2,3.0);
+
+    solver->SetDiffusion(&dicoef);
+    solver->SetLoad(&incoef);
+    solver->SetReaction(&mucoef);
+    solver->SetVelocity(&vecoef);
+
+    solver->FSolve();
+
+
+    {
+        mfem::ParGridFunction& fprim=solver->GetPrimField();
+        mfem::ParGridFunction& fflux=solver->GetFluxField();
+        mfem::ParGridFunction& fmult=solver->GetMultField();
+
+        mfem::ParaViewDataCollection paraview_dc("AdvectionDiffusion", &pmesh);
+        paraview_dc.SetPrefixPath("ParaView");
+        paraview_dc.SetLevelsOfDetail(order);
+        paraview_dc.SetDataFormat(mfem::VTKFormat::BINARY);
+        paraview_dc.SetHighOrderOutput(true);
+        paraview_dc.SetCycle(0);
+        paraview_dc.SetTime(0.0);
+        paraview_dc.RegisterField("prim",&fprim);
+        paraview_dc.RegisterField("flux",&fflux);
+        paraview_dc.RegisterField("mult",&fmult);
+        paraview_dc.Save();
+
+    }
+
+
+    delete solver;
+    mfem::MFEMFinalizePetsc();
+    MPI_Finalize();
+    return 0;
+}

miniapps/autodiff/example.hpp

@@ -0,0 +1,844 @@
+// Shared implementation ex71p/ex71 for the AD integrands and the manually
+// implemented integrators
+
+#ifndef ADEXAMPLE_HPP
+#define ADEXAMPLE_HPP
+
+#include "mfem.hpp"
+#include "admfem.hpp"
+#include <memory>
+#include <iostream>
+#include <fstream>
+
+namespace mfem
+{
+
+///Example: Implementation of the energy and the residual
+/// for p-Laplacian problem. Both, the energy and the residual
+/// are evaluated at the integration points for PDE parameters
+/// vparam and state fields (derivatives with respect to x,y,z
+/// and primal field) stored in vector uu.
+template<typename TDataType, typename TParamVector, typename TStateVector,
+         int residual_size, int state_size, int param_size>
+class MyVFunctor
+{
+public:
+   ///The operator returns the first derivative of the energy with respect
+   /// to all state variables. These are set in vector uu and consist of the
+   /// derivatives with respect to x,y,z and the primal field. The derivative
+   /// is stored in vector rr with length equal to the length of vector uu.
+   void operator()(TParamVector &vparam, TStateVector &uu, TStateVector &rr)
+   {
+      MFEM_ASSERT(residual_size==4,"PLaplacianResidual residual_size should be equal to 4!")
+      double pp = vparam[0];
+      double ee = vparam[1];
+      double ff = vparam[2];
+
+      TDataType norm2 = uu[0] * uu[0] + uu[1] * uu[1] + uu[2] * uu[2];
+      TDataType 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;
+   }
+};
+
+///Defines template class (functor) for evaluating the energy
+/// of the p-Laplacian problem. The input parameters vparam are:
+/// vparam[0] - the p-Laplacian power, vparam[1] small value
+/// ensuring exsitance of an unique solution, and vparam[2] -
+/// the distributed extenal input to the PDE.
+template<typename TDataType, typename TParamVector, typename TStateVector
+         , int state_size, int param_size>
+class MyQFunctor
+{
+public:
+   ///Returns the energy of a  p-Laplacian for state field input
+   /// provided in vector uu and parameters provided in vector
+   /// vparam.
+   TDataType operator()(TParamVector &vparam, TStateVector &uu)
+   {
+      MFEM_ASSERT(state_size==4,"MyQFunctor state_size should be equal to 4!");
+      MFEM_ASSERT(param_size==3,"MyQFunctor param_size should be equal to 3!");
+      double pp = vparam[0];
+      double ee = vparam[1];
+      double ff = vparam[2];
+
+      TDataType u = uu[3];
+      TDataType norm2 = uu[0] * uu[0] + uu[1] * uu[1] + uu[2] * uu[2];
+
+      TDataType rez = pow(ee * ee + norm2, pp / 2.0) / pp - ff * u;
+      return rez;
+   }
+};
+
+
+///Implements integrator for a p-Laplacian problem.
+/// The integrator is based on a class QFunction utilized for
+/// evaluating the energy, the first derivative (residual) and
+/// the Hessian of the energy (the Jacobian of the residual).
+template<class CQVectAutoDiff>
+class pLaplaceAD : public NonlinearFormIntegrator
+{
+protected:
+   Coefficient *pp;
+   Coefficient *coeff;
+   Coefficient *load;
+
+   CQVectAutoDiff rdf;
+
+public:
+   pLaplaceAD()
+   {
+      coeff = nullptr;
+      pp = nullptr;
+   }
+
+   pLaplaceAD(Coefficient &pp_) : pp(&pp_), coeff(nullptr), load(nullptr) {}
+
+   pLaplaceAD(Coefficient &pp_, Coefficient &q, Coefficient &ld_)
+      : pp(&pp_), coeff(&q), load(&ld_)
+   {}
+
+   virtual ~pLaplaceAD() {}
+
+   virtual double GetElementEnergy(const FiniteElement &el,
+                                   ElementTransformation &trans,
+                                   const Vector &elfun) override
+   {
+      double energy = 0.0;
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      bool square = (ndim == spaceDim);
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      Vector shapef(ndof);
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      Vector grad(spaceDim);
+
+      Vector vparam(3); //[power, epsilon, load]
+      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 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
+         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+         // dshape_xyz should be divided 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 positiveness 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;
+         // the energy is taken directly the the templated function
+         MyQFunctor<double,Vector,Vector,4,3> qfunc;
+         energy = energy + w * qfunc(vparam,uu);
+      }
+      return energy;
+   }
+
+   virtual void AssembleElementVector(const FiniteElement &el,
+                                      ElementTransformation &trans,
+                                      const Vector &elfun,
+                                      Vector &elvect) override
+   {
+MFEM_PERF_BEGIN("AssembleElementVector");
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      Vector shapef(ndof);
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      Vector lvec(ndof);
+      elvect.SetSize(ndof);
+      elvect = 0.0;
+
+      DenseMatrix B(ndof, 4); //[diff_x,diff_y,diff_z, shape]
+      Vector vparam(3);       //[power, epsilon, load]
+      Vector uu(4);           //[diff_x,diff_y,diff_z,u]
+      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;
+
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const 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);
+         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 positiveness 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
+         rdf.QVectorFunc(vparam,uu,du);
+         B.Mult(du, lvec);
+         elvect.Add(w, lvec);
+      } // end integration loop
+MFEM_PERF_END("AssembleElementVector");
+   }
+
+   virtual void AssembleElementGrad(const FiniteElement &el,
+                                    ElementTransformation &trans,
+                                    const Vector &elfun,
+                                    DenseMatrix &elmat) override
+   {
+MFEM_PERF_BEGIN("AssembleElementGrad");
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      Vector shapef(ndof);
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      elmat.SetSize(ndof, ndof);
+      elmat = 0.0;
+
+      DenseMatrix B(ndof, 4); // [diff_x,diff_y,diff_z, shape]
+      DenseMatrix A(ndof, 4);
+      Vector vparam(3); // [power, epsilon, load]
+      Vector uu(4);     // [diff_x,diff_y,diff_z,u]
+      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;
+
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         trans.SetIntPoint(&ip);
+         w = trans.Weight();
+         w = ip.weight * w;
+
+         el.CalcDShape(ip, dshape_iso);
+         el.CalcShape(ip, shapef);
+         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 positiveness 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
+         rdf.QJacobian(vparam,uu,duu);
+         Mult(B, duu, A);
+         AddMult_a_ABt(w, A, B, elmat);
+
+      } // end integration loop
+MFEM_PERF_END("AssembleElementGrad");
+   }
+};
+
+///Implements hand-coded integrator for a p-Laplacian problem.
+/// Utilized as alternative for the  pLaplaceAD class based on
+/// automatic differentiation.
+class pLaplace : public NonlinearFormIntegrator
+{
+protected:
+   Coefficient *pp;
+   Coefficient *coeff;
+   Coefficient *load;
+
+public:
+   pLaplace()
+   {
+      coeff = nullptr;
+      pp = nullptr;
+   }
+
+   pLaplace(Coefficient &pp_) : pp(&pp_), coeff(nullptr), load(nullptr) {}
+
+   pLaplace(Coefficient &pp_, Coefficient &q, Coefficient &ld_)
+      : pp(&pp_), coeff(&q), load(&ld_)
+   {}
+
+   virtual ~pLaplace() {}
+
+   virtual double GetElementEnergy(const FiniteElement &el,
+                                   ElementTransformation &trans,
+                                   const Vector &elfun) override
+   {
+      double energy = 0.0;
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      bool square = (ndim == spaceDim);
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      Vector shapef(ndof);
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      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 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
+         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+         // dshape_xyz should be divided 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 positiveness 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 FiniteElement &el,
+                                      ElementTransformation &trans,
+                                      const Vector &elfun,
+                                      Vector &elvect) override
+   {
+
+MFEM_PERF_BEGIN("AssembleElementVector");
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      bool square = (ndim == spaceDim);
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      Vector shapef(ndof);
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      Vector grad(spaceDim);
+      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 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
+         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+         // dshape_xyz should be divided 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 positiveness 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
+MFEM_PERF_END("AssembleElementVector");
+   }
+
+   virtual void AssembleElementGrad(const FiniteElement &el,
+                                    ElementTransformation &trans,
+                                    const Vector &elfun,
+                                    DenseMatrix &elmat) override
+   {
+MFEM_PERF_BEGIN("AssembleElementGrad");
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      bool square = (ndim == spaceDim);
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      Vector grad(spaceDim);
+      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 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
+         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+         // dshape_xyz should be divided 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 positiveness 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);
+         AddMult_a_VVt(w * aa0 / (detJ * detJ), lvec, elmat);
+         AddMult_a_AAt(w * aa1, dshape_xyz, elmat);
+
+      } // end integration loop
+MFEM_PERF_END("AssembleElementGrad");
+   }
+};
+
+/// Implements AD enabled integrator for a p-Laplacian problem.
+/// The tangent matrix is computed using the residual of the
+/// element.
+
+template<int sizeres=10>
+class pLaplaceSL : public NonlinearFormIntegrator
+{
+protected:
+   Coefficient *pp;
+   Coefficient *coeff;
+   Coefficient *load;
+
+public:
+   pLaplaceSL()
+   {
+      coeff = nullptr;
+      pp = nullptr;
+   }
+
+   pLaplaceSL(Coefficient &pp_) : pp(&pp_), coeff(nullptr), load(nullptr) {}
+
+   pLaplaceSL(Coefficient &pp_, Coefficient &q, Coefficient &ld_)
+      : pp(&pp_), coeff(&q), load(&ld_)
+   {}
+
+   virtual ~pLaplaceSL() {}
+
+   virtual double GetElementEnergy(const FiniteElement &el,
+                                   ElementTransformation &trans,
+                                   const Vector &elfun) override
+   {
+      double energy = 0.0;
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      bool square = (ndim == spaceDim);
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      Vector shapef(ndof);
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      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 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
+         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+         // dshape_xyz should be divided 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 positiveness 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 FiniteElement &el,
+                                      ElementTransformation &trans,
+                                      const Vector &elfun,
+                                      Vector &elvect) override
+   {
+MFEM_PERF_BEGIN("AssembleElementVector");
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      bool square = (ndim == spaceDim);
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      Vector shapef(ndof);
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      Vector grad(spaceDim);
+      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 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
+         Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+         // dshape_xyz should be divided 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 positiveness 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
+MFEM_PERF_END("AssembleElementVector");
+   }
+
+   virtual void AssembleElementGrad(const FiniteElement &el,
+                                    ElementTransformation &trans,
+                                    const Vector &elfun,
+                                    DenseMatrix &elmat) override
+   {
+MFEM_PERF_BEGIN("AssembleElementGrad");
+      int ndof = el.GetDof();
+      int ndim = el.GetDim();
+      int spaceDim = trans.GetSpaceDim();
+      bool square = (ndim == spaceDim);
+      const IntegrationRule *ir = NULL;
+      int order = 2 * el.GetOrder() + trans.OrderGrad(&el);
+      ir = &IntRules.Get(el.GetGeomType(), order);
+
+      DenseMatrix dshape_iso(ndof, ndim);
+      DenseMatrix dshape_xyz(ndof, spaceDim);
+      elmat.SetSize(ndof, ndof);
+      elmat = 0.0;
+
+      double w;
+      double detJ;
+      double ppp = 2.0;
+      double eee = 0.0;
+
+      mfem::Vector param(3); param=0.0;
+
+      // Computes the residual at an integration point.
+      // The implementation is a copy of the integration loop
+      // in AssembleElementVector.
+      auto resfun = [&](mfem::Vector& vparam, mfem::ad::ADVectorType& uu, mfem::ad::ADVectorType& vres){
+
+         vres.SetSize(uu.Size()); vres=0.0;
+         mfem::ad::ADVectorType grad(spaceDim);
+         mfem::ad::ADFloatType nrgrad;
+         mfem::ad::ADFloatType aa;
+         mfem::ad::ADVectorType lvec(ndof);
+
+         for (int i = 0; i < ir->GetNPoints(); i++)
+         {
+             lvec=0.0;
+
+             const 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
+             Mult(dshape_iso, trans.AdjugateJacobian(), dshape_xyz);
+             // dshape_xyz should be divided 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 positiveness of the tangent matrix
+             if (coeff != nullptr)
+             {
+                eee = coeff->Eval(trans, ip);
+             }
+
+             grad=0.0;
+             //calculate the gradient
+             for(int i=0;i<spaceDim;i++){
+             for(int j=0;j<ndof;j++){
+                grad[i]= grad[i]+ dshape_xyz(j,i)*uu[j];
+             }}
+
+
+             nrgrad= (grad*grad)/(detJ*detJ);
+
+             aa = nrgrad + eee * eee;
+             aa = pow(aa, (ppp - 2.0) / 2.0);
+
+             for(int i=0;i<spaceDim;i++){
+             for(int j=0;j<ndof;j++){
+                lvec[j] = lvec[j] + dshape_xyz(j,i) * grad[i];
+             }}
+
+             for(int j=0;j<ndof;j++)
+             {
+                vres[j]=vres[j] + lvec[j] * (w*aa/(detJ*detJ));
+             }
+         }
+      };
+
+      mfem::Vector bla(elfun);
+      //calculate the gradient - only for a fixed ndof
+      mfem::VectorFuncAutoDiff<sizeres,sizeres,3> fdr(resfun);
+      fdr.QJacobian(param, bla, elmat);
+MFEM_PERF_END("AssembleElementGrad");
+   }
+
+};
+
+} // namespace mfem
+#endif

miniapps/autodiff/fdual.hpp

@@ -0,0 +1,548 @@
+// 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
+{
+/** The FDual template class provides forward automatic differentiation
+   implementation based on dual numbers. The derivative of an arbitrary
+ function double f(double a)  can be obtained by replacing the double
+ type for the return value and the argument a  with FDual<double>, i.e.,
+ FDual<double> f(FDual<double> a). The derivative is evaluated automatically
+ by calling the function r=f(a). The value of the function is stored in r.pr
+ and the derivative in r.du. These can be extracted by the corresponding
+  methods real()/prim() and dual(). Internally, the function f can be
+ composed of standard functions predefined for FDual type. These consist
+ of a large set of functions replicating the functionality of the standard
+ math library, i.e., sin, cos, exp, log, ...
+
+New functions (non-member) can be eaily added to the class. Example:
+
+template<typename tbase>
+inline FDual<tbase> cos(const FDual<tbase> &f)
+{
+   return FDual<tbase>(cos(f.real()), -f.dual() * sin(f.real()));
+}
+
+The real part of the return value consists of the standard real value
+ of the function, i.e., cos(f.real()).
+
+The dual part of the return value consists of the first derivative of
+ the function with respect to the real part of the argument -sin(f.reaf)
+ multiplied with the dual part of the argument f.dual().
+*/
+template<typename tbase>
+class FDual
+{
+private:
+   /// Real value
+   tbase pr;
+   /// Dual value holding derivative information
+   tbase du;
+
+public:
+   /// Standard contructor - both values are set to zero.
+   FDual() : pr(0), du(0) {}
+
+   /// The constructor utilized in nested definition of dual numbers.
+   /// It is used for second and higher order derivatives.
+   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)
+   {}
+
+   /// Standard constructor with user supplied inpur for both parts of the dual number.
+   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) {}
+
+   /// Return the real value of the dual number.
+   tbase prim() const { return pr; }
+
+   /// Same as prim(). Return the real value of the dual number.
+   tbase real() const { return pr; }
+
+   /// Return the dual value of the dual number.
+   tbase dual() const { return du; }
+
+   /// Set the primal and the dual values.
+   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));
+}
+
+} // namespace ad
+
+} // namespace mfem
+
+#endif

miniapps/autodiff/par_example.cpp

@@ -0,0 +1,524 @@
+//                       MFEM Example 71 - Parallel Version
+//
+// Compile with: make ex71p
+//
+// Sample runs:
+//   mpirun -np 2 ex71p -m ../data/beam-quad.mesh -pp 3.8
+//   mpirun -np 2 ex71p -m ../data/beam-tri.mesh  -pp 7.2
+//   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
+//               p-Laplacian 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 the manually implemented integrator.  Selecting
+//           integrator=1,2 will utilize the AD integrator.
+//
+//        The AD integrators are implemented in ex71.hpp (pLaplaceAD).
+//           The integrand qint is a function which is evaluated at every
+//           integration point. For implementations utilizing ADQFunctionTJ,
+//           the user has to implement the function and the residual
+//           evaluation. The Jacobian of the residual is evaluated using AD
+//           For implementations utilizing ADQFunctionTH, the user has to
+//           implement only the function evaluation (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 "example.hpp"
+
+using namespace mfem;
+
+///Non-linear solver for the p-Laplacian problem.
+class ParNLSolverPLaplacian
+{
+public:
+   ///Constructor Input: imesh - FE mesh, finite element space,
+   /// power for the p-Laplacian, external load (source, input),
+   /// regularization parameter
+   ParNLSolverPLaplacian(MPI_Comm comm, ParMesh& imesh,
+                         ParFiniteElementSpace& ifespace,
+                         double powerp=2,
+                         Coefficient* load=nullptr,
+                         double regularizationp=1e-7)
+   {
+      lcomm = comm;
+
+      //default parameters for
+      //the Newton solver
+      newton_rtol = 1e-4;
+      newton_atol = 1e-8;
+      newton_iter = 10;
+
+      //linear solver
+      linear_rtol = 1e-7;
+      linear_atol = 1e-15;
+      linear_iter = 500;
+
+      print_level = 0;
+
+      //set the mesh
+      mesh=&imesh;
+
+      //set the fespace
+      fespace=&ifespace;
+
+      //set the parameters
+      plap_epsilon=new ConstantCoefficient(regularizationp);
+      plap_power=new ConstantCoefficient(powerp);
+      if (load==nullptr)
+      {
+         plap_input=new ConstantCoefficient(1.0);
+         input_ownership=true;
+      }
+      else
+      {
+         plap_input=load;
+         input_ownership=false;
+      }
+
+      nf=nullptr;
+      ns=nullptr;
+      gmres=nullptr;
+      prec=nullptr;
+
+      //set the default integrator
+      integ=0; //hand coded
+   }
+
+   ~ParNLSolverPLaplacian()
+   {
+      if (nf!=nullptr) { delete nf;}
+      if (ns!=nullptr) { delete ns;}
+      if (prec!=nullptr) { delete prec;}
+      if (gmres!=nullptr) { delete gmres;}
+      if (input_ownership) { delete plap_input;}
+      delete plap_epsilon;
+      delete plap_power;
+   }
+
+   ///Set the integrator.
+   /// 0 - hand coded, 1 - AD based (compute only Heassian by AD),
+   /// 2 - AD based (compute residual and Hessian by AD)
+   void SetIntegrator(int intr)
+   {
+      integ=intr;
+   }
+
+   //set relative tolerance for the Newton solver
+   void SetNRRTol(double rtol)
+   {
+      newton_rtol=rtol;
+   }
+
+   //set absolute tolerance for the Newton solver
+   void SetNRATol(double atol)
+   {
+      newton_atol=atol;
+   }
+
+   //set max iterations for the NR solver
+   void SetMaxNRIter(int miter)
+   {
+      newton_iter=miter;
+   }
+
+   void SetLSRTol(double rtol)
+   {
+      linear_rtol=rtol;
+   }
+
+   void SetLSATol(double atol)
+   {
+      linear_atol=atol;
+   }
+
+   //set max iterations for the linear solver
+   void SetMaxLSIter(int miter)
+   {
+      linear_iter=miter;
+   }
+
+   //set the print level
+   void SetPrintLevel(int plev)
+   {
+      print_level=plev;
+   }
+
+   ///The state vector is used as initial condition for the NR solver.
+   /// On return the statev holds the solution to the problem.
+   void Solve(Vector& statev)
+   {
+      if (nf==nullptr)
+      {
+         AllocSolvers();
+      }
+      Vector b; //RHS is zero
+      ns->Mult(b, statev);
+   }
+
+   ///Compute the energy
+   double GetEnergy(Vector& statev)
+   {
+      if (nf==nullptr)
+      {
+         //allocate the solvers
+         AllocSolvers();
+      }
+      return nf->GetEnergy(statev);
+   }
+
+private:
+   void AllocSolvers()
+   {
+      if (nf!=nullptr) { delete nf;}
+      if (ns!=nullptr) { delete ns;}
+      if (gmres!=nullptr) { delete gmres;}
+      if (prec!=nullptr) { delete prec;}
+
+      // Define the essential boundary attributes
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+
+      nf = new ParNonlinearForm(fespace);
+      if (integ==0)
+      {
+         nf->AddDomainIntegrator(new pLaplace(*plap_power,*plap_epsilon,*plap_input));
+      }
+      else if (integ==1)
+      {
+         nf->AddDomainIntegrator(new pLaplaceAD<mfem::QVectorFuncAutoDiff<MyVFunctor,4,4,3>>(*plap_power,*plap_epsilon,*plap_input));
+      }
+      else if (integ==2)
+      {
+         nf->AddDomainIntegrator(new pLaplaceAD<mfem::QFunctionAutoDiff<MyQFunctor,4,3>>(*plap_power,*plap_epsilon,*plap_input));
+      }
+      else{
+         nf->AddDomainIntegrator(new pLaplaceSL<56>(*plap_power,*plap_epsilon,*plap_input));
+      }
+
+
+      nf->SetEssentialBC(ess_bdr);
+
+      prec = new HypreBoomerAMG();
+      prec->SetPrintLevel(print_level);
+
+      gmres = new GMRESSolver(lcomm);
+      gmres->SetAbsTol(linear_atol);
+      gmres->SetRelTol(linear_rtol);
+      gmres->SetMaxIter(linear_iter);
+      gmres->SetPrintLevel(print_level);
+      gmres->SetPreconditioner(*prec);
+
+      ns = new NewtonSolver(lcomm);
+
+      ns->iterative_mode = true;
+      ns->SetSolver(*gmres);
+      ns->SetOperator(*nf);
+      ns->SetPrintLevel(print_level);
+      ns->SetRelTol(newton_rtol);
+      ns->SetAbsTol(newton_atol);
+      ns->SetMaxIter(newton_iter);
+   }
+
+   double newton_rtol;
+   double newton_atol;
+   int newton_iter;
+
+   double linear_rtol;
+   double linear_atol;
+   int linear_iter;
+
+   int print_level;
+
+   //power of the p-laplacian
+   Coefficient* plap_power;
+   //regularization parammeter
+   Coefficient* plap_epsilon;
+   //load(input) paramater
+   Coefficient* plap_input;
+   bool input_ownership;
+
+   MPI_Comm lcomm;
+
+   ParMesh *mesh;
+   ParFiniteElementSpace *fespace;
+
+   ParNonlinearForm *nf;
+
+   HypreBoomerAMG *prec;
+   GMRESSolver *gmres;
+   NewtonSolver *ns;
+   int integ;
+
+};
+
+
+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);
+   // Define Caliper ConfigManager
+#ifdef MFEM_USE_CALIPER
+   cali::ConfigManager mgr;
+#endif
+   // Caliper instrumentation
+   MFEM_PERF_FUNCTION;
+
+   // 2. Parse command-line options
+   const char *mesh_file = "../../data/beam-tet.mesh";
+   int ser_ref_levels = 3;
+   int par_ref_levels = 1;
+   int order = 1;
+   bool visualization = true;
+   double newton_rel_tol = 1e-4;
+   double newton_abs_tol = 1e-6;
+   int newton_iter = 10;
+   int print_level = 0;
+   double pp = 2.0;
+   int integrator = 2; //use AD
+   const char* cali_config = "runtime-report";
+
+   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 for Hessian; 2: AD for residual and Hessian");
+   args.AddOption(&cali_config, "-p", "--caliper",
+                  "Caliper configuration string.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myrank == 0)
+      {
+         args.PrintUsage(std::cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myrank == 0)
+   {
+      args.PrintOptions(std::cout);
+   }
+
+   StopWatch *timer = new StopWatch();
+
+   // Caliper configuration
+#ifdef MFEM_USE_CALIPER
+   mgr.add(cali_config);
+   mgr.start();
+#endif
+   // 3. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral and hexahedral meshes
+   //    with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter.
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // 6. Define the load for the p-Laplacian
+   ConstantCoefficient load(1.00);
+
+   // 7. Define the finite element spaces for the solution
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace fespace(pmesh, &fec, 1, 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 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.
+   ParGridFunction x(&fespace);
+   x = 0.0;
+   HypreParVector *sv = x.GetTrueDofs();
+
+   // 9. Define ParaView DataCollection
+   ParaViewDataCollection *dacol = new ParaViewDataCollection("Example71",
+                                                              pmesh);
+   dacol->SetLevelsOfDetail(order);
+   dacol->RegisterField("sol", &x);
+
+   // 10. Define the NR solver
+   ParNLSolverPLaplacian* nr;
+
+   // 11. Start with linear diffusion - solvable for any initial guess
+   nr=new ParNLSolverPLaplacian(MPI_COMM_WORLD,*pmesh, fespace, 2.0, &load);
+   nr->SetIntegrator(integrator);
+   nr->SetMaxNRIter(newton_iter);
+   nr->SetNRATol(newton_abs_tol);
+   nr->SetNRRTol(newton_rel_tol);
+   nr->SetPrintLevel(print_level);
+   timer->Clear();
+   timer->Start();
+   nr->Solve(*sv);
+   timer->Stop();
+   if (myrank==0)
+   {
+      std::cout << "[pp=2] The solution time is: " << timer->RealTime()
+                << std::endl;
+   }
+   // Compute the energy
+   double energy = nr->GetEnergy(*sv);
+   if (myrank==0)
+   {
+      std::cout << "[pp=2] The total energy of the system is E=" << energy
+                << std::endl;
+   }
+   delete nr;
+   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++)
+   {
+      nr=new ParNLSolverPLaplacian(MPI_COMM_WORLD,*pmesh, fespace, (double)i, &load);
+      nr->SetIntegrator(integrator);
+      nr->SetMaxNRIter(newton_iter);
+      nr->SetNRATol(newton_abs_tol);
+      nr->SetNRRTol(newton_rel_tol);
+      nr->SetPrintLevel(print_level);
+      timer->Clear();
+      timer->Start();
+      nr->Solve(*sv);
+      timer->Stop();
+      if (myrank==0)
+      {
+         std::cout << "[pp="<<i<<"] The solution time is: " << timer->RealTime()
+                   << std::endl;
+      }
+      // Compute the energy
+      double energy = nr->GetEnergy(*sv);
+      if (myrank==0)
+      {
+         std::cout << "[pp="<<i<<"] The total energy of the system is E=" << energy
+                   << std::endl;
+      }
+      delete nr;
+      x.SetFromTrueDofs(*sv);
+      dacol->SetTime((double)i);
+      dacol->SetCycle(i);
+      dacol->Save();
+   }
+
+   // 13. Continue with the final power
+   if (std::abs(pp - 2.0) > std::numeric_limits<double>::epsilon())
+   {
+      nr=new ParNLSolverPLaplacian(MPI_COMM_WORLD,*pmesh, fespace, pp, &load);
+      nr->SetIntegrator(integrator);
+      nr->SetMaxNRIter(newton_iter);
+      nr->SetNRATol(newton_abs_tol);
+      nr->SetNRRTol(newton_rel_tol);
+      nr->SetPrintLevel(print_level);
+      timer->Clear();
+      timer->Start();
+      nr->Solve(*sv);
+      timer->Stop();
+      if (myrank==0)
+      {
+         std::cout << "[pp="<<pp<<"] The solution time is: " << timer->RealTime()
+                   << std::endl;
+      }
+      // Compute the energy
+      double energy = nr->GetEnergy(*sv);
+      if (myrank==0)
+      {
+         std::cout << "[pp="<<pp<<"] The total energy of the system is E=" << energy
+                   << std::endl;
+      }
+      delete nr;
+      x.SetFromTrueDofs(*sv);
+      dacol->SetTime(pp);
+      if (pp < 2.0)
+      {
+         dacol->SetCycle(std::floor(pp));
+      }
+      else
+      {
+         dacol->SetCycle(std::ceil(pp));
+      }
+      dacol->Save();
+   }
+
+   // 14. Free the used memory
+   delete dacol;
+   delete sv;
+   delete pmesh;
+   delete timer;
+
+   // Flush output before MPI_finalize
+#ifdef MFEM_USE_CALIPER
+   mgr.flush();
+#endif
+   MPI_Finalize();
+
+   return 0;
+}

miniapps/autodiff/seq_example.cpp

@@ -0,0 +1,454 @@
+//                       MFEM Example 71 - Serial Version
+//
+// Compile with: make ex71
+//
+// Sample runs:
+//    ex71 -m ../data/beam-quad.mesh -pp 3.5
+//    ex71 -m ../data/beam-tri.mesh  -pp 4.6
+//    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
+//               p-Laplacian 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 the manually implemented integrator. Selecting
+//           integrator=1,2 will utilize one of the AD integrators.
+//
+//           The AD integrators are implemented in ex71.hpp (pLaplaceAD).
+//           The integrand qint is a function which is evaluated at every
+//           integration point. For implementations utilizing ADQFunctionTJ,
+//           the user has to implement the function and the residual
+//           evaluation. The Jacobian of the residual is evaluated using AD
+//
+//           For implementations utilizing ADQFunctionTH, the user has to
+//           implement only the function evaluation (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 "example.hpp"
+
+using namespace mfem;
+
+///Non-linear solver for the p-Laplacian problem.
+class NLSolverPLaplacian
+{
+public:
+   ///Constructor Input: imesh - FE mesh, finite element space,
+   /// power for the p-Laplacian, external load (source, input),
+   /// regularization parameter
+   NLSolverPLaplacian(Mesh& imesh, FiniteElementSpace& ifespace,
+                      double powerp=2,
+                      Coefficient* load=nullptr,
+                      double regularizationp=1e-7)
+   {
+      //default parameters for
+      //the Newton solver
+      newton_rtol = 1e-4;
+      newton_atol = 1e-8;
+      newton_iter = 10;
+
+      //linear solver
+      linear_rtol = 1e-7;
+      linear_atol = 1e-15;
+      linear_iter = 500;
+
+      print_level = 0;
+
+      //set the mesh
+      mesh=&imesh;
+
+      //set the fespace
+      fespace=&ifespace;
+
+      //set the parameters
+      plap_epsilon=new ConstantCoefficient(regularizationp);
+      plap_power=new ConstantCoefficient(powerp);
+      if (load==nullptr)
+      {
+         plap_input=new ConstantCoefficient(1.0);
+         input_ownership=true;
+      }
+      else
+      {
+         plap_input=load;
+         input_ownership=false;
+      }
+
+      //set the nonlinear form
+      nf=nullptr;
+      lsolver=nullptr;
+      prec=nullptr;
+      ns=nullptr;
+
+      //set the default integrator
+      integ=0; //hand coded
+
+   }
+
+   ~NLSolverPLaplacian()
+   {
+      if (nf!=nullptr) { delete nf;}
+      if (ns!=nullptr) { delete ns;}
+      if (prec!=nullptr) { delete prec;}
+      if (lsolver!=nullptr) { delete lsolver;}
+      if (input_ownership) { delete plap_input;}
+      delete plap_epsilon;
+      delete plap_power;
+   }
+
+   ///Set the integrator.
+   /// 0 - hand coded, 1 - AD based (compute only Heassian by AD),
+   /// 2 - AD based (compute residual and Hessian by AD)
+   void SetIntegrator(int intr)
+   {
+      integ=intr;
+   }
+
+
+   //set relative tolerance for the Newton solver
+   void SetNRRTol(double rtol)
+   {
+      newton_rtol=rtol;
+   }
+
+   //set absolute tolerance for the Newton solver
+   void SetNRATol(double atol)
+   {
+      newton_atol=atol;
+   }
+
+   //set max iterations for the NR solver
+   void SetMaxNRIter(int miter)
+   {
+      newton_iter=miter;
+   }
+
+   void SetLSRTol(double rtol)
+   {
+      linear_rtol=rtol;
+   }
+
+   void SetLSATol(double atol)
+   {
+      linear_atol=atol;
+   }
+
+   //set max iterations for the linear solver
+   void SetMaxLSIter(int miter)
+   {
+      linear_iter=miter;
+   }
+
+   //set the print level
+   void SetPrintLevel(int plev)
+   {
+      print_level=plev;
+   }
+
+   ///The state vector is used as initial condition for the NR solver.
+   /// On return the statev holds the solution to the problem.
+   void Solve(Vector& statev)
+   {
+      if (nf==nullptr)
+      {
+         AllocSolvers();
+      }
+      Vector b; //RHS is zero
+      ns->Mult(b, statev);
+   }
+
+   ///Compute the energy
+   double GetEnergy(Vector& statev)
+   {
+      if (nf==nullptr)
+      {
+         //allocate the solvers
+         AllocSolvers();
+      }
+      return nf->GetEnergy(statev);
+   }
+
+
+private:
+
+   void AllocSolvers()
+   {
+      if (nf!=nullptr) { delete nf;}
+      if (ns!=nullptr) {delete ns;}
+      if (prec!=nullptr) {delete prec;}
+      if (lsolver!=nullptr) { delete lsolver;}
+
+      // Define the essential boundary attributes
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+
+      nf = new NonlinearForm(fespace);
+
+      if (integ==0)
+      {
+         nf->AddDomainIntegrator(new pLaplace(*plap_power,*plap_epsilon,*plap_input));
+      }
+      else if (integ==1)
+      {
+         nf->AddDomainIntegrator(new pLaplaceAD<mfem::QVectorFuncAutoDiff<MyVFunctor,4,4,3>>(*plap_power,*plap_epsilon,*plap_input));
+      }
+      else
+      {
+         nf->AddDomainIntegrator(new pLaplaceAD<mfem::QFunctionAutoDiff<MyQFunctor,4,3>>(*plap_power,*plap_epsilon,*plap_input));
+      }
+
+      nf->SetEssentialBC(ess_bdr);
+
+#ifdef MFEM_USE_SUITESPARSE
+      prec = new UMFPackSolver();
+#else
+      prec = new GSSmoother();
+#endif
+
+      //allocate the linear solver
+      lsolver=new CGSolver();
+      lsolver->SetRelTol(linear_rtol);
+      lsolver->SetAbsTol(linear_atol);
+      lsolver->SetMaxIter(linear_iter);
+      lsolver->SetPrintLevel(print_level);
+      lsolver->SetPreconditioner(*prec);
+
+      //allocate the NR solver
+      ns = new NewtonSolver();
+      ns->iterative_mode = true;
+      ns->SetSolver(*lsolver);
+      ns->SetOperator(*nf);
+      ns->SetPrintLevel(print_level);
+      ns->SetRelTol(newton_rtol);
+      ns->SetAbsTol(newton_atol);
+      ns->SetMaxIter(newton_iter);
+   }
+
+   double newton_rtol;
+   double newton_atol;
+   int newton_iter;
+
+   double linear_rtol;
+   double linear_atol;
+   int linear_iter;
+
+   int print_level;
+
+
+   //reference to the mesh
+   Mesh* mesh;
+   //reference to the fespace
+   FiniteElementSpace *fespace;
+
+   //nonlinear form for the p-laplacian
+   NonlinearForm *nf;
+   CGSolver *lsolver; //linear solver
+   Solver *prec; //preconditioner for the linear solver
+   NewtonSolver *ns; //NR solver
+   int integ;
+
+   //power of the p-laplacian
+   Coefficient* plap_power;
+   //regularization parammeter
+   Coefficient* plap_epsilon;
+   //load(input) paramater
+   Coefficient* plap_input;
+   bool input_ownership;
+};
+
+
+
+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 = 10;
+   int print_level = 0;
+
+   double pp = 2.0; // p-Lapalacian power
+
+   int integrator = 2; // 2 - use AD for Residual and Hessian
+   // 1 - use AD for Hessian only
+   // 0 - do not use AD (hand coded)
+   StopWatch *timer = new StopWatch();
+   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 for Hessian; 2: AD for residual and Hessian");
+
+   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.
+   Mesh *mesh = new 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 load parameter for the p-Laplacian
+   ConstantCoefficient load(1.00);
+
+   // 5. Define the finite element spaces for the solution
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace fespace(mesh, &fec, 1, Ordering::byVDIM);
+   int glob_size = fespace.GetTrueVSize();
+
+   std::cout << "Number of finite element unknowns: " << glob_size << std::endl;
+
+   // 6. Define the solution grid function
+   GridFunction x(&fespace);
+   x = 0.0;
+
+   // 7. Define the solution true vector
+   Vector sv(fespace.GetTrueVSize());
+   sv = 0.0;
+
+   // 8. Define ParaView DataCollection
+   ParaViewDataCollection *dacol = new ParaViewDataCollection("Example71", mesh);
+   dacol->SetLevelsOfDetail(order);
+   dacol->RegisterField("sol", &x);
+
+   // 9. Define the NR solver
+   NLSolverPLaplacian* nr;
+
+   // 10. Start with linear diffusion - solvable for any initial guess
+   nr=new NLSolverPLaplacian(*mesh, fespace, 2.0, &load);
+   nr->SetIntegrator(integrator);
+   nr->SetMaxNRIter(newton_iter);
+   nr->SetNRATol(newton_abs_tol);
+   nr->SetNRRTol(newton_rel_tol);
+   timer->Clear();
+   timer->Start();
+   nr->Solve(sv);
+   timer->Stop();
+   std::cout << "[pp=2] The solution time is: " << timer->RealTime()
+             << std::endl;
+   // Compute the energy
+   double energy = nr->GetEnergy(sv);
+   std::cout << "[pp=2] The total energy of the system is E=" << energy
+             << std::endl;
+   delete nr;
+   x.SetFromTrueDofs(sv);
+   dacol->SetTime(2.0);
+   dacol->SetCycle(2);
+   dacol->Save();
+
+
+   // 11. Continue with powers higher than 2
+   for (int i = 3; i < pp; i++)
+   {
+      nr=new NLSolverPLaplacian(*mesh, fespace, (double)i, &load);
+      nr->SetIntegrator(integrator);
+      nr->SetMaxNRIter(newton_iter);
+      nr->SetNRATol(newton_abs_tol);
+      nr->SetNRRTol(newton_rel_tol);
+      timer->Clear();
+      timer->Start();
+      nr->Solve(sv);
+      timer->Stop();
+      std::cout << "[pp=" << i
+                << "] The solution time is: " << timer->RealTime() << std::endl;
+      energy = nr->GetEnergy(sv);
+      std::cout << "[pp="<< i<<"] The total energy of the system is E=" << energy
+                << std::endl;
+      delete nr;
+      x.SetFromTrueDofs(sv);
+      dacol->SetTime(i);
+      dacol->SetCycle(i);
+      dacol->Save();
+   }
+
+   // 12. Continue with the final power
+   if (std::abs(pp - 2.0) > std::numeric_limits<double>::epsilon())
+   {
+      nr=new NLSolverPLaplacian(*mesh, fespace, pp, &load);
+      nr->SetIntegrator(integrator);
+      nr->SetMaxNRIter(newton_iter);
+      nr->SetNRATol(newton_abs_tol);
+      nr->SetNRRTol(newton_rel_tol);
+      timer->Clear();
+      timer->Start();
+      nr->Solve(sv);
+      timer->Stop();
+      std::cout << "[pp=" << pp
+                << "] The solution time is: " << timer->RealTime() << std::endl;
+      energy = nr->GetEnergy(sv);
+      std::cout << "[pp="<<pp<<"] The total energy of the system is E=" << energy
+                << std::endl;
+      delete nr;
+      x.SetFromTrueDofs(sv);
+      dacol->SetTime(pp);
+      if (pp < 2.0)
+      {
+         dacol->SetCycle(std::floor(pp));
+      }
+      else
+      {
+         dacol->SetCycle(std::ceil(pp));
+      }
+      dacol->Save();
+   }
+
+   // 13. Free the memory
+   delete dacol;
+   delete mesh;
+   delete timer;
+   return 0;
+}

miniapps/autodiff/seq_test.cpp

@@ -0,0 +1,141 @@
+#include "admfem.hpp"
+#include "mfem.hpp"
+
+template<typename TDataType, typename TParamVector, typename TStateVector
+         , int state_size, int param_size>
+class DiffusionFunctional
+{
+public:
+    TDataType operator() (TParamVector& vparam, TStateVector& uu)
+    {
+        MFEM_ASSERT(state_size==4,"ExampleFunctor state_size should be equal to 4!");
+        MFEM_ASSERT(param_size==2,"ExampleFunctor param_size should be equal to 2!");
+        auto kapa = vparam[0]; //diffusion coefficient
+        auto load = vparam[1]; //volumetric influx
+        TDataType rez = kapa*(uu[0]*uu[0]+uu[1]*uu[1]+uu[2]*uu[2])/2.0 - load*uu[3];
+        return rez;
+    }
+
+};
+
+template<typename TDataType, typename TParamVector, typename TStateVector,
+         int residual_size, int state_size, int param_size>
+class DiffusionResidual
+{
+public:
+    void operator ()(TParamVector& vparam, TStateVector& uu, TStateVector& rr)
+    {
+        MFEM_ASSERT(residual_size==4,"DiffusionResidual residual_size should be equal to 4!")
+        MFEM_ASSERT(state_size==4,"ExampleFunctor state_size should be equal to 4!");
+        MFEM_ASSERT(param_size==2,"ExampleFunctor param_size should be equal to 2!");
+        auto kapa = vparam[0]; //diffusion coefficient
+        auto load = vparam[1]; //volumetric influx
+
+        rr[0] = kapa * uu[0];
+        rr[1] = kapa * uu[1];
+        rr[2] = kapa * uu[2];
+        rr[3] = -load;
+    }
+};
+
+
+
+int main(int argc, char *argv[])
+{
+
+#ifdef MFEM_USE_ADFORWARD
+    std::cout<<"MFEM_USE_ADFORWARD == true"<<std::endl;
+#else
+    std::cout<<"MFEM_USE_ADFORWARD == false"<<std::endl;
+#endif
+
+#ifdef MFEM_USE_CALIPER
+    cali::ConfigManager mgr;
+#endif
+    // Caliper instrumentation
+    MFEM_PERF_FUNCTION;
+    const char* cali_config = "runtime-report";
+#ifdef MFEM_USE_CALIPER
+    mgr.add(cali_config);
+    mgr.start();
+#endif
+    mfem::Vector param(2);
+    param[0]=3.0; //diffusion coefficient
+    param[1]=2.0; //volumetric influx
+
+    mfem::Vector state(4);
+    state[0]=1.0; // grad_x
+    state[1]=2.0; // grad_y
+    state[2]=3.0; // grad_z
+    state[3]=4.0; // state value
+
+    mfem::QFunctionAutoDiff<DiffusionFunctional,4,2> adf;
+    mfem::QVectorFuncAutoDiff<DiffusionResidual,4,4,2> rdf;
+
+    mfem::Vector rr0(4);
+    mfem::DenseMatrix hh0(4,4);
+
+    mfem::Vector rr1(4);
+    mfem::DenseMatrix hh1(4,4);
+MFEM_PERF_BEGIN("QGrad");
+    adf.QGrad(param,state,rr0);
+MFEM_PERF_END("QGrad");
+MFEM_PERF_BEGIN("QHessian");
+    adf.QHessian(param, state, hh0);
+MFEM_PERF_END("QHessian");
+    // dump out the results
+    std::cout<<"FunctionAutoDiff"<<std::endl;
+    std::cout<< adf.QEval(param,state)<<std::endl;
+    rr0.Print(std::cout);
+    hh0.Print(std::cout);
+
+MFEM_PERF_BEGIN("QJacobian");
+    rdf.QJacobian(param, state, hh1);
+MFEM_PERF_END("QJacobian");
+
+    std::cout<<"ResidualAutoDiff"<<std::endl;
+    hh1.Print(std::cout);
+
+    //using lambda expression
+    auto func = [](mfem::Vector& vparam, mfem::ad::ADVectorType& uu, mfem::ad::ADVectorType& vres) {
+    //auto func = [](auto& vparam, auto& uu, auto& vres) { //c++14
+        auto kappa = vparam[0]; //diffusion coefficient
+        auto load = vparam[1]; //volumetric influx
+
+        vres[0] = kappa * uu[0];
+        vres[1] = kappa * uu[1];
+        vres[2] = kappa * uu[2];
+        vres[3] = -load;
+    };
+
+    mfem::VectorFuncAutoDiff<4,4,2> fdr(func);
+MFEM_PERF_BEGIN("QJacobianV");
+    fdr.QJacobian(param,state,hh1); //computes the gradient of func and stores the result in hh1
+MFEM_PERF_END("QJacobianV");
+    std::cout<<"LambdaAutoDiff"<<std::endl;
+    hh1.Print(std::cout);
+
+
+    double kappa = param[0];
+    double load =  param[1];
+    //using lambda expression
+    auto func01 = [&kappa,&load](mfem::Vector& vparam, mfem::ad::ADVectorType& uu, mfem::ad::ADVectorType& vres) {
+    //auto func = [](auto& vparam, auto& uu, auto& vres) { //c++14
+
+        vres[0] = kappa * uu[0];
+        vres[1] = kappa * uu[1];
+        vres[2] = kappa * uu[2];
+        vres[3] = -load;
+    };
+
+    mfem::VectorFuncAutoDiff<4,4,2> fdr01(func01);
+MFEM_PERF_BEGIN("QJacobian1");
+    fdr01.QJacobian(param,state,hh1);
+MFEM_PERF_END("QJacobian1");
+    std::cout<<"LambdaAutoDiff 01"<<std::endl;
+    hh1.Print(std::cout);
+
+#ifdef MFEM_USE_CALIPER
+    mgr.flush();
+#endif
+}

miniapps/autodiff/stokes.cpp

@@ -0,0 +1,208 @@
+#include "stokes.hpp"
+#include "petsc.h"
+
+namespace mfem {
+
+void StokesSolver::FSolve()
+{
+
+    if(pmat!=nullptr)
+    {
+        delete psol;
+        delete prec;
+        delete pmat;
+        ess_tdofv.DeleteAll();
+    }
+
+
+    sol=0.0;
+    // Set the BC
+    ess_tdofv.DeleteAll();
+    Array<int> ess_tdofx;
+    Array<int> ess_tdofy;
+    Array<int> ess_tdofz;
+    Array<int> ess_tdofp;
+    int dim=pmesh->Dimension();
+    {
+        for(auto it=bccx.begin();it!=bccx.end();it++)
+        {
+            mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+            ess_bdr=0;
+            ess_bdr[it->first -1]=1;
+            mfem::Array<int> ess_tdof_list;
+            vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,0);
+            ess_tdofx.Append(ess_tdof_list);
+
+            mfem::VectorArrayCoefficient pcoeff(dim);
+            pcoeff.Set(0, it->second, false);
+            fvelocity.ProjectBdrCoefficient(pcoeff, ess_bdr);
+        }
+        //copy tdofsx from velocity grid function
+        {
+            fvelocity.GetTrueDofs(rhs.GetBlock(0)); // use the rhs vector as a tmp vector
+            for(int ii=0;ii<ess_tdofx.Size();ii++)
+            {
+                sol.GetBlock(0)[ess_tdofx[ii]]=rhs.GetBlock(0)[ess_tdofx[ii]];
+            }
+        }
+        ess_tdofv.Append(ess_tdofx);
+
+
+        for(auto it=bccy.begin();it!=bccy.end();it++)
+        {
+            mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+            ess_bdr=0;
+            ess_bdr[it->first -1]=1;
+            mfem::Array<int> ess_tdof_list;
+            vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,1);
+            ess_tdofy.Append(ess_tdof_list);
+
+            mfem::VectorArrayCoefficient pcoeff(dim);
+            pcoeff.Set(1, it->second, false);
+            fvelocity.ProjectBdrCoefficient(pcoeff, ess_bdr);
+        }
+        //copy tdofsy from velocity grid function
+        {
+            fvelocity.GetTrueDofs(rhs.GetBlock(0)); // use the rhs vector as a tmp vector
+            for(int ii=0;ii<ess_tdofy.Size();ii++)
+            {
+                sol.GetBlock(0)[ess_tdofy[ii]]=rhs.GetBlock(0)[ess_tdofy[ii]];
+            }
+        }
+        ess_tdofv.Append(ess_tdofy);
+
+        for(auto it=bccz.begin();it!=bccz.end();it++)
+        {
+            mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+            ess_bdr=0;
+            ess_bdr[it->first -1]=1;
+            mfem::Array<int> ess_tdof_list;
+            vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,2);
+            ess_tdofz.Append(ess_tdof_list);
+
+            mfem::VectorArrayCoefficient pcoeff(dim);
+
+            pcoeff.Set(2, it->second, false);
+            fvelocity.ProjectBdrCoefficient(pcoeff, ess_bdr);
+        }
+        //copy tdofsz from velocity grid function
+        {
+            fvelocity.GetTrueDofs(rhs.GetBlock(0)); // use the rhs vector as a tmp vector
+            for(int ii=0;ii<ess_tdofz.Size();ii++)
+            {
+                sol.GetBlock(0)[ess_tdofz[ii]]=rhs.GetBlock(0)[ess_tdofz[ii]];
+            }
+        }
+        ess_tdofv.Append(ess_tdofz);
+    }
+
+    if(nfin==nullptr)
+    {
+        nfin=new StokesIntegratorTH(viscosity,bpenal,load);
+        nf->AddDomainIntegrator(nfin);
+        nf->SetGradientType(mfem::Operator::Type::Hypre_ParCSR);
+    }
+    // set the RHS
+    nf->Mult(sol,rhs);
+    rhs.Neg();
+
+    mfem::BlockOperator& A=nf->GetGradient(sol);
+
+    mfem::HypreParMatrix* A00=static_cast<mfem::HypreParMatrix*>(&(A.GetBlock(0,0)));
+    mfem::HypreParMatrix* A01=static_cast<mfem::HypreParMatrix*>(&(A.GetBlock(0,1)));
+    mfem::HypreParMatrix* A10=static_cast<mfem::HypreParMatrix*>(&(A.GetBlock(1,0)));
+    mfem::HypreParMatrix* A11=static_cast<mfem::HypreParMatrix*>(&(A.GetBlock(1,1)));
+
+    mfem::HypreParMatrix* A00elim=A00->EliminateRowsCols(ess_tdofv);
+    mfem::HypreParMatrix* A11elim=A11->EliminateRowsCols(ess_tdofp);
+    mfem::HypreParMatrix* A01elim=A01->EliminateCols(ess_tdofp); A01->EliminateRows(ess_tdofv);
+    mfem::HypreParMatrix* A10elim=A10->EliminateCols(ess_tdofv); A10->EliminateRows(ess_tdofp);
+
+    //copy BC to RHS
+
+    for(int ii=0;ii<ess_tdofv.Size();ii++)
+    {
+        rhs.GetBlock(0)[ess_tdofv[ii]]=sol.GetBlock(0)[ess_tdofv[ii]];
+    }
+
+    delete psol;
+    delete prec;
+    delete pmat;
+    pmat= new mfem::PetscParMatrix(pmesh->GetComm(),&A, mfem::Operator::PETSC_MATAIJ);
+    //set the local block size of the matrix
+    Mat sub;
+    MatNestGetSubMat(pmat->operator petsc::Mat(),0,0,&sub);
+    MatSetBlockSize(sub,dim);
+    // construct the preconditioner
+    prec = new mfem::PetscFieldSplitSolver(pmesh->GetComm(),*pmat,"prec_");
+    // construct the linear solver
+    psol = new mfem::PetscLinearSolver(pmesh->GetComm());
+
+    /*
+    {
+        std::fstream out("full.mat",std::ios::out);
+        pmat->PrintMatlab(out);
+        out.close();
+    }
+    */
+
+
+    psol->SetOperator(*pmat);
+    psol->SetPreconditioner(*prec);
+    psol->SetAbsTol(abs_tol);
+    psol->SetRelTol(rel_tol);
+    psol->SetMaxIter(max_iter);
+    psol->SetPrintLevel(print_level);
+    psol->Mult(rhs, sol);
+
+    //psol->GetConverged();
+
+    delete A11elim;
+    delete A01elim;
+    delete A10elim;
+    delete A00elim;
+
+}
+
+
+void StokesSolver::ASolve(BlockVector& arhs)
+{
+    if(pmat==nullptr)
+    {
+        MFEM_ABORT("StokesSolve::Adjoint - The forward solver should be called first!!!")
+    }
+
+    //set BC
+    rhs=arhs;
+    for(int ii=0;ii<ess_tdofv.Size();ii++)
+    {
+        adj.GetBlock(0)[ess_tdofv[ii]]=0.0;
+        rhs.GetBlock(0)[ess_tdofv[ii]]=0.0;
+    }
+
+    psol->Mult(rhs, adj);
+}
+
+void StokesSolver::GradD(Vector &grad)
+{
+    if(dfes==nullptr)
+    {
+        MFEM_ABORT("StokesSolve::GradD - The design space in not set!!!")
+    }
+
+    //set vector size
+    grad.SetSize(dfes->GetTrueVSize());
+    grad=0.0;
+
+    fvelocity.SetFromTrueDofs(sol.GetBlock(0));
+    avelocity.SetFromTrueDofs(adj.GetBlock(0));
+
+    mfem::ParLinearForm lf(dfes);
+    lf.AddDomainIntegrator(new StokesGradIntergrator(fvelocity, avelocity
+                                                     ,*(ltopopt.fcoef),
+                                                     vfes->GetOrder(0)));
+    lf.ParallelAssemble(grad);
+}
+
+
+}

miniapps/autodiff/stokes.hpp

@@ -0,0 +1,813 @@
+#ifndef STOKES_HPP
+#define STOKES_HPP
+
+#include "mfem.hpp"
+
+
+namespace mfem {
+
+namespace PointwiseTrans
+{
+
+/*  Standrd "Heaviside" projection in topology optimization with threhold eta
+ * and steepness of the projection beta.
+ * */
+inline
+double HProject(double rho, double eta, double beta)
+{
+    // tanh projection - Wang&Lazarov&Sigmund2011
+    double a=std::tanh(eta*beta);
+    double b=std::tanh(beta*(1.0-eta));
+    double c=std::tanh(beta*(rho-eta));
+    double rez=(a+c)/(a+b);
+    return rez;
+}
+
+/// Gradient of the "Heaviside" projection with respect to rho.
+inline
+double HGrad(double rho, double eta, double beta)
+
+{
+    double c=std::tanh(beta*(rho-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;
+}
+
+/// Second derivative of the "Heaviside" projection with respect to rho.
+inline
+double HHess(double rho,double eta, double beta)
+{
+    double c=std::tanh(beta*(rho-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;
+}
+
+
+inline
+double FluidInterpolation(double rho,double q)
+{
+    return q*(1.0-rho)/(q+rho);
+}
+
+inline
+double GradFluidInterpolation(double rho, double q)
+{
+    double tt=q+rho;
+    return -q/tt-q*(1.0-rho)/(tt*tt);
+}
+
+
+}
+
+
+
+
+// Taylor-Hood finite elements
+class StokesIntegratorTH:public mfem::BlockNonlinearFormIntegrator
+{
+public:
+    StokesIntegratorTH()
+    {
+        mu=nullptr;
+        bc=nullptr;
+        ff=nullptr;
+
+        ss.SetSize(13);
+        rr.SetSize(13);
+        mm.SetSize(13);
+    }
+
+    StokesIntegratorTH(mfem::Coefficient* mu_, mfem::Coefficient* bc_, mfem::VectorCoefficient* ff_)
+    {
+        mu=mu_;
+        bc=bc_;
+        ff=ff_;
+
+        ss.SetSize(13);
+        rr.SetSize(13);
+        mm.SetSize(13);
+    }
+
+    virtual
+    double GetElementEnergy(const Array<const FiniteElement *> &el,
+                            ElementTransformation &Tr,
+                            const Array<const Vector *> &elfun)
+    {
+        return 0.0;
+    }
+
+    virtual
+    void AssembleElementVector(const Array<const FiniteElement *> &el,
+                               ElementTransformation &Tr,
+                               const Array<const Vector *> &elfun,
+                               const Array<Vector *> &elvec)
+    {
+        int dof_u = el[0]->GetDof();
+        int dof_p = el[1]->GetDof();
+
+        int dim =  el[0]->GetDim();
+
+        elvec[0]->SetSize(dim*dof_u);
+        elvec[1]->SetSize(dof_p);
+
+        int spaceDim = Tr.GetDimension();
+        if (dim != spaceDim)
+        {
+           mfem::mfem_error("StokesIntegrator::AssembleElementVector"
+                            " is not defined on manifold meshes");
+        }
+
+        // gradients
+        bsu.SetSize(dof_u,4);
+        bsp.SetSize(dof_p,1);
+
+        Vector uu(elfun[0]->GetData()+0*dof_u, dof_u);
+        Vector vv(elfun[0]->GetData()+1*dof_u, dof_u);
+
+        Vector ru(elvec[0]->GetData()+0*dof_u, dof_u); ru=0.0;
+        Vector rv(elvec[0]->GetData()+1*dof_u, dof_u); rv=0.0;
+
+        Vector ww;
+        Vector rw;
+        if(dim==2){
+            ww.SetSize(dof_u); ww=0.0;
+        }
+        else{
+            ww.SetDataAndSize(elfun[0]->GetData()+2*dof_u, dof_u);
+            rw.SetDataAndSize(elvec[0]->GetData()+2*dof_u, dof_u); rw=0.0;
+        }
+
+        Vector pp(elfun[1]->GetData(), dof_p);
+        Vector rp(elvec[1]->GetData(), dof_p); rp=0.0;
+
+        // temp storages for vectors and matrices
+        Vector sh;
+        DenseMatrix dh;
+
+        const IntegrationRule *ir = nullptr;
+        int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0]);
+        ir=&IntRules.Get(Tr.GetGeometryType(),order);
+
+        double bpenal; // Brinkmann penalization
+        double mmu;
+        Vector fv(3); fv=0.0;
+        double w;
+
+        for (int i = 0; i < ir->GetNPoints(); i++)
+        {
+           const IntegrationPoint &ip = ir->IntPoint(i);
+           Tr.SetIntPoint(&ip);
+           w=Tr.Weight();
+           w = ip.weight * w;
+
+           sh.SetDataAndSize(bsu.GetData(),dof_u);
+           el[0]->CalcPhysShape(Tr,sh);
+           sh.SetDataAndSize(bsp.GetData(),dof_p);
+           el[1]->CalcPhysShape(Tr,sh);
+
+           dh.UseExternalData(bsu.GetData()+dof_u,dof_u,dim);
+           el[0]->CalcPhysDShape(Tr,dh);
+           if(dim=2){
+               sh.SetDataAndSize(bsu.GetData()+3*dof_u,dof_u);
+               sh=0.0;
+           }
+
+           sh.SetDataAndSize(ss.GetData(),4);
+           bsu.MultTranspose(uu,sh);
+           sh.SetDataAndSize(ss.GetData()+4,4);
+           bsu.MultTranspose(vv,sh);
+           sh.SetDataAndSize(ss.GetData()+8,4);
+           if(dim==3){
+               bsu.MultTranspose(ww,sh);}
+           else{
+               sh=0.0;}
+           sh.SetDataAndSize(ss.GetData()+12,1);
+           bsp.MultTranspose(pp,sh);
+
+           mmu=1.0;
+           if(mu!=nullptr){
+               mmu=mu->Eval(Tr,ip);}
+
+           bpenal=0.0;
+           if(bc!=nullptr){
+               bpenal=bc->Eval(Tr,ip);}
+
+           fv=0.0;
+           if(ff!=nullptr){
+               ff->Eval(fv,Tr,ip);}
+
+           EvalQres(mmu,bpenal,fv[0],fv[1],fv[2],ss.GetData(),rr.GetData());
+
+           sh.SetDataAndSize(rr.GetData(),4);
+           bsu.AddMult_a(w,sh,ru);
+           sh.SetDataAndSize(rr.GetData()+4,4);
+           bsu.AddMult_a(w,sh,rv);
+           if(dim==3){
+               sh.SetDataAndSize(rr.GetData()+8,4);
+               bsu.AddMult_a(w,sh,rw);
+           }
+           sh.SetDataAndSize(rr.GetData()+12,1);
+           bsp.AddMult_a(w,sh,rp);
+        }
+
+    }
+
+    virtual
+    void AssembleElementGrad(const Array<const FiniteElement *> &el,
+                             ElementTransformation &Tr,
+                             const Array<const Vector *> &elfun,
+                             const Array2D<DenseMatrix *> &elmats)
+    {
+        int dof_u = el[0]->GetDof();
+        int dof_p = el[1]->GetDof();
+
+        int dim =  el[0]->GetDim();
+
+        int spaceDim = Tr.GetDimension();
+        if (dim != spaceDim)
+        {
+           mfem::mfem_error("StokesIntegrator::AssembleElementVector"
+                            " is not defined on manifold meshes");
+        }
+
+        elmats(0,0)->SetSize(dof_u*dim);
+        elmats(0,1)->SetSize(dof_u*dim,dof_p);
+        elmats(1,0)->SetSize(dof_p,dim*dof_u);
+        elmats(1,1)->SetSize(dof_p,dof_p);
+
+        (*elmats(0,0))=0.0;
+        (*elmats(0,1))=0.0;
+        (*elmats(1,0))=0.0;
+        (*elmats(1,1))=0.0;
+
+        // gradients
+        DenseMatrix bsu(dof_u,4);
+        DenseMatrix bsp(dof_p,1);
+
+        mm.SetSize(13,13); // state matrix
+
+        // temp storages for vectors and matrices
+        Vector sh;
+        DenseMatrix mh;
+        DenseMatrix th;
+        DenseMatrix rh;
+        DenseMatrix dh;
+
+        const IntegrationRule *ir = nullptr;
+        int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0]);
+        ir=&IntRules.Get(Tr.GetGeometryType(),order);
+
+
+        double mmu;
+        double bpenal; // Brinkmann penalization
+        double w;
+
+        for (int i = 0; i < ir->GetNPoints(); i++)
+        {
+           const IntegrationPoint &ip = ir->IntPoint(i);
+           Tr.SetIntPoint(&ip);
+           w=Tr.Weight();
+
+           // Primal
+           sh.SetDataAndSize(bsu.GetData(),dof_u);
+           el[0]->CalcPhysShape(Tr,sh);
+           sh.SetDataAndSize(bsp.GetData(),dof_p);
+           el[1]->CalcPhysShape(Tr,sh);
+
+
+           // Gradients
+           dh.UseExternalData(bsu.GetData()+dof_u,dof_u,dim);
+           el[0]->CalcPhysDShape(Tr,dh);
+
+           if(dim=2){
+               sh.SetDataAndSize(bsu.GetData()+3*dof_u,dof_u);
+               sh=0.0;
+           }
+
+           mmu=1.0;
+           if(mu!=nullptr)
+           {
+               mmu=mu->Eval(Tr,ip);
+           }
+
+           bpenal=0.0;
+           if(bc!=nullptr){
+               bpenal=bc->Eval(Tr,ip);}
+
+           EvalQMat(mmu,bpenal,mm.GetData());
+
+           w = ip.weight * w;
+
+           th.SetSize(dof_u,4);
+           rh.SetSize(dof_u);
+           mh.SetSize(4,4);
+           for(int ii=0;ii<dim;ii++){
+           for(int jj=0;jj<dim;jj++){
+               mh.CopyMN(mm,4,4,ii*4,jj*4);
+               mh.Transpose();
+               MultABt(bsu,mh,th);
+               MultABt(th,bsu,rh);
+               elmats(0,0)->AddMatrix(w,rh,ii*dof_u,jj*dof_u);
+           }}
+
+           th.SetSize(dof_u,1);
+           rh.SetSize(dof_u,dof_p);
+           mh.SetSize(4,1);
+           for(int jj=0;jj<dim;jj++){
+               mh.CopyMN(mm,4,1,jj*4,12);
+               mh.Transpose();
+               MultABt(bsu,mh,th);
+               MultABt(th,bsp,rh);
+               elmats(0,1)->AddMatrix(w,rh,jj*dof_u,0);
+           }
+
+           th.SetSize(dof_p,1);
+           rh.SetSize(dof_p,dof_p);
+           mh.SetSize(1,1);
+           mh.CopyMN(mm,1,1,12,12);
+           mh.Transpose();
+           MultABt(bsp,mh,th);
+           MultABt(th,bsp,rh);
+           elmats(1,1)->AddMatrix(w,rh,0,0);
+
+        }
+
+        elmats(1,0)->CopyMNt(*elmats(0,1),0,0);
+    }
+
+
+
+private:
+    mfem::Coefficient* mu;
+    mfem::Coefficient* bc;
+    mfem::VectorCoefficient* ff;
+
+    mfem::DenseMatrix bsu;
+    mfem::DenseMatrix bsp;
+    DenseMatrix mm;
+    Vector rr;
+    Vector ss;
+
+    void EvalQres(double mmu, double bpenal,
+                  double fx, double fy, double fz,
+                  double* uu, double* rr)
+    {
+        double t7,t9,t16;
+        t7 = mmu*(uu[2]+uu[5]);
+        t9 = mmu*(uu[3]+uu[9]);
+        t16 = mmu*(uu[7]+uu[10]);
+        rr[0] = bpenal*uu[0]-fx;
+        rr[1] = 2.0*mmu*uu[1]-uu[12];
+        rr[2] = t7;
+        rr[3] = t9;
+        rr[4] = bpenal*uu[4]-fy;
+        rr[5] = t7;
+        rr[6] = 2.0*mmu*uu[6]-uu[12];
+        rr[7] = t16;
+        rr[8] = bpenal*uu[8]-fz;
+        rr[9] = t9;
+        rr[10] = t16;
+        rr[11] = 2.0*mmu*uu[11]-uu[12];
+        rr[12] = -uu[1]-uu[6]-uu[11];
+    }
+
+    void EvalQMat(double mmu, double bpenal, double* kmat)
+    {
+        for(int i=0;i<169;i++){
+            kmat[i]=0.0;
+        }
+        double t1 = 2.0*mmu;
+        kmat[0] = bpenal;
+        kmat[14] = t1;
+        kmat[25] = -1.0;
+        kmat[28] = mmu;
+        kmat[31] = mmu;
+        kmat[42] = mmu;
+        kmat[48] = mmu;
+        kmat[56] = bpenal;
+        kmat[67] = mmu;
+        kmat[70] = mmu;
+        kmat[84] = t1;
+        kmat[90] = -1.0;
+        kmat[98] = mmu;
+        kmat[101] = mmu;
+        kmat[112] = bpenal;
+        kmat[120] = mmu;
+        kmat[126] = mmu;
+        kmat[137] = mmu;
+        kmat[140] = mmu;
+        kmat[154] = t1;
+        kmat[155] = -1.0;
+        kmat[157] = -1.0;
+        kmat[162] = -1.0;
+        kmat[167] = -1.0;
+    }
+
+};
+
+
+class FluidInterpolationCoefficient: public mfem::Coefficient
+{
+public:
+    FluidInterpolationCoefficient()
+    {
+        eta=0.5;
+        beta=8.0;
+        q=1.0;
+        lambda=100;
+    }
+
+    FluidInterpolationCoefficient(double eta_, double beta_,
+                                  double q_, double lambda_, mfem::GridFunction& gf):gfc(&gf)
+    {
+        eta=eta_;
+        beta=beta_;
+        q=q_;
+        lambda=lambda_;
+    }
+
+    virtual
+    double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+    {
+        double rhop=PointwiseTrans::HProject(gfc.Eval(T,ip),eta,beta);
+        return lambda*PointwiseTrans::FluidInterpolation(rhop,q);
+
+    }
+
+    virtual
+    double GradEval(ElementTransformation &T, const IntegrationPoint &ip)
+    {
+        double rhoo=gfc.Eval(T,ip);
+        double rhop=PointwiseTrans::HProject(rhoo,eta,beta);
+        double g1=lambda*PointwiseTrans::GradFluidInterpolation(rhop,q);
+        double g2=PointwiseTrans::HGrad(rhoo,eta,beta);
+        return g1*g2;
+
+    }
+
+    void SetGridFunction(mfem::GridFunction* gf)
+    {
+        gfc.SetGridFunction(gf);
+    }
+
+    void SetGridFunction(mfem::GridFunction& gf)
+    {
+        gfc.SetGridFunction(&gf);
+    }
+
+
+private:
+    mfem::GridFunctionCoefficient gfc;
+    double eta;  // threshold level
+    double beta; // steepness of the projection
+    double q; // interpolation parameter
+    double lambda; // penalization
+
+};
+
+class StokesGradIntergrator: public mfem::LinearFormIntegrator
+{
+public:
+    StokesGradIntergrator(mfem::GridFunction& velocity,
+                          mfem::GridFunction& adjoint,
+                          mfem::FluidInterpolationCoefficient& bcoef,
+                          int intorder_)
+        :vel(velocity), adj(adjoint), brm(bcoef)
+    {
+        intorder=intorder_; // should be equal to the order of the velocity field
+    }
+
+    virtual
+    void AssembleRHSElementVect(const mfem::FiniteElement &el, mfem::ElementTransformation &Tr, mfem::Vector &elvect)
+    {
+         int dof = el.GetDof();
+         int dim = el.GetDim();
+         elvect.SetSize(dof);
+         elvect=0.0;
+
+         const mfem::IntegrationRule *ir = nullptr;
+         int order= 2 * intorder + Tr.OrderGrad(&el);
+         ir=&IntRules.Get(Tr.GetGeometryType(),order);
+
+         mfem::Vector sh(dof);
+         mfem::Vector vv(dim);
+         mfem::Vector aa(dim);
+
+         double gbpenal; //gradient of the Brinkmann penalization
+         double dp;
+         double w;
+
+
+         for (int i = 0; i < ir->GetNPoints(); i++)
+         {
+            const IntegrationPoint &ip = ir->IntPoint(i);
+            Tr.SetIntPoint(&ip);
+            w=Tr.Weight();
+            w = ip.weight * w;
+
+            gbpenal=brm.GradEval(Tr,ip);
+
+            vel.GetVectorValue(Tr,ip,vv);
+            adj.GetVectorValue(Tr,ip,aa);
+
+            el.CalcPhysShape(Tr,sh);
+
+            for(int i=0;i<dof;i++)
+            {
+                dp=vv*aa;
+                elvect[i]=elvect[i]-w*sh[i]*gbpenal*dp;
+            }
+
+         }
+    }
+
+private:
+    mfem::GridFunction& vel;
+    mfem::GridFunction& adj;
+    mfem::FluidInterpolationCoefficient& brm;
+    int intorder;
+
+};
+
+
+
+class StokesSolver
+{
+public:
+    StokesSolver(mfem::ParMesh* pmesh_, int vorder=2)
+    {
+        if(vorder<2){vorder=2;}
+        int porder=vorder-1;
+        pmesh=pmesh_;
+
+        int dim=pmesh->Dimension();
+
+        vfec=new H1_FECollection(vorder,dim);
+        pfec=new H1_FECollection(porder,dim);
+
+
+        vfes=new mfem::ParFiniteElementSpace(pmesh,vfec,dim, Ordering::byVDIM);
+        pfes=new mfem::ParFiniteElementSpace(pmesh,pfec);
+
+        sfes.Append(vfes);
+        sfes.Append(pfes);
+
+        bpenal=nullptr;
+        load=nullptr;
+        viscosity=nullptr;
+
+        mfem::Array<mfem::ParFiniteElementSpace*> pf;
+        pf.Append(vfes);
+        pf.Append(pfes);
+
+        nf=new mfem::ParBlockNonlinearForm(pf);
+        nfin=nullptr;
+
+        rhs.Update(nf->GetBlockTrueOffsets()); rhs=0.0;
+        sol.Update(nf->GetBlockTrueOffsets()); sol=0.0;
+        adj.Update(nf->GetBlockTrueOffsets()); adj=0.0;
+
+        fvelocity.SetSpace(vfes); fvelocity=0.0;
+        fpressure.SetSpace(pfes); fpressure=0.0;
+        avelocity.SetSpace(vfes); avelocity=0.0;
+
+
+        pmat=nullptr;
+        prec=nullptr;
+        psol=nullptr;
+
+        dfes=nullptr;
+        ltopopt.fcoef=nullptr;
+        SetDesignParameters();
+
+        SetSolver();
+
+    }
+
+    ~StokesSolver()
+    {
+        delete ltopopt.fcoef;
+
+        delete psol;
+        delete prec;
+        delete pmat;
+
+        delete nf;
+
+        delete vfes;
+        delete pfes;
+        delete vfec;
+        delete pfec;
+
+    }
+
+    void SetViscosity(mfem::Coefficient& coef_)
+    {
+        viscosity=&coef_;
+    }
+
+    // Set the penalization field to coef_.
+    void SetBrinkmanPenal(mfem::Coefficient& coef_)
+    {
+        bpenal=&coef_;
+    }
+
+    void SetVolForces(mfem::VectorCoefficient& load_)
+    {
+        load=&load_;
+    }
+
+    void SetSolver(double rtol=1e-8, double atol=1e-12,int miter=1000, int prt_level=1)
+    {
+        rel_tol=rtol;
+        abs_tol=atol;
+        max_iter=miter;
+        print_level=prt_level;
+    }
+
+
+    /// Solves the forward problem.
+    void FSolve();
+
+    /// Solves the adjoint with the provided rhs.
+    void ASolve(mfem::BlockVector& rhs);
+
+    /// Return adj*d(residual(sol))/d(design). The dimension
+    /// of the vector grad is the save as the dimension of the
+    /// true design vector.
+    void GradD(mfem::Vector& grad);
+
+    mfem::ParGridFunction& GetVelocity()
+    {
+        fvelocity.SetFromTrueDofs(sol.GetBlock(0));
+        return fvelocity;
+    }
+
+    mfem::ParGridFunction& GetPressure()
+    {
+        fpressure.SetFromTrueDofs(sol.GetBlock(1));
+        return fpressure;
+    }
+
+    void AddVelocityBC(int id, int dir, double val)
+    {
+        if(dir==0){
+            bcx[id]=mfem::ConstantCoefficient(val);
+            AddVelocityBC(id,dir,bcx[id]);
+        }
+        if(dir==1){
+            bcy[id]=mfem::ConstantCoefficient(val);
+            AddVelocityBC(id,dir,bcy[id]);
+
+        }
+        if(dir==2){
+            bcz[id]=mfem::ConstantCoefficient(val);
+            AddVelocityBC(id,dir,bcz[id]);
+        }
+        if(dir==4){
+            bcx[id]=mfem::ConstantCoefficient(val);
+            bcy[id]=mfem::ConstantCoefficient(val);
+            bcz[id]=mfem::ConstantCoefficient(val);
+            AddVelocityBC(id,0,bcx[id]);
+            AddVelocityBC(id,1,bcy[id]);
+            AddVelocityBC(id,2,bcz[id]);
+        }
+    }
+
+    void AddVelocityBC(int id, int dir, mfem::Coefficient& val)
+    {
+        if(dir==0){ bccx[id]=&val; }
+        if(dir==1){ bccy[id]=&val; }
+        if(dir==2){ bccz[id]=&val; }
+        if(dir==4){ bccx[id]=&val; bccy[id]=&val; bccz[id]=&val;}
+        if(pmesh->Dimension()==2)
+        {
+            bccz.clear();
+        }
+    }
+
+
+    mfem::BlockVector& GetSol(){return sol;}
+    mfem::BlockVector& GetAdj(){return adj;}
+
+    mfem::ParFiniteElementSpace* GetVelocityFES(){return vfes;}
+    mfem::ParFiniteElementSpace* GetPressureFES(){return pfes;}
+
+    ///Get state spaces
+    mfem::Array<mfem::ParFiniteElementSpace*>& GetStateFES()
+    {
+        return sfes;
+    }
+
+    void SetDesignSpace(mfem::ParFiniteElementSpace* dfes_)
+    {
+        dfes=dfes_;
+        density.SetSpace(dfes);
+        density=0.0;
+
+        if(ltopopt.fcoef!=nullptr)
+        {
+            delete ltopopt.fcoef;
+        }
+
+        ltopopt.fcoef=new FluidInterpolationCoefficient(ltopopt.eta,ltopopt.beta,
+                                                        ltopopt.q,ltopopt.lambda,density);
+        bpenal=ltopopt.fcoef;
+    }
+
+    // Sets the design field using true vector designvec
+    // and set the penalization field
+    void SetDesign(mfem::Vector& designvec)
+    {
+        density.SetFromTrueDofs(designvec);
+        if(ltopopt.fcoef!=nullptr)
+        {
+            delete ltopopt.fcoef;
+        }
+
+        ltopopt.fcoef=new FluidInterpolationCoefficient(ltopopt.eta,ltopopt.beta,
+                                                        ltopopt.q,ltopopt.lambda,density);
+        bpenal=ltopopt.fcoef;
+    }
+
+    void SetDesignParameters(double eta_=0.5, double beta_=8.0, double q_=1, double lambda_=100)
+    {
+        ltopopt.eta=eta_;
+        ltopopt.beta=beta_;
+        ltopopt.lambda=lambda_;
+        ltopopt.q=q_;
+    }
+
+
+private:
+    double mu;
+    double alpha;
+    mfem::Coefficient* viscosity;
+    mfem::Coefficient* bpenal;
+    mfem::VectorCoefficient* load;
+
+    mfem::ParMesh* pmesh;
+
+    mfem::ParFiniteElementSpace* vfes;
+    mfem::ParFiniteElementSpace* pfes;
+    mfem::FiniteElementCollection* vfec;
+    mfem::FiniteElementCollection* pfec;
+    mfem::Array<mfem::ParFiniteElementSpace*> sfes;
+
+    // boundary conditions
+    std::map<int, mfem::ConstantCoefficient> bcx;
+    std::map<int, mfem::ConstantCoefficient> bcy;
+    std::map<int, mfem::ConstantCoefficient> bcz;
+
+    std::map<int, mfem::Coefficient*> bccx;
+    std::map<int, mfem::Coefficient*> bccy;
+    std::map<int, mfem::Coefficient*> bccz;
+
+    mfem::Array<int> ess_tdofv;
+
+
+    mfem::BlockNonlinearFormIntegrator* nfin;
+    mfem::ParBlockNonlinearForm* nf;
+
+    mfem::BlockVector rhs;
+    mfem::BlockVector sol;
+    mfem::BlockVector adj;
+
+    mfem::ParGridFunction fvelocity;
+    mfem::ParGridFunction avelocity;
+    mfem::ParGridFunction fpressure;
+
+    /// The PETSc objects are allocated once the problem is
+    /// assembled. They are utilized in computing the adjoint
+    /// solutions.
+    mfem::PetscParMatrix* pmat;
+    mfem::PetscPreconditioner* prec;
+    mfem::PetscLinearSolver*   psol;
+
+    double abs_tol;
+    double rel_tol;
+    int print_level;
+    int max_iter;
+
+    mfem::ParFiniteElementSpace* dfes; //design space
+    mfem::ParGridFunction density;
+
+    struct{
+      double eta;
+      double beta;
+      double lambda;
+      double q;
+      mfem::FluidInterpolationCoefficient* fcoef;
+    } ltopopt;
+
+};
+
+
+}
+
+
+#endif

miniapps/autodiff/stokes_test.cpp

@@ -0,0 +1,179 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+#include "stokes.hpp"
+
+// example runs
+// mpirun -np 4 ./stokes -m ./ball2D.msh -petscopts ./stokes_fieldsplit
+// mpirun -np 4 ./stokes -m ./ball2D.msh -petscopts ./stokes_fieldsplit_01
+
+double bpenal(const mfem::Vector &x)
+{
+    double nx=(x[0]-1.5)*(x[0]-1.5)+(x[1]-0.5)*(x[1]-0.5);
+    if(std::sqrt(nx)<0.1){ return 1e6;}
+    return 0.0;
+}
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI.
+   int nprocs, myrank;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+
+
+
+   // Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int order = 1;
+   bool static_cond = false;
+   int ser_ref_levels = 1;
+   int par_ref_levels = 1;
+   double newton_rel_tol = 1e-7;
+   double newton_abs_tol = 1e-12;
+   int newton_iter = 10;
+   int print_level = 1;
+   bool visualization = false;
+
+   const char *petscrc_file = "stokes_fieldsplit";
+
+   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",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&visualization,
+                  "-vis",
+                  "--visualization",
+                  "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   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(&petscrc_file, "-petscopts", "--petscopts",
+                                 "PetscOptions file to use.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myrank == 0)
+      {
+         args.PrintUsage(std::cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+
+   if (myrank == 0)
+   {
+      args.PrintOptions(std::cout);
+   }
+
+   mfem::MFEMInitializePetsc(NULL,NULL,petscrc_file,NULL);
+
+   // Read the (serial) mesh from the given mesh file on all processors.  We
+   // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   // and volume meshes with the same code.
+   mfem::Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // Refine the serial mesh on all processors 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 10,000 elements.
+   {
+      int ref_levels =
+         (int)floor(log(100./mesh.GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   // 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(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   /*
+   {
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh.UniformRefinement();
+      }
+   }*/
+
+
+   mfem::StokesSolver* solver=new mfem::StokesSolver(&pmesh,2);
+
+   mfem::ConstantCoefficient viscosity(1);
+   solver->SetViscosity(viscosity);
+   solver->AddVelocityBC(2,0,1.0);
+   solver->AddVelocityBC(1,4,0.0);
+
+   mfem::Vector vload(2); vload(0)=1.0; vload(1)=0.0; //vload(2)=0.0;
+   mfem::VectorConstantCoefficient load(vload);
+   solver->SetVolForces(load);
+
+   mfem::FunctionCoefficient brink(bpenal);
+   solver->SetBrinkmanPenal(brink);
+
+   //solver->AddVelocityBC(3,4,0.0);
+   //solver->AddVelocityBC(1,4,0.0);
+
+   //mfem::ConstantCoefficient bcvelx(2.0);
+   //mfem::ConstantCoefficient bcvely(1.5);
+   //solver->AddVelocityBC(1,0,bcvelx);
+   //solver->AddVelocityBC(2,1,bcvely);
+
+   solver->FSolve();
+
+
+   {
+       mfem::ParGridFunction& veloc=solver->GetVelocity();
+       mfem::ParGridFunction& press=solver->GetPressure();
+
+       mfem::ParaViewDataCollection paraview_dc("Stokes", &pmesh);
+       paraview_dc.SetPrefixPath("ParaView");
+       paraview_dc.SetLevelsOfDetail(order);
+       paraview_dc.SetDataFormat(mfem::VTKFormat::BINARY);
+       paraview_dc.SetHighOrderOutput(true);
+       paraview_dc.SetCycle(0);
+       paraview_dc.SetTime(0.0);
+       paraview_dc.RegisterField("velocity",&veloc);
+       paraview_dc.RegisterField("pressure",&press);
+       paraview_dc.Save();
+   }
+
+
+
+   delete solver;
+
+   mfem::MFEMFinalizePetsc();
+   MPI_Finalize();
+   return 0;
+}
+
+

miniapps/autodiff/taddensemat.hpp

@@ -0,0 +1,513 @@
+// 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 "mfem.hpp"
+#include "tadvector.hpp"
+
+namespace mfem
+{
+/// Templated dense matrix data type.
+/** The main goal of the TADDenseMatrix class is to serve as a data
+  container for representing dense matrices in classes, methods, and
+  functions utilized with automatic differentiation (AD). The
+  functionality/interface is copied from the standard MFEM dense
+  matrix mfem::DenseMatrix. The basic idea is to utilize the templated
+  vector class in combination with AD during the development phase.
+  The AD parts can be replaced with optimized code once the initial
+  development of the application is complete.  The common interface
+  between TADDenseMatrix and DenseMatrix will ease the transition
+  from AD to hand-optimized code as it does not require a change
+  in the interface or the code structure. TADDenseMatrix is intended
+  to be utilized for dense serial matrices. The objects can be combined
+  with TADVector or standard Vector.*/
+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;
+      }
+   }
+   /// Copy constructor using standard DenseMatrix
+   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);
+   }
+}
+
+} // namespace mfem
+
+#endif

miniapps/autodiff/tadvector.hpp

@@ -0,0 +1,699 @@
+// 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 "mfem.hpp"
+
+#include <cmath>
+#include <iostream>
+#include <limits>
+#if defined(_MSC_VER) && (_MSC_VER < 1800)
+#include <float.h>
+#define isfinite _finite
+#endif
+
+namespace mfem
+{
+/// Templated vector data type.
+/** The main goal of the TADVector class is to serve as a data
+  container for representing vectors in classes, methods, and
+  functions utilized with automatic differentiation (AD). The
+  functionality/interface is copied from the standard MFEM dense
+  vector mfem::Vector.  The basic idea is to utilize the templated
+  vector class in combination with AD during the development phase.
+  The AD parts can be replaced with optimized code once the initial
+  development of the application is complete.  The common interface
+  between TADVector and Vector will ease the transition from AD to
+  hand-optimized code as it does not require a change in the
+  interface or the code structure. TADVector is intended to be
+  utilized for dense serial vectors. */
+template<typename dtype>
+class TADVector
+{
+protected:
+   dtype *data;
+   int size;
+   int capacity;
+
+public:
+   /// Default constructor for Vector. Sets size = 0 and data = NULL.
+   TADVector()
+   {
+      data = nullptr;
+      size = 0;
+      capacity = 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 dtype[s];
+         capacity = s;
+         for (int i = 0; i < s; i++)
+         {
+            data[i] = v[i];
+         }
+      }
+      else
+      {
+         size = 0;
+         capacity = 0;
+         data = nullptr;
+      }
+   }
+
+   TADVector(const Vector &v)
+   {
+      const int s = v.Size();
+      if (s > 0)
+      {
+         size = s;
+         capacity = s;
+         data = new dtype[s];
+         for (int i = 0; i < s; i++)
+         {
+            data[i] = v[i];
+         }
+      }
+      else
+      {
+         size = 0;
+         capacity = 0;
+         data = nullptr;
+      }
+   }
+
+   /// @brief Creates vector of size s.
+   /// @warning Entries are not initialized to zero!
+   explicit TADVector(int s)
+   {
+      if (s > 0)
+      {
+         size = s;
+         capacity = s;
+         data = new dtype[size];
+      }
+      else
+      {
+         size = 0;
+         capacity = 0;
+         data = nullptr;
+      }
+   }
+
+   /// 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)
+   {
+      if (capacity > 0)
+      {
+         delete[] data;
+         capacity = 0;
+      }
+      size = _size;
+      data = _data;
+   }
+
+   /// 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 <= capacity)
+      {
+         size = s;
+         return;
+      }
+
+      delete[] data;
+      data = new dtype[s];
+      size = s;
+      capacity = s;
+   }
+
+   /// Set the Vector data and size.
+   /** The Vector does not assume ownership of the new data. The new size is
+       @warning This method should be called only when OwnsData() is false.
+       @sa NewDataAndSize(). */
+   void SetDataAndSize(dtype *d, int s)
+   {
+      if (OwnsData())
+      {
+         delete[] data;
+         capacity = 0;
+      }
+      data = d;
+      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) { SetDataAndSize(d, s); }
+
+   /// Reset the Vector to be a reference to a sub-vector of @a base.
+   inline void MakeRef(TADVector<dtype> &base, int offset, int size_)
+   {
+      NewDataAndSize(base.GetData() + 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)
+   {
+      int tsiz = size;
+      NewDataAndSize(base.GetData() + offset, tsiz);
+   }
+
+   /// Destroy a vector
+   void Destroy()
+   {
+      size = 0;
+      capacity = 0;
+      delete[] data;
+   }
+
+   /// 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 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; }
+
+   /// Read the Vector data (host pointer) ownership flag.
+   inline bool OwnsData() const { return (capacity > 0); }
+
+   /// Changes the ownership of the data; after the call the Vector is empty
+   inline void StealData(dtype **p)
+   {
+      *p = data;
+      delete[] data;
+      size = 0;
+      capacity = 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 dtype &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 dtype &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 dtype &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;
+      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] = (dtype) 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;
+   }
+
+   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;
+   }
+
+   template<typename ivtype>
+   TADVector &operator-=(ivtype v)
+   {
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = data[i] - v;
+      }
+      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;
+   }
+
+   template<typename ivtype>
+   TADVector &operator+=(ivtype v)
+   {
+      for (int i = 0; i < size; i++)
+      {
+         data[i] = data[i] + v;
+      }
+      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<dtype> &other)
+   {
+      Swap(data, other.data);
+      Swap(size, other.size);
+      Swap(capacity, other.capacity);
+   }
+
+   /// 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 vtype2>
+   friend void add(const vtype1 &v1,
+                   dtype 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];
+      }
+   }
+
+   template<typename vtype1, typename vtype2>
+   friend void add(const dtype a,
+                   const vtype1 &x,
+                   const dtype b,
+                   const vtype2 &y,
+                   TADVector<dtype> &z)
+   {
+      MFEM_ASSERT(x.Size() == y.Size() && x.Size() == z.Size(),
+                  "incompatible Vectors!");
+
+      for (int i = 0; i < z.Size(); i++)
+      {
+         z[i] = a * x[i] + b * y[i];
+      }
+   }
+
+   template<typename vtype1, typename vtype2>
+   friend void add(const dtype a,
+                   const vtype1 &x,
+                   const vtype2 &y,
+                   TADVector<dtype> &z)
+   {
+      MFEM_ASSERT(x.Size() == y.Size() && x.Size() == z.Size(),
+                  "incompatible Vectors!");
+
+      for (int i = 0; i < z.Size(); i++)
+      {
+         z[i] = a * x[i] + y[i];
+      }
+   }
+
+   template<typename vtype1, typename vtype2>
+   friend void subtract(const vtype1 &x, const vtype2 &y, TADVector<dtype> &z)
+   {
+      MFEM_ASSERT(x.Size() == y.Size() && x.Size() == z.Size(),
+                  "incompatible Vectors!");
+      for (int i = 0; i < z.Size(); i++)
+      {
+         z[i] = x[i] - y[i];
+      }
+   }
+
+   template<typename ivtype, typename vtype1, typename vtype2>
+   friend void subtract(const ivtype a,
+                        const vtype1 &x,
+                        const vtype2 &y,
+                        TADVector<dtype> &z)
+   {
+      MFEM_ASSERT(x.Size() == y.Size() && x.Size() == z.Size(),
+                  "incompatible Vectors!");
+      for (int i = 0; i < z.Size(); i++)
+      {
+         z[i] = a * (x[i] - y[i]);
+      }
+   }
+
+   /// Destroys vector.
+   ~TADVector() { delete[] data; }
+
+   /// Prints vector to stream out.
+   void Print(std::ostream &out = mfem::out, int width = 8) const
+   {
+      if (!size)
+      {
+         return;
+      }
+      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 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