minor changes

CHANGELOG

@@ -33,6 +33,11 @@ Linear and nonlinear solvers
 
 New and updated examples and miniapps
 -------------------------------------
+- Added a new example code, Example 34/34p, for solving the obstacle problem.
+
+- Added a new example code, Example 35/35p, for topology optimization with
+  MFEM.
+
 - The mesh-explorer miniapp can now save mesh files in the VisIt or ParaView
   formats using the corresponding DataCollection objects. See option 'D' in the
   main menu.

data/disk.mesh

@@ -0,0 +1,109 @@
+MFEM NURBS mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# SEGMENT     = 1
+# SQUARE      = 3
+# CUBE        = 5
+#
+
+dimension
+2
+
+elements
+5
+1 3 4 5 6 7
+1 3 0 1 5 4
+1 3 1 2 6 5
+1 3 3 7 6 2
+1 3 0 4 7 3
+
+boundary
+4
+1 1 0 1
+1 1 2 3
+1 1 1 2
+1 1 3 0
+
+edges
+12
+0 0 1
+0 4 5
+0 7 6
+0 3 2
+1 1 2
+1 5 6
+1 4 7
+1 0 3
+2 0 4
+2 1 5
+2 2 6
+2 3 7
+
+vertices
+8
+
+knotvectors
+3
+2 3 0 0 0 1 1 1
+2 3 0 0 0 1 1 1
+2 3 0 0 0 1 1 1
+
+weights
+1
+1
+1
+1
+1
+1
+1
+1
+0.70710678118655
+1
+1
+0.70710678118655
+0.70710678118655
+1
+1
+0.70710678118655
+1
+1
+1
+1
+1
+0.85355339059327
+0.85355339059327
+0.85355339059327
+0.85355339059327
+
+FiniteElementSpace
+FiniteElementCollection: NURBS2
+VDim: 2
+Ordering: 1
+
+-0.70710678118 -0.70710678118
+0.70710678118 -0.70710678118
+0.70710678118 0.70710678118
+-0.70710678118 0.70710678118
+-0.35355339059 -0.35355339059
+0.35355339059 -0.35355339059
+0.35355339059 0.35355339059
+-0.35355339059 0.35355339059
+0 -1.41421356236
+0 -0.35355339059
+0 0.35355339059
+0 1.41421356236
+1.41421356236 0
+0.35355339059 0
+-0.35355339059 0
+-1.41421356236 0
+-0.530330085885 -0.530330085885
+0.530330085885 -0.530330085885
+0.530330085885 0.530330085885
+-0.530330085885 0.530330085885
+0 0
+0 -0.883883476475
+0.883883476475 0
+0 0.883883476475
+-0.883883476475 0

examples/CMakeLists.txt

@@ -40,6 +40,9 @@ list(APPEND ALL_EXE_SRCS
   ex30.cpp
   ex31.cpp
   ex33.cpp
+  ex34.cpp
+  ex35.cpp
+  efem_thermal.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -77,6 +80,8 @@ if (MFEM_USE_MPI)
     ex31p.cpp
     ex32p.cpp
     ex33p.cpp
+    ex34p.cpp
+    ex35p.cpp
     )
 endif()
 

examples/efem.hpp

@@ -0,0 +1,636 @@
+#ifndef MFEM_EFEM_HPP
+#define MFEM_EFEM_HPP
+
+#include "mfem.hpp"
+
+
+namespace mfem
+{
+/**
+ * @brief Inverse sigmoid, log(x/(1-x))
+ *
+ * @param x -
+ * @param tol tolerance to force x ∈ (tol, 1 - tol)
+ * @return double log(x/(1-x))
+ */
+double invsigmoid(const double x, const double tol=1e-12)
+{
+   // forcing x to be in (0, 1)
+   const double clipped_x = std::min(std::max(tol,x),1.0-tol);
+   return std::log(clipped_x/(1.0-clipped_x));
+}
+
+// Sigmoid function
+double sigmoid(const double x)
+{
+   return x >= 0 ? 1.0 / (1.0 + std::exp(-x)) : std::exp(x) / (1.0 + std::exp(x));
+}
+
+// Derivative of sigmoid function d(sigmoid)/dx
+double dsigdx(const double x)
+{
+   double tmp = sigmoid(-x);
+   return tmp - std::pow(tmp,2);
+}
+
+
+/**
+ * @brief A coefficient that maps u to f(u).
+ *
+ */
+class MappedGridFunctionCoefficient : public GridFunctionCoefficient
+{
+   // lambda function maps double to double
+   typedef std::__1::function<double(const double)> __LambdaFunction;
+private:
+   __LambdaFunction fun; // a lambda function f(u(x))
+protected:
+   std::string name = "NONE";
+public:
+   /**
+    * @brief Construct a mapped grid function coefficient with given gridfunction and lambda function
+    *
+    * @param[in] gf u
+    * @param[in] double_to_double lambda function, f(x)
+    * @param[in] comp (Optional) index of a vector if u is a vector
+    */
+   MappedGridFunctionCoefficient(const GridFunction *gf,
+                                 __LambdaFunction double_to_double, int comp = 1): GridFunctionCoefficient(gf,
+                                          comp), fun(double_to_double) {}
+
+   /// Evaluate the coefficient at @a ip.
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+   {
+      const double value = GridFunctionCoefficient::Eval(T, ip);
+      return fun(value);
+   }
+};
+
+/**
+ * @brief GridFunctionCoefficient that returns exp(u)
+ *
+ */
+class ExponentialGridFunctionCoefficient : public MappedGridFunctionCoefficient
+{
+public:
+   ExponentialGridFunctionCoefficient(const GridFunction *gf,
+                                      int comp=1):MappedGridFunctionCoefficient(gf, [](const double x) {return std::exp(x);},
+   comp) {name = "EXP";}
+};
+
+
+/**
+ * @brief GridFunctionCoefficient that returns log(u)
+ *
+ */
+class LogarithmicGridFunctionCoefficient : public MappedGridFunctionCoefficient
+{
+public:
+   LogarithmicGridFunctionCoefficient(const GridFunction *gf,
+                                      int comp=1):MappedGridFunctionCoefficient(gf, [](const double x) {return std::log(x);},
+   comp) {name = "LOG";}
+};
+/**
+ * @brief GridFunctionCoefficient that returns log(max(u, tolerance))
+ *
+ */
+class SafeLogarithmicGridFunctionCoefficient : public
+   MappedGridFunctionCoefficient
+{
+public:
+   SafeLogarithmicGridFunctionCoefficient(const GridFunction *gf,
+                                          const double tolerance,
+                                          int comp=1):MappedGridFunctionCoefficient(gf, [tolerance](
+                                                      const double x) {return std::log(std::max(x, tolerance));},
+   comp) {name = "SAFE LOG";}
+};
+
+
+/**
+ * @brief GridFunctionCoefficient that returns sigmoid(u)
+ *
+ */
+class SigmoidGridFunctionCoefficient : public MappedGridFunctionCoefficient
+{
+public:
+   SigmoidGridFunctionCoefficient(const GridFunction *gf,
+                                  int comp=1):MappedGridFunctionCoefficient(gf, [](const double x) {return sigmoid(x);},
+   comp) {name = "SIGMOID";}
+};
+
+
+/**
+ * @brief GridFunctionCoefficient that returns dsigdx(u) = sigmoid'(u)
+ *
+ */
+class DerSigmoidGridFunctionCoefficient : public MappedGridFunctionCoefficient
+{
+public:
+   DerSigmoidGridFunctionCoefficient(const GridFunction *gf,
+                                     int comp=1):MappedGridFunctionCoefficient(gf, [](const double x) {return dsigdx(x);},
+   comp) {name = "D(SIGMOID)/DX";}
+};
+
+
+/**
+ * @brief GridFunctionCoefficient that returns invsigmoid(u)
+ *
+ */
+class InvSigmoidGridFunctionCoefficient : public MappedGridFunctionCoefficient
+{
+public:
+   InvSigmoidGridFunctionCoefficient(const GridFunction *gf,
+                                     int comp=1):MappedGridFunctionCoefficient(gf, [](const double x) {return invsigmoid(x);},
+   comp) {name = "INVSIGMOID";}
+};
+
+
+/**
+ * @brief GridFunctionCoefficient that returns pow(u, exponent)
+ *
+ */
+class PowerGridFunctionCoefficient : public MappedGridFunctionCoefficient
+{
+public:
+   PowerGridFunctionCoefficient(const GridFunction *gf, int exponent, int comp=1)
+      : MappedGridFunctionCoefficient(gf, [exponent](double x) {return std::pow(x, exponent);},
+   comp) {name = "POWER";}
+};
+
+
+/**
+ * @brief GridFunctionCoefficient that returns u^2
+ *
+ */
+class SquaredGridFunctionCoefficient : public MappedGridFunctionCoefficient
+{
+public:
+   SquaredGridFunctionCoefficient(const GridFunction *gf, int exponent, int comp=1)
+      : MappedGridFunctionCoefficient(gf, [](const double x) {return x*x;},
+   comp) {name = "SQUARE";}
+};
+
+/**
+ * @brief SIMP Rule, r(ρ) = ρ_0 + (1-ρ_0)ρ^p
+ *
+ */
+class SIMPCoefficient : public MappedGridFunctionCoefficient
+{
+public:
+   /**
+    * @brief Make a GridFunctionCoefficient that computes r(ρ) = ρ_0 + (1-ρ_0)ρ^p
+    *
+    * @param gf Density, ρ
+    * @param exponent Exponent, p
+    * @param rho_min minimum density, ρ_0
+    */
+   SIMPCoefficient(const GridFunction *gf, const double exponent,
+                   const double rho_min=1e-12)
+      : MappedGridFunctionCoefficient(gf, [rho_min, exponent](const double x) {return rho_min + (1-rho_min)*std::pow(x, exponent);}) {name = "SIMP";}
+};
+
+/**
+ * @brief Derivative of SIMP Rule, r'(ρ) = p(1-ρ_0)ρ^(p-1). Used when computing RHS
+ *
+ */
+class SIMPDerCoefficient : public MappedGridFunctionCoefficient
+{
+public:
+   /**
+    * @brief Make a GridFunctionCoefficient that computes r'(ρ) = p(1-ρ_0)ρ^(p-1)
+    *
+    * @param gf Density, ρ
+    * @param exponent Exponent, p
+    * @param rho_min minimum density, ρ_0
+    */
+   SIMPDerCoefficient(const GridFunction *gf, const double exponent,
+                      const double rho_min=1e-12)
+      : MappedGridFunctionCoefficient(gf, [rho_min, exponent](const double x) {return exponent*(1-rho_min)*std::pow(x, exponent - 1.0);}) {name = "SIMPDER";}
+};
+
+/**
+ * @brief Projector Π : ψ → ψ + c so that ∫ ρ = θ|Ω| where ρ = sigmoid(ψ + c)
+ *
+ */
+class SigmoidDensityProjector
+{
+private:
+   FiniteElementSpace *fes;
+   Mesh *mesh;
+   const double target_volume;
+   SigmoidGridFunctionCoefficient *rho = nullptr; // ρ = sigmoid(ψ)
+   DerSigmoidGridFunctionCoefficient *dsigPsi = nullptr; // d(sigmoid(ψ))/dψ
+   LinearForm *intRho = nullptr; // ∫ ρ = ∫ sigmoid(ψ)
+   LinearForm *intDerSigPsi = nullptr; // ∫ d(sigmoid(ψ))/dψ
+   bool isParallel = false;
+
+public:
+   /**
+    * @brief Projector Π : ψ → ψ + c so that ∫ ρ = θ|Ω| where ρ = sigmoid(ψ + c)
+    *
+    * @param fespace Finite element space for ψ
+    * @param volume_fraction Volume fraction, θ
+    * @param volume Total volume of the domain, |Ω|
+    */
+   SigmoidDensityProjector(FiniteElementSpace *fespace,
+                           const double volume_fraction,
+                           const double volume)
+      :fes(fespace),
+       mesh(fespace->GetMesh()),
+       target_volume(volume_fraction*volume) {}
+
+   /**
+    * @brief Update ψ ↦ ψ + c so that ∫ ρ = θ |Ω|.
+    *
+    * Using Newton's method, find c such that
+    *
+    * ∫ sigmoid(ψ + c) = θ |Ω|
+    *
+    * @param psi ρ = sigmoid(ψ)
+    * @param max_iteration Maximum iteration for Newton iteration
+    * @param tolerance Newton update tolerance
+    */
+   double Apply(GridFunction &psi, const int max_iteration,
+                const double tolerance=1e-12)
+   {
+      // 0. Make or Update Helper objects
+
+      if (rho) // if helper objects are already created,
+      {
+         // update with the current GridFunction
+         rho->SetGridFunction(&psi);
+         dsigPsi->SetGridFunction(&psi);
+      }
+      else // if Apply is not called at all
+      {
+         // Create MappedGridFunctionCoefficients
+         rho = new SigmoidGridFunctionCoefficient(&psi);
+         dsigPsi = new DerSigmoidGridFunctionCoefficient(&psi);
+
+         // Create ∫ sigmoid(ψ) and ∫ sigmoid'(ψ)
+#ifdef MFEM_USE_MPI // if Using MPI,
+         // try convert it to parallel version
+         ParFiniteElementSpace * pfes = dynamic_cast<ParFiniteElementSpace *>(fes);
+         if (pfes)
+         {
+            isParallel = true;
+            // make parallel linear forms
+            intRho = new ParLinearForm(pfes);
+            intDerSigPsi = new ParLinearForm(pfes);
+         }
+         else
+         {
+            // make serial linear forms
+            intRho = new LinearForm(fes);
+            intDerSigPsi = new LinearForm(fes);
+         }
+#else
+         intRho = new LinearForm(fes);
+         intDerSigPsi = new LinearForm(fes);
+#endif
+         intRho->AddDomainIntegrator(new DomainLFIntegrator(*rho, 2, 0));
+         intDerSigPsi->AddDomainIntegrator(new DomainLFIntegrator(*dsigPsi, 2, 0));
+      }
+
+      // Newton Method
+      for (int i=0; i<max_iteration; i++)
+      {
+         // Compute ∫ sigmoid(ψ + c)
+         intRho->Assemble(); // necessary whenever ψ is updated
+         double f = intRho->Sum();
+         // Compute ∫ sigmoid'(ψ + c)
+         intDerSigPsi->Assemble();
+         double df = intDerSigPsi->Sum();
+
+#ifdef MFEM_USE_MPI
+         if (isParallel)
+         {
+            MPI_Allreduce(MPI_IN_PLACE, &f, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+            MPI_Allreduce(MPI_IN_PLACE, &df, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+         }
+#endif
+         f -= target_volume;
+
+
+         // Newton increment
+         const double dc = - f / df;
+         // Update ψ
+         psi += dc;
+         out << "Iteration: " << i << " (θ|Ω|, ∫ρ - θ|Ω|, Δc) = (" <<
+             target_volume << ", " <<
+             f << ", " << dc << ")" << std::endl;
+
+         if (abs(dc) < tolerance)
+         {
+            break;
+         }
+         MFEM_VERIFY(std::isfinite(dc), "Projection failed");
+      }
+      intRho->Assemble();
+      return intRho->Sum();
+   }
+};
+
+class EllipticSolver
+{
+private:
+   FiniteElementSpace *fes; // finite element space
+   BilinearForm *bilinForm; // main bilinear form
+   Array<int> ess_tdof_list; // essential boundary dof list
+   bool isParallel = false; // whether input fespace is parallel or not
+#ifdef MFEM_USE_MPI
+   ParFiniteElementSpace *pfes; // parallel
+   ParMesh *pmesh;
+#endif
+   bool pa; // partial assembly flag
+public:
+   EllipticSolver(FiniteElementSpace *fespace,
+                  BilinearForm *bilinearForm, Array<int> ess_bdr)
+      :fes(fespace),
+       bilinForm(bilinearForm)
+   {
+      fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      out << "Essential dof size: " << ess_tdof_list.Size() << std::endl;
+#ifdef MFEM_USE_MPI
+      {
+         pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
+         if (pfes) { isParallel = true; pmesh = pfes->GetParMesh();}
+      }
+#endif
+      pa = bilinForm->GetAssemblyLevel() == AssemblyLevel::PARTIAL;
+   }
+
+   void Solve(LinearForm *b, GridFunction *sol)
+   {
+      OperatorPtr A;
+      Vector B, X;
+      bilinForm->FormLinearSystem(ess_tdof_list, *sol, *b, A, X, B);
+
+      // 11. Solve the linear system A X = B.
+      CGSolver * cg = nullptr;
+      Solver * M = nullptr;
+#ifdef MFEM_USE_MPI
+      if (isParallel)
+      {
+         M = new HypreBoomerAMG;
+         dynamic_cast<HypreBoomerAMG*>(M)->SetPrintLevel(0);
+         cg = new CGSolver(pmesh->GetComm());
+      }
+      else
+      {
+         M = new GSSmoother((SparseMatrix&)(*A));
+         cg = new CGSolver;
+      }
+#else
+      M = new GSSmoother((SparseMatrix&)(*A));
+      cg = new CGSolver;
+#endif
+      cg->SetRelTol(1e-12);
+      cg->SetMaxIter(10000);
+      cg->SetPrintLevel(0);
+      cg->SetPreconditioner(*M);
+      cg->SetOperator(*A);
+      cg->Mult(B, X);
+      delete M;
+      delete cg;
+      bilinForm->RecoverFEMSolution(X, *b, *sol);
+   }
+};
+
+//  Class for solving Poisson's equation:
+//
+//       - ∇ ⋅(κ ∇ u) = f  in Ω
+//
+class DiffusionSolver
+{
+private:
+   Mesh * mesh = nullptr;
+   // diffusion coefficient
+   Coefficient * diffcf = nullptr;
+   // mass coefficient
+   Coefficient * masscf = nullptr;
+   Coefficient * rhscf = nullptr;
+   Coefficient * essbdr_cf = nullptr;
+   Coefficient * neumann_cf = nullptr;
+   VectorCoefficient * gradient_cf = nullptr;
+
+   // FEM solver
+   int dim;
+   FiniteElementCollection * fec = nullptr;
+   FiniteElementSpace * fes = nullptr;
+   Array<int> ess_bdr;
+   Array<int> neumann_bdr;
+   LinearForm * b = nullptr;
+   bool parallel = false;
+#ifdef MFEM_USE_MPI
+   ParMesh * pmesh = nullptr;
+   ParFiniteElementSpace * pfes = nullptr;
+#endif
+
+public:
+   DiffusionSolver() { }
+
+   void SetMesh(Mesh * mesh_)
+   {
+      mesh = mesh_;
+#ifdef MFEM_USE_MPI
+      pmesh = dynamic_cast<ParMesh *>(mesh);
+      if (pmesh) { parallel = true; }
+#endif
+   }
+   void SetDiffusionCoefficient(Coefficient * diffcf_) { diffcf = diffcf_; }
+   void SetMassCoefficient(Coefficient * masscf_) { masscf = masscf_; }
+   void SetRHSCoefficient(Coefficient * rhscf_) { rhscf = rhscf_; }
+   void SetEssentialBoundary(const Array<int> & ess_bdr_) { ess_bdr = ess_bdr_;};
+   void SetNeumannBoundary(const Array<int> & neumann_bdr_) { neumann_bdr = neumann_bdr_;};
+   void SetNeumannData(Coefficient * neumann_cf_) {neumann_cf = neumann_cf_;}
+   void SetEssBdrData(Coefficient * essbdr_cf_) {essbdr_cf = essbdr_cf_;}
+   void SetGradientData(VectorCoefficient * gradient_cf_) {gradient_cf = gradient_cf_;}
+   void SetFESpace(FiniteElementSpace * fespace)
+   {
+      fes = fespace;
+#ifdef MFEM_USE_MPI
+      pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
+      if (pmesh) {parallel = true;};
+#endif
+   }
+
+   void ResetFEM();
+   void SetupFEM();
+
+   void Solve(GridFunction *u);
+   LinearForm * GetLinearForm() {return b;}
+#ifdef MFEM_USE_MPI
+   ParLinearForm * GetParLinearForm()
+   {
+      if (parallel)
+      {
+         return dynamic_cast<ParLinearForm *>(b);
+      }
+      else
+      {
+         MFEM_ABORT("Wrong code path. Call GetLinearForm");
+         return nullptr;
+      }
+   }
+#endif
+
+   ~DiffusionSolver();
+
+};
+
+
+void DiffusionSolver::SetupFEM()
+{
+   dim = mesh->Dimension();
+
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      b = new ParLinearForm(pfes);
+   }
+   else
+   {
+      b = new LinearForm(fes);
+   }
+#else
+   fes = new FiniteElementSpace(mesh, fec);
+   b = new LinearForm(fes);
+#endif
+
+   if (!ess_bdr.Size())
+   {
+      if (mesh->bdr_attributes.Size())
+      {
+         ess_bdr.SetSize(mesh->bdr_attributes.Max());
+         ess_bdr = 1;
+      }
+   }
+}
+
+void DiffusionSolver::Solve(GridFunction *u)
+{
+   OperatorPtr A;
+   Vector B, X;
+   Array<int> ess_tdof_list;
+
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      pfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
+   }
+   else
+   {
+      fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
+   }
+#else
+   fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
+#endif
+   if (b)
+   {
+      delete b;
+#ifdef MFEM_USE_MPI
+      if (parallel)
+      {
+         b = new ParLinearForm(pfes);
+      }
+      else
+      {
+         b = new LinearForm(fes);
+      }
+#else
+      b = new LinearForm(fes);
+#endif
+   }
+   if (rhscf)
+   {
+      b->AddDomainIntegrator(new DomainLFIntegrator(*rhscf));
+   }
+   if (neumann_cf)
+   {
+      MFEM_VERIFY(neumann_bdr.Size(), "neumann_bdr attributes not provided");
+      b->AddBoundaryIntegrator(new BoundaryLFIntegrator(*neumann_cf),neumann_bdr);
+   }
+   else if (gradient_cf)
+   {
+      MFEM_VERIFY(neumann_bdr.Size(), "neumann_bdr attributes not provided");
+      b->AddBoundaryIntegrator(new BoundaryNormalLFIntegrator(*gradient_cf),
+                               neumann_bdr);
+   }
+
+   b->Assemble();
+
+   BilinearForm * a = nullptr;
+
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      a = new ParBilinearForm(pfes);
+   }
+   else
+   {
+      a = new BilinearForm(fes);
+   }
+#else
+   a = new BilinearForm(fes);
+#endif
+   a->AddDomainIntegrator(new DiffusionIntegrator(*diffcf));
+   if (masscf)
+   {
+      a->AddDomainIntegrator(new MassIntegrator(*masscf));
+   }
+   a->Assemble();
+   if (essbdr_cf)
+   {
+      u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
+   }
+   a->FormLinearSystem(ess_tdof_list, *u, *b, A, X, B, 1);
+
+   CGSolver * cg = nullptr;
+   Solver * M = nullptr;
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      M = new HypreBoomerAMG;
+      dynamic_cast<HypreBoomerAMG*>(M)->SetPrintLevel(0);
+      cg = new CGSolver(pmesh->GetComm());
+   }
+   else
+   {
+      M = new GSSmoother((SparseMatrix&)(*A));
+      cg = new CGSolver;
+   }
+#else
+   M = new GSSmoother((SparseMatrix&)(*A));
+   cg = new CGSolver;
+#endif
+   cg->SetRelTol(1e-12);
+   cg->SetMaxIter(10000);
+   cg->SetPrintLevel(0);
+   cg->SetPreconditioner(*M);
+   cg->SetOperator(*A);
+   cg->Mult(B, X);
+   delete M;
+   delete cg;
+   a->RecoverFEMSolution(X, *b, *u);
+   delete a;
+}
+
+void DiffusionSolver::ResetFEM()
+{
+   delete fes; fes = nullptr;
+   delete fec; fec = nullptr;
+   delete b;
+}
+
+
+DiffusionSolver::~DiffusionSolver()
+{
+   ResetFEM();
+}
+
+}
+
+
+#endif

examples/efem_thermal.cpp

@@ -0,0 +1,552 @@
+//                              MFEM Example 35
+//
+//
+// Compile with: make ex35
+//
+// Sample runs:
+//     ex35 -alpha 10
+//     ex35 -lambda 0.1 -mu 0.1
+//     ex35 -r 5 -o 2 -alpha 5.0 -epsilon 0.01 -mi 50 -mf 0.5 -tol 1e-5
+//     ex35 -r 6 -o 1 -alpha 10.0 -epsilon 0.01 -mi 50 -mf 0.5 -tol 1e-5
+//
+//
+// Description: This example code demonstrates the use of MFEM to solve a
+//              density-filtered [3] topology optimization problem. The
+//              objective is to minimize the thermal compliance
+//
+//                  minimize ∫_Ω f u dx over u ∈ H¹(Ω) and ρ ∈ L²(Ω)
+//
+//                  subject to
+//
+//                    -∇⋅(r(ρ̃)∇ u) = f       in Ω + BCs
+//                    -ϵ²Δρ̃ + ρ̃ = ρ          in Ω + Neumann BCs
+//                    0 ≤ ρ ≤ 1              in Ω
+//                    u ≤ 1                  in Ω
+//                    ∫_Ω ρ dx = θ vol(Ω)
+//
+//              Here, r(ρ̃) = ρ₀ + ρ̃³ (1-ρ₀) is the solid isotropic material
+//              penalization (SIMP) law, ϵ > 0 is the design length scale,
+//              and 0 < θ < 1 is the volume fraction. Note that we have
+//
+//              More specifically, we have f = 1 in an insulated rectagular
+//              domain Ω = (0, 1) x (0, 1) where the left middle section
+//              {x = 0} x (0.4, 0.6) is held at temperature 0.
+//
+//                                INSULATED
+//                       ---------------------------  1
+//                       |                         |
+//                       |                         |
+//                       * -                       |
+//                 u = 0 * |  0.2                  |
+//                       * -                       |
+//                       |                         |
+//                       |                         |
+//                       ---------------------------   0
+//                       0                         1
+//
+//              The problem is discretized and gradients are computing using
+//              finite elements [1]. The design is optimized using an entropic
+//              mirror descent algorithm introduced by Keith and Surowiec [2]
+//              that is tailored to the bound constraint 0 ≤ ρ ≤ 1.
+//
+//              This example highlights the ability of MFEM to deliver high-
+//              order solutions to inverse design problems and showcases how
+//              to set up and solve PDE-constrained optimization problems
+//              using the so-called reduced space approach.
+//
+//
+// [1] Andreassen, E., Clausen, A., Schevenels, M., Lazarov, B. S., & Sigmund, O.
+//    (2011). Efficient topology optimization in MATLAB using 88 lines of
+//    code. Structural and Multidisciplinary Optimization, 43(1), 1-16.
+// [2] Keith, B. and Surowiec, T. (2023) The entropic finite element method
+//     (in preparation).
+// [3] Lazarov, B. S., & Sigmund, O. (2011). Filters in topology optimization
+//     based on Helmholtz‐type differential equations. International Journal
+//     for Numerical Methods in Engineering, 86(6), 765-781.
+
+#include "mfem.hpp"
+#include <iostream>
+#include <fstream>
+#include "efem.hpp"
+
+using namespace std;
+using namespace mfem;
+/**
+ * ---------------------------------------------------------------
+ *                      ALGORITHM PREAMBLE
+ * ---------------------------------------------------------------
+ *
+ *  The Lagrangian for this problem is
+ *
+ *          L(u,ρ,ρ̃,w,w̃) = (f,u) + (r(ρ̃)∇u, ∇w) - (f,w) + ϵ^2(∇ρ̃, ∇w̃) + (ρ̃ - ρ, w̃)
+ *                       + α⁻¹D≤(u, uk) + α⁻¹(D≥(ρ, ρk) + D≤(ρ, ρk))
+ *
+ *  where
+ *
+ *    r(ρ̃) = ρ₀ + ρ̃³ (1 - ρ₀)            (SIMP rule)
+ *
+ *    D≥(x, y) = ∫ xlog(x/y) - (x - y)   (Lower Bound, away from 0)
+ *
+ *    D≤(x, y) = D≥(1 - x, 1 - y)        (Upper Bound, away from 1)
+ *
+ * ---------------------------------------------------------------
+ *
+ *  Discretization choices:
+ *
+ *     u ∈ Vh ⊂ H¹ (order p)
+ *     w ∈ Vh ⊂ H¹ (order p)
+ *     ρ̃ ∈ Vl ⊂ H¹ (order p)
+ *     w̃ ∈ Vl ⊂ H¹ (order p)
+ *     ψ ∈ Wl ⊂ L² (order p - 1)
+ *
+ *  where ρ = sigmoid(ψ) so that 0≤ρ≤1 is strongly enforced
+ *
+ * ---------------------------------------------------------------
+ *                          ALGORITHM
+ * ---------------------------------------------------------------
+ *
+ *  Update ψ with projected mirror descent via the following algorithm.
+ *
+ *  0. Initialize density field ψ = sigmoid⁻¹(θ) so that ∫ρ = ∫sigmoid(ψ) = θ|Ω|
+ *
+ *  While not converged:
+ *
+ *     1. Solve filter equation ∂_w̃ L = 0; i.e.,
+ *
+ *           (ϵ² ∇ ρ̃, ∇ v ) + (ρ̃,v) = (ρ,v)   ∀ v ∈ Vl.
+ *
+ *     2. Solve primal problem ∂_w L = 0; i.e.,
+ *
+ *                    (r(ρ̃) ∇u, ∇v) = (f,v)   ∀ v ∈ Vh.
+ *
+ *     3. Solve dual problem ∂_u L = 0; i.e.,
+ *
+ *                    (r(ρ̃) ∇w, ∇v) = (f,v) + α⁻¹(log(u/uk), v)    ∀ v ∈ Vh.
+ *
+ *        NOTE: Currently, log(u/uk) is not implemented here.
+ *
+ *        NOTE: When there is no constraint u≤1, then w = u.
+ *              In that case, we do not have to solve the dual problem.
+ *
+ *     4. Solve for filtered gradient ∂_ρ̃ L = 0; i.e.,
+ *
+ *      (ϵ² ∇ w̃ , ∇ v ) + (w̃ ,v) = ( r'(ρ̃) (∇ u ⋅ ∇ w), v)   ∀ v ∈ Vl.
+ *
+ *     5. Set intermediate variable ψ⋆ = ψ - α⁻¹ w̃.
+ *
+ *     6. Update ψ by ψ = proj(ψ⋆) = ψ⋆ + c where c is chosen to be
+ *
+ *                ∫ sigmoid(ψ⋆ + c) = θ|Ω|.
+ *
+ *  end
+ *
+ */
+
+
+/**
+ * @brief alpha*(log(max(a, tol)) - log(max(b, tol)))
+ *
+ */
+class SafeLogDiffGridFunctionCoefficient : public
+   SafeLogarithmicGridFunctionCoefficient
+{
+private:
+   SafeLogarithmicGridFunctionCoefficient
+   *gf_other; // gridfunction log(b) to be subtracted
+   double a = 1.0;
+
+public:
+
+   /**
+    * @brief log(max(a, tol)) - log(max(b, tol))
+    *
+    */
+   SafeLogDiffGridFunctionCoefficient(GridFunction *self_gf,
+                                      GridFunction *other_gf, const double tolerance):
+      SafeLogarithmicGridFunctionCoefficient(self_gf, tolerance),
+      gf_other(new SafeLogarithmicGridFunctionCoefficient(other_gf, tolerance)) {}
+
+   /// Evaluate the coefficient at @a ip.
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+   {
+      return a*(SafeLogarithmicGridFunctionCoefficient::Eval(T, ip)
+                - gf_other->Eval(T, ip));
+   }
+   void SetAlpha(const double alpha) { a = alpha; }
+};
+
+
+/**
+ * @brief sigmoid(u) - sigmoid(w), used for computing successive difference
+ *
+ */
+class SigmoidDiffGridFunctionCoefficient : public
+   SigmoidGridFunctionCoefficient
+{
+private:
+   SigmoidGridFunctionCoefficient
+   *gf_other; // gridfunction log(b) to be subtracted
+
+public:
+
+   /**
+    * @brief log(max(a, tol)) - log(max(b, tol))
+    *
+    */
+   SigmoidDiffGridFunctionCoefficient(GridFunction *a,
+                                      GridFunction *b):
+      SigmoidGridFunctionCoefficient(a),
+      gf_other(new SigmoidGridFunctionCoefficient(b)) {}
+
+   /// Evaluate the coefficient at @a ip.
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+   {
+      return SigmoidGridFunctionCoefficient::Eval(T, ip)
+             - gf_other->Eval(T, ip);
+   }
+};
+
+/**
+ * @brief -r'(ρ̃)(∇ u ⋅ ∇ w)
+ *
+ */
+class SIMPDerEnergyCoefficient : public GridFunctionCoefficient
+{
+private:
+   SIMPDerCoefficient *r_prime_rho;
+   GradientGridFunctionCoefficient *gradu;
+   GradientGridFunctionCoefficient *gradw;
+   Vector gradu_val, gradw_val;
+public:
+
+   SIMPDerEnergyCoefficient(GridFunction *rho_filter, const double exponent,
+                            const double rho_min, GridFunction *u,
+                            GridFunction *w):GridFunctionCoefficient()
+   {
+      r_prime_rho = new SIMPDerCoefficient(rho_filter, exponent, rho_min);
+      gradu = new GradientGridFunctionCoefficient(u);
+      gradw = new GradientGridFunctionCoefficient(w);
+      const int dim = u->FESpace()->GetMesh()->Dimension();
+      gradu_val = Vector(dim);
+      gradw_val = Vector(dim);
+   }
+
+   /// Evaluate the coefficient at @a ip.
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+   {
+      gradu->Eval(gradu_val, T, ip);
+      gradw->Eval(gradw_val, T, ip);
+      return -r_prime_rho->Eval(T, ip)*(gradu_val*gradw_val);
+   }
+};
+
+inline void clip(GridFunction &gf, const double lower, const double upper)
+{
+   for (auto &x : gf)
+   {
+      x = min(max(x, lower), upper);
+   }
+}
+
+int main(int argc, char *argv[])
+{
+
+   // 0 - 1. Parse command-line options.
+   int ref_levels = 4; // The number of initial mesh refinement
+   int order = 2; // Polynomial order p. State - p, Design - p - 1, Filter - p
+   bool visualization = true;
+   double alpha0 = 1.0; // Update rule
+   double epsilon = 0.01; // Design parameter, ϵ.
+   double mass_fraction = 0.3; // mass fraction, θ.
+   int max_it = 1e2; // projected mirror gradient maximum iteration
+   double tol = 1e-4; // Projected mirror gradient tolerance
+   double rho_min = 1e-6; // SIMP ρ0
+   double exponent = 3; // SIMP exponent
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&alpha0, "-alpha", "--alpha-step-length",
+                  "Step length for gradient descent.");
+   args.AddOption(&epsilon, "-epsilon", "--epsilon-thickness",
+                  "epsilon phase field thickness");
+   args.AddOption(&max_it, "-mi", "--max-it",
+                  "Maximum number of gradient descent iterations.");
+   args.AddOption(&tol, "-tol", "--tol",
+                  "Exit tolerance for ρ ");
+   args.AddOption(&mass_fraction, "-mf", "--mass-fraction",
+                  "Mass fraction for diffusion coefficient.");
+   args.AddOption(&rho_min, "-rmin", "--rho-min",
+                  "Minimum of density coefficient.");
+   args.AddOption(&exponent, "-exp", "--exponent",
+                  "SIMP exponent.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   Mesh mesh = Mesh::MakeCartesian2D(10,10,mfem::Element::Type::QUADRILATERAL,true,
+                                     1.0,1.0);
+
+   int dim = mesh.Dimension();
+
+   // 2. Set BCs.
+   for (int i = 0; i<mesh.GetNBE(); i++)
+   {
+      Element * be = mesh.GetBdrElement(i);
+      Array<int> vertices;
+      be->GetVertices(vertices);
+
+      double * coords1 = mesh.GetVertex(vertices[0]);
+      double * coords2 = mesh.GetVertex(vertices[1]);
+
+      Vector center(2);
+      center(0) = 0.5*(coords1[0] + coords2[0]);
+      center(1) = 0.5*(coords1[1] + coords2[1]);
+
+      // if (abs(center(1) - 0.5) < 0.1 && center(0) < 1e-12)
+      if (center(0) < 1e-12 // left
+          //  && std::abs(center(1) - 0.5) < 0.1 // middle
+         )
+      {
+         // the left center
+         be->SetAttribute(1);
+      }
+      else
+      {
+         // all other boundaries
+         be->SetAttribute(2);
+      }
+   }
+   mesh.SetAttributes();
+
+   // 3. Refine the mesh.
+   for (int lev = 0; lev < ref_levels; lev++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   const int filter_order = order - 1;
+
+
+   // 4. Define the necessary finite element spaces on the mesh.
+   H1_FECollection state_fec(order, dim); // FE collection for u
+   H1_FECollection filter_fec(order-1, dim); // FE collection for ρ̃
+   L2_FECollection control_fec(order-1, dim); // FE collection for ρ
+   FiniteElementSpace state_fes(&mesh, &state_fec); // Space for u
+   FiniteElementSpace filter_fes(&mesh, &filter_fec); // space for ρ̃
+   FiniteElementSpace control_fes(&mesh, &control_fec); // space for ρ
+
+   int state_size = state_fes.GetTrueVSize();
+   int control_size = control_fes.GetTrueVSize();
+   int filter_size = filter_fes.GetTrueVSize();
+   cout << "Number of state unknowns: " << state_size << endl;
+   cout << "Number of filter unknowns: " << filter_size << endl;
+   cout << "Number of control unknowns: " << control_size << endl;
+
+   // 5. Set the initial guess for ρ.
+   GridFunction u(&state_fes);
+   GridFunction w(&state_fes);
+   GridFunction psi(&control_fes);
+   GridFunction rho_filter(&filter_fes);
+   GridFunction w_filter(&filter_fes);
+
+   SigmoidGridFunctionCoefficient rho(&psi); // sigmoid(ρ)
+
+   u = 0.0;
+   w = 0.0;
+   w_filter = 0.0;
+   psi = invsigmoid(mass_fraction);
+   rho_filter.ProjectCoefficient(rho);
+
+   GridFunction u_old(u);
+   GridFunction psi_old(psi);
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient zero(0.0);
+
+   // 6. Set-up the physics solver.
+
+   // 6 - 1. State problem LHS
+   Array<int> ess_bdr_state(mesh.bdr_attributes.Max()); // Dirichlet at bdr == 1
+   // Only the first component is essential bdr
+   ess_bdr_state = 0;
+   ess_bdr_state[0] = 1;
+   // r(ρ̃) = ρ0 + (1-ρ0)ρ̃^p
+   SIMPCoefficient r_rho_filter(&rho_filter, exponent, rho_min);
+   // heat source
+   ConstantCoefficient f(1.0);
+   // (r(ρ̃)∇ u, ∇ v)
+   DiffusionSolver *state_solver = new DiffusionSolver();
+   state_solver->SetMesh(&mesh);
+   state_solver->SetFESpace(&state_fes);
+   state_solver->SetEssentialBoundary(ess_bdr_state);
+   state_solver->SetDiffusionCoefficient(&r_rho_filter);
+   state_solver->SetupFEM();
+
+   // 6 - 3. Filter problem LHS
+   Array<int> ess_bdr_filter(mesh.bdr_attributes.Max()); // Pure Neumann
+   ess_bdr_filter = 0;
+   // ϵ^2, filter diffusion coeff
+   ConstantCoefficient eps_squared(epsilon*epsilon);
+   // (ϵ∇ ρ̃, ∇ v) + (ρ̃, v)
+   DiffusionSolver *filter_solver = new DiffusionSolver();
+   filter_solver->SetMesh(&mesh);
+   filter_solver->SetFESpace(&filter_fes);
+   filter_solver->SetEssentialBoundary(ess_bdr_filter);
+   filter_solver->SetDiffusionCoefficient(&eps_squared);
+   filter_solver->SetMassCoefficient(&one);
+   filter_solver->SetupFEM();
+   SIMPDerEnergyCoefficient r_energy(&rho_filter, exponent, rho_min, &u, &w);
+
+   // 6 - 5. Prepare for Projection
+   LinearForm volForm(&control_fes);
+   volForm.AddDomainIntegrator(new DomainLFIntegrator(one, 0, 0));
+   volForm.Assemble();
+   const double vol = volForm.Sum(); // domain volume
+   out << "|Ω| = " << vol << std::endl;
+   SigmoidDensityProjector volProj(&control_fes, mass_fraction, vol);
+
+   // 6 - 6. M⁻¹: Vl -> Wl
+   BilinearForm invMass(&control_fes);
+   invMass.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator()));
+   invMass.Assemble();
+   GridFunctionCoefficient w_filter_cf(&w_filter);
+   LinearForm w_filter_load(&control_fes);
+   w_filter_load.AddDomainIntegrator(new DomainLFIntegrator(w_filter_cf));
+
+   // 10. Connect to GLVis. Prepare for VisIt output.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sout_u,sout_r,sout_rho;
+   if (visualization)
+   {
+      sout_u.open(vishost, visport);
+      sout_rho.open(vishost, visport);
+      sout_r.open(vishost, visport);
+      sout_u.precision(8);
+      sout_rho.precision(8);
+      sout_r.precision(8);
+
+      sout_u << "solution\n" << mesh << u;
+      sout_u << "view 0 0\n";  // view from top
+      sout_u << "keys jl********\n";  // turn off perspective and light
+      sout_u << "window_title 'Temperature u'";
+      sout_u.flush();
+
+      GridFunction rho_gf(&control_fes);
+      rho_gf.ProjectCoefficient(rho);
+
+      sout_rho << "solution\n" << mesh << rho_gf;
+      sout_rho << "view 0 0\n";  // view from top
+      sout_rho << "keys jl********\n";  // turn off perspective and light
+      sout_rho << "window_title 'Density ρ'";
+      sout_rho.flush();
+
+      sout_r << "solution\n" << mesh << rho_filter;
+      sout_r << "view 0 0\n";  // view from top
+      sout_r << "keys jl********\n";  // turn off perspective and light
+      sout_r << "window_title 'Filtered density ρ̃'";
+      sout_r.flush();
+   }
+
+   // mfem::ParaViewDataCollection paraview_dc("Elastic_compliance", &mesh);
+   // paraview_dc.SetPrefixPath("ParaView");
+   // paraview_dc.SetLevelsOfDetail(order);
+   // paraview_dc.SetCycle(0);
+   // paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   // paraview_dc.SetHighOrderOutput(true);
+   // paraview_dc.SetTime(0.0);
+   // paraview_dc.RegisterField("displacement",&u);
+   // paraview_dc.RegisterField("density",&rho);
+   // paraview_dc.RegisterField("filtered_density",&rho_filter);
+
+   // 11. Iterate
+   double c0 = 0.0;
+   SigmoidDiffGridFunctionCoefficient succ_err(&psi, &psi_old);
+   GridFunction zero_gf(&control_fes);
+   zero_gf = 0.0;
+   for (int k = 1; k < max_it; k++)
+   {
+      const double alpha = alpha0*k;
+
+      cout << "\nStep = " << k << endl;
+
+      // Step 1 - Filter solve
+      mfem::out << "(ϵ^2 ∇ ρ̃, ∇ v) + (ρ̃,v) = (ρ,v)" << std::endl;
+      filter_solver->SetRHSCoefficient(&rho);
+      filter_solver->Solve(&rho_filter);
+
+      // Step 2 - Primal solve
+      mfem::out << "(r(ρ̃) ∇ u, ∇ v) = (f, v)" << std::endl;
+      state_solver->SetRHSCoefficient(&f);
+      state_solver->Solve(&u);
+
+      // Step 3 - Dual solve
+      // @note w is actually -w as we do not negate the RHS.
+      mfem::out << "(r(ρ̃) ∇ w, ∇ v) = (f, v) + α⁻¹(log(u/uk), v)" <<
+                std::endl;
+      state_solver->SetRHSCoefficient(&f);
+      w = u;
+      state_solver->Solve(&w);
+
+      // Step 4 - Dual filter solve
+      // @note Because of Step 3, we also solving -w̃ instead of w̃.
+      mfem::out <<
+                "(ϵ^2 ∇ w̃, ∇ v) + (w̃, v) = (r'(ρ̃)(∇ u ⋅ ∇ w), v)" <<
+                std::endl;
+      filter_solver->SetRHSCoefficient(&r_energy);
+      filter_solver->Solve(&w_filter);
+
+      // Step 5 - Get ψ⋆ = ψ - α⁻¹ w̃
+      w_filter_load.Assemble();
+      psi_old = psi;
+      invMass.AddMult(w_filter_load, psi, 1/alpha);
+
+      // Step 6 - ψ = proj(ψ⋆)
+      // bound psi so that 0≈sigmoid(-100) < rho < sigmoid(100)≈1
+      // project
+      clip(psi, -100.0, 100.0);
+      const double currVol = volProj.Apply(psi, 20);
+
+      if (visualization)
+      {
+
+         GridFunction rho_gf(&state_fes); // use continuous fes for visualization
+         rho_gf.ProjectCoefficient(rho);
+         sout_rho << "solution\n" << mesh << rho_gf;
+         sout_rho.flush();
+
+         sout_r << "solution\n" << mesh << rho_filter;
+         sout_r.flush();
+
+         sout_u << "solution\n" << mesh << u;
+         sout_u.flush();
+
+         // paraview_dc.SetCycle(k);
+         // paraview_dc.SetTime((double)k);
+         // paraview_dc.Save();
+      }
+      const double norm_reduced_gradient = zero_gf.ComputeL2Error(succ_err);
+
+      mfem::out << "||ψ-ψk||: " << norm_reduced_gradient << std::endl;
+      mfem::out << "Volume Fraction: " << currVol / vol << std::endl;
+
+      if (norm_reduced_gradient < tol)
+      {
+         break;
+      }
+   }
+
+   return 0;
+}

examples/ex34.cpp

@@ -0,0 +1,490 @@
+//                                MFEM Example 34
+//
+//
+// Compile with: make ex34
+//
+// Sample runs: ex34p -o 2
+//              ex34p -o 2 -r 4
+//
+//
+// Description: This example code demonstrates the use of MFEM to solve the
+//              bound-constrained energy minimization problem
+//
+//                      minimize ||∇u||² subject to u ≥ ϕ in H¹₀.
+//
+//              This is known as the obstacle problem, and it is a simple
+//              mathematical model for contact mechanics.
+//
+//              In this example, the obstacle ϕ is a half-sphere centered
+//              at the origin of a circular domain Ω. After solving to a
+//              specified tolerance, the numerical solution is compared to
+//              a closed-form exact solution to assess accuracy.
+//
+//              The problem is discretized and solved using the entropic
+//              finite element method (EFEM) introduced by Keith and
+//              Surowiec [1].
+//
+//              This example highlights the ability of MFEM to deliver high-
+//              order solutions to variation inequality problems and
+//              showcases how to set up and solve nonlinear mixed methods.
+//
+//
+// [1] Keith, B. and Surowiec, T. (2023) The entropic finite element method
+//     (in preparation).
+
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+double spherical_obstacle(const Vector &pt);
+void spherical_obstacle_gradient(const Vector &pt, Vector &grad);
+double exact_solution_obstacle(const Vector &pt);
+void exact_solution_gradient_obstacle(const Vector &pt, Vector &grad);
+
+class LogarithmGridFunctionCoefficient : public Coefficient
+{
+protected:
+   GridFunction *u; // grid function
+   Coefficient *obstacle;
+   double min_val;
+
+public:
+   LogarithmGridFunctionCoefficient(GridFunction &u_, Coefficient &obst_,
+                                    double min_val_=-36)
+      : u(&u_), obstacle(&obst_), min_val(min_val_) { }
+
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+};
+
+class ExponentialGridFunctionCoefficient : public Coefficient
+{
+protected:
+   GridFunction *u; // grid function
+   Coefficient *obstacle;
+   double min_val;
+   double max_val;
+
+public:
+   ExponentialGridFunctionCoefficient(GridFunction &u_, Coefficient &obst_,
+                                      double min_val_=0.0, double max_val_=1e6)
+      : u(&u_), obstacle(&obst_), min_val(min_val_), max_val(max_val_) { }
+
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+};
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/disk.mesh";
+   int order = 1;
+   bool visualization = true;
+   int max_it = 10;
+   double tol = 1e-5;
+   int ref_levels = 3;
+   double alpha0 = 1.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  "isoparametric space.");
+   args.AddOption(&ref_levels, "-r", "--refs",
+                  "Number of h-refinements.");
+   args.AddOption(&max_it, "-mi", "--max-it",
+                  "Maximum number of iterations");
+   args.AddOption(&tol, "-tol", "--tol",
+                  "Stopping criteria based on the difference between"
+                  "successive solution updates");
+   args.AddOption(&alpha0, "-step", "--step",
+                  "Initial step size alpha");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the given mesh file.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 3. Refine the mesh to increase the resolution.
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   int curvature_order = max(order,2);
+   mesh.SetCurvature(curvature_order);
+   mesh.EnsureNCMesh();
+
+   // 4. Define the necessary finite element spaces on the mesh.
+   H1_FECollection H1fec(order, dim);
+   FiniteElementSpace H1fes(&mesh, &H1fec);
+
+   L2_FECollection L2fec(order-1, dim);
+   FiniteElementSpace L2fes(&mesh, &L2fec);
+
+   cout << "Number of finite element unknowns: "
+        << H1fes.GetTrueVSize()
+        << " "
+        << L2fes.GetTrueVSize() << endl;
+
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = H1fes.GetVSize();
+   offsets[2] = L2fes.GetVSize();
+   offsets.PartialSum();
+
+   BlockVector x(offsets), rhs(offsets);
+   x = 0.0; rhs = 0.0;
+
+   // 5. Determine the list of true (i.e., conforming) essential boundary dofs.
+   Array<int> empty;
+   Array<int> ess_tdof_list;
+   Array<int> ess_bdr(mesh.bdr_attributes.Max());
+   if (mesh.bdr_attributes.Size())
+   {
+      ess_bdr = 1;
+      H1fes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 6. Define an initial guess for the solution.
+   auto IC_func = [](const Vector &x)
+   {
+      double r0 = 1.0;
+      double rr = 0.0;
+      for (int i=0; i<x.Size(); i++)
+      {
+         rr += x(i)*x(i);
+      }
+      return r0*r0 - rr;
+   };
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient zero(0.0);
+
+   // 7. Define the solution vectors as a finite element grid functions
+   //    corresponding to the fespaces.
+   GridFunction u_gf, delta_psi_gf;
+   u_gf.MakeRef(&H1fes,x.GetBlock(0).GetData());
+   delta_psi_gf.MakeRef(&L2fes,x.GetBlock(1).GetData());
+   delta_psi_gf = 0.0;
+
+   GridFunction u_old_gf(&H1fes);
+   GridFunction psi_old_gf(&L2fes);
+   GridFunction psi_gf(&L2fes);
+   u_old_gf = 0.0;
+   psi_old_gf = 0.0;
+
+   // 8. Define the function coefficients for the solution and use them to
+   //    initialize the initial guess
+   FunctionCoefficient exact_coef(exact_solution_obstacle);
+   VectorFunctionCoefficient exact_grad_coef(dim,exact_solution_gradient_obstacle);
+   FunctionCoefficient IC_coef(IC_func);
+   ConstantCoefficient f(0.0);
+   FunctionCoefficient obstacle(spherical_obstacle);
+   u_gf.ProjectCoefficient(IC_coef);
+   u_old_gf = u_gf;
+
+   // 9. Initialize the slack variable ψₕ = exp(uₕ)
+   LogarithmGridFunctionCoefficient ln_u(u_gf, obstacle);
+   psi_gf.ProjectCoefficient(ln_u);
+   psi_old_gf = psi_gf;
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+
+   GridFunction u_alt_gf(&L2fes);
+   GridFunction error_gf(&L2fes);
+
+   ExponentialGridFunctionCoefficient exp_psi(psi_gf,obstacle);
+   u_alt_gf.ProjectCoefficient(exp_psi);
+
+   if (visualization)
+   {
+      sol_sock.open(vishost,visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << mesh << u_alt_gf <<
+               "window_title 'Discrete solution'" << flush;
+   }
+
+   // 10. Iterate
+   int k;
+   int total_iterations = 0;
+   double increment_u = 0.1;
+   for (k = 0; k < max_it; k++)
+   {
+      double alpha = alpha0 * (k+1);
+
+      GridFunction u_tmp(&H1fes);
+      u_tmp = u_old_gf;
+
+      mfem::out << "\nOUTER ITERATION " << k+1 << endl;
+
+      int j;
+      for ( j = 0; j < 15; j++)
+      {
+         total_iterations++;
+
+         ConstantCoefficient alpha_cf(alpha);
+
+         LinearForm b0,b1;
+         b0.Update(&H1fes,rhs.GetBlock(0),0);
+         b1.Update(&L2fes,rhs.GetBlock(1),0);
+
+         ExponentialGridFunctionCoefficient exp_psi(psi_gf, zero);
+         ProductCoefficient neg_exp_psi(-1.0,exp_psi);
+         GradientGridFunctionCoefficient grad_u_old(&u_old_gf);
+         ProductCoefficient alpha_f(alpha, f);
+         GridFunctionCoefficient psi_cf(&psi_gf);
+         GridFunctionCoefficient psi_old_cf(&psi_old_gf);
+         SumCoefficient psi_old_minus_psi(psi_old_cf, psi_cf, 1.0, -1.0);
+
+         b0.AddDomainIntegrator(new DomainLFIntegrator(alpha_f));
+         b0.AddDomainIntegrator(new DomainLFIntegrator(psi_old_minus_psi));
+         b0.Assemble();
+
+         b1.AddDomainIntegrator(new DomainLFIntegrator(exp_psi));
+         b1.AddDomainIntegrator(new DomainLFIntegrator(obstacle));
+         b1.Assemble();
+
+         BilinearForm a00(&H1fes);
+         a00.SetDiagonalPolicy(mfem::Operator::DIAG_ONE);
+         a00.AddDomainIntegrator(new DiffusionIntegrator(alpha_cf));
+         a00.Assemble();
+         a00.EliminateEssentialBC(ess_bdr,x.GetBlock(0),rhs.GetBlock(0),
+                                  mfem::Operator::DIAG_ONE);
+         a00.Finalize();
+         SparseMatrix &A00 = a00.SpMat();
+
+         MixedBilinearForm a10(&H1fes,&L2fes);
+         a10.AddDomainIntegrator(new MixedScalarMassIntegrator());
+         a10.Assemble();
+         a10.EliminateTrialDofs(ess_bdr, x.GetBlock(0), rhs.GetBlock(1));
+         a10.Finalize();
+         SparseMatrix &A10 = a10.SpMat();
+
+         SparseMatrix &A01 = *Transpose(A10);
+
+         BilinearForm a11(&L2fes);
+         a11.AddDomainIntegrator(new MassIntegrator(neg_exp_psi));
+         ConstantCoefficient eps_cf(-1e-6);
+         if (order == 1)
+         {
+            a11.AddDomainIntegrator(new MassIntegrator(eps_cf));
+         }
+         else
+         {
+            a11.AddDomainIntegrator(new DiffusionIntegrator(eps_cf));
+         }
+         a11.Assemble();
+         a11.Finalize();
+         SparseMatrix &A11 = a11.SpMat();
+
+         BlockOperator A(offsets);
+         A.SetBlock(0,0,&A00);
+         A.SetBlock(1,0,&A10);
+         A.SetBlock(0,1,&A01);
+         A.SetBlock(1,1,&A11);
+
+         BlockDiagonalPreconditioner prec(offsets);
+         prec.SetDiagonalBlock(0,new GSSmoother(A00));
+         prec.SetDiagonalBlock(1,new GSSmoother(A11));
+
+         GMRES(A,prec,rhs,x,0,200, 50, 1e-12,0.0);
+
+         u_gf.MakeRef(&H1fes, x.GetBlock(0), 0);
+         delta_psi_gf.MakeRef(&L2fes, x.GetBlock(1), 0);
+
+         u_tmp -= u_gf;
+         double Newton_update_size = u_tmp.ComputeL2Error(zero);
+         u_tmp = u_gf;
+
+         double gamma = 1.0;
+         delta_psi_gf *= gamma;
+         psi_gf += delta_psi_gf;
+
+         if (visualization)
+         {
+            sol_sock << "solution\n" << mesh << u_gf << "window_title 'Discrete solution'"
+                     << flush;
+            mfem::out << "Newton_update_size = " << Newton_update_size << endl;
+         }
+
+         if (Newton_update_size < increment_u)
+         {
+            break;
+         }
+      }
+
+      u_tmp = u_gf;
+      u_tmp -= u_old_gf;
+      increment_u = u_tmp.ComputeL2Error(zero);
+
+      mfem::out << "Number of Newton iterations = " << j+1 << endl;
+      mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << increment_u << endl;
+
+      u_old_gf = u_gf;
+      psi_old_gf = psi_gf;
+
+      if (increment_u < tol || k == max_it-1)
+      {
+         break;
+      }
+
+      double L2_error = u_gf.ComputeL2Error(exact_coef);
+      mfem::out << "L2-error  (|| u - uₕ||)       = " << L2_error << endl;
+
+   }
+
+   mfem::out << "\n Outer iterations: " << k+1
+             << "\n Total iterations: " << total_iterations
+             << "\n dofs:             " << H1fes.GetTrueVSize() + L2fes.GetTrueVSize()
+             << endl;
+
+   // 11. Exact solution.
+   if (visualization)
+   {
+      socketstream err_sock(vishost, visport);
+      err_sock.precision(8);
+
+      GridFunction error(&H1fes);
+      error = 0.0;
+      error.ProjectCoefficient(exact_coef);
+      error -= u_gf;
+
+      err_sock << "solution\n" << mesh << error << "window_title 'Error'"  << flush;
+   }
+
+   {
+      ExponentialGridFunctionCoefficient exp_psi(psi_gf,obstacle);
+      u_alt_gf.ProjectCoefficient(exp_psi);
+      error_gf = 0.0;
+      error_gf.ProjectCoefficient(exact_coef);
+      error_gf -= u_alt_gf;
+      error_gf *= -1.0;
+
+      double L2_error = u_gf.ComputeL2Error(exact_coef);
+      double H1_error = u_gf.ComputeH1Error(&exact_coef,&exact_grad_coef);
+      double L2_error_alt = u_alt_gf.ComputeL2Error(exact_coef);
+
+      mfem::out << "\n Final L2-error (|| u - uₕ||)          = " << L2_error <<
+                endl;
+      mfem::out << " Final H1-error (|| u - uₕ||)          = " << H1_error << endl;
+      mfem::out << " Final L2-error (|| u - ϕ - exp(ψₕ)||) = " << L2_error_alt <<
+                endl;
+   }
+
+   return 0;
+}
+
+double LogarithmGridFunctionCoefficient::Eval(ElementTransformation &T,
+                                              const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(u != NULL, "grid function is not set");
+
+   double val = u->GetValue(T, ip) - obstacle->Eval(T, ip);
+   return max(min_val, log(val));
+}
+
+double ExponentialGridFunctionCoefficient::Eval(ElementTransformation &T,
+                                                const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(u != NULL, "grid function is not set");
+
+   double val = u->GetValue(T, ip);
+   return min(max_val, max(min_val, exp(val) + obstacle->Eval(T, ip)));
+}
+
+double spherical_obstacle(const Vector &pt)
+{
+   double x = pt(0), y = pt(1);
+   double r = sqrt(x*x + y*y);
+   double r0 = 0.5;
+   double beta = 0.9;
+
+   double b = r0*beta;
+   double tmp = sqrt(r0*r0 - b*b);
+   double B = tmp + b*b/tmp;
+   double C = -b/tmp;
+
+   if (r > b)
+   {
+      return B + r * C;
+   }
+   else
+   {
+      return sqrt(r0*r0 - r*r);
+   }
+}
+
+void spherical_obstacle_gradient(const Vector &pt, Vector &grad)
+{
+   double x = pt(0), y = pt(1);
+   double r = sqrt(x*x + y*y);
+   double r0 = 0.5;
+   double beta = 0.9;
+
+   double b = r0*beta;
+   double tmp = sqrt(r0*r0-b*b);
+   double C = -b/tmp;
+
+   if (r > b)
+   {
+      grad(0) = C * x / r;
+      grad(1) = C * y / r;
+   }
+   else
+   {
+      grad(0) = - x / sqrt( r0*r0 - r*r );
+      grad(1) = - y / sqrt( r0*r0 - r*r );
+   }
+}
+
+double exact_solution_obstacle(const Vector &pt)
+{
+   double x = pt(0), y = pt(1);
+   double r = sqrt(x*x + y*y);
+   double r0 = 0.5;
+   double a =  0.348982574111686;
+   double A = -0.340129705945858;
+
+   if (r > a)
+   {
+      return A * log(r);
+   }
+   else
+   {
+      return sqrt(r0*r0-r*r);
+   }
+}
+
+void exact_solution_gradient_obstacle(const Vector &pt, Vector &grad)
+{
+   double x = pt(0), y = pt(1);
+   double r = sqrt(x*x + y*y);
+   double r0 = 0.5;
+   double a =  0.348982574111686;
+   double A = -0.340129705945858;
+
+   if (r > a)
+   {
+      grad(0) =  A * x / (r*r);
+      grad(1) =  A * y / (r*r);
+   }
+   else
+   {
+      grad(0) = - x / sqrt( r0*r0 - r*r );
+      grad(1) = - y / sqrt( r0*r0 - r*r );
+   }
+}

examples/ex34p.cpp

@@ -0,0 +1,551 @@
+//                         MFEM Example 34 - Parallel Version
+//
+//
+// Compile with: make ex34p
+//
+// Sample runs: mpirun -np 2 ex34p -o 2
+//              mpirun -np 2 ex34p -o 2 -r 4
+//
+//
+// Description: This example code demonstrates the use of MFEM to solve the
+//              bound-constrained energy minimization problem
+//
+//                      minimize ||∇u||² subject to u ≥ ϕ in H¹₀.
+//
+//              This is known as the obstacle problem, and it is a simple
+//              mathematical model for contact mechanics.
+//
+//              In this example, the obstacle ϕ is a half-sphere centered
+//              at the origin of a circular domain Ω. After solving to a
+//              specified tolerance, the numerical solution is compared to
+//              a closed-form exact solution to assess accuracy.
+//
+//              The problem is discretized and solved using the entropic
+//              finite element method (EFEM) introduced by Keith and
+//              Surowiec [1].
+//
+//              This example highlights the ability of MFEM to deliver high-
+//              order solutions to variation inequality problems and
+//              showcases how to set up and solve nonlinear mixed methods.
+//
+//
+// [1] Keith, B. and Surowiec, T. (2023) The entropic finite element method
+//     (in preparation).
+
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+double spherical_obstacle(const Vector &pt);
+void spherical_obstacle_gradient(const Vector &pt, Vector &grad);
+double exact_solution_obstacle(const Vector &pt);
+void exact_solution_gradient_obstacle(const Vector &pt, Vector &grad);
+
+class LogarithmGridFunctionCoefficient : public Coefficient
+{
+protected:
+   GridFunction *u; // grid function
+   Coefficient *obstacle;
+   double min_val;
+
+public:
+   LogarithmGridFunctionCoefficient(GridFunction &u_, Coefficient &obst_,
+                                    double min_val_=-36)
+      : u(&u_), obstacle(&obst_), min_val(min_val_) { }
+
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+};
+
+class ExponentialGridFunctionCoefficient : public Coefficient
+{
+protected:
+   GridFunction *u; // grid function
+   Coefficient *obstacle;
+   double min_val;
+   double max_val;
+
+public:
+   ExponentialGridFunctionCoefficient(GridFunction &u_, Coefficient &obst_,
+                                      double min_val_=0.0, double max_val_=1e6)
+      : u(&u_), obstacle(&obst_), min_val(min_val_), max_val(max_val_) { }
+
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+};
+
+int main(int argc, char *argv[])
+{
+   // 0. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/disk.mesh";
+   int order = 1;
+   bool visualization = true;
+   int max_it = 10;
+   double tol = 1e-5;
+   int ref_levels = 3;
+   double alpha0 = 1.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  "isoparametric space.");
+   args.AddOption(&ref_levels, "-r", "--refs",
+                  "Number of h-refinements.");
+   args.AddOption(&max_it, "-mi", "--max-it",
+                  "Maximum number of iterations");
+   args.AddOption(&tol, "-tol", "--tol",
+                  "Stopping criteria based on the difference between"
+                  "successive solution updates");
+   args.AddOption(&alpha0, "-step", "--step",
+                  "Initial step size alpha");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 2. Read the mesh from the given mesh file.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 3. Refine the mesh to increase the resolution.
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   int curvature_order = max(order,2);
+   mesh.SetCurvature(curvature_order);
+   mesh.EnsureNCMesh();
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   // 4. Define the necessary finite element spaces on the mesh.
+   H1_FECollection H1fec(order, dim);
+   ParFiniteElementSpace H1fes(&pmesh, &H1fec);
+
+   L2_FECollection L2fec(order-1, dim);
+   ParFiniteElementSpace L2fes(&pmesh, &L2fec);
+
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: "
+           << H1fes.GetTrueVSize()
+           << " "
+           << L2fes.GetTrueVSize() << endl;
+   }
+
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = H1fes.GetVSize();
+   offsets[2] = L2fes.GetVSize();
+   offsets.PartialSum();
+
+   Array<int> toffsets(3);
+   toffsets[0] = 0;
+   toffsets[1] = H1fes.GetTrueVSize();
+   toffsets[2] = L2fes.GetTrueVSize();
+   toffsets.PartialSum();
+
+   BlockVector x(offsets), rhs(offsets);
+   x = 0.0; rhs = 0.0;
+
+   BlockVector tx(toffsets), trhs(toffsets);
+   tx = 0.0; trhs = 0.0;
+
+   // 5. Determine the list of true (i.e. conforming) essential boundary dofs.
+   Array<int> empty;
+   Array<int> ess_tdof_list;
+   if (pmesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      H1fes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 6. Define an initial guess for the solution.
+   auto IC_func = [](const Vector &x)
+   {
+      double r0 = 1.0;
+      double rr = 0.0;
+      for (int i=0; i<x.Size(); i++)
+      {
+         rr += x(i)*x(i);
+      }
+      return r0*r0 - rr;
+   };
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient zero(0.0);
+
+   // 7. Define the solution vectors as a finite element grid functions
+   //    corresponding to the fespaces.
+   ParGridFunction u_gf, delta_psi_gf;
+   u_gf.MakeRef(&H1fes,x.GetBlock(0).GetData());
+   delta_psi_gf.MakeRef(&L2fes,x.GetBlock(1).GetData());
+   delta_psi_gf = 0.0;
+
+   ParGridFunction u_old_gf(&H1fes);
+   ParGridFunction psi_old_gf(&L2fes);
+   ParGridFunction psi_gf(&L2fes);
+   u_old_gf = 0.0;
+   psi_old_gf = 0.0;
+
+   // 8. Define the function coefficients for the solution and use them to
+   //    initialize the initial guess
+   FunctionCoefficient exact_coef(exact_solution_obstacle);
+   VectorFunctionCoefficient exact_grad_coef(dim,exact_solution_gradient_obstacle);
+   FunctionCoefficient IC_coef(IC_func);
+   ConstantCoefficient f(0.0);
+   FunctionCoefficient obstacle(spherical_obstacle);
+   u_gf.ProjectCoefficient(IC_coef);
+   u_old_gf = u_gf;
+
+   // 9. Initialize the slack variable ψₕ = exp(uₕ)
+   LogarithmGridFunctionCoefficient ln_u(u_gf, obstacle);
+   psi_gf.ProjectCoefficient(ln_u);
+   psi_old_gf = psi_gf;
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+
+   ParGridFunction u_alt_gf(&L2fes);
+   ParGridFunction error_gf(&L2fes);
+
+   ExponentialGridFunctionCoefficient exp_psi(psi_gf,obstacle);
+   u_alt_gf.ProjectCoefficient(exp_psi);
+
+   if (visualization)
+   {
+      sol_sock.open(vishost,visport);
+      sol_sock.precision(8);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock << "solution\n" << pmesh << u_alt_gf <<
+               "window_title 'Discrete solution'" << flush;
+   }
+
+   // 10. Iterate
+   int k;
+   int total_iterations = 0;
+   double increment_u = 0.1;
+   for (k = 0; k < max_it; k++)
+   {
+      double alpha = alpha0 * (k+1);
+
+      ParGridFunction u_tmp(&H1fes);
+      u_tmp = u_old_gf;
+
+      if (myid == 0)
+      {
+         mfem::out << "\nOUTER ITERATION " << k+1 << endl;
+      }
+
+      int j;
+      for ( j = 0; j < 15; j++)
+      {
+         total_iterations++;
+
+         ConstantCoefficient alpha_cf(alpha);
+
+         ParLinearForm b0,b1;
+         b0.Update(&H1fes,rhs.GetBlock(0),0);
+         b1.Update(&L2fes,rhs.GetBlock(1),0);
+
+         ExponentialGridFunctionCoefficient exp_psi(psi_gf, zero);
+         ProductCoefficient neg_exp_psi(-1.0,exp_psi);
+         GradientGridFunctionCoefficient grad_u_old(&u_old_gf);
+         ProductCoefficient alpha_f(alpha, f);
+         GridFunctionCoefficient psi_cf(&psi_gf);
+         GridFunctionCoefficient psi_old_cf(&psi_old_gf);
+         SumCoefficient psi_old_minus_psi(psi_old_cf, psi_cf, 1.0, -1.0);
+
+         b0.AddDomainIntegrator(new DomainLFIntegrator(alpha_f));
+         b0.AddDomainIntegrator(new DomainLFIntegrator(psi_old_minus_psi));
+         b0.Assemble();
+
+         b1.AddDomainIntegrator(new DomainLFIntegrator(exp_psi));
+         b1.AddDomainIntegrator(new DomainLFIntegrator(obstacle));
+         b1.Assemble();
+
+         ParBilinearForm a00(&H1fes);
+         a00.SetDiagonalPolicy(mfem::Operator::DIAG_ONE);
+         a00.AddDomainIntegrator(new DiffusionIntegrator(alpha_cf));
+         a00.Assemble();
+         HypreParMatrix A00;
+         a00.FormLinearSystem(ess_tdof_list, x.GetBlock(0), rhs.GetBlock(0),
+                              A00, tx.GetBlock(0), trhs.GetBlock(0));
+
+
+         ParMixedBilinearForm a10(&H1fes,&L2fes);
+         a10.AddDomainIntegrator(new MixedScalarMassIntegrator());
+         a10.Assemble();
+         HypreParMatrix A10;
+         a10.FormRectangularLinearSystem(ess_tdof_list, empty, x.GetBlock(0),
+                                         rhs.GetBlock(1),
+                                         A10, tx.GetBlock(0), trhs.GetBlock(1));
+
+         HypreParMatrix &A01 = *A10.Transpose();
+
+         ParBilinearForm a11(&L2fes);
+         a11.AddDomainIntegrator(new MassIntegrator(neg_exp_psi));
+         ConstantCoefficient eps_cf(-1e-6);
+         if (order == 1)
+         {
+            a11.AddDomainIntegrator(new MassIntegrator(eps_cf));
+         }
+         else
+         {
+            a11.AddDomainIntegrator(new DiffusionIntegrator(eps_cf));
+         }
+         a11.Assemble();
+         a11.Finalize();
+         HypreParMatrix A11;
+         a11.FormSystemMatrix(empty, A11);
+
+         BlockOperator A(toffsets);
+         A.SetBlock(0,0,&A00);
+         A.SetBlock(1,0,&A10);
+         A.SetBlock(0,1,&A01);
+         A.SetBlock(1,1,&A11);
+
+         BlockDiagonalPreconditioner prec(toffsets);
+         HypreBoomerAMG P00(A00);
+         P00.SetPrintLevel(0);
+         HypreSmoother P11(A11);
+         prec.SetDiagonalBlock(0,&P00);
+         prec.SetDiagonalBlock(1,new HypreSmoother(A11));
+
+         GMRESSolver gmres(MPI_COMM_WORLD);
+         gmres.SetPrintLevel(-1);
+         gmres.SetRelTol(1e-8);
+         gmres.SetMaxIter(20000);
+         gmres.SetKDim(500);
+         gmres.SetOperator(A);
+         gmres.SetPreconditioner(prec);
+         gmres.Mult(trhs,tx);
+
+         u_gf.SetFromTrueDofs(tx.GetBlock(0));
+         delta_psi_gf.SetFromTrueDofs(tx.GetBlock(1));
+
+         u_tmp -= u_gf;
+         double Newton_update_size = u_tmp.ComputeL2Error(zero);
+         u_tmp = u_gf;
+
+         double gamma = 1.0;
+         delta_psi_gf *= gamma;
+         psi_gf += delta_psi_gf;
+
+         if (visualization)
+         {
+            sol_sock << "parallel " << num_procs << " " << myid << "\n";
+            sol_sock << "solution\n" << pmesh << u_gf << "window_title 'Discrete solution'"
+                     << flush;
+         }
+
+         if (myid == 0)
+         {
+            mfem::out << "Newton_update_size = " << Newton_update_size << endl;
+         }
+
+         if (Newton_update_size < increment_u)
+         {
+            break;
+         }
+      }
+
+      u_tmp = u_gf;
+      u_tmp -= u_old_gf;
+      increment_u = u_tmp.ComputeL2Error(zero);
+
+      if (myid == 0)
+      {
+         mfem::out << "Number of Newton iterations = " << j+1 << endl;
+         mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << increment_u << endl;
+      }
+
+      u_old_gf = u_gf;
+      psi_old_gf = psi_gf;
+
+      if (increment_u < tol || k == max_it-1)
+      {
+         break;
+      }
+
+      double L2_error = u_gf.ComputeL2Error(exact_coef);
+      if (myid == 0)
+      {
+         mfem::out << "L2-error  (|| u - uₕ||)       = " << L2_error << endl;
+      }
+
+   }
+
+   if (myid == 0)
+   {
+      mfem::out << "\n Outer iterations: " << k+1
+                << "\n Total iterations: " << total_iterations
+                << "\n dofs:             " << H1fes.GetTrueVSize() + L2fes.GetTrueVSize()
+                << endl;
+   }
+
+   // 11. Exact solution.
+   if (visualization)
+   {
+      socketstream err_sock(vishost, visport);
+      err_sock.precision(8);
+
+      ParGridFunction error(&H1fes);
+      error = 0.0;
+      error.ProjectCoefficient(exact_coef);
+      error -= u_gf;
+
+      err_sock << "parallel " << num_procs << " " << myid << "\n";
+      err_sock << "solution\n" << pmesh << error << "window_title 'Error'"  << flush;
+   }
+
+   {
+      ExponentialGridFunctionCoefficient exp_psi(psi_gf,obstacle);
+      u_alt_gf.ProjectCoefficient(exp_psi);
+      error_gf = 0.0;
+      error_gf.ProjectCoefficient(exact_coef);
+      error_gf -= u_alt_gf;
+      error_gf *= -1.0;
+
+      double L2_error = u_gf.ComputeL2Error(exact_coef);
+      double H1_error = u_gf.ComputeH1Error(&exact_coef,&exact_grad_coef);
+      double L2_error_alt = u_alt_gf.ComputeL2Error(exact_coef);
+
+      if (myid == 0)
+      {
+         mfem::out << "\n Final L2-error (|| u - uₕ||)          = " << L2_error <<
+                   endl;
+         mfem::out << " Final H1-error (|| u - uₕ||)          = " << H1_error << endl;
+         mfem::out << " Final L2-error (|| u - ϕ - exp(ψₕ)||) = " << L2_error_alt <<
+                   endl;
+      }
+   }
+
+   return 0;
+}
+
+double LogarithmGridFunctionCoefficient::Eval(ElementTransformation &T,
+                                              const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(u != NULL, "grid function is not set");
+
+   double val = u->GetValue(T, ip) - obstacle->Eval(T, ip);
+   return max(min_val, log(val));
+}
+
+double ExponentialGridFunctionCoefficient::Eval(ElementTransformation &T,
+                                                const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(u != NULL, "grid function is not set");
+
+   double val = u->GetValue(T, ip);
+   return min(max_val, max(min_val, exp(val) + obstacle->Eval(T, ip)));
+}
+
+double spherical_obstacle(const Vector &pt)
+{
+   double x = pt(0), y = pt(1);
+   double r = sqrt(x*x + y*y);
+   double r0 = 0.5;
+   double beta = 0.9;
+
+   double b = r0*beta;
+   double tmp = sqrt(r0*r0 - b*b);
+   double B = tmp + b*b/tmp;
+   double C = -b/tmp;
+
+   if (r > b)
+   {
+      return B + r * C;
+   }
+   else
+   {
+      return sqrt(r0*r0 - r*r);
+   }
+}
+
+void spherical_obstacle_gradient(const Vector &pt, Vector &grad)
+{
+   double x = pt(0), y = pt(1);
+   double r = sqrt(x*x + y*y);
+   double r0 = 0.5;
+   double beta = 0.9;
+
+   double b = r0*beta;
+   double tmp = sqrt(r0*r0-b*b);
+   double C = -b/tmp;
+
+   if (r > b)
+   {
+      grad(0) = C * x / r;
+      grad(1) = C * y / r;
+   }
+   else
+   {
+      grad(0) = - x / sqrt( r0*r0 - r*r );
+      grad(1) = - y / sqrt( r0*r0 - r*r );
+   }
+}
+
+double exact_solution_obstacle(const Vector &pt)
+{
+   double x = pt(0), y = pt(1);
+   double r = sqrt(x*x + y*y);
+   double r0 = 0.5;
+   double a =  0.348982574111686;
+   double A = -0.340129705945858;
+
+   if (r > a)
+   {
+      return A * log(r);
+   }
+   else
+   {
+      return sqrt(r0*r0-r*r);
+   }
+}
+
+void exact_solution_gradient_obstacle(const Vector &pt, Vector &grad)
+{
+   double x = pt(0), y = pt(1);
+   double r = sqrt(x*x + y*y);
+   double r0 = 0.5;
+   double a =  0.348982574111686;
+   double A = -0.340129705945858;
+
+   if (r > a)
+   {
+      grad(0) =  A * x / (r*r);
+      grad(1) =  A * y / (r*r);
+   }
+   else
+   {
+      grad(0) = - x / sqrt( r0*r0 - r*r );
+      grad(1) = - y / sqrt( r0*r0 - r*r );
+   }
+}

examples/ex35.cpp

@@ -0,0 +1,473 @@
+//                              MFEM Example 35
+//
+//
+// Compile with: make ex35
+//
+// Sample runs:
+//     ex35 -alpha 10
+//     ex35 -lambda 0.1 -mu 0.1
+//     ex35 -r 5 -o 2 -alpha 5.0 -epsilon 0.01 -mi 50 -mf 0.5 -tol 1e-5
+//     ex35 -r 6 -o 1 -alpha 10.0 -epsilon 0.01 -mi 50 -mf 0.5 -tol 1e-5
+//
+//
+// Description: This example code demonstrates the use of MFEM to solve a
+//              density-filtered [3] topology optimization problem. The
+//              objective is to minimize the compliance
+//
+//                  minimize ∫_Ω f⋅u dx over u ∈ [H¹(Ω)]² and ρ ∈ L²(Ω)
+//
+//                  subject to
+//
+//                    -Div(r(ρ̃)Cε(u)) = f       in Ω + BCs
+//                    -ϵ²Δρ̃ + ρ̃ = ρ             in Ω + Neumann BCs
+//                    0 ≤ ρ ≤ 1                 in Ω
+//                    ∫_Ω ρ dx = θ vol(Ω)
+//
+//              Here, r(ρ̃) = ρ₀ + ρ̃³ (1-ρ₀) is the solid isotropic material
+//              penalization (SIMP) law, C is the elasticity tensor for an
+//              isotropic linearly elastic material, ϵ > 0 is the design
+//              length scale, and 0 < θ < 1 is the volume fraction.
+//
+//              The problem is discretized and gradients are computing using
+//              finite elements [1]. The design is optimized using an entropic
+//              mirror descent algorithm introduced by Keith and Surowiec [2]
+//              that is tailored to the bound constraint 0 ≤ ρ ≤ 1.
+//
+//              This example highlights the ability of MFEM to deliver high-
+//              order solutions to inverse design problems and showcases how
+//              to set up and solve PDE-constrained optimization problems
+//              using the so-called reduced space approach.
+//
+//
+// [1] Andreassen, E., Clausen, A., Schevenels, M., Lazarov, B. S., & Sigmund, O.
+//    (2011). Efficient topology optimization in MATLAB using 88 lines of
+//    code. Structural and Multidisciplinary Optimization, 43(1), 1-16.
+// [2] Keith, B. and Surowiec, T. (2023) The entropic finite element method
+//     (in preparation).
+// [3] Lazarov, B. S., & Sigmund, O. (2011). Filters in topology optimization
+//     based on Helmholtz‐type differential equations. International Journal
+//     for Numerical Methods in Engineering, 86(6), 765-781.
+
+#include "mfem.hpp"
+#include <iostream>
+#include <fstream>
+#include "ex35.hpp"
+
+/**
+ * @brief Nonlinear projection of 0 < τ < 1 onto the subspace
+ *        ∫_Ω τ dx = θ vol(Ω) as follows.
+ *
+ *        1. Compute the root of the R → R function
+ *            f(c) = ∫_Ω expit(lnit(τ) + c) dx - θ vol(Ω)
+ *        2. Set τ ← expit(lnit(τ) + c).
+ *
+ */
+void projit(GridFunction &tau, double &c, LinearForm &vol_form,
+            double volume_fraction, double tol=1e-12, int max_its=10)
+{
+   GridFunction ftmp(tau.FESpace());
+   GridFunction dftmp(tau.FESpace());
+   for (int k=0; k<max_its; k++)
+   {
+      // Compute f(c) and dfdc(c)
+      for (int i=0; i<tau.Size(); i++)
+      {
+         ftmp[i]  = expit(lnit(tau[i]) + c) - volume_fraction;
+         dftmp[i] = dexpitdx(lnit(tau[i]) + c);
+      }
+      double f = vol_form(ftmp);
+      double df = vol_form(dftmp);
+
+      double dc = -f/df;
+      c += dc;
+      if (abs(dc) < tol) { break; }
+   }
+   tau = ftmp;
+   tau += volume_fraction;
+}
+
+using namespace std;
+using namespace mfem;
+/**
+ * ---------------------------------------------------------------
+ *                      ALGORITHM PREAMBLE
+ * ---------------------------------------------------------------
+ *
+ *  The Lagrangian for this problem is
+ *
+ *          L(u,ρ,ρ̃,w,w̃) = (f,u) - (r(ρ̃) C ε(u),ε(w)) + (f,w)
+ *                       - (ϵ² ∇ρ̃,∇w̃) - (ρ̃,w̃) + (ρ,w̃)
+ *
+ *  where
+ *
+ *    r(ρ̃) = ρ₀ + ρ̃³ (1 - ρ₀)       (SIMP rule)
+ *
+ *    ε(u) = (∇u + ∇uᵀ)/2           (symmetric gradient)
+ *
+ *    C e = λtr(e)I + 2μe           (isotropic material)
+ *
+ *  NOTE: The Lame parameters can be computed from Young's modulus E
+ *        and Poisson's ratio ν as follows:
+ *
+ *             λ = E ν/((1+ν)(1-2ν)),      μ = E/(2(1+ν))
+ *
+ * ---------------------------------------------------------------
+ *
+ *  Discretization choices:
+ *
+ *     u ∈ V ⊂ (H¹)ᵈ (order p)
+ *     ρ ∈ L² (order p - 1)
+ *     ρ̃ ∈ H¹ (order p - 1)
+ *     w ∈ V  (order p)
+ *     w̃ ∈ H¹ (order p - 1)
+ *
+ * ---------------------------------------------------------------
+ *                          ALGORITHM
+ * ---------------------------------------------------------------
+ *
+ *  Update ρ with projected mirror descent via the following algorithm.
+ *
+ *  1. Initialize density field 0 < ρ(x) < 1.
+ *
+ *  While not converged:
+ *
+ *     2. Solve filter equation ∂_w̃ L = 0; i.e.,
+ *
+ *           (ϵ² ∇ ρ̃, ∇ v ) + (ρ̃,v) = (ρ,v)   ∀ v ∈ H¹.
+ *
+ *     3. Solve primal problem ∂_w L = 0; i.e.,
+ *
+ *      (λ(ρ̃) ∇⋅u, ∇⋅v) + (2 μ(ρ̃) ε(u), ε(v)) = (f,v)   ∀ v ∈ V,
+ *
+ *     where λ(ρ̃) := λ r(ρ̃) and  μ(ρ̃) := μ r(ρ̃).
+ *
+ *     NB. The dual problem ∂_u L = 0 is the same as the primal problem due to symmetry.
+ *
+ *     4. Solve for filtered gradient ∂_ρ̃ L = 0; i.e.,
+ *
+ *      (ϵ² ∇ w̃ , ∇ v ) + (w̃ ,v) = (-r'(ρ̃) ( λ(ρ̃) |∇⋅u|² + 2 μ(ρ̃) |ε(u)|²),v)   ∀ v ∈ H¹.
+ *
+ *     5. Construct gradient G ∈ L²; i.e.,
+ *
+ *                         (G,v) = (w̃,v)   ∀ v ∈ L².
+ *
+ *     6. Mirror descent update until convergence; i.e.,
+ *
+ *                      ρ ← projit(expit(linit(ρ) - αG)),
+ *
+ *     where
+ *
+ *          α > 0                            (step size parameter)
+ *
+ *          expit(x) = eˣ/(1+eˣ)             (sigmoid)
+ *
+ *          linit(y) = ln(y) - ln(1-y)       (inverse of sigmoid)
+ *
+ *     and projit is a (compatible) projection operator enforcing ∫_Ω ρ dx = θ vol(Ω).
+ *
+ *  end
+ *
+ */
+
+int main(int argc, char *argv[])
+{
+
+   // 1. Parse command-line options.
+   int ref_levels = 4;
+   int order = 2;
+   bool visualization = true;
+   double alpha0 = 1.0;
+   double epsilon = 0.01;
+   double mass_fraction = 0.5;
+   int max_it = 1e2;
+   double tol = 1e-4;
+   double rho_min = 1e-6;
+   double lambda = 1.0;
+   double mu = 1.0;
+   double exponent = 3;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&alpha0, "-alpha", "--alpha-step-length",
+                  "Step length for gradient descent.");
+   args.AddOption(&epsilon, "-epsilon", "--epsilon-thickness",
+                  "epsilon phase field thickness");
+   args.AddOption(&max_it, "-mi", "--max-it",
+                  "Maximum number of gradient descent iterations.");
+   args.AddOption(&tol, "-tol", "--tol",
+                  "Exit tolerance for ρ ");
+   args.AddOption(&mass_fraction, "-mf", "--mass-fraction",
+                  "Mass fraction for diffusion coefficient.");
+   args.AddOption(&lambda, "-lambda", "--lambda",
+                  "Lame constant λ");
+   args.AddOption(&mu, "-mu", "--mu",
+                  "Lame constant μ");
+   args.AddOption(&rho_min, "-rmin", "--rho-min",
+                  "Minimum of density coefficient.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   Mesh mesh = Mesh::MakeCartesian2D(3,1,mfem::Element::Type::QUADRILATERAL,true,
+                                     3.0,1.0);
+
+   int dim = mesh.Dimension();
+
+   // 2. Set BCs.
+   for (int i = 0; i<mesh.GetNBE(); i++)
+   {
+      Element * be = mesh.GetBdrElement(i);
+      Array<int> vertices;
+      be->GetVertices(vertices);
+
+      double * coords1 = mesh.GetVertex(vertices[0]);
+      double * coords2 = mesh.GetVertex(vertices[1]);
+
+      Vector center(2);
+      center(0) = 0.5*(coords1[0] + coords2[0]);
+      center(1) = 0.5*(coords1[1] + coords2[1]);
+
+      if (abs(center(0) - 0.0) < 1e-10)
+      {
+         // the left edge
+         be->SetAttribute(1);
+      }
+      else
+      {
+         // all other boundaries
+         be->SetAttribute(2);
+      }
+   }
+   mesh.SetAttributes();
+
+   // 3. Refine the mesh.
+   for (int lev = 0; lev < ref_levels; lev++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   // 4. Define the necessary finite element spaces on the mesh.
+   H1_FECollection state_fec(order, dim); // space for u
+   H1_FECollection filter_fec(order, dim); // space for ρ̃
+   L2_FECollection control_fec(order-1, dim,
+                               BasisType::Positive); // space for ρ
+   FiniteElementSpace state_fes(&mesh, &state_fec,dim);
+   FiniteElementSpace filter_fes(&mesh, &filter_fec);
+   FiniteElementSpace control_fes(&mesh, &control_fec);
+
+   int state_size = state_fes.GetTrueVSize();
+   int control_size = control_fes.GetTrueVSize();
+   int filter_size = filter_fes.GetTrueVSize();
+   cout << "Number of state unknowns: " << state_size << endl;
+   cout << "Number of filter unknowns: " << filter_size << endl;
+   cout << "Number of control unknowns: " << control_size << endl;
+
+   // 5. Set the initial guess for ρ.
+   GridFunction u(&state_fes);
+   GridFunction rho(&control_fes);
+   GridFunction rho_old(&control_fes);
+   GridFunction rho_filter(&filter_fes);
+   u = 0.0;
+   rho_filter = 0.0;
+   rho = 0.5;
+   rho_old = 0.5;
+
+   // 6. Set-up the physics solver.
+   int maxat = mesh.bdr_attributes.Max();
+   Array<int> ess_bdr(maxat);
+   ess_bdr = 0;
+   ess_bdr[0] = 1;
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient lambda_cf(lambda);
+   ConstantCoefficient mu_cf(mu);
+   LinearElasticitySolver * ElasticitySolver = new LinearElasticitySolver();
+   ElasticitySolver->SetMesh(&mesh);
+   ElasticitySolver->SetOrder(state_fec.GetOrder());
+   ElasticitySolver->SetupFEM();
+   Vector center(2); center(0) = 2.9; center(1) = 0.5;
+   Vector force(2); force(0) = 0.0; force(1) = -1.0;
+   double r = 0.05;
+   VolumeForceCoefficient vforce_cf(r,center,force);
+   ElasticitySolver->SetRHSCoefficient(&vforce_cf);
+   ElasticitySolver->SetEssentialBoundary(ess_bdr);
+
+   // 7. Set-up the filter solver.
+   ConstantCoefficient eps2_cf(epsilon*epsilon);
+   DiffusionSolver * FilterSolver = new DiffusionSolver();
+   FilterSolver->SetMesh(&mesh);
+   FilterSolver->SetOrder(filter_fec.GetOrder());
+   FilterSolver->SetDiffusionCoefficient(&eps2_cf);
+   FilterSolver->SetMassCoefficient(&one);
+   Array<int> ess_bdr_filter;
+   if (mesh.bdr_attributes.Size())
+   {
+      ess_bdr_filter.SetSize(mesh.bdr_attributes.Max());
+      ess_bdr_filter = 0;
+   }
+   FilterSolver->SetEssentialBoundary(ess_bdr_filter);
+   FilterSolver->SetupFEM();
+
+   BilinearForm mass(&control_fes);
+   mass.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator(one)));
+   mass.Assemble();
+   SparseMatrix M;
+   Array<int> empty;
+   mass.FormSystemMatrix(empty,M);
+
+   // 8. Define the Lagrange multiplier and gradient functions
+   GridFunction grad(&control_fes);
+   GridFunction w_filter(&filter_fes);
+
+   // 9. Define some tools for later
+   ConstantCoefficient zero(0.0);
+   GridFunction onegf(&control_fes);
+   onegf = 1.0;
+   LinearForm vol_form(&control_fes);
+   vol_form.AddDomainIntegrator(new DomainLFIntegrator(one));
+   vol_form.Assemble();
+   double domain_volume = vol_form(onegf);
+
+   // 10. Connect to GLVis. Prepare for VisIt output.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sout_u,sout_r,sout_rho;
+   if (visualization)
+   {
+      sout_u.open(vishost, visport);
+      sout_rho.open(vishost, visport);
+      sout_r.open(vishost, visport);
+      sout_u.precision(8);
+      sout_rho.precision(8);
+      sout_r.precision(8);
+
+      sout_u << "solution\n" << mesh << u;
+      sout_u << "view 0 0\n";  // view from top
+      sout_u << "keys jl********\n";  // turn off perspective and light
+      sout_u << "window_title 'Temperature u'";
+      sout_u.flush();
+
+      sout_rho << "solution\n" << mesh << rho;
+      sout_rho << "view 0 0\n";  // view from top
+      sout_rho << "keys jl********\n";  // turn off perspective and light
+      sout_rho << "window_title 'Density ρ'";
+      sout_rho.flush();
+
+      sout_r << "solution\n" << mesh << rho_filter;
+      sout_r << "view 0 0\n";  // view from top
+      sout_r << "keys jl********\n";  // turn off perspective and light
+      sout_r << "window_title 'Filtered density ρ̃'";
+      sout_r.flush();
+   }
+
+   mfem::ParaViewDataCollection paraview_dc("Elastic_compliance", &mesh);
+   paraview_dc.SetPrefixPath("ParaView");
+   paraview_dc.SetLevelsOfDetail(order);
+   paraview_dc.SetCycle(0);
+   paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   paraview_dc.SetHighOrderOutput(true);
+   paraview_dc.SetTime(0.0);
+   paraview_dc.RegisterField("displacement",&u);
+   paraview_dc.RegisterField("density",&rho);
+   paraview_dc.RegisterField("filtered_density",&rho_filter);
+
+   // 11. Iterate
+   int step = 0;
+   double c0 = 0.0;
+   for (int k = 1; k < max_it; k++)
+   {
+      const double alpha = alpha0 * k;
+
+      cout << "\nStep = " << k << endl;
+
+      // Step 1 - Filter solve
+      // Solve (ϵ^2 ∇ ρ̃, ∇ v ) + (ρ̃,v) = (ρ,v)
+      GridFunctionCoefficient rho_cf(&rho);
+      FilterSolver->SetRHSCoefficient(&rho_cf);
+      FilterSolver->Solve();
+      rho_filter = *FilterSolver->GetFEMSolution();
+
+      // Step 2 - State solve
+      // Solve (λ(ρ̃) ∇⋅u, ∇⋅v) + (2 μ(ρ̃) ε(u), ε(v)) = (f,v)
+      SIMPInterpolationCoefficient SIMP_cf(&rho_filter,rho_min, 1.0);
+      ProductCoefficient lambda_SIMP_cf(lambda_cf,SIMP_cf);
+      ProductCoefficient mu_SIMP_cf(mu_cf,SIMP_cf);
+      ElasticitySolver->SetLameCoefficients(&lambda_SIMP_cf,&mu_SIMP_cf);
+      ElasticitySolver->Solve();
+      u = *ElasticitySolver->GetFEMSolution();
+
+      // Step 3 - Adjoint filter solve
+      // Solve (ϵ² ∇ w̃, ∇ v) + (w̃ ,v) = (-r'(ρ̃) ( λ(ρ̃) |∇⋅u|² + 2 μ(ρ̃) |ε(u)|²),v)
+      StrainEnergyDensityCoefficient rhs_cf(&lambda_cf,&mu_cf,&u, &rho_filter,
+                                            rho_min);
+      FilterSolver->SetRHSCoefficient(&rhs_cf);
+      FilterSolver->Solve();
+      w_filter = *FilterSolver->GetFEMSolution();
+
+      // Step 4 - Compute gradient
+      // Solve G = M⁻¹w̃
+      GridFunctionCoefficient w_cf(&w_filter);
+      LinearForm w_rhs(&control_fes);
+      w_rhs.AddDomainIntegrator(new DomainLFIntegrator(w_cf));
+      w_rhs.Assemble();
+      M.Mult(w_rhs,grad);
+
+      // Step 5 - Update design variable ρ ← projit(expit(linit(ρ) - αG))
+      // Note: The update here is performed on the coefficients of the linear
+      //       representation of the design variable in the Berstein basis. It would
+      //       be more mathematically sound to update the design variable so that
+      //       the values at the integration points follow this update rule.
+      for (int i = 0; i < rho.Size(); i++)
+      {
+         rho[i] = expit(lnit(rho[i]) - alpha*grad[i]);
+      }
+      projit(rho, c0, vol_form, mass_fraction);
+
+      GridFunctionCoefficient tmp(&rho_old);
+      double norm_reduced_gradient = rho.ComputeL2Error(tmp)/alpha;
+      rho_old = rho;
+
+      double compliance = (*(ElasticitySolver->GetLinearForm()))(u);
+      double material_volume = vol_form(rho);
+      mfem::out << "norm of reduced gradient = " << norm_reduced_gradient << endl;
+      mfem::out << "compliance = " << compliance << endl;
+      mfem::out << "mass_fraction = " << material_volume / domain_volume << endl;
+
+      if (visualization)
+      {
+         sout_u << "solution\n" << mesh << u
+                << "window_title 'Displacement u'" << flush;
+
+         sout_rho << "solution\n" << mesh << rho
+                  << "window_title 'Control variable ρ'" << flush;
+
+         GridFunction r_gf(&filter_fes);
+         r_gf.ProjectCoefficient(SIMP_cf);
+         sout_r << "solution\n" << mesh << r_gf
+                << "window_title 'Design density r(ρ̃)'" << flush;
+
+         paraview_dc.SetCycle(k);
+         paraview_dc.SetTime((double)k);
+         paraview_dc.Save();
+      }
+
+      if (norm_reduced_gradient < tol)
+      {
+         break;
+      }
+   }
+
+   delete ElasticitySolver;
+   delete FilterSolver;
+
+   return 0;
+}

examples/ex35.hpp

@@ -0,0 +1,729 @@
+//                  MFEM Example 35 - Serial/Parallel Shared Code
+//
+//
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Inverse sigmoid function
+double lnit(double x)
+{
+   double tol = 1e-12;
+   x = min(max(tol,x),1.0-tol);
+   return log(x/(1.0-x));
+}
+
+// Sigmoid function
+double expit(double x)
+{
+   if (x >= 0)
+   {
+      return 1.0/(1.0+exp(-x));
+   }
+   else
+   {
+      return exp(x)/(1.0+exp(x));
+   }
+}
+
+// Derivative of sigmoid function
+double dexpitdx(double x)
+{
+   double tmp = expit(-x);
+   return tmp - pow(tmp,2);
+}
+
+//  Solid isotropic material penalization (SIMP) coefficient
+class SIMPInterpolationCoefficient : public Coefficient
+{
+protected:
+   GridFunction *rho_filter; // grid function
+   double min_val;
+   double max_val;
+   double exponent;
+
+public:
+   SIMPInterpolationCoefficient(GridFunction *rho_filter_, double min_val_= 1e-6,
+                                double max_val_ = 1.0, double exponent_ = 3)
+      : rho_filter(rho_filter_), min_val(min_val_), max_val(max_val_),
+        exponent(exponent_) { }
+
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      double val = rho_filter->GetValue(T, ip);
+      double coeff = min_val + pow(val,exponent)*(max_val-min_val);
+      return coeff;
+   }
+};
+
+
+// Strain energy density coefficient
+class StrainEnergyDensityCoefficient : public Coefficient
+{
+protected:
+   Coefficient * lambda=nullptr;
+   Coefficient * mu=nullptr;
+   GridFunction *u = nullptr; // displacement
+   GridFunction *rho_filter = nullptr; // filter density
+   DenseMatrix grad; // auxiliary matrix, used in Eval
+   double exponent;
+   double rho_min;
+
+public:
+   StrainEnergyDensityCoefficient(Coefficient *lambda_, Coefficient *mu_,
+                                  GridFunction * u_, GridFunction * rho_filter_, double rho_min_=1e-6,
+                                  double exponent_ = 3.0)
+      : lambda(lambda_), mu(mu_),  u(u_), rho_filter(rho_filter_),
+        exponent(exponent_), rho_min(rho_min_)
+   {
+      MFEM_ASSERT(rho_min_ >= 0.0, "rho_min must be >= 0");
+      MFEM_ASSERT(rho_min_ < 1.0,  "rho_min must be > 1");
+      MFEM_ASSERT(u, "displacement field is not set");
+      MFEM_ASSERT(rho_filter, "density field is not set");
+   }
+
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      double L = lambda->Eval(T, ip);
+      double M = mu->Eval(T, ip);
+      u->GetVectorGradient(T, grad);
+      double div_u = grad.Trace();
+      double density = L*div_u*div_u;
+      int dim = T.GetSpaceDim();
+      for (int i=0; i<dim; i++)
+      {
+         for (int j=0; j<dim; j++)
+         {
+            density += M*grad(i,j)*(grad(i,j)+grad(j,i));
+         }
+      }
+      double val = rho_filter->GetValue(T,ip);
+
+      return -exponent * pow(val, exponent-1.0) * (1-rho_min) * density;
+   }
+};
+
+
+// Strain energy density coefficient
+class DiffusionEnergyCoefficient : public Coefficient
+{
+protected:
+   Coefficient * K=nullptr;
+   GridFunction *u = nullptr; // displacement
+   GridFunction *rho_filter = nullptr; // filter density
+   Vector grad; // auxiliary matrix, used in Eval
+   double exponent;
+   double rho_min;
+
+public:
+   DiffusionEnergyCoefficient(Coefficient *K_,
+                                  GridFunction * u_, GridFunction * rho_filter_, double rho_min_=1e-6,
+                                  double exponent_ = 3.0)
+      : K(K_),  u(u_), rho_filter(rho_filter_),
+        exponent(exponent_), rho_min(rho_min_)
+   {
+      MFEM_ASSERT(rho_min_ >= 0.0, "rho_min must be >= 0");
+      MFEM_ASSERT(rho_min_ < 1.0,  "rho_min must be > 1");
+      MFEM_ASSERT(u, "displacement field is not set");
+      MFEM_ASSERT(rho_filter, "density field is not set");
+   }
+
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      double Kval = K->Eval(T, ip);
+      u->GetGradient(T, grad);
+      double density = Kval*(grad*grad);
+      double val = rho_filter->GetValue(T,ip);
+
+      return -exponent * pow(val, exponent-1.0) * (1-rho_min) * density;
+   }
+};
+
+// Volumetric force for linear elasticity
+class VolumeForceCoefficient : public VectorCoefficient
+{
+private:
+   double r;
+   Vector center;
+   Vector force;
+public:
+   VolumeForceCoefficient(double r_,Vector &  center_, Vector & force_) :
+      VectorCoefficient(center_.Size()), r(r_), center(center_), force(force_) { }
+
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      Vector xx; xx.SetSize(T.GetDimension());
+      T.Transform(ip,xx);
+      for (int i=0; i<xx.Size(); i++)
+      {
+         xx[i]=xx[i]-center[i];
+      }
+
+      double cr=xx.Norml2();
+      V.SetSize(T.GetDimension());
+      if (cr <= r)
+      {
+         V = force;
+      }
+      else
+      {
+         V = 0.0;
+      }
+   }
+
+   void Set(double r_,Vector & center_, Vector & force_)
+   {
+      r=r_;
+      center = center_;
+      force = force_;
+   }
+};
+
+//  Class for solving Poisson's equation:
+//
+//       - ∇ ⋅(κ ∇ u) = f  in Ω
+//
+class DiffusionSolver
+{
+private:
+   Mesh * mesh = nullptr;
+   int order = 1;
+   // diffusion coefficient
+   Coefficient * diffcf = nullptr;
+   // mass coefficient
+   Coefficient * masscf = nullptr;
+   Coefficient * rhscf = nullptr;
+   Coefficient * essbdr_cf = nullptr;
+   Coefficient * neumann_cf = nullptr;
+   VectorCoefficient * gradient_cf = nullptr;
+
+   // FEM solver
+   int dim;
+   FiniteElementCollection * fec = nullptr;
+   FiniteElementSpace * fes = nullptr;
+   Array<int> ess_bdr;
+   Array<int> neumann_bdr;
+   GridFunction * u = nullptr;
+   LinearForm * b = nullptr;
+   bool parallel;
+#ifdef MFEM_USE_MPI
+   ParMesh * pmesh = nullptr;
+   ParFiniteElementSpace * pfes = nullptr;
+#endif
+
+public:
+   DiffusionSolver() { }
+   DiffusionSolver(Mesh * mesh_, int order_, Coefficient * diffcf_,
+                   Coefficient * cf_);
+
+   void SetMesh(Mesh * mesh_)
+   {
+      mesh = mesh_;
+      parallel = false;
+#ifdef MFEM_USE_MPI
+      pmesh = dynamic_cast<ParMesh *>(mesh);
+      if (pmesh) { parallel = true; }
+#endif
+   }
+   void SetOrder(int order_) { order = order_ ; }
+   void SetDiffusionCoefficient(Coefficient * diffcf_) { diffcf = diffcf_; }
+   void SetMassCoefficient(Coefficient * masscf_) { masscf = masscf_; }
+   void SetRHSCoefficient(Coefficient * rhscf_) { rhscf = rhscf_; }
+   void SetEssentialBoundary(const Array<int> & ess_bdr_) { ess_bdr = ess_bdr_;};
+   void SetNeumannBoundary(const Array<int> & neumann_bdr_) { neumann_bdr = neumann_bdr_;};
+   void SetNeumannData(Coefficient * neumann_cf_) {neumann_cf = neumann_cf_;}
+   void SetEssBdrData(Coefficient * essbdr_cf_) {essbdr_cf = essbdr_cf_;}
+   void SetGradientData(VectorCoefficient * gradient_cf_) {gradient_cf = gradient_cf_;}
+
+   void ResetFEM();
+   void SetupFEM();
+
+   void Solve();
+   GridFunction * GetFEMSolution();
+   LinearForm * GetLinearForm() {return b;}
+#ifdef MFEM_USE_MPI
+   ParGridFunction * GetParFEMSolution();
+   ParLinearForm * GetParLinearForm()
+   {
+      if (parallel)
+      {
+         return dynamic_cast<ParLinearForm *>(b);
+      }
+      else
+      {
+         MFEM_ABORT("Wrong code path. Call GetLinearForm");
+         return nullptr;
+      }
+   }
+#endif
+
+   ~DiffusionSolver();
+
+};
+
+
+//  Class for solving linear elasticity:
+//
+//        -∇ ⋅ σ(u) = f  in Ω  + BCs
+//
+//  where
+//
+//        σ(u) = λ ∇⋅u I + μ (∇ u + ∇uᵀ)
+//
+class LinearElasticitySolver
+{
+private:
+   Mesh * mesh = nullptr;
+   int order = 1;
+   Coefficient * lambda_cf = nullptr;
+   Coefficient * mu_cf = nullptr;
+   VectorCoefficient * essbdr_cf = nullptr;
+   VectorCoefficient * rhs_cf = nullptr;
+
+   // FEM solver
+   int dim;
+   FiniteElementCollection * fec = nullptr;
+   FiniteElementSpace * fes = nullptr;
+   Array<int> ess_bdr;
+   Array<int> neumann_bdr;
+   GridFunction * u = nullptr;
+   LinearForm * b = nullptr;
+   bool parallel;
+#ifdef MFEM_USE_MPI
+   ParMesh * pmesh = nullptr;
+   ParFiniteElementSpace * pfes = nullptr;
+#endif
+
+public:
+   LinearElasticitySolver() { }
+   LinearElasticitySolver(Mesh * mesh_, int order_,
+                          Coefficient * lambda_cf_, Coefficient * mu_cf_);
+
+   void SetMesh(Mesh * mesh_)
+   {
+      mesh = mesh_;
+      parallel = false;
+#ifdef MFEM_USE_MPI
+      pmesh = dynamic_cast<ParMesh *>(mesh);
+      if (pmesh) { parallel = true; }
+#endif
+   }
+   void SetOrder(int order_) { order = order_ ; }
+   void SetLameCoefficients(Coefficient * lambda_cf_, Coefficient * mu_cf_) { lambda_cf = lambda_cf_; mu_cf = mu_cf_;  }
+   void SetRHSCoefficient(VectorCoefficient * rhs_cf_) { rhs_cf = rhs_cf_; }
+   void SetEssentialBoundary(const Array<int> & ess_bdr_) { ess_bdr = ess_bdr_;};
+   void SetNeumannBoundary(const Array<int> & neumann_bdr_) { neumann_bdr = neumann_bdr_;};
+   void SetEssBdrData(VectorCoefficient * essbdr_cf_) {essbdr_cf = essbdr_cf_;}
+
+   void ResetFEM();
+   void SetupFEM();
+
+   void Solve();
+   GridFunction * GetFEMSolution();
+   LinearForm * GetLinearForm() {return b;}
+#ifdef MFEM_USE_MPI
+   ParGridFunction * GetParFEMSolution();
+   ParLinearForm * GetParLinearForm()
+   {
+      if (parallel)
+      {
+         return dynamic_cast<ParLinearForm *>(b);
+      }
+      else
+      {
+         MFEM_ABORT("Wrong code path. Call GetLinearForm");
+         return nullptr;
+      }
+   }
+#endif
+
+   ~LinearElasticitySolver();
+
+};
+
+// -----------------------------------------------------------------------
+// --------------------      Poisson solver     --------------------------
+// -----------------------------------------------------------------------
+
+DiffusionSolver::DiffusionSolver(Mesh * mesh_, int order_,
+                                 Coefficient * diffcf_, Coefficient * rhscf_)
+   : mesh(mesh_), order(order_), diffcf(diffcf_), rhscf(rhscf_)
+{
+
+#ifdef MFEM_USE_MPI
+   pmesh = dynamic_cast<ParMesh *>(mesh);
+   if (pmesh) { parallel = true; }
+#endif
+
+   SetupFEM();
+}
+
+void DiffusionSolver::SetupFEM()
+{
+   dim = mesh->Dimension();
+   fec = new H1_FECollection(order, dim);
+
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      pfes = new ParFiniteElementSpace(pmesh, fec);
+      u = new ParGridFunction(pfes);
+      b = new ParLinearForm(pfes);
+   }
+   else
+   {
+      fes = new FiniteElementSpace(mesh, fec);
+      u = new GridFunction(fes);
+      b = new LinearForm(fes);
+   }
+#else
+   fes = new FiniteElementSpace(mesh, fec);
+   u = new GridFunction(fes);
+   b = new LinearForm(fes);
+#endif
+   *u=0.0;
+
+   if (!ess_bdr.Size())
+   {
+      if (mesh->bdr_attributes.Size())
+      {
+         ess_bdr.SetSize(mesh->bdr_attributes.Max());
+         ess_bdr = 1;
+      }
+   }
+}
+
+void DiffusionSolver::Solve()
+{
+   OperatorPtr A;
+   Vector B, X;
+   Array<int> ess_tdof_list;
+
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      pfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
+   }
+   else
+   {
+      fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
+   }
+#else
+   fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
+#endif
+   *u=0.0;
+   if (b)
+   {
+      delete b;
+#ifdef MFEM_USE_MPI
+      if (parallel)
+      {
+         b = new ParLinearForm(pfes);
+      }
+      else
+      {
+         b = new LinearForm(fes);
+      }
+#else
+      b = new LinearForm(fes);
+#endif
+   }
+   if (rhscf)
+   {
+      b->AddDomainIntegrator(new DomainLFIntegrator(*rhscf));
+   }
+   if (neumann_cf)
+   {
+      MFEM_VERIFY(neumann_bdr.Size(), "neumann_bdr attributes not provided");
+      b->AddBoundaryIntegrator(new BoundaryLFIntegrator(*neumann_cf),neumann_bdr);
+   }
+   else if (gradient_cf)
+   {
+      MFEM_VERIFY(neumann_bdr.Size(), "neumann_bdr attributes not provided");
+      b->AddBoundaryIntegrator(new BoundaryNormalLFIntegrator(*gradient_cf),
+                               neumann_bdr);
+   }
+
+   b->Assemble();
+
+   BilinearForm * a = nullptr;
+
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      a = new ParBilinearForm(pfes);
+   }
+   else
+   {
+      a = new BilinearForm(fes);
+   }
+#else
+   a = new BilinearForm(fes);
+#endif
+   a->AddDomainIntegrator(new DiffusionIntegrator(*diffcf));
+   if (masscf)
+   {
+      a->AddDomainIntegrator(new MassIntegrator(*masscf));
+   }
+   a->Assemble();
+   if (essbdr_cf)
+   {
+      u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
+   }
+   a->FormLinearSystem(ess_tdof_list, *u, *b, A, X, B);
+
+   CGSolver * cg = nullptr;
+   Solver * M = nullptr;
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      M = new HypreBoomerAMG;
+      dynamic_cast<HypreBoomerAMG*>(M)->SetPrintLevel(0);
+      cg = new CGSolver(pmesh->GetComm());
+   }
+   else
+   {
+      M = new GSSmoother((SparseMatrix&)(*A));
+      cg = new CGSolver;
+   }
+#else
+   M = new GSSmoother((SparseMatrix&)(*A));
+   cg = new CGSolver;
+#endif
+   cg->SetRelTol(1e-12);
+   cg->SetMaxIter(10000);
+   cg->SetPrintLevel(0);
+   cg->SetPreconditioner(*M);
+   cg->SetOperator(*A);
+   cg->Mult(B, X);
+   delete M;
+   delete cg;
+   a->RecoverFEMSolution(X, *b, *u);
+   delete a;
+}
+
+GridFunction * DiffusionSolver::GetFEMSolution()
+{
+   return u;
+}
+
+#ifdef MFEM_USE_MPI
+ParGridFunction * DiffusionSolver::GetParFEMSolution()
+{
+   if (parallel)
+   {
+      return dynamic_cast<ParGridFunction*>(u);
+   }
+   else
+   {
+      MFEM_ABORT("Wrong code path. Call GetFEMSolution");
+      return nullptr;
+   }
+}
+#endif
+
+void DiffusionSolver::ResetFEM()
+{
+   delete u; u = nullptr;
+   delete fes; fes = nullptr;
+   delete fec; fec = nullptr;
+   delete b;
+}
+
+
+DiffusionSolver::~DiffusionSolver()
+{
+   ResetFEM();
+}
+
+
+
+// -----------------------------------------------------------------------
+// ------------------      Elasticity solver     -------------------------
+// -----------------------------------------------------------------------
+
+LinearElasticitySolver::LinearElasticitySolver(Mesh * mesh_, int order_,
+                                               Coefficient * lambda_cf_, Coefficient * mu_cf_)
+   : mesh(mesh_), order(order_), lambda_cf(lambda_cf_), mu_cf(mu_cf_)
+{
+#ifdef MFEM_USE_MPI
+   pmesh = dynamic_cast<ParMesh *>(mesh);
+   if (pmesh) { parallel = true; }
+#endif
+   SetupFEM();
+}
+
+void LinearElasticitySolver::SetupFEM()
+{
+   dim = mesh->Dimension();
+   fec = new H1_FECollection(order, dim,BasisType::Positive);
+
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      pfes = new ParFiniteElementSpace(pmesh, fec, dim);
+      u = new ParGridFunction(pfes);
+      b = new ParLinearForm(pfes);
+   }
+   else
+   {
+      fes = new FiniteElementSpace(mesh, fec,dim);
+      u = new GridFunction(fes);
+      b = new LinearForm(fes);
+   }
+#else
+   fes = new FiniteElementSpace(mesh, fec, dim);
+   u = new GridFunction(fes);
+   b = new LinearForm(fes);
+#endif
+   *u=0.0;
+
+   if (!ess_bdr.Size())
+   {
+      if (mesh->bdr_attributes.Size())
+      {
+         ess_bdr.SetSize(mesh->bdr_attributes.Max());
+         ess_bdr = 1;
+      }
+   }
+}
+
+void LinearElasticitySolver::Solve()
+{
+   GridFunction * x = nullptr;
+   OperatorPtr A;
+   Vector B, X;
+   Array<int> ess_tdof_list;
+
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      x = new ParGridFunction(pfes);
+      pfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
+   }
+   else
+   {
+      x = new GridFunction(fes);
+      fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
+   }
+#else
+   x = new GridFunction(fes);
+   fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
+#endif
+   *u=0.0;
+   if (b)
+   {
+      delete b;
+#ifdef MFEM_USE_MPI
+      if (parallel)
+      {
+         b = new ParLinearForm(pfes);
+      }
+      else
+      {
+         b = new LinearForm(fes);
+      }
+#else
+      b = new LinearForm(fes);
+#endif
+   }
+   if (rhs_cf)
+   {
+      b->AddDomainIntegrator(new VectorDomainLFIntegrator(*rhs_cf));
+   }
+
+   b->Assemble();
+
+   *x = 0.0;
+
+   BilinearForm * a = nullptr;
+
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      a = new ParBilinearForm(pfes);
+   }
+   else
+   {
+      a = new BilinearForm(fes);
+   }
+#else
+   a = new BilinearForm(fes);
+#endif
+   a->AddDomainIntegrator(new ElasticityIntegrator(*lambda_cf, *mu_cf));
+   a->Assemble();
+   if (essbdr_cf)
+   {
+      u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
+   }
+   a->FormLinearSystem(ess_tdof_list, *x, *b, A, X, B);
+
+   CGSolver * cg = nullptr;
+   Solver * M = nullptr;
+#ifdef MFEM_USE_MPI
+   if (parallel)
+   {
+      M = new HypreBoomerAMG;
+      dynamic_cast<HypreBoomerAMG*>(M)->SetPrintLevel(0);
+      cg = new CGSolver(pmesh->GetComm());
+   }
+   else
+   {
+      M = new GSSmoother((SparseMatrix&)(*A));
+      cg = new CGSolver;
+   }
+#else
+   M = new GSSmoother((SparseMatrix&)(*A));
+   cg = new CGSolver;
+#endif
+   cg->SetRelTol(1e-10);
+   cg->SetMaxIter(10000);
+   cg->SetPrintLevel(0);
+   cg->SetPreconditioner(*M);
+   cg->SetOperator(*A);
+   cg->Mult(B, X);
+   delete M;
+   delete cg;
+   a->RecoverFEMSolution(X, *b, *x);
+   *u+=*x;
+   delete a;
+   delete x;
+}
+
+GridFunction * LinearElasticitySolver::GetFEMSolution()
+{
+   return u;
+}
+
+#ifdef MFEM_USE_MPI
+ParGridFunction * LinearElasticitySolver::GetParFEMSolution()
+{
+   if (parallel)
+   {
+      return dynamic_cast<ParGridFunction*>(u);
+   }
+   else
+   {
+      MFEM_ABORT("Wrong code path. Call GetFEMSolution");
+      return nullptr;
+   }
+}
+#endif
+
+void LinearElasticitySolver::ResetFEM()
+{
+   delete u; u = nullptr;
+   delete fes; fes = nullptr;
+   delete fec; fec = nullptr;
+   delete b;
+}
+
+LinearElasticitySolver::~LinearElasticitySolver()
+{
+   ResetFEM();
+}

examples/ex35p.cpp

@@ -0,0 +1,487 @@
+//                            MFEM Example 35 - Parallel Version
+//
+//
+// Compile with: make ex35p
+//
+// Sample runs:
+// mpirun -np 4 ex35p -lambda 0.1 -mu 0.1
+// mpirun -np 6 ex35p -r 5 -o 2 -alpha 5.0 -epsilon 0.01 -mi 50 -mf 0.5 -tol 1e-5
+// mpirun -np 8 ex35p -r 6 -o 2 -alpha 10.0 -epsilon 0.02 -mi 50 -mf 0.5 -tol 1e-5
+// mpirun -np 8 ex35p -r 7 -o 1 -alpha 10.0 -epsilon 0.01 -mi 50 -mf 0.5 -tol 1e-5
+//
+//
+// Description: This example code demonstrates the use of MFEM to solve a
+//              density-filtered [3] topology optimization problem. The
+//              objective is to minimize the compliance
+//
+//                  minimize ∫_Ω f⋅u dx over u ∈ [H¹(Ω)]² and ρ ∈ L²(Ω)
+//
+//                  subject to
+//
+//                    -Div(r(ρ̃)Cε(u)) = f       in Ω + BCs
+//                    -ϵ²Δρ̃ + ρ̃ = ρ             in Ω + Neumann BCs
+//                    0 ≤ ρ ≤ 1                 in Ω
+//                    ∫_Ω ρ dx = θ vol(Ω)
+//
+//              Here, r(ρ̃) = ρ₀ + ρ̃³ (1-ρ₀) is the solid isotropic material
+//              penalization (SIMP) law, C is the elasticity tensor for an
+//              isotropic linearly elastic material, ϵ > 0 is the design
+//              length scale, and 0 < θ < 1 is the volume fraction.
+//
+//              The problem is discretized and gradients are computing using
+//              finite elements [1]. The design is optimized using an entropic
+//              mirror descent algorithm introduced by Keith and Surowiec [2]
+//              that is tailored to the bound constraint 0 ≤ ρ ≤ 1.
+//
+//              This example highlights the ability of MFEM to deliver high-
+//              order solutions to inverse design problems and showcases how
+//              to set up and solve PDE-constrained optimization problems
+//              using the so-called reduced space approach.
+//
+//
+// [1] Andreassen, E., Clausen, A., Schevenels, M., Lazarov, B. S., & Sigmund, O.
+//    (2011). Efficient topology optimization in MATLAB using 88 lines of
+//    code. Structural and Multidisciplinary Optimization, 43(1), 1-16.
+// [2] Keith, B. and Surowiec, T. (2023) The entropic finite element method
+//     (in preparation).
+// [3] Lazarov, B. S., & Sigmund, O. (2011). Filters in topology optimization
+//     based on Helmholtz‐type differential equations. International Journal
+//     for Numerical Methods in Engineering, 86(6), 765-781.
+
+#include "mfem.hpp"
+#include <iostream>
+#include <fstream>
+#include "ex35.hpp"
+
+/**
+ * @brief Nonlinear projection of 0 < τ < 1 onto the subspace
+ *        ∫_Ω τ dx = θ vol(Ω) as follows.
+ *
+ *        1. Compute the root of the R → R function
+ *            f(c) = ∫_Ω expit(lnit(τ) + c) dx - θ vol(Ω)
+ *        2. Set τ ← expit(lnit(τ) + c).
+ *
+ */
+void projit(GridFunction &tau, double &c, LinearForm &vol_form,
+            double volume_fraction, double tol=1e-12, int max_its=10)
+{
+   GridFunction ftmp(tau.FESpace());
+   GridFunction dftmp(tau.FESpace());
+   for (int k=0; k<max_its; k++)
+   {
+      // Compute f(c) and dfdc(c)
+      for (int i=0; i<tau.Size(); i++)
+      {
+         ftmp[i]  = expit(lnit(tau[i]) + c) - volume_fraction;
+         dftmp[i] = dexpitdx(lnit(tau[i]) + c);
+      }
+      double f = vol_form(ftmp);
+      double df = vol_form(dftmp);
+
+      MPI_Allreduce(MPI_IN_PLACE,&f,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+      MPI_Allreduce(MPI_IN_PLACE,&df,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+
+      double dc = -f/df;
+      c += dc;
+      if (abs(dc) < tol) { break; }
+   }
+   tau = ftmp;
+   tau += volume_fraction;
+}
+
+using namespace std;
+using namespace mfem;
+/**
+ * ---------------------------------------------------------------
+ *                      ALGORITHM PREAMBLE
+ * ---------------------------------------------------------------
+ *
+ *  The Lagrangian for this problem is
+ *
+ *          L(u,ρ,ρ̃,w,w̃) = (f,u) - (r(ρ̃) C ε(u),ε(w)) + (f,w)
+ *                       - (ϵ² ∇ρ̃,∇w̃) - (ρ̃,w̃) + (ρ,w̃)
+ *
+ *  where
+ *
+ *    r(ρ̃) = ρ₀ + ρ̃³ (1 - ρ₀)       (SIMP rule)
+ *
+ *    ε(u) = (∇u + ∇uᵀ)/2           (symmetric gradient)
+ *
+ *    C e = λtr(e)I + 2μe           (isotropic material)
+ *
+ *  NOTE: The Lame parameters can be computed from Young's modulus E
+ *        and Poisson's ratio ν as follows:
+ *
+ *             λ = E ν/((1+ν)(1-2ν)),      μ = E/(2(1+ν))
+ *
+ * ---------------------------------------------------------------
+ *
+ *  Discretization choices:
+ *
+ *     u ∈ V ⊂ (H¹)ᵈ (order p)
+ *     ρ ∈ L² (order p - 1)
+ *     ρ̃ ∈ H¹ (order p - 1)
+ *     w ∈ V  (order p)
+ *     w̃ ∈ H¹ (order p - 1)
+ *
+ * ---------------------------------------------------------------
+ *                          ALGORITHM
+ * ---------------------------------------------------------------
+ *
+ *  Update ρ with projected mirror descent via the following algorithm.
+ *
+ *  1. Initialize density field 0 < ρ(x) < 1.
+ *
+ *  While not converged:
+ *
+ *     2. Solve filter equation ∂_w̃ L = 0; i.e.,
+ *
+ *           (ϵ² ∇ ρ̃, ∇ v ) + (ρ̃,v) = (ρ,v)   ∀ v ∈ H¹.
+ *
+ *     3. Solve primal problem ∂_w L = 0; i.e.,
+ *
+ *      (λ(ρ̃) ∇⋅u, ∇⋅v) + (2 μ(ρ̃) ε(u), ε(v)) = (f,v)   ∀ v ∈ V,
+ *
+ *     where λ(ρ̃) := λ r(ρ̃) and  μ(ρ̃) := μ r(ρ̃).
+ *
+ *     NB. The dual problem ∂_u L = 0 is the same as the primal problem due to symmetry.
+ *
+ *     4. Solve for filtered gradient ∂_ρ̃ L = 0; i.e.,
+ *
+ *      (ϵ² ∇ w̃ , ∇ v ) + (w̃ ,v) = (-r'(ρ̃) ( λ(ρ̃) |∇⋅u|² + 2 μ(ρ̃) |ε(u)|²),v)   ∀ v ∈ H¹.
+ *
+ *     5. Construct gradient G ∈ L²; i.e.,
+ *
+ *                         (G,v) = (w̃,v)   ∀ v ∈ L².
+ *
+ *     6. Mirror descent update until convergence; i.e.,
+ *
+ *                      ρ ← projit(expit(linit(ρ) - αG)),
+ *
+ *     where
+ *
+ *          α > 0                            (step size parameter)
+ *
+ *          expit(x) = eˣ/(1+eˣ)             (sigmoid)
+ *
+ *          linit(y) = ln(y) - ln(1-y)       (inverse of sigmoid)
+ *
+ *     and projit is a (compatible) projection operator enforcing ∫_Ω ρ dx = θ vol(Ω).
+ *
+ *  end
+ *
+ */
+
+int main(int argc, char *argv[])
+{
+   // 0. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 1. Parse command-line options.
+   int ref_levels = 4;
+   int order = 2;
+   bool visualization = true;
+   double alpha = 1.0;
+   double epsilon = 0.01;
+   double mass_fraction = 0.5;
+   int max_it = 1e2;
+   double tol = 1e-4;
+   double rho_min = 1e-6;
+   double lambda = 1.0;
+   double mu = 1.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&alpha, "-alpha", "--alpha-step-length",
+                  "Step length for gradient descent.");
+   args.AddOption(&epsilon, "-epsilon", "--epsilon-thickness",
+                  "epsilon phase field thickness");
+   args.AddOption(&max_it, "-mi", "--max-it",
+                  "Maximum number of gradient descent iterations.");
+   args.AddOption(&tol, "-tol", "--tol",
+                  "Exit tolerance for ρ ");
+   args.AddOption(&mass_fraction, "-mf", "--mass-fraction",
+                  "Mass fraction for diffusion coefficient.");
+   args.AddOption(&lambda, "-lambda", "--lambda",
+                  "Lame constant λ");
+   args.AddOption(&mu, "-mu", "--mu",
+                  "Lame constant μ");
+   args.AddOption(&rho_min, "-rmin", "--rho-min",
+                  "Minimum of density coefficient.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   Mesh mesh = Mesh::MakeCartesian2D(3,1,mfem::Element::Type::QUADRILATERAL,true,
+                                     3.0,1.0);
+
+   int dim = mesh.Dimension();
+
+   // 2. Set BCs.
+   for (int i = 0; i<mesh.GetNBE(); i++)
+   {
+      Element * be = mesh.GetBdrElement(i);
+      Array<int> vertices;
+      be->GetVertices(vertices);
+
+      double * coords1 = mesh.GetVertex(vertices[0]);
+      double * coords2 = mesh.GetVertex(vertices[1]);
+
+      Vector center(2);
+      center(0) = 0.5*(coords1[0] + coords2[0]);
+      center(1) = 0.5*(coords1[1] + coords2[1]);
+
+      if (abs(center(0) - 0.0) < 1e-10)
+      {
+         // the left edge
+         be->SetAttribute(1);
+      }
+      else
+      {
+         // all other boundaries
+         be->SetAttribute(2);
+      }
+   }
+   mesh.SetAttributes();
+
+   // 3. Refine the mesh.
+   for (int lev = 0; lev < ref_levels; lev++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   // 4. Define the necessary finite element spaces on the mesh.
+   H1_FECollection state_fec(order, dim); // space for u
+   H1_FECollection filter_fec(order, dim); // space for ρ̃
+   L2_FECollection control_fec(order-1, dim,
+                               BasisType::Positive); // space for ρ
+   ParFiniteElementSpace state_fes(&pmesh, &state_fec,dim);
+   ParFiniteElementSpace filter_fes(&pmesh, &filter_fec);
+   ParFiniteElementSpace control_fes(&pmesh, &control_fec);
+
+   HYPRE_BigInt state_size = state_fes.GlobalTrueVSize();
+   HYPRE_BigInt control_size = control_fes.GlobalTrueVSize();
+   HYPRE_BigInt filter_size = filter_fes.GlobalTrueVSize();
+   if (myid==0)
+   {
+      cout << "Number of state unknowns: " << state_size << endl;
+      cout << "Number of filter unknowns: " << filter_size << endl;
+      cout << "Number of control unknowns: " << control_size << endl;
+   }
+
+   // 5. Set the initial guess for ρ.
+   ParGridFunction u(&state_fes);
+   ParGridFunction rho(&control_fes);
+   ParGridFunction rho_old(&control_fes);
+   ParGridFunction rho_filter(&filter_fes);
+   u = 0.0;
+   rho_filter = 0.0;
+   rho = 0.5;
+   rho_old = 0.5;
+
+   // 6. Set-up the physics solver.
+   int maxat = pmesh.bdr_attributes.Max();
+   Array<int> ess_bdr(maxat);
+   ess_bdr = 0;
+   ess_bdr[0] = 1;
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient lambda_cf(lambda);
+   ConstantCoefficient mu_cf(mu);
+   LinearElasticitySolver * ElasticitySolver = new LinearElasticitySolver();
+   ElasticitySolver->SetMesh(&pmesh);
+   ElasticitySolver->SetOrder(state_fec.GetOrder());
+   ElasticitySolver->SetupFEM();
+   Vector center(2); center(0) = 2.9; center(1) = 0.5;
+   Vector force(2); force(0) = 0.0; force(1) = -1.0;
+   double r = 0.05;
+   VolumeForceCoefficient vforce_cf(r,center,force);
+   ElasticitySolver->SetRHSCoefficient(&vforce_cf);
+   ElasticitySolver->SetEssentialBoundary(ess_bdr);
+
+   // 7. Set-up the filter solver.
+   ConstantCoefficient eps2_cf(epsilon*epsilon);
+   DiffusionSolver * FilterSolver = new DiffusionSolver();
+   FilterSolver->SetMesh(&pmesh);
+   FilterSolver->SetOrder(filter_fec.GetOrder());
+   FilterSolver->SetDiffusionCoefficient(&eps2_cf);
+   FilterSolver->SetMassCoefficient(&one);
+   Array<int> ess_bdr_filter;
+   if (pmesh.bdr_attributes.Size())
+   {
+      ess_bdr_filter.SetSize(pmesh.bdr_attributes.Max());
+      ess_bdr_filter = 0;
+   }
+   FilterSolver->SetEssentialBoundary(ess_bdr_filter);
+   FilterSolver->SetupFEM();
+
+   ParBilinearForm mass(&control_fes);
+   mass.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator(one)));
+   mass.Assemble();
+   HypreParMatrix M;
+   Array<int> empty;
+   mass.FormSystemMatrix(empty,M);
+
+   // 8. Define the Lagrange multiplier and gradient functions
+   ParGridFunction grad(&control_fes);
+   ParGridFunction w_filter(&filter_fes);
+
+   // 9. Define some tools for later
+   ConstantCoefficient zero(0.0);
+   ParGridFunction onegf(&control_fes);
+   onegf = 1.0;
+   ParLinearForm vol_form(&control_fes);
+   vol_form.AddDomainIntegrator(new DomainLFIntegrator(one));
+   vol_form.Assemble();
+   double domain_volume = vol_form(onegf);
+
+   // 10. Connect to GLVis. Prepare for VisIt output.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sout_u,sout_r,sout_rho;
+   if (visualization)
+   {
+      sout_u.open(vishost, visport);
+      sout_rho.open(vishost, visport);
+      sout_r.open(vishost, visport);
+      sout_u.precision(8);
+      sout_rho.precision(8);
+      sout_r.precision(8);
+   }
+
+   mfem::ParaViewDataCollection paraview_dc("Elastic_compliance", &pmesh);
+   paraview_dc.SetPrefixPath("ParaView");
+   paraview_dc.SetLevelsOfDetail(order);
+   paraview_dc.SetCycle(0);
+   paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   paraview_dc.SetHighOrderOutput(true);
+   paraview_dc.SetTime(0.0);
+   paraview_dc.RegisterField("displacement",&u);
+   paraview_dc.RegisterField("density",&rho);
+   paraview_dc.RegisterField("filtered_density",&rho_filter);
+
+   // 11. Iterate
+   int step = 0;
+   double c0 = 0.0;
+   for (int k = 1; k < max_it; k++)
+   {
+      if (k > 1) { alpha *= ((double) k) / ((double) k-1); }
+      step++;
+
+      if (myid == 0)
+      {
+         cout << "\nStep = " << k << endl;
+      }
+
+      // Step 1 - Filter solve
+      // Solve (ϵ^2 ∇ ρ̃, ∇ v ) + (ρ̃,v) = (ρ,v)
+      GridFunctionCoefficient rho_cf(&rho);
+      FilterSolver->SetRHSCoefficient(&rho_cf);
+      FilterSolver->Solve();
+      rho_filter = *FilterSolver->GetFEMSolution();
+
+      // Step 2 - State solve
+      // Solve (λ(ρ̃) ∇⋅u, ∇⋅v) + (2 μ(ρ̃) ε(u), ε(v)) = (f,v)
+      SIMPInterpolationCoefficient SIMP_cf(&rho_filter,rho_min, 1.0);
+      ProductCoefficient lambda_SIMP_cf(lambda_cf,SIMP_cf);
+      ProductCoefficient mu_SIMP_cf(mu_cf,SIMP_cf);
+      ElasticitySolver->SetLameCoefficients(&lambda_SIMP_cf,&mu_SIMP_cf);
+      ElasticitySolver->Solve();
+      u = *ElasticitySolver->GetFEMSolution();
+
+      // Step 3 - Adjoint filter solve
+      // Solve (ϵ² ∇ w̃, ∇ v) + (w̃ ,v) = (-r'(ρ̃) ( λ(ρ̃) |∇⋅u|² + 2 μ(ρ̃) |ε(u)|²),v)
+      StrainEnergyDensityCoefficient rhs_cf(&lambda_cf,&mu_cf,&u, &rho_filter,
+                                            rho_min);
+      FilterSolver->SetRHSCoefficient(&rhs_cf);
+      FilterSolver->Solve();
+      w_filter = *FilterSolver->GetFEMSolution();
+
+      // Step 4 - Compute gradient
+      // Solve G = M⁻¹w̃
+      GridFunctionCoefficient w_cf(&w_filter);
+      ParLinearForm w_rhs(&control_fes);
+      w_rhs.AddDomainIntegrator(new DomainLFIntegrator(w_cf));
+      w_rhs.Assemble();
+      M.Mult(w_rhs,grad);
+
+      // Step 5 - Update design variable ρ ← projit(expit(linit(ρ) - αG))
+      // Note: The update here is performed on the coefficients of the linear
+      //       representation of the design variable in the Berstein basis. It would
+      //       be more mathematically sound to update the design variable so that
+      //       the values at the integration points follow this update rule.
+      for (int i = 0; i < rho.Size(); i++)
+      {
+         rho[i] = expit(lnit(rho[i]) - alpha*grad[i]);
+      }
+      projit(rho, c0, vol_form, mass_fraction);
+
+      GridFunctionCoefficient tmp(&rho_old);
+      double norm_reduced_gradient = rho.ComputeL2Error(tmp)/alpha;
+      rho_old = rho;
+
+      double compliance = (*(ElasticitySolver->GetLinearForm()))(u);
+      MPI_Allreduce(MPI_IN_PLACE,&compliance,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+      double material_volume = vol_form(rho);
+      if (myid == 0)
+      {
+         mfem::out << "norm of reduced gradient = " << norm_reduced_gradient << endl;
+         mfem::out << "compliance = " << compliance << endl;
+         mfem::out << "mass_fraction = " << material_volume / domain_volume << endl;
+      }
+
+      if (visualization)
+      {
+         sout_u << "parallel " << num_procs << " " << myid << "\n";
+         sout_u << "solution\n" << pmesh << u
+                << "window_title 'Displacement u'" << flush;
+
+         sout_rho << "parallel " << num_procs << " " << myid << "\n";
+         sout_rho << "solution\n" << pmesh << rho
+                  << "window_title 'Control variable ρ'" << flush;
+
+         ParGridFunction r_gf(&filter_fes);
+         r_gf.ProjectCoefficient(SIMP_cf);
+         sout_r << "parallel " << num_procs << " " << myid << "\n";
+         sout_r << "solution\n" << pmesh << r_gf
+                << "window_title 'Design density r(ρ̃)'" << flush;
+
+         paraview_dc.SetCycle(k);
+         paraview_dc.SetTime((double)k);
+         paraview_dc.Save();
+      }
+
+      if (norm_reduced_gradient < tol)
+      {
+         break;
+      }
+   }
+
+   delete ElasticitySolver;
+   delete FilterSolver;
+
+   return 0;
+}

examples/ex36.cpp

@@ -0,0 +1,540 @@
+//                              MFEM Example 35
+//
+//
+// Compile with: make ex35
+//
+// Sample runs:
+//     ex35 -alpha 10
+//     ex35 -lambda 0.1 -mu 0.1
+//     ex35 -r 5 -o 2 -alpha 5.0 -epsilon 0.01 -mi 50 -mf 0.5 -tol 1e-5
+//     ex35 -r 6 -o 1 -alpha 10.0 -epsilon 0.01 -mi 50 -mf 0.5 -tol 1e-5
+//
+//
+// Description: This example code demonstrates the use of MFEM to solve a
+//              density-filtered [3] topology optimization problem. The
+//              objective is to minimize the compliance
+//
+//                  minimize ∫_Ω f⋅u dx over u ∈ [H¹(Ω)]² and ρ ∈ L²(Ω)
+//
+//                  subject to
+//
+//                    -Div(r(ρ̃)Cε(u)) = f       in Ω + BCs
+//                    -ϵ²Δρ̃ + ρ̃ = ρ             in Ω + Neumann BCs
+//                    0 ≤ ρ ≤ 1                 in Ω
+//                    ∫_Ω ρ dx = θ vol(Ω)
+//
+//              Here, r(ρ̃) = ρ₀ + ρ̃³ (1-ρ₀) is the solid isotropic material
+//              penalization (SIMP) law, C is the elasticity tensor for an
+//              isotropic linearly elastic material, ϵ > 0 is the design
+//              length scale, and 0 < θ < 1 is the volume fraction.
+//
+//              The problem is discretized and gradients are computing using
+//              finite elements [1]. The design is optimized using an entropic
+//              mirror descent algorithm introduced by Keith and Surowiec [2]
+//              that is tailored to the bound constraint 0 ≤ ρ ≤ 1.
+//
+//              This example highlights the ability of MFEM to deliver high-
+//              order solutions to inverse design problems and showcases how
+//              to set up and solve PDE-constrained optimization problems
+//              using the so-called reduced space approach.
+//
+//
+// [1] Andreassen, E., Clausen, A., Schevenels, M., Lazarov, B. S., & Sigmund, O.
+//    (2011). Efficient topology optimization in MATLAB using 88 lines of
+//    code. Structural and Multidisciplinary Optimization, 43(1), 1-16.
+// [2] Keith, B. and Surowiec, T. (2023) The entropic finite element method
+//     (in preparation).
+// [3] Lazarov, B. S., & Sigmund, O. (2011). Filters in topology optimization
+//     based on Helmholtz‐type differential equations. International Journal
+//     for Numerical Methods in Engineering, 86(6), 765-781.
+
+#include "mfem.hpp"
+#include <iostream>
+#include <fstream>
+#include "ex35.hpp"
+
+/**
+ * @brief Nonlinear projection of 0 < τ < 1 onto the subspace
+ *        ∫_Ω τ dx = θ vol(Ω) as follows.
+ *
+ *        1. Compute the root of the R → R function
+ *            f(c) = ∫_Ω expit(lnit(τ) + c) dx - θ vol(Ω)
+ *        2. Set τ ← expit(lnit(τ) + c).
+ *
+ */
+// void projit(GridFunction &tau, double &c, LinearForm &vol_form,
+//             double volume_fraction, double tol=1e-12, int max_its=10)
+// {
+//    GridFunction ftmp(tau.FESpace());
+//    GridFunction dftmp(tau.FESpace());
+//    for (int k=0; k<max_its; k++)
+//    {
+//       // Compute f(c) and dfdc(c)
+//       for (int i=0; i<tau.Size(); i++)
+//       {
+//          ftmp[i]  = expit(lnit(tau[i]) + c) - volume_fraction;
+//          dftmp[i] = dexpitdx(lnit(tau[i]) + c);
+//       }
+//       double f = vol_form(ftmp);
+//       double df = vol_form(dftmp);
+
+//       double dc = -f/df;
+//       c += dc;
+//       if (abs(dc) < tol) { break; }
+//    }
+//    tau = ftmp;
+//    tau += volume_fraction;
+// }
+class MappedGridFunctionCoefficient : public GridFunctionCoefficient
+{
+private:
+   std::__1::function<double(const double)> fun;
+public:
+   MappedGridFunctionCoefficient(GridFunction *gf,
+                                 std::__1::function<double(const double)> fun_):GridFunctionCoefficient(gf),
+      fun(fun_) {}
+
+   /// Evaluate the coefficient at @a ip.
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+   {
+      return fun(GridFunctionCoefficient::Eval(T, ip));
+   }
+};
+class GridFunctionPlusCoefficient : public GridFunctionCoefficient
+{
+private:
+   Coefficient *cf;
+public:
+   GridFunctionPlusCoefficient(GridFunction *gf,
+                               Coefficient *coeff): GridFunctionCoefficient(gf), cf(coeff) {}
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      const double gf_val = GridFunctionCoefficient::Eval(T, ip);
+      const double cf_val = cf->Eval(T, ip);
+      return gf_val + cf_val;
+   }
+};
+
+/**
+ * @brief Volume compliant projection operator for ψ
+ * 
+ * @param psi ψ, where ρ = sigmoid(ψ) 
+ * @param target_volume Target ∫ρ
+ * @param tol Tolerance for Newton method
+ * @param max_its Maximum iteration for Newton method
+ */
+void projit(GridFunction &psi, double target_volume,
+            const double tol=1e-12, const int max_its=30)
+{  
+   // ρ = sigmoid(ψ).
+   MappedGridFunctionCoefficient rho(&psi, [](const double x) {return expit(x);});
+   LinearForm rho_form(psi.FESpace()); // ∫ ρ
+   rho_form.AddDomainIntegrator(new DomainLFIntegrator(rho));
+
+   // dsigmoid = d(sigmoid(ψ))/dψ
+   MappedGridFunctionCoefficient dsigmoid(&psi, [](const double x) {return dexpitdx(x);});
+   LinearForm dsigmoid_form(psi.FESpace()); // ∫ dsigmoid(ψ)
+   dsigmoid_form.AddDomainIntegrator(new DomainLFIntegrator(dsigmoid));
+
+   // Newton method with respect to f(c) = ∫sigmoid(ψ + c) - θ|Ω|
+   //
+   //   ψ_new = ψ_old - f(ψ_old) / f'(ψ_old)
+   //         = ψ_old - (∫sigmoid(ψ) - θ|Ω|) / ∫dsigmoid(ψ)
+   for (int i=0; i<max_its; i++)
+   {
+      // Compute ∫ρ with updated ψ
+      rho_form.Assemble();
+      const double f = rho_form.Sum() - target_volume; // ∫sigmoid(ψ) - θ|Ω|
+      // Compute ∫dsigmoid(ψ) with updated ψ
+      dsigmoid_form.Assemble();
+      const double df = dsigmoid_form.Sum(); // ∫dsigmoid(ψ)
+      // Newton increment
+      const double dc = - f / df;
+
+      // For debugging, put assert here.
+      MFEM_ASSERT(isfinite(dc), "Newton increment is not finite.");
+      psi += dc;
+      // tolerance check
+      if (abs(dc) < tol)
+      {
+         return;
+      }
+   }
+   // If you are here, then Newton method failed to converge
+   mfem_error("Projection failed to converge");
+}
+
+using namespace std;
+using namespace mfem;
+/**
+ * ---------------------------------------------------------------
+ *                      ALGORITHM PREAMBLE
+ * ---------------------------------------------------------------
+ *
+ *  The Lagrangian for this problem is
+ *
+ *          L(u,ρ,ρ̃,w,w̃) = (f,u) - (r(ρ̃) C ε(u),ε(w)) + (f,w)
+ *                       - (ϵ² ∇ρ̃,∇w̃) - (ρ̃,w̃) + (ρ,w̃)
+ *
+ *  where
+ *
+ *    r(ρ̃) = ρ₀ + ρ̃³ (1 - ρ₀)       (SIMP rule)
+ *
+ *    ε(u) = (∇u + ∇uᵀ)/2           (symmetric gradient)
+ *
+ *    C e = λtr(e)I + 2μe           (isotropic material)
+ *
+ *  NOTE: The Lame parameters can be computed from Young's modulus E
+ *        and Poisson's ratio ν as follows:
+ *
+ *             λ = E ν/((1+ν)(1-2ν)),      μ = E/(2(1+ν))
+ *
+ * ---------------------------------------------------------------
+ *
+ *  Discretization choices:
+ *
+ *     u ∈ V ⊂ (H¹)ᵈ (order p)
+ *     ρ ∈ L² (order p - 1)
+ *     ρ̃ ∈ H¹ (order p - 1)
+ *     w ∈ V  (order p)
+ *     w̃ ∈ H¹ (order p - 1)
+ *
+ * ---------------------------------------------------------------
+ *                          ALGORITHM
+ * ---------------------------------------------------------------
+ *
+ *  Update ρ with projected mirror descent via the following algorithm.
+ *
+ *  1. Initialize density field 0 < ρ(x) < 1.
+ *
+ *  While not converged:
+ *
+ *     2. Solve filter equation ∂_w̃ L = 0; i.e.,
+ *
+ *           (ϵ² ∇ ρ̃, ∇ v ) + (ρ̃,v) = (ρ,v)   ∀ v ∈ H¹.
+ *
+ *     3. Solve primal problem ∂_w L = 0; i.e.,
+ *
+ *      (λ(ρ̃) ∇⋅u, ∇⋅v) + (2 μ(ρ̃) ε(u), ε(v)) = (f,v)   ∀ v ∈ V,
+ *
+ *     where λ(ρ̃) := λ r(ρ̃) and  μ(ρ̃) := μ r(ρ̃).
+ *
+ *     NB. The dual problem ∂_u L = 0 is the same as the primal problem due to symmetry.
+ *
+ *     4. Solve for filtered gradient ∂_ρ̃ L = 0; i.e.,
+ *
+ *      (ϵ² ∇ w̃ , ∇ v ) + (w̃ ,v) = (-r'(ρ̃) ( λ(ρ̃) |∇⋅u|² + 2 μ(ρ̃) |ε(u)|²),v)   ∀ v ∈ H¹.
+ *
+ *     5. Construct gradient G ∈ L²; i.e.,
+ *
+ *                         (G,v) = (w̃,v)   ∀ v ∈ L².
+ *
+ *     6. Mirror descent update until convergence; i.e.,
+ *
+ *                      ρ ← projit(expit(linit(ρ) - αG)),
+ *
+ *     where
+ *
+ *          α > 0                            (step size parameter)
+ *
+ *          expit(x) = eˣ/(1+eˣ)             (sigmoid)
+ *
+ *          linit(y) = ln(y) - ln(1-y)       (inverse of sigmoid)
+ *
+ *     and projit is a (compatible) projection operator enforcing ∫_Ω ρ dx = θ vol(Ω).
+ *
+ *  end
+ *
+ */
+
+int main(int argc, char *argv[])
+{
+
+   // 1. Parse command-line options.
+   int ref_levels = 5;
+   int order = 2;
+   bool visualization = true;
+   double alpha0 = 1.0;
+   double epsilon = 0.01;
+   double mass_fraction = 0.3;
+   int max_it = 1e3;
+   double tol = 1e-4;
+   double rho_min = 1e-3;
+   double exponent = 3;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&alpha0, "-alpha", "--alpha-step-length",
+                  "Step length for gradient descent.");
+   args.AddOption(&epsilon, "-epsilon", "--epsilon-thickness",
+                  "epsilon phase field thickness");
+   args.AddOption(&max_it, "-mi", "--max-it",
+                  "Maximum number of gradient descent iterations.");
+   args.AddOption(&tol, "-tol", "--tol",
+                  "Exit tolerance for ρ ");
+   args.AddOption(&mass_fraction, "-mf", "--mass-fraction",
+                  "Mass fraction for diffusion coefficient.");
+   args.AddOption(&rho_min, "-rmin", "--rho-min",
+                  "Minimum of density coefficient.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   Mesh mesh = Mesh::MakeCartesian2D(1,1,mfem::Element::Type::QUADRILATERAL,true,
+                                     1.0,1.0);
+
+   int dim = mesh.Dimension();
+
+   // 2. Set BCs.
+   for (int i = 0; i<mesh.GetNBE(); i++)
+   {
+      Element * be = mesh.GetBdrElement(i);
+      Array<int> vertices;
+      be->GetVertices(vertices);
+
+      double * coords1 = mesh.GetVertex(vertices[0]);
+      double * coords2 = mesh.GetVertex(vertices[1]);
+
+      Vector center(2);
+      center(0) = 0.5*(coords1[0] + coords2[0]);
+      center(1) = 0.5*(coords1[1] + coords2[1]);
+
+      if (abs(center(0) - 0.0) < 1e-10)
+      {
+         // the left edge
+         be->SetAttribute(1);
+      }
+      else
+      {
+         // all other boundaries
+         be->SetAttribute(2);
+      }
+   }
+   mesh.SetAttributes();
+
+   // 3. Refine the mesh.
+   for (int lev = 0; lev < ref_levels; lev++)
+   {
+      mesh.UniformRefinement();
+   }
+   ConstantCoefficient zero(0.0);
+   ConstantCoefficient one(1.0);
+
+   // 4. Define the necessary finite element spaces on the mesh.
+   H1_FECollection state_fec(order, dim); // space for u
+   H1_FECollection filter_fec(order, dim); // space for ρ̃
+   L2_FECollection control_fec(order-1, dim,
+                               BasisType::GaussLobatto); // space for ρ
+   FiniteElementSpace state_fes(&mesh, &state_fec);
+   FiniteElementSpace filter_fes(&mesh, &filter_fec);
+   FiniteElementSpace control_fes(&mesh, &control_fec);
+
+   int state_size = state_fes.GetTrueVSize();
+   int control_size = control_fes.GetTrueVSize();
+   int filter_size = filter_fes.GetTrueVSize();
+   cout << "Number of state unknowns: " << state_size << endl;
+   cout << "Number of filter unknowns: " << filter_size << endl;
+   cout << "Number of control unknowns: " << control_size << endl;
+
+   // 5. Set the initial guess for ρ.
+   GridFunction u(&state_fes);
+   GridFunction psi(&control_fes);
+   GridFunction rho_filter_old(&filter_fes);
+   GridFunction rho_filter(&filter_fes);
+   u = 0.0;
+   rho_filter = mass_fraction;
+   psi = lnit(mass_fraction);
+   rho_filter_old = mass_fraction;
+
+   MappedGridFunctionCoefficient rho(&psi, [](const double x) {return expit(x);});
+
+   // 6. Set-up the physics solver.
+   int maxat = mesh.bdr_attributes.Max();
+   Array<int> ess_bdr(maxat);
+   ess_bdr = 0;
+   ess_bdr[0] = 1;
+   DiffusionSolver * diffusionSolver = new DiffusionSolver();
+   diffusionSolver->SetMesh(&mesh);
+   diffusionSolver->SetOrder(state_fec.GetOrder());
+   diffusionSolver->SetRHSCoefficient(&one);
+   diffusionSolver->SetDiffusionCoefficient(&one);
+   diffusionSolver->SetMassCoefficient(&one);
+   diffusionSolver->SetEssentialBoundary(ess_bdr);
+   diffusionSolver->SetupFEM();
+
+   // 7. Set-up the filter solver.
+   Array<int> ess_bdr_filter(maxat);
+   ess_bdr_filter = 0;
+   ConstantCoefficient eps2_cf(epsilon*epsilon);
+   DiffusionSolver * filterSolver = new DiffusionSolver();
+   filterSolver->SetMesh(&mesh);
+   filterSolver->SetOrder(filter_fec.GetOrder());
+   filterSolver->SetDiffusionCoefficient(&eps2_cf);
+   filterSolver->SetMassCoefficient(&one);
+   filterSolver->SetEssentialBoundary(ess_bdr_filter);
+   filterSolver->SetupFEM();
+
+   BilinearForm mass(&control_fes);
+   mass.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator(one)));
+   mass.Assemble();
+   SparseMatrix M;
+   Array<int> empty;
+   mass.FormSystemMatrix(empty,M);
+
+   // 8. Define the Lagrange multiplier and gradient functions
+   GridFunction grad(&control_fes);
+   GridFunction w_filter(&filter_fes);
+
+   // 9. Define some tools for later
+   GridFunction onegf(&control_fes);
+   onegf = 1.0;
+   LinearForm vol_form(&control_fes);
+   vol_form.AddDomainIntegrator(new DomainLFIntegrator(one));
+   vol_form.Assemble();
+   double domain_volume = vol_form(onegf);
+
+   // 10. Connect to GLVis. Prepare for VisIt output.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sout_u,sout_r,sout_rho;
+   if (visualization)
+   {
+      sout_u.open(vishost, visport);
+      sout_rho.open(vishost, visport);
+      sout_r.open(vishost, visport);
+      sout_u.precision(8);
+      sout_rho.precision(8);
+      sout_r.precision(8);
+
+      sout_u << "solution\n" << mesh << u;
+      sout_u << "view 0 0\n";  // view from top
+      sout_u << "keys jl********\n";  // turn off perspective and light
+      sout_u << "window_title 'Temperature u'";
+      sout_u.flush();
+
+      GridFunction rho_gf(&control_fes);
+      rho_gf.ProjectCoefficient(rho);
+
+      sout_rho << "solution\n" << mesh << rho_gf;
+      sout_rho << "view 0 0\n";  // view from top
+      sout_rho << "keys jl********\n";  // turn off perspective and light
+      sout_rho << "window_title 'Density ρ'";
+      sout_rho.flush();
+
+      sout_r << "solution\n" << mesh << rho_filter;
+      sout_r << "view 0 0\n";  // view from top
+      sout_r << "keys jl********\n";  // turn off perspective and light
+      sout_r << "window_title 'Filtered density ρ̃'";
+      sout_r.flush();
+   }
+
+   // mfem::ParaViewDataCollection paraview_dc("Elastic_compliance", &mesh);
+   // paraview_dc.SetPrefixPath("ParaView");
+   // paraview_dc.SetLevelsOfDetail(order);
+   // paraview_dc.SetCycle(0);
+   // paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   // paraview_dc.SetHighOrderOutput(true);
+   // paraview_dc.SetTime(0.0);
+   // paraview_dc.RegisterField("displacement",&u);
+   // paraview_dc.RegisterField("density",&rho);
+   // paraview_dc.RegisterField("filtered_density",&rho_filter);
+
+   // 11. Iterate
+   int step = 0;
+   double c0 = 0.0;
+   GridFunction zero_gf(&control_fes);
+   for (int k = 1; k < max_it; k++)
+   {
+      const double alpha = alpha0 * k;
+
+      cout << "\nStep = " << k << endl;
+
+      // Step 1 - Filter solve
+      // Solve (ϵ^2 ∇ ρ̃, ∇ v ) + (ρ̃,v) = (ρ,v)
+      // GridFunctionCoefficient rho_cf(&rho);
+      rho_filter_old = rho_filter;
+      filterSolver->SetRHSCoefficient(&rho);
+      filterSolver->Solve();
+      rho_filter = *filterSolver->GetFEMSolution();
+
+      // Step 2 - State solve
+      // Solve (r(ρ̃) ∇u, ∇v) = (f,v)
+      SIMPInterpolationCoefficient SIMP_cf(&rho_filter, rho_min, 1.0);
+      GridFunctionPlusCoefficient u_plus_SIMP(&u, &SIMP_cf);
+      diffusionSolver->SetDiffusionCoefficient(&u_plus_SIMP);
+      diffusionSolver->Solve();
+      u = *diffusionSolver->GetFEMSolution();
+
+      // Step 3 - Adjoint filter solve
+      // Solve (ϵ² ∇ w̃, ∇ v) + (w̃ ,v) = (-r'(ρ̃) (|∇ u|²),v)
+      DiffusionEnergyCoefficient rhs_cf(&one, &u, &rho_filter,
+                                        rho_min);
+      filterSolver->SetRHSCoefficient(&rhs_cf);
+      filterSolver->Solve();
+      w_filter = *filterSolver->GetFEMSolution();
+
+      // Step 4 - Compute gradient
+      // Solve G = M⁻¹w̃
+      GridFunctionCoefficient w_cf(&w_filter);
+      LinearForm w_rhs(&control_fes);
+      w_rhs.AddDomainIntegrator(new DomainLFIntegrator(w_cf));
+      w_rhs.Assemble();
+      M.Mult(w_rhs,grad);
+
+      // Step 5 - Update design variable ψ ← projit(ψ - αG + c)
+      //          where c is a constant so that
+      //
+      //               ∫ρ = ∫sigmoid(ψ) = θ|Ω|
+      grad *= alpha;
+      psi -= grad;
+      projit(psi, mass_fraction*domain_volume);
+
+
+      // Step 6 - Compute other quantities
+      GridFunctionCoefficient rho_filter_cf(&rho_filter);
+      const double norm_reduced_gradient = rho_filter_old.ComputeL2Error(
+                                              rho_filter_cf)/alpha;
+      const double compliance = (*(diffusionSolver->GetLinearForm()))(u);
+      const double material_volume = zero_gf.ComputeL2Error(rho);
+      mfem::out << "norm of reduced gradient = " << norm_reduced_gradient << endl;
+      mfem::out << "compliance = " << compliance << endl;
+      mfem::out << "mass_fraction = " << material_volume / domain_volume << endl;
+
+      if (visualization)
+      {
+         sout_u << "solution\n" << mesh << u
+                << "window_title 'Displacement u'" << flush;
+
+         GridFunction rho_gf(&control_fes);
+         rho_gf.ProjectCoefficient(rho);
+         sout_rho << "solution\n" << mesh << rho_gf
+                  << "window_title 'Control variable ρ'" << flush;
+
+         GridFunction r_gf(&filter_fes);
+         r_gf.ProjectCoefficient(SIMP_cf);
+         sout_r << "solution\n" << mesh << r_gf
+                << "window_title 'Design density r(ρ̃)'" << flush;
+      }
+
+      if (norm_reduced_gradient < tol && k > 1)
+      {
+         break;
+      }
+   }
+
+   delete diffusionSolver;
+   delete filterSolver;
+
+   return 0;
+}

examples/makefile

@@ -23,10 +23,10 @@ MFEM_LIB_FILE = mfem_is_not_built
 
 SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
   ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29 ex30 \
-  ex31 ex33
+  ex31 ex33 ex34 ex35 efem_thermal
 PAR_EXAMPLES = ex0p 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 ex28p ex29p ex30p ex31p ex32p ex33p
+  ex25p ex26p ex27p ex28p ex29p ex30p ex31p ex32p ex33p ex34p ex35p
 SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex22 ex24 ex25 ex26
 PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex22p \
   ex24p ex25p ex26p
@@ -92,10 +92,13 @@ $(SUBDIRS_TPRINT):
 # Additional dependencies
 ex18: $(SRC)ex18.hpp
 ex33: $(SRC)ex33.hpp
+ex35: $(SRC)ex35.hpp
+efem_thermal: $(SRC)efem.hpp
 
 ifeq ($(MFEM_USE_MPI),YES)
 ex18p: $(SRC)ex18.hpp
 ex33p: $(SRC)ex33.hpp
+ex35p: $(SRC)ex35.hpp
 endif
 
 MFEM_TESTS = EXAMPLES