minor changes

examples/CMakeLists.txt

@@ -42,6 +42,7 @@ list(APPEND ALL_EXE_SRCS
   ex33.cpp
   ex34.cpp
   ex35.cpp
+  efem_thermal.cpp
   )
 
 if (MFEM_USE_MPI)

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/ex35.cpp

@@ -176,7 +176,7 @@ int main(int argc, char *argv[])
    int ref_levels = 4;
    int order = 2;
    bool visualization = true;
-   double alpha = 1.0;
+   double alpha0 = 1.0;
    double epsilon = 0.01;
    double mass_fraction = 0.5;
    int max_it = 1e2;
@@ -184,13 +184,14 @@ int main(int argc, char *argv[])
    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(&alpha, "-alpha", "--alpha-step-length",
+   args.AddOption(&alpha0, "-alpha", "--alpha-step-length",
                   "Step length for gradient descent.");
    args.AddOption(&epsilon, "-epsilon", "--epsilon-thickness",
                   "epsilon phase field thickness");
@@ -348,6 +349,24 @@ int main(int argc, char *argv[])
       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);
@@ -366,8 +385,7 @@ int main(int argc, char *argv[])
    double c0 = 0.0;
    for (int k = 1; k < max_it; k++)
    {
-      if (k > 1) { alpha *= ((double) k) / ((double) k-1); }
-      step++;
+      const double alpha = alpha0 * k;
 
       cout << "\nStep = " << k << endl;
 

examples/ex35.hpp

@@ -107,6 +107,42 @@ public:
    }
 };
 
+
+// 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
 {

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,7 +23,7 @@ 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 ex34 ex35
+  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 ex34p ex35p
@@ -93,6 +93,7 @@ $(SUBDIRS_TPRINT):
 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