filter implementation and test

*AD integrator for pLaplacian
*Select between AD integrator and hond coded integrator

.gitignore

@@ -147,6 +147,11 @@ examples/hiop/ex9-mesh.*
 examples/hiop/ex9-init.*
 examples/hiop/ex9-final.*
 
+
+examples/ex71
+examples/ex71p
+examples/Example71*
+
 examples/pumi/refined.mesh
 examples/pumi/sol.gf
 examples/pumi/mesh.*

CMakeLists.txt

@@ -296,6 +296,18 @@ 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()
+
 # CUDA
 if (MFEM_USE_CUDA)
    set(CMAKE_CUDA_STANDARD 11)
@@ -357,7 +369,8 @@ endif()
 #    be before SuiteSparse.
 set(MFEM_TPLS MPI_CXX OPENMP BLAS LAPACK METIS HYPRE SuiteSparse SUNDIALS PETSC
     SLEPC MESQUITE SuperLUDist STRUMPACK AXOM CONDUIT Ginkgo GNUTLS GSLIB NETCDF
-    MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE ADIOS2)
+    ADEPT FADBADPP MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA 
+    UMPIRE ADIOS2)
 # Add all *_FOUND libraries in the variable TPL_LIBRARIES.
 set(TPL_LIBRARIES "")
 set(TPL_INCLUDE_DIRS "")

INSTALL

@@ -448,6 +448,19 @@ 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_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
@@ -635,6 +648,16 @@ The specific libraries and their options are:
   URL: https://github.com/LLNL/hiop
   Options: HIOP_OPT, HIOP_LIB.
 
+- 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
+
 - 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.5, untar it at
   the same level as MFEM and create a symbolic link: "ln -s gslib-1.0.5 gslib".
@@ -813,6 +836,9 @@ MFEM_USE_MPFR
 MFEM_USE_ZLIB
 MFEM_USE_PUMI
 MFEM_USE_HIOP
+MFEM_USE_ADEPT
+MFEM_USE_FADBADPP
+MFEM_USE_ADFORWARD
 MFEM_USE_CUDA
 MFEM_USE_OCCA
 MFEM_USE_CEED
@@ -867,6 +893,8 @@ The CMake build system adds auto-detection for the following packages/libraries:
  - POSIXCLOCKS
  - PUMI
  - HIOP
+ - ADEPT
+ - FADBAD++
  - OCCA
  - RAJA
  - UMPIRE

config/cmake/MFEMConfig.cmake.in

@@ -50,6 +50,9 @@ 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_ADFORWARD @MFEM_USE_ADFORWARD@)
 
 set(MFEM_CXX_COMPILER "@CMAKE_CXX_COMPILER@")
 set(MFEM_CXX_FLAGS "@CMAKE_CXX_FLAGS@")

config/cmake/config.hpp.in

@@ -156,4 +156,13 @@
 // library.
 #cmakedefine MFEM_USE_SIMMETRIX
 
+// use ADEPT library for AD
+#cmakedefine MFEM_USE_ADEPT
+
+// use FADBAD++ library for AD
+#cmakedefine MFEM_USE_FADBADPP
+
+// use forward mode for automatic differentiation
+#cmakedefine MFEM_USE_ADFORWARD
+
 #endif // MFEM_CONFIG_HEADER

config/cmake/modules/FindADEPT.cmake

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

config/cmake/modules/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/cmake/modules/MfemCmakeUtilities.cmake

@@ -733,7 +733,8 @@ function(mfem_export_mk_files)
       MFEM_USE_SUPERLU MFEM_USE_STRUMPACK MFEM_USE_GNUTLS
       MFEM_USE_GSLIB MFEM_USE_NETCDF MFEM_USE_PETSC MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE
       MFEM_USE_CONDUIT MFEM_USE_PUMI MFEM_USE_CUDA MFEM_USE_OCCA MFEM_USE_RAJA
-      MFEM_USE_UMPIRE MFEM_USE_SIMD MFEM_USE_ADIOS2)
+      MFEM_USE_UMPIRE MFEM_USE_SIMD MFEM_USE_ADEPT MFEM_USE_FADBADPP 
+      MFEM_USE_ADFORWARD MFEM_USE_ADIOS2)
   foreach(var ${CONFIG_MK_BOOL_VARS})
     if (${var})
       set(${var} YES)

config/config.hpp.in

@@ -163,4 +163,14 @@
 // library.
 // #define MFEM_USE_SIMMETRIX
 
+// use ADEPT library for AD
+// #define MFEM_USE_ADEPT
+
+// use FADBAD++ library for AD
+// #define MFEM_USE_FADBADPP
+
+// use forward mode for automatic differentiation
+// #define MFEM_USE_ADFORWARD
+
+
 #endif // MFEM_CONFIG_HEADER

config/config.mk.in

@@ -43,6 +43,9 @@ 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_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

@@ -52,6 +52,9 @@ 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_ADFORWARD "Enable forward mode for AD" OFF)
 
 set(MFEM_MPI_NP 4 CACHE STRING "Number of processes used for MPI tests")
 
@@ -190,6 +193,13 @@ 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 "")
+
+
 # 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

@@ -140,6 +140,9 @@ MFEM_USE_CEED          = NO
 MFEM_USE_UMPIRE        = NO
 MFEM_USE_SIMD          = NO
 MFEM_USE_ADIOS2        = NO
+MFEM_USE_ADEPT         = NO
+MFEM_USE_FADBADPP      = NO
+MFEM_USE_ADFORWARD     = NO
 
 # Compile and link options for zlib.
 ZLIB_DIR =
@@ -336,6 +339,16 @@ 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.
+
 # GSLIB library
 GSLIB_DIR = @MFEM_DIR@/../gslib/build
 GSLIB_OPT = -I$(GSLIB_DIR)/include

examples/CMakeLists.txt

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

examples/ex71.cpp

@@ -0,0 +1,348 @@
+//                       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 handcoded integrator.
+//           Selecting integrator=1 will utilize the AD integrator.
+//           The AD integrator can be modifief to use ADQFunctionTJ.
+//
+//           qint (the integrand) 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 "ex71.hpp"
+
+#undef MFEM_USE_SUITESPARSE
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options
+   const char *mesh_file = "../data/beam-tet.mesh";
+   int ser_ref_levels = 3;
+   int order = 1;
+   bool visualization = true;
+   double newton_rel_tol = 1e-4;
+   double newton_abs_tol = 1e-6;
+   int newton_iter = 500;
+   int print_level = 0;
+   double pp = 2.0;
+   int integrator=1; //use AD
+   mfem::StopWatch* timer=new mfem::StopWatch();
+
+   mfem::OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&newton_rel_tol, "-rel", "--relative-tolerance",
+                  "Relative tolerance for the Newton solve.");
+   args.AddOption(&newton_abs_tol, "-abs", "--absolute-tolerance",
+                  "Absolute tolerance for the Newton solve.");
+   args.AddOption(&newton_iter, "-it", "--newton-iterations",
+                  "Maximum iterations for the Newton solve.");
+   args.AddOption(&pp, "-pp", "--power-parameter",
+                  "Power parameter (>=2.0) for the p-Laplacian.");
+   args.AddOption((&print_level),"-prt","--print-level",
+                  "Print level.");
+   args.AddOption(&integrator, "-int","--integrator",
+                  "Integrator 0: standard; 1: AD;");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(std::cout);
+      return 1;
+   }
+   args.PrintOptions(std::cout);
+
+
+   // 2. Read the (serial) mesh from the given mesh file.
+   mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 3. Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter.
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 4. Define the power parameter for the p-Laplacian and all other
+   //    coefficients
+   mfem::ConstantCoefficient c_pp(pp);
+   mfem::ConstantCoefficient load(1.000000000);
+   mfem::ConstantCoefficient c_ee(0.000000001);
+
+   // 5. Define the finite element spaces for the solution
+   mfem::H1_FECollection fec(order,dim);
+   mfem::FiniteElementSpace fespace(mesh,&fec,1,mfem::Ordering::byVDIM);
+   int glob_size=fespace.GetTrueVSize();
+   std::cout << "Number of finite element unknowns: " << glob_size << std::endl;
+
+   // 6. Define the Dirichlet conditions
+   mfem::Array<int> ess_bdr(mesh->bdr_attributes.Max());
+   ess_bdr = 1;
+
+   // 7. Define the nonlinear form
+   mfem::NonlinearForm* nf=new mfem::NonlinearForm(&fespace);
+
+   // 8. Define the solution vector x
+   mfem::GridFunction x(&fespace);
+   x = 0.0;
+   mfem::Vector tv(fespace.GetTrueVSize());
+   mfem::Vector sv(fespace.GetTrueVSize());
+   tv=0.0;
+   sv=0.0;
+
+   // 9. Define ParaView DataCollection
+   mfem::ParaViewDataCollection *dacol=new
+   mfem::ParaViewDataCollection("Example71",
+                                mesh);
+   dacol->SetLevelsOfDetail(order);
+   dacol->RegisterField("sol",&x);
+
+   // 11. Set domain integrators - start with linear diffusion
+   {
+      // the default power coefficient is 2.0
+      mfem::ConstantCoefficient lpp(2.0);
+      if (integrator==0)
+      {
+         nf->AddDomainIntegrator(new mfem::pLaplace(lpp,c_ee,load));
+      }
+      else if (integrator==1)
+      {
+         nf->AddDomainIntegrator(new mfem::pLaplaceAD(lpp,c_ee,load));
+      }
+      nf->SetEssentialBC(ess_bdr);
+      // compute the energy
+      double energy=nf->GetEnergy(tv);
+      std::cout<<"[2] The total energy of the system is E="<<energy<<std::endl;
+      // time the assembly
+      timer->Clear();
+      timer->Start();
+      nf->GetGradient(sv);
+      timer->Stop();
+      std::cout<<"[2] The assembly time is: "<<timer->RealTime()<<std::endl;
+      mfem::Solver *prec;
+#ifdef MFEM_USE_SUITESPARSE
+      prec=new mfem::UMFPackSolver();
+#else
+      prec=new mfem::GSSmoother();
+#endif
+      mfem::CGSolver *j_pcg = new mfem::CGSolver();
+      j_pcg->SetRelTol(1e-7);
+      j_pcg->SetAbsTol(1e-15);
+      j_pcg->SetMaxIter(500);
+      j_pcg->SetPrintLevel(print_level);
+      j_pcg->SetPreconditioner(*prec);
+
+      mfem::NewtonSolver* ns;
+      ns=new mfem::NewtonSolver();
+      ns->iterative_mode = true;
+      ns->SetSolver(*j_pcg);
+      ns->SetOperator(*nf);
+      ns->SetPrintLevel(print_level);
+      ns->SetRelTol(1e-6);
+      ns->SetAbsTol(1e-12);
+      ns->SetMaxIter(10);
+      //solve the problem
+      timer->Clear();
+      timer->Start();
+      ns->Mult(tv, sv);
+      timer->Stop();
+      std::cout<<"Time for the NewtonSolver: "<<timer->RealTime()<<std::endl;
+
+      energy=nf->GetEnergy(sv);
+      std::cout<<"[pp=2] The total energy of the system is E="<<energy<<std::endl;
+
+      delete ns;
+      delete j_pcg;
+      delete prec;
+
+      x.SetFromTrueDofs(sv);
+      dacol->SetTime(2.0);
+      dacol->SetCycle(2);
+      dacol->Save();
+   }
+
+   // 12. Continue with powers higher than 2
+   for (int i=3; i<pp; i++)
+   {
+      delete nf;
+      nf=new mfem::NonlinearForm(&fespace);
+      mfem::ConstantCoefficient lpp((double)i);
+      if (integrator==0)
+      {
+         nf->AddDomainIntegrator(new mfem::pLaplace(lpp,c_ee,load));
+      }
+      else if (integrator==1)
+      {
+         nf->AddDomainIntegrator(new mfem::pLaplaceAD(lpp,c_ee,load));
+      }
+      nf->SetEssentialBC(ess_bdr);
+      // compute the energy
+      double energy=nf->GetEnergy(sv);
+      std::cout<<"[pp="<<i<<"] The total energy of the system is E="<<energy<<std::endl;
+      // time the assembly
+      timer->Clear();
+      timer->Start();
+      nf->GetGradient(sv);
+      timer->Stop();
+      std::cout<<"[pp="<<i<<"] The assembly time is: "<<timer->RealTime()<<std::endl;
+      mfem::Solver *prec;
+#ifdef MFEM_USE_SUITESPARSE
+      prec=new mfem::UMFPackSolver();
+#else
+      prec=new mfem::GSSmoother();
+#endif
+      mfem::CGSolver *j_pcg = new mfem::CGSolver();
+      j_pcg->SetRelTol(1e-7);
+      j_pcg->SetAbsTol(1e-15);
+      j_pcg->SetMaxIter(500);
+      j_pcg->SetPrintLevel(print_level);
+      j_pcg->SetPreconditioner(*prec);
+
+      mfem::NewtonSolver* ns;
+      ns=new mfem::NewtonSolver();
+      ns->iterative_mode = true;
+      ns->SetSolver(*j_pcg);
+      ns->SetOperator(*nf);
+      ns->SetPrintLevel(print_level);
+      ns->SetRelTol(1e-6);
+      ns->SetAbsTol(1e-12);
+      ns->SetMaxIter(10);
+      //solve the problem
+      timer->Clear();
+      timer->Start();
+      ns->Mult(tv, sv);
+      timer->Stop();
+      std::cout<<"Time for the NewtonSolver: "<<timer->RealTime()<<std::endl;
+
+      energy=nf->GetEnergy(sv);
+      std::cout<<"[pp="<<i<<"] The total energy of the system is E="<<energy<<std::endl;
+
+      delete ns;
+      delete j_pcg;
+      delete prec;
+
+      x.SetFromTrueDofs(sv);
+      dacol->SetTime(i);
+      dacol->SetCycle(i);
+      dacol->Save();
+   }
+
+   // 13. Continue with the final power
+   if ( std::abs(pp-2.0) > std::numeric_limits<double>::epsilon())
+   {
+      delete nf;
+      nf=new mfem::NonlinearForm(&fespace);
+      if (integrator==0)
+      {
+         nf->AddDomainIntegrator(new mfem::pLaplace(c_pp,c_ee,load));
+      }
+      else if (integrator==1)
+      {
+         nf->AddDomainIntegrator(new mfem::pLaplaceAD(c_pp,c_ee,load));
+      }
+      nf->SetEssentialBC(ess_bdr);
+      // compute the energy
+      double energy=nf->GetEnergy(sv);
+      std::cout<<"[pp="<<pp<<"] The total energy of the system is E="<<energy<<std::endl;
+      // time the assembly
+      timer->Clear();
+      timer->Start();
+      nf->GetGradient(sv);
+      timer->Stop();
+      std::cout<<"[pp="<<pp<<"] The assembly time is: "<<timer->RealTime()<<std::endl;
+      mfem::Solver *prec;
+#ifdef MFEM_USE_SUITESPARSE
+      prec=new mfem::UMFPackSolver();
+#else
+      prec=new mfem::GSSmoother();
+#endif
+      mfem::CGSolver *j_pcg = new mfem::CGSolver();
+      j_pcg->SetRelTol(1e-7);
+      j_pcg->SetAbsTol(1e-15);
+      j_pcg->SetMaxIter(500);
+      j_pcg->SetPrintLevel(print_level);
+      j_pcg->SetPreconditioner(*prec);
+
+      mfem::NewtonSolver* ns;
+      ns=new mfem::NewtonSolver();
+      ns->iterative_mode = true;
+      ns->SetSolver(*j_pcg);
+      ns->SetOperator(*nf);
+      ns->SetPrintLevel(print_level);
+      ns->SetRelTol(1e-6);
+      ns->SetAbsTol(1e-12);
+      ns->SetMaxIter(10);
+      //solve the problem
+      timer->Clear();
+      timer->Start();
+      ns->Mult(tv, sv);
+      timer->Stop();
+      std::cout<<"Time for the NewtonSolver: "<<timer->RealTime()<<std::endl;
+
+      energy=nf->GetEnergy(sv);
+      std::cout<<"[pp="<<pp<<"] The total energy of the system is E="<<energy<<std::endl;
+
+      delete ns;
+      delete j_pcg;
+      delete prec;
+
+      x.SetFromTrueDofs(sv);
+      dacol->SetTime(pp);
+      if (pp<2.0)
+      {
+         dacol->SetCycle(std::floor(pp));
+      }
+      else
+      {
+         dacol->SetCycle(std::ceil(pp));
+      }
+      dacol->Save();
+   }
+
+
+
+   // 19. Free the used memory
+   delete dacol;
+   delete nf;
+   delete mesh;
+   delete timer;
+
+
+   return 0;
+}
+
+

examples/ex71.hpp

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

examples/ex71p.cpp

@@ -0,0 +1,401 @@
+//                       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 handcoded integrator.
+//           Selecting integrator=1 will utilize the AD integrator.
+//           The AD integrator can be modifief to use ADQFunctionTJ.
+//
+//           qint (the integrand) 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 "ex71.hpp"
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI
+   int num_procs, myrank;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+
+   // 2. Parse command-line options
+   const char *mesh_file = "../data/beam-tet.mesh";
+   int ser_ref_levels = 0;
+   int par_ref_levels = 0;
+   int order = 2;
+   bool visualization = true;
+   double newton_rel_tol = 1e-4;
+   double newton_abs_tol = 1e-6;
+   int newton_iter = 500;
+   int print_level = 0;
+   double pp = 2.0;
+   int integrator=1; //use AD
+   mfem::StopWatch* timer=new mfem::StopWatch();
+
+   mfem::OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&newton_rel_tol, "-rel", "--relative-tolerance",
+                  "Relative tolerance for the Newton solve.");
+   args.AddOption(&newton_abs_tol, "-abs", "--absolute-tolerance",
+                  "Absolute tolerance for the Newton solve.");
+   args.AddOption(&newton_iter, "-it", "--newton-iterations",
+                  "Maximum iterations for the Newton solve.");
+   args.AddOption(&pp, "-pp", "--power-parameter",
+                  "Power parameter (>=2.0) for the p-Laplacian.");
+   args.AddOption((&print_level),"-prt","--print-level",
+                  "Print level.");
+   args.AddOption(&integrator, "-int","--integrator",
+                  "Integrator 0: standard; 1: AD");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myrank == 0)
+      {
+         args.PrintUsage(std::cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myrank == 0)
+   {
+      args.PrintOptions(std::cout);
+   }
+
+   // 3. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral and hexahedral meshes
+   //    with the same code.
+   mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter.
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   mfem::ParMesh *pmesh = new mfem::ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // 6. Define the power parameter for the p-Laplacian and all other
+   //    coefficients
+   mfem::ConstantCoefficient c_pp(pp);
+   mfem::ConstantCoefficient load(1.000000000);
+   mfem::ConstantCoefficient c_ee(0.000000001);
+
+   // 7. Define the finite element spaces for the solution
+   mfem::H1_FECollection fec(order,dim);
+   mfem::ParFiniteElementSpace fespace(pmesh,&fec,1,mfem::Ordering::byVDIM);
+   HYPRE_Int glob_size=fespace.GlobalTrueVSize();
+   if (myrank == 0)
+   {
+      std::cout << "Number of finite element unknowns: " << glob_size << std::endl;
+   }
+
+   // 8. Define the Dirichlet conditions
+   mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+   ess_bdr = 1;
+
+   // 9. Define the nonlinear form
+   mfem::ParNonlinearForm* nf=new mfem::ParNonlinearForm(&fespace);
+
+   // 10. Define the solution vector x as a parallel finite element grid function
+   //     corresponding to fespace. Initialize x with initial guess of zero,
+   //     which satisfies the boundary conditions.
+   mfem::ParGridFunction x(&fespace);
+   x = 0.0;
+   mfem::HypreParVector* tv=x.GetTrueDofs();
+   mfem::HypreParVector* sv=x.GetTrueDofs();
+
+   // 11. Define ParaView DataCollection
+   mfem::ParaViewDataCollection *dacol=new
+   mfem::ParaViewDataCollection("Example71",
+                                pmesh);
+   dacol->SetLevelsOfDetail(order);
+   dacol->RegisterField("sol",&x);
+
+
+   // 11. Set domain integrators - start with linear diffusion
+   {
+      // the default power coefficient is 2.0
+      mfem::ConstantCoefficient lpp(2.0);
+      if (integrator==0)
+      {
+         nf->AddDomainIntegrator(new mfem::pLaplace(lpp,c_ee,load));
+      }
+      else if (integrator==1)
+      {
+         nf->AddDomainIntegrator(new mfem::pLaplaceAD(lpp,c_ee,load));
+      }
+      nf->SetEssentialBC(ess_bdr);
+      // compute the energy
+      double energy=nf->GetEnergy(*tv);
+      if (myrank==0)
+      {
+         std::cout<<"[2] The total energy of the system is E="<<energy<<std::endl;
+      }
+      // time the assembly
+      timer->Clear();
+      timer->Start();
+      nf->GetGradient(*sv);
+      timer->Stop();
+      if (myrank==0)
+      {
+         std::cout<<"[2] The assembly time is: "<<timer->RealTime()<<std::endl;
+      }
+      mfem::Solver *prec=new mfem::HypreBoomerAMG();
+      mfem::GMRESSolver *j_gmres = new mfem::GMRESSolver(MPI_COMM_WORLD);
+      j_gmres->SetRelTol(1e-7);
+      j_gmres->SetAbsTol(1e-15);
+      j_gmres->SetMaxIter(300);
+      j_gmres->SetPrintLevel(print_level);
+      j_gmres->SetPreconditioner(*prec);
+
+      mfem::NewtonSolver* ns;
+      ns=new mfem::NewtonSolver(MPI_COMM_WORLD);
+      ns->iterative_mode = true;
+      ns->SetSolver(*j_gmres);
+      ns->SetOperator(*nf);
+      ns->SetPrintLevel(print_level);
+      ns->SetRelTol(1e-6);
+      ns->SetAbsTol(1e-12);
+      ns->SetMaxIter(3);
+      //solve the problem
+      timer->Clear();
+      timer->Start();
+      ns->Mult(*tv, *sv);
+      timer->Stop();
+      if (myrank==0)
+      {
+         std::cout<<"Time for the NewtonSolver: "<<timer->RealTime()<<std::endl;
+      }
+
+      energy=nf->GetEnergy(*sv);
+      if (myrank==0)
+      {
+         std::cout<<"[pp=2] The total energy of the system is E="<<energy<<std::endl;
+      }
+
+      delete ns;
+      delete j_gmres;
+      delete prec;
+
+      x.SetFromTrueDofs(*sv);
+      dacol->SetTime(2.0);
+      dacol->SetCycle(2);
+      dacol->Save();
+   }
+
+   // 12. Continue with powers higher than 2
+   for (int i=3; i<pp; i++)
+   {
+      delete nf;
+      nf=new mfem::ParNonlinearForm(&fespace);
+      mfem::ConstantCoefficient lpp((double)i);
+      if (integrator==0)
+      {
+         nf->AddDomainIntegrator(new mfem::pLaplace(lpp,c_ee,load));
+      }
+      else if (integrator==1)
+      {
+         nf->AddDomainIntegrator(new mfem::pLaplaceAD(lpp,c_ee,load));
+      }
+      nf->SetEssentialBC(ess_bdr);
+      // compute the energy
+      double energy=nf->GetEnergy(*sv);
+      if (myrank==0)
+      {
+         std::cout<<"[pp="<<i<<"] The total energy of the system is E="<<energy<<std::endl;
+      }
+      // time the assembly
+      timer->Clear();
+      timer->Start();
+      nf->GetGradient(*sv);
+      timer->Stop();
+      if (myrank==0)
+      {
+         std::cout<<"[pp="<<i<<"] The assembly time is: "<<timer->RealTime()<<std::endl;
+      }
+      mfem::Solver *prec=new mfem::HypreBoomerAMG();
+      mfem::GMRESSolver *j_gmres = new mfem::GMRESSolver(MPI_COMM_WORLD);
+      j_gmres->SetRelTol(1e-7);
+      j_gmres->SetAbsTol(1e-15);
+      j_gmres->SetMaxIter(300);
+      j_gmres->SetPrintLevel(print_level);
+      j_gmres->SetPreconditioner(*prec);
+
+      mfem::NewtonSolver* ns;
+      ns=new mfem::NewtonSolver(MPI_COMM_WORLD);
+      ns->iterative_mode = true;
+      ns->SetSolver(*j_gmres);
+      ns->SetOperator(*nf);
+      ns->SetPrintLevel(print_level);
+      ns->SetRelTol(1e-6);
+      ns->SetAbsTol(1e-12);
+      ns->SetMaxIter(3);
+      //solve the problem
+      timer->Clear();
+      timer->Start();
+      ns->Mult(*tv, *sv);
+      timer->Stop();
+      if (myrank==0)
+      {
+         std::cout<<"Time for the NewtonSolver: "<<timer->RealTime()<<std::endl;
+      }
+
+      energy=nf->GetEnergy(*sv);
+      if (myrank==0)
+      {
+         std::cout<<"[pp="<<i<<"] The total energy of the system is E="<<energy<<std::endl;
+      }
+
+      delete ns;
+      delete j_gmres;
+      delete prec;
+
+      x.SetFromTrueDofs(*sv);
+      dacol->SetTime(i);
+      dacol->SetCycle(i);
+      dacol->Save();
+   }
+
+   // 13. Continue with the final power
+   if ( std::abs(pp-2.0) > std::numeric_limits<double>::epsilon())
+   {
+      delete nf;
+      nf=new mfem::ParNonlinearForm(&fespace);
+      if (integrator==0)
+      {
+         nf->AddDomainIntegrator(new mfem::pLaplace(c_pp,c_ee,load));
+      }
+      else if (integrator==1)
+      {
+         nf->AddDomainIntegrator(new mfem::pLaplaceAD(c_pp,c_ee,load));
+      }
+      nf->SetEssentialBC(ess_bdr);
+      // compute the energy
+      double energy=nf->GetEnergy(*sv);
+      if (myrank==0)
+      {
+         std::cout<<"[pp="<<pp<<"] The total energy of the system is E="<<energy<<std::endl;
+      }
+      // time the assembly
+      timer->Clear();
+      timer->Start();
+      nf->GetGradient(*sv);
+      timer->Stop();
+      if (myrank==0)
+      {
+         std::cout<<"[pp="<<pp<<"] The assembly time is: "<<timer->RealTime()<<std::endl;
+      }
+      mfem::Solver *prec=new mfem::HypreBoomerAMG();
+      mfem::GMRESSolver *j_gmres = new mfem::GMRESSolver(MPI_COMM_WORLD);
+      j_gmres->SetRelTol(1e-8);
+      j_gmres->SetAbsTol(1e-15);
+      j_gmres->SetMaxIter(300);
+      j_gmres->SetPrintLevel(print_level);
+      j_gmres->SetPreconditioner(*prec);
+
+      mfem::NewtonSolver* ns;
+      ns=new mfem::NewtonSolver(MPI_COMM_WORLD);
+      ns->iterative_mode = true;
+      ns->SetSolver(*j_gmres);
+      ns->SetOperator(*nf);
+      ns->SetPrintLevel(print_level);
+      ns->SetRelTol(1e-6);
+      ns->SetAbsTol(1e-12);
+      ns->SetMaxIter(3);
+      //solve the problem
+      timer->Clear();
+      timer->Start();
+      ns->Mult(*tv, *sv);
+      timer->Stop();
+      if (myrank==0)
+      {
+         std::cout<<"Time for the NewtonSolver: "<<timer->RealTime()<<std::endl;
+      }
+
+      energy=nf->GetEnergy(*sv);
+      if (myrank==0)
+      {
+         std::cout<<"[pp="<<pp<<"] The total energy of the system is E="<<energy<<std::endl;
+      }
+
+      delete ns;
+      delete j_gmres;
+      delete prec;
+
+      x.SetFromTrueDofs(*sv);
+      dacol->SetTime(pp);
+      if (pp<2.0)
+      {
+         dacol->SetCycle(std::floor(pp));
+      }
+      else
+      {
+         dacol->SetCycle(std::ceil(pp));
+      }
+      dacol->Save();
+   }
+
+
+
+   // 19. Free the used memory
+   delete dacol;
+   delete sv;
+   delete tv;
+   delete nf;
+   delete pmesh;
+   delete timer;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+

examples/makefile

@@ -22,10 +22,10 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_EXAMPLES = ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 ex17\
-  ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27
+  ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex71
 PAR_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p ex12p\
   ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p ex25p\
-  ex26p ex27p
+  ex26p ex27p ex71p
 
 ifeq ($(MFEM_USE_MPI),NO)
    EXAMPLES = $(SEQ_EXAMPLES)
@@ -146,3 +146,4 @@ clean-exec:
 	@rm -f ex21*.mesh ex21*.sol ex21p_*.*
 	@rm -f ex23.mesh ex23-*.gf
 	@rm -f ex25.mesh ex25-*.gf ex25p-*.*
+	@rm -rf Example71

fem/CMakeLists.txt

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

fem/adnonlininteg.hpp

@@ -0,0 +1,402 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+
+
+#ifndef MFEM_ADNONLININTEG
+#define MFEM_ADNONLININTEG
+
+#include "../config/config.hpp"
+#include "fe.hpp"
+#include "coefficient.hpp"
+#include "fespace.hpp"
+#include "nonlininteg.hpp"
+#include "../linalg/tadvector.hpp"
+#include "../linalg/taddensemat.hpp"
+#include "../linalg/fdual.hpp"
+
+#if defined MFEM_USE_ADEPT
+#include <adept.h>
+#elif defined  MFEM_USE_FADBADPP
+#include <fadiff.h>
+#include <badiff.h>
+#endif
+
+//define Forward AD mode
+//#define MFEM_USE_ADFORWARD
+
+namespace mfem
+{
+
+// m - dimension of the residual vector
+// the Jacobian will have dimensions [m,length(uu)]
+template<template <typename, typename> class CTD, int m>
+class ADQFunctionTJ
+{
+protected:
+#ifdef MFEM_USE_ADEPT
+   adept::Stack  m_stack;
+#endif
+
+public:
+#if defined MFEM_USE_ADEPT
+   typedef adept::adouble             ADFType;
+   typedef TADVector<ADFType>         ADFVector;
+   typedef TADDenseMatrix<ADFType>    ADFDenseMatrix;
+#elif defined MFEM_USE_FADBADPP
+#ifdef MFEM_USE_ADFORWARD
+   typedef fadbad::F<double>         ADFType;
+   typedef TADVector<ADFType>         ADFVector;
+   typedef TADDenseMatrix<ADFType>    ADFDenseMatrix;
+#else
+   typedef fadbad::B<double>         ADFType;
+   typedef TADVector<ADFType>         ADFVector;
+   typedef TADDenseMatrix<ADFType>    ADFDenseMatrix;
+#endif
+#else
+   typedef mfem::ad::FDual<double>    ADFType;
+   typedef TADVector<ADFType>         ADFVector;
+   typedef TADDenseMatrix<ADFType>    ADFDenseMatrix;
+#endif
+
+#ifdef MFEM_USE_ADEPT
+   ADQFunctionTJ():m_stack(false)  {}
+#else
+   ADQFunctionTJ() {}
+#endif
+
+   ~ADQFunctionTJ() {}
+
+   double QFunction(const mfem::Vector& vparam, mfem::Vector& uu)
+   {
+      CTD<double,mfem::Vector> func;
+      return func(vparam,uu);
+   }
+
+   void   QFunctionDU(const mfem::Vector& vparam, ADFVector& uu,
+                      ADFVector& rr)
+   {
+      CTD<ADFType,ADFVector> func;
+      func(vparam,uu,rr);
+   }
+
+   void  QFunctionAU(const Vector &vparam, mfem::Vector &uu,
+                     mfem::Vector &rr)
+   {
+      //the result is computed automaticaly 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
+
+   }
+
+   void QFunctionDU(const mfem::Vector& vparam, mfem::Vector& uu,
+                    mfem::Vector& rr)
+   {
+      CTD<double,mfem::Vector> func;
+      func(vparam,uu,rr);
+   }
+
+   void QFunctionDD(const mfem::Vector& vparam, mfem::Vector& uu,
+                    mfem::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
+template<template <typename, typename> class CTD>
+class ADQFunctionTH
+{
+public:
+#if defined MFEM_USE_FADBADPP
+   typedef fadbad::B<double>             ADFType;
+   typedef TADVector<ADFType>            ADFVector;
+   typedef TADDenseMatrix<ADFType>       ADFDenseMatrix;
+
+   typedef fadbad::B<fadbad::F<double>>    ADSType;
+   typedef TADVector<ADSType>              ADSVector;
+   typedef TADDenseMatrix<ADSType>         ADSDenseMatrix;
+#else
+   typedef mfem::ad::FDual<double>    ADFType;
+   typedef TADVector<ADFType>         ADFVector;
+   typedef TADDenseMatrix<ADFType>    ADFDenseMatrix;
+
+   typedef mfem::ad::FDual<ADFType>   ADSType;
+   typedef TADVector<ADSType>         ADSVector;
+   typedef TADDenseMatrix<ADSType>    ADSDenseMatrix;
+#endif
+
+   ADQFunctionTH() {}
+
+   ~ADQFunctionTH() {}
+
+   double QFunction(const mfem::Vector& vparam, mfem::Vector& uu)
+   {
+      CTD<double, mfem::Vector> tf;
+      return tf(vparam, uu);
+   }
+
+   ADFType QFunction(const mfem::Vector& vparam, ADFVector& uu)
+   {
+      CTD<ADFType,ADFVector> tf;
+      return tf(vparam, uu);
+   }
+
+   ADSType QFunction(const mfem::Vector &vparam, ADSVector& uu)
+   {
+      CTD<ADSType,ADSVector> tf;
+      return tf(vparam, uu);
+   }
+
+   void QFunctionDU(const mfem::Vector& vparam, mfem::Vector& uu,
+                    mfem::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
+   }
+
+   void QFunctionDD(const mfem::Vector& vparam, const mfem::Vector& uu,
+                    mfem::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
+   }
+
+
+};// end template ADFunctionTH
+
+
+
+}
+
+#endif
+

fem/fem.hpp

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

fem/nonlinearform.cpp

@@ -581,6 +581,13 @@ double BlockNonlinearForm::GetEnergyBlocked(const BlockVector &bx) const
          }
       }
 
+   //free the allocated memory
+   for (int i=0; i<fes.Size(); ++i)
+   {
+      delete el_x[i];
+      delete vdofs[i];
+   }
+
    if (fnfi.Size())
    {
       MFEM_ABORT("TODO: add energy contribution from interior face terms");

linalg/fdual.hpp

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

linalg/linalg.hpp

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

linalg/taddensemat.hpp

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

linalg/tadvector.hpp

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

makefile

@@ -261,7 +261,8 @@ endif
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
 MFEM_REQ_LIB_DEPS = SUPERLU METIS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
  SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP GSLIB\
- OCCA CEED RAJA UMPIRE
+ ADEPT FADBADPP OCCA CEED RAJA UMPIRE
+
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
 SLEPC_ERROR_MSG = $(if $(SLEPC_FOUND),,. SLEPC config not found: $(SLEPC_VARS))
 
@@ -325,7 +326,8 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_NETCDF MFEM_USE_PETSC MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_CONDUIT\
  MFEM_USE_PUMI 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_SOURCE_DIR MFEM_INSTALL_DIR
+ MFEM_USE_ADEPT MFEM_USE_FADBADPP MFEM_USE_ADFORWARD\
+ MFEM_USE_ADIOS2 MFEM_SOURCE_DIR MFEM_INSTALL_DIR 
 
 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -647,6 +649,9 @@ 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_ADFORWARD     = $(MFEM_USE_ADFORWARD))
 	$(info MFEM_CXX               = $(value MFEM_CXX))
 	$(info MFEM_HOST_CXX          = $(value MFEM_HOST_CXX))
 	$(info MFEM_CPPFLAGS          = $(value MFEM_CPPFLAGS))

miniapps/CMakeLists.txt

@@ -20,6 +20,7 @@ endif()
 add_subdirectory(common)
 add_subdirectory(electromagnetics)
 add_subdirectory(navier)
+add_subdirectory(mtop)
 add_subdirectory(meshing)
 add_subdirectory(performance)
 add_subdirectory(tools)

miniapps/mtop/CMakeLists.txt

@@ -0,0 +1,18 @@
+# Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+
+#PDEFilter
+add_mfem_miniapp(pdefilt
+    MAIN pdefilt.cpp
+    EXTRA_SOURCES pdenssolver.cpp pphyssolvers.cpp
+    EXTRA_HEADERS pdenssolver.hpp pphyssolvers.hpp
+    LIBRARIES mfem)

miniapps/mtop/pdefilt.cpp

@@ -0,0 +1,168 @@
+#include<mfem.hpp>
+#include <fstream>
+#include <iostream>
+#include <cmath>
+#include "pdenssolver.hpp"
+
+
+double DensFunc(const mfem::Vector& a){
+    double sca=(4.0*M_PI);
+    double rez=(std::sin(sca*a[0])*std::sin(sca*a[1])*std::sin(sca*a[2]));
+    if(rez>0.0){ rez=1.0;} else {rez=0.0;}
+    return rez;
+}
+
+
+int main(int argc, char *argv[])
+{
+    // 1. Initialize MPI.
+    int nprocs, myrank;
+    MPI_Init(&argc, &argv);
+    MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+
+    // 2. Parse command-line options.
+    const char *mesh_file = "";
+    int  element_order = 1;
+    int  input_order = 2;
+    bool static_cond = false;
+    bool visualization = true;
+    double len_scale=0.1;
+
+    mfem::OptionsParser args(argc, argv);
+    args.AddOption(&mesh_file, "-m", "--mesh",
+                          "Mesh file to use.");
+    args.AddOption(&element_order, "-o", "--order",
+                          "Finite element order (filtered field - polynomial degree) or -1 for"
+                          " isoparametric space.");
+    args.AddOption(&input_order, "-io", "--iorder",
+                          "Finite element order (input field - polynomial degree) or -1 for"
+                          " isoparametric space.");
+    args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                          "--no-static-condensation", "Enable static condensation.");
+    args.AddOption(&len_scale, "-ls","--lscale",
+                          "Length scale for the PDE filter.");
+    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                          "--no-visualization",
+                          "Enable or disable GLVis visualization.");
+    args.Parse();
+
+    if (!args.Good())
+    {
+        if (myrank == 0)
+        {
+           args.PrintUsage(std::cout);
+        }
+        MPI_Finalize();
+        return 1;
+    }
+
+    if (myrank == 0)
+    {
+        args.PrintOptions(std::cout);
+    }
+
+    //generate parallel mesh
+    mfem::ParMesh *pmesh;
+    if(strlen(mesh_file)==0)
+    {
+        //generate the mesh
+        int nx=10;
+        int ny=10;
+        int nz=10;
+
+        double sx=1.0;
+        double sy=1.0;
+        double sz=1.0;
+        //alternative
+        //mfem::Element::Type::HEXAHEDRON
+        mfem::Mesh *mesh=new mfem::Mesh(nx,ny,nz, mfem::Element::Type::TETRAHEDRON,sx,sy,sz);
+        int dim = mesh->Dimension();
+        {
+            int ref_levels = (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
+            for (int l = 0; l < ref_levels; l++)
+            {
+                mesh->UniformRefinement();
+            }
+        }
+        //create the parallel mesh
+        pmesh = new mfem::ParMesh(MPI_COMM_WORLD, *mesh);
+        delete mesh;
+    }
+    else{
+        int generate_edges=0;
+        int refine=1;
+        bool fix_orientation=true;
+        mfem::Mesh *mesh = new mfem::Mesh(mesh_file, generate_edges, refine, fix_orientation);
+        int dim = mesh->Dimension();
+        {
+            int ref_levels = (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
+            for (int l = 0; l < ref_levels; l++)
+            {
+                mesh->UniformRefinement();
+            }
+        }
+        //create the parallel mesh
+        pmesh = new mfem::ParMesh(MPI_COMM_WORLD, *mesh);
+        delete mesh;
+    }
+
+    mfem::FunctionCoefficient fco(DensFunc);
+
+
+    int dim = pmesh->Dimension();
+    mfem::FiniteElementCollection* ffec=new mfem::H1_FECollection(element_order ,dim);
+    mfem::FiniteElementCollection* ifec=new mfem::L2_FECollection(input_order, dim,mfem::BasisType::Positive);
+    mfem::ParFiniteElementSpace* ffs= new mfem::ParFiniteElementSpace(pmesh,ffec,1,mfem::Ordering::byNODES);
+    mfem::ParFiniteElementSpace* ifs= new mfem::ParFiniteElementSpace(pmesh,ifec,1,mfem::Ordering::byNODES);
+
+    mfem::ParaViewDataCollection *dacol=new mfem::ParaViewDataCollection("filt",pmesh);
+    dacol->SetLevelsOfDetail(2);
+
+    mfem::PDEFilter* filt=new mfem::PDEFilter(pmesh,ifs,ffs);
+    filt->SetLenScale(len_scale);
+
+    //define the grid functions
+    mfem::ParGridFunction* gfin=new mfem::ParGridFunction(ifs); //input field
+    mfem::ParGridFunction* gfft=new mfem::ParGridFunction(ffs); //filtered filed
+    //true-dof vectors
+    mfem::HypreParVector* vin=gfin->GetTrueDofs();
+    mfem::HypreParVector* vft=gfft->GetTrueDofs();
+    *vft=0.0;
+
+    gfin->ProjectCoefficient(fco);
+    filt->FFilter(fco,*vft);
+    gfft->SetFromTrueDofs(*vft);
+    filt->FFilter(fco,*vft);
+    gfft->SetFromTrueDofs(*vft);
+
+    dacol->RegisterField("inp",gfin);
+    dacol->RegisterField("flt",gfft);
+
+    dacol->SetTime(0.0);
+    dacol->SetCycle(0);
+    dacol->Save();
+
+    gfin->GetTrueDofs(*vin);
+    filt->FFilter(*vin,*vft);
+    gfft->SetFromTrueDofs(*vft);
+    dacol->SetTime(1.0);
+    dacol->SetCycle(1);
+    dacol->Save();
+
+    delete dacol;
+
+    delete vft;
+    delete vin;
+    delete gfft;
+    delete gfin;
+    delete filt;
+    delete ifs;
+    delete ffs;
+    delete ifec;
+    delete ffec;
+    delete pmesh;
+
+    MPI_Finalize();
+    return 0;
+}

miniapps/mtop/pdenssolver.cpp

@@ -0,0 +1,253 @@
+#include "pdenssolver.hpp"
+
+namespace mfem {
+
+PDEFilter::PDEFilter(mfem::ParMesh *mesh,
+                     mfem::ParFiniteElementSpace *pfin_,
+                     mfem::ParFiniteElementSpace *pfout_, double r)
+{
+    mfem_solver.mesh=mesh;
+    mfem_solver.pfin=pfin_;
+    mfem_solver.pfout=pfout_;
+
+    SetLenScale(r);
+
+    mfem_solver.a=nullptr;
+    mfem_solver.bl=nullptr;
+    mfem_solver.rl=nullptr;
+
+    mfem_solver.mc=new mfem::ConstantCoefficient(1.0);
+
+    mfem_solver.prec=nullptr;
+    mfem_solver.solv=nullptr;
+
+    mfem_solver.A=nullptr;
+
+    mfem_solver.gfin.SetSpace(mfem_solver.pfin);
+    mfem_solver.gfft.SetSpace(mfem_solver.pfout);
+
+    mfem_solver.B.SetSize(mfem_solver.pfout->GetTrueVSize());
+
+    realloc_required=true;
+}
+
+
+PDEFilter::~PDEFilter()
+{
+
+    if(mfem_solver.A)
+    {
+        delete mfem_solver.A;
+    }
+    if(mfem_solver.prec)
+    {
+        delete mfem_solver.prec;
+    }
+    if(mfem_solver.solv)
+    {
+        delete mfem_solver.solv;
+    }
+    if(mfem_solver.a)
+    {
+        delete mfem_solver.a;
+        delete mfem_solver.dc;
+    }
+    if(mfem_solver.bl)
+    {
+        delete mfem_solver.bl;
+    }
+    if(mfem_solver.rl)
+    {
+        delete mfem_solver.rl;
+    }
+
+    delete mfem_solver.mc;
+}
+
+void PDEFilter::ClearLenScale()
+{
+    mcmap.clear();
+    default_diffusion=0.0;
+    realloc_required=true;
+}
+
+void PDEFilter::SetDiffusion(double a)
+{//set directly the default diffusion parameter
+    default_diffusion=a;
+    realloc_required=true;
+}
+void PDEFilter::SetDiffusion(int mark, double a)
+{
+    mcmap[mark]=a;
+    realloc_required=true;
+}
+void PDEFilter::SetLenScale(double r)
+{//set the default length scale
+    default_diffusion=r*r/12.0;
+    realloc_required=true;
+}
+void PDEFilter::SetLenScale(int mark, double r)
+{//set length scale for region with a specified mark
+    mcmap[mark]=r*r/12.0;
+    realloc_required=true;
+}
+
+void PDEFilter::Allocate()
+{
+    if(mfem_solver.solv)
+    {
+        delete mfem_solver.solv;
+        mfem_solver.solv=nullptr;
+    }
+
+    if(mfem_solver.prec)
+    {
+        delete mfem_solver.prec;
+        mfem_solver.prec=nullptr;
+    }
+
+    if(mfem_solver.bl)
+    {
+        delete mfem_solver.bl;
+        mfem_solver.bl=nullptr;
+    }
+
+    if(mfem_solver.rl)
+    {
+        delete mfem_solver.rl;
+        mfem_solver.rl=nullptr;
+    }
+
+    if(mfem_solver.a)
+    {
+        delete mfem_solver.a;
+        delete mfem_solver.dc;
+    }
+
+    if(mfem_solver.A)
+    {
+        delete mfem_solver.A;
+    }
+
+    mfem_solver.a=new mfem::ParBilinearForm(mfem_solver.pfout);
+    //allocate the diffusion coefficicent
+    {
+         mfem::Vector vv(mfem_solver.mesh->attributes.Max());
+         vv=default_diffusion;
+         for(auto it=mcmap.begin();it!=mcmap.end();it++)
+         {
+             vv(it->first-1)=it->second;
+         }
+         mfem_solver.dc=new mfem::PWConstCoefficient(vv);
+    }
+
+    //add integrators
+    mfem_solver.a->AddDomainIntegrator(new mfem::DiffusionIntegrator(*mfem_solver.dc));
+    mfem_solver.a->AddDomainIntegrator(new mfem::MassIntegrator(*mfem_solver.mc));
+    mfem_solver.a->Assemble();
+    mfem_solver.a->Finalize();
+    mfem_solver.A=mfem_solver.a->ParallelAssemble();
+    realloc_required=false;
+}
+
+
+void PDEFilter::FFilter(mfem::Coefficient& in, mfem::Vector& out)
+{
+    if(realloc_required)
+    {
+       Allocate();
+    }
+
+    if(mfem_solver.bl==nullptr)
+    {
+        //allocate the linear form
+        int io=mfem_solver.pfin->GetOrder(0);
+        int fo=mfem_solver.pfout->GetOrder(0);
+        mfem_solver.bl=new mfem::ParLinearForm(mfem_solver.pfout);
+        mfem_solver.bl->AddDomainIntegrator(new mfem::DomainLFIntegrator(in,0,io+fo+1));
+    }else{
+        //change only the integrator
+        Array<LinearFormIntegrator*>* ints = mfem_solver.bl->GetDLFI();
+        delete (*ints)[0];
+        int io=mfem_solver.pfin->GetOrder(0);
+        int fo=mfem_solver.pfout->GetOrder(0);
+        (*ints)[0]=new mfem::DomainLFIntegrator(in,0,io+fo+1);
+    }
+    (*mfem_solver.bl)=0.0;
+    mfem_solver.bl->Assemble();
+    mfem_solver.bl->ParallelAssemble(mfem_solver.B);
+
+    //set the prec
+    if(mfem_solver.prec==nullptr)
+    {
+        mfem_solver.prec=new mfem::HypreBoomerAMG(*mfem_solver.A);
+    }
+    //set the solver
+    if(mfem_solver.solv==nullptr)
+    {
+        mfem_solver.solv=new mfem::HyprePCG(mfem_solver.mesh->GetComm());
+    }
+
+    mfem_solver.solv->SetOperator(*mfem_solver.A);
+    mfem_solver.solv->SetTol(1e-8);
+    mfem_solver.solv->SetMaxIter(500);
+    mfem_solver.solv->SetPrintLevel(2);
+    mfem_solver.solv->SetPreconditioner(*mfem_solver.prec);
+    mfem_solver.solv->Mult(mfem_solver.B, out);
+}
+
+void PDEFilter::FFilter(mfem::Vector &in, mfem::Vector &out)
+{
+    mfem_solver.gfin.SetFromTrueDofs(in);
+    mfem::GridFunctionCoefficient inco(&mfem_solver.gfin);
+    FFilter(inco,out);
+}
+
+void PDEFilter::RFilter(mfem::Vector &in, mfem::Vector &out)
+{
+    if(realloc_required)
+    {
+       Allocate();
+    }
+
+    //set the prec
+    if(mfem_solver.prec==nullptr)
+    {
+        mfem_solver.prec=new mfem::HypreBoomerAMG(*mfem_solver.A);
+    }
+    //set the solver
+    if(mfem_solver.solv==nullptr)
+    {
+        mfem_solver.solv=new mfem::HyprePCG(mfem_solver.mesh->GetComm());
+    }
+
+    mfem_solver.solv->SetOperator(*mfem_solver.A);
+    mfem_solver.solv->SetTol(1e-8);
+    mfem_solver.solv->SetMaxIter(500);
+    mfem_solver.solv->SetPrintLevel(2);
+    mfem_solver.solv->SetPreconditioner(*mfem_solver.prec);
+    mfem_solver.solv->Mult(in,mfem_solver.B);
+
+    mfem_solver.gfft.SetFromTrueDofs(mfem_solver.B);
+    mfem::GridFunctionCoefficient inco(&mfem_solver.gfft);
+    if(mfem_solver.rl==nullptr)
+    {
+        //allocate the linear form
+        int io=mfem_solver.pfin->GetOrder(0);
+        int fo=mfem_solver.pfout->GetOrder(0);
+        mfem_solver.rl=new mfem::ParLinearForm(mfem_solver.pfin);
+        mfem_solver.rl->AddDomainIntegrator(new mfem::DomainLFIntegrator(inco,0,io+fo+1));
+    }else{
+        //change only the integrator
+        Array<LinearFormIntegrator*>* ints = mfem_solver.bl->GetDLFI();
+        delete (*ints)[0];
+        int io=mfem_solver.pfin->GetOrder(0);
+        int fo=mfem_solver.pfout->GetOrder(0);
+        (*ints)[0]=new mfem::DomainLFIntegrator(inco,0,io+fo+1);
+    }
+    (*mfem_solver.rl)=0.0;
+    mfem_solver.rl->Assemble();
+    mfem_solver.rl->ParallelAssemble(out);
+}
+
+}

miniapps/mtop/pdenssolver.hpp

@@ -0,0 +1,98 @@
+#ifndef PDENSSOLVER_H
+#define PDENSSOLVER_H
+
+#include <mfem.hpp>
+#include <map>
+#include <vector>
+#include <tuple>
+
+namespace mfem {
+
+class ParFilter{
+public:
+    ParFilter(){}
+    virtual ~ParFilter(){}
+
+    //input and output must be a true-dof vector.
+    virtual void FFilter(mfem::Vector& in, mfem::Vector& out)=0;
+    //output must be a true-dof vector.
+    virtual void FFilter(mfem::Coefficient& in, mfem::Vector& out)=0;
+    //input and output must be a true-dof vector.
+    virtual void RFilter(mfem::Vector& in, mfem::Vector& out)=0;
+};
+
+
+class PDEFilter: public ParFilter{
+public:
+    //The input parameter r is the support radius of a cone filter.
+    //The diffusion parameter is obtained as r^2/((2*sqrt(3))^2).
+    //For details see:
+    //Lazarov, B. S. & Sigmund, O.
+    //Filters in topology optimization based on Helmholtz-type differential equations
+    //International Journal for Numerical Methods in Engineering, 2011, 86, 765-781
+    //int order is utilized for the RHS of the filter
+    PDEFilter(mfem::ParMesh* mesh,  mfem::ParFiniteElementSpace *pfin,  //input field
+                                    mfem::ParFiniteElementSpace *pfout, //filtered field
+                                    double r=0.0);
+
+
+    virtual ~PDEFilter();
+
+    //in -true-dof vector derived from pfin
+    //out -true-dof vector derived from pfout
+    virtual void FFilter(mfem::Vector& in, mfem::Vector& out);
+    virtual void FFilter(mfem::Coefficient& in, mfem::Vector& out);
+
+    //in  -gradients true-dof vector derived from pfout
+    //out -gradients true-dof vector derived from pfin
+    virtual void RFilter(mfem::Vector& in, mfem::Vector& out);
+
+    void ClearLenScale(); //clear all length scales,i.e., set them to zero
+    void SetDiffusion(double a); //set directly the default diffusion parameter
+    void SetDiffusion(int mark, double a);
+    void SetLenScale(double r); //set the default length scale
+    void SetLenScale(int mark, double r); //set length scale for region with a specified mark
+
+private:
+    //define coefficients
+    double default_diffusion;
+    std::map<int,double> mcmap; //<mark,diffusion coefficient>
+
+    struct{
+       mfem::ParMesh* mesh;
+       mfem::ParFiniteElementSpace *pfin;
+       mfem::ParFiniteElementSpace *pfout;
+
+       mfem::Coefficient* dc; //diffusion coefficient
+       mfem::Coefficient* mc; //mass coefficient
+
+       mfem::ParBilinearForm *a;
+       mfem::ParLinearForm *bl;
+       mfem::ParLinearForm *rl;
+
+       mfem::HypreParMatrix *A;//assembled matrix
+       mfem::Vector B;
+
+       mfem::ParGridFunction gfin;//input density field
+       mfem::ParGridFunction gfft;//filtered density field
+
+       mfem::HypreSolver *prec;
+       mfem::HyprePCG    *solv;
+    }mfem_solver;
+
+    bool realloc_required;
+
+    void Allocate();
+
+};
+
+
+
+
+}
+
+
+
+
+
+#endif

miniapps/mtop/pphyssolvers.hpp

@@ -0,0 +1,122 @@
+#ifndef PPHYSSOLVERS_H
+#define PPHYSSOLVERS_H
+
+#include <mfem.hpp>
+#include <map>
+#include <vector>
+#include <tuple>
+
+namespace mfem {
+
+//block form of the PPhysSolvers
+class BPPhysSolvers
+{
+protected:
+
+
+public:
+    virtual void UpdateDesign(mfem::BlockVector& desf)=0;
+    //solve for the the state field
+    //for non-linear problems the initial guess is solf
+    virtual void FSolve(mfem::BlockVector& solf)=0;
+    //solve the adjoint problem
+    //the method should be called always after
+    virtual void ASolve(const mfem::BlockVector& solf, const mfem::BlockVector& arhs,
+                        mfem::BlockVector& adjf)=0;
+    //return adjf'*dr/ddesign
+    virtual void GradD(const mfem::BlockVector& solf, const mfem::BlockVector& adjf,
+                      mfem::BlockVector& grad)=0;
+
+    const mfem::ParFiniteElementSpace* GetSFES(int k)=0; //return solver FES
+    const mfem::ParFiniteElementSpace* GetDFES(int k)=0; //return design FES
+
+
+};
+
+class PPhysSolvers
+{
+protected:
+
+
+public:
+    virtual void UpdateDesign(mfem::Vector& desf)=0;
+    //solve for the the state field
+    //for non-linear problems the initial guess is solf
+    virtual void FSolve(mfem::Vector& solf)=0;
+    //solve the adjoint problem
+    //the method should be called always after
+    virtual void ASolve(const mfem::Vector& solf, const mfem::Vector& arhs,
+                        mfem::Vector& adjf)=0;
+    //return adjf'*dr/ddesign
+    virtual void GradD(const mfem::Vector& solf, const mfem::Vector& adjf,
+                      mfem::Vector& grad)=0;
+
+    const mfem::ParFiniteElementSpace* GetSFES()=0; //return solver FES
+    const mfem::ParFiniteElementSpace* GetDFES()=0; //return design FES
+
+
+};
+
+
+class ElastSolver3D:public PPhysSolvers
+{
+  public:
+    ElastSolver(mfem::ParMesh* mesh,
+                mfem::ParFiniteElementSpace* desfes);
+
+    virtual ~ElastSolver() override;
+
+    virtual void UpdateDesign(mfem::Vector& desf) override;
+
+    virtual void FSolve(mfem::Vector& solf) override;
+
+    virtual void ASolve(const mfem::Vector& solf,
+                        const mfem::Vector& arhs,
+                        mfem::Vector& adjf) override;
+
+    virtual void GradD(const mfem::Vector& solf,
+                       const mfem::Vector& adjf,
+                       mfem::Vector& grad) override;
+
+    const mfem::ParFiniteElementSpace* GetSFES() override;
+    const mfem::ParFiniteElementSpace* GetDFES() override;
+
+    void SetOrder(int order){ mfem_solv.order=order; }
+
+    void SetMaterial(double lam, double mu) {
+        mfem_solv.lam=lam;
+        mfem_solv.mu=mu;
+    }
+    //solver BC
+
+
+    //solver loads and BC
+private:
+    //BC map <mark, dof, val>
+    std::vector< std::tuple<int,int,double> > bcmap;
+    //load map <mark, pressure>
+    std::map<int, double> lcmap;
+    struct{
+        int order; //order of the elements
+
+        //Lame parameters
+        double lam;
+        double mu;
+
+        mfem::Vector* pdesf;//pointer to the design field
+        mfem::ParFiniteElementSpace* solfes;
+        mfem::ParFiniteElementSpace* desfes;
+    } mfem_solv;
+
+};
+
+}
+
+
+
+
+
+
+
+
+#endif

tests/unit/CMakeLists.txt

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

tests/unit/linalg/test_fdual.cpp

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