stfc · nmnobre · Jun 29, 2026 · Jun 29, 2026 · Jun 29, 2026 · Jun 29, 2026
diff --git a/coefficients/EPGradientLogCoefficient.hpp b/coefficients/EPGradientLogCoefficient.hpp
@@ -0,0 +1,37 @@
+#include "mfem.hpp"
+using namespace mfem;
+
+class EPGradientLogCoefficient : public VectorCoefficient
+{
+private:
+  GridFunctionCoefficient * _electrolyte_concentration_gfc = nullptr;
+  GradientGridFunctionCoefficient _electrolyte_concentration_ggfc;
+
+public:
+  EPGradientLogCoefficient(GridFunctionCoefficient & ec)
+    : VectorCoefficient(1),
+      _electrolyte_concentration_gfc(&ec),
+      _electrolyte_concentration_ggfc(ec.GetGridFunction())
+  {
+  }
+
+  virtual void Eval(Vector &V, ElementTransformation & T, const IntegrationPoint & ip) override
+  {
+    _electrolyte_concentration_ggfc.Eval(V, T, ip);
+    const real_t ec = _electrolyte_concentration_gfc->Eval(T, ip);
+    V[0] *= 2 * constants::T * (1 - TPLUS) * Kappa(ec) / ec;
+
+    switch (T.Attribute)
+    {
+      case NE:
+        V[0] *= BNE / (LNE / NNE * NX);
+        break;
+      case SEP:
+        V[0] *= BSEP / (LSEP / NSEP * NX);
+        break;
+      case PE:
+        V[0] *= BPE / (LPE / NPE * NX);
+        break;
+    }
+  }
+};
diff --git a/coefficients/ElectrodeReactionCurrentCoefficient.hpp b/coefficients/ElectrodeReactionCurrentCoefficient.hpp
@@ -7,26 +7,26 @@ class ElectrodeReactionCurrentCoefficient : public Coefficient
   ExchangeCurrentCoefficient * _jex = nullptr;
   OverPotentialCoefficient * _op = nullptr;
 
-  const real_t _a;
-  const Region _r;
-  const int _sign;
+  real_t _a;
+  Region _r;
+  int _sign;
 
 public:
   /// P2D
-  ElectrodeReactionCurrentCoefficient(const Region & r,
-                                      const int & sign,
-                                      ExchangeCurrentCoefficient & jex,
+  ElectrodeReactionCurrentCoefficient(ExchangeCurrentCoefficient & jex,
                                       OverPotentialCoefficient & op)
     : _jex(&jex),
-      _op(&op),
-      _a(r == NE   ? AN * LNE / NNE * NX
-         : r == PE ? AP * LPE / NPE * NX
-                   : 0),
-      _r(r),
-      _sign(sign)
+      _op(&op)
   {
   }
 
+  void SetRegionSign(const Region &r, const int & sign)
+  {
+    _a = r == NE ? AN * LNE / NNE * NX : r == PE ? AP * LPE / NPE * NX : 0;
+    _r = r;
+    _sign = sign;
+  }
+
   /// P2D
   virtual real_t Eval(ElementTransformation & T, const IntegrationPoint & ip) override
   {

diff --git a/equations/ElectrolyteConcentration.cpp b/equations/ElectrolyteConcentration.cpp
@@ -2,7 +2,7 @@
 #include "equations/ElectrolyteConcentration.hpp"
 
 void
-ElectrolyteConcentration::Update(const GridFunctionCoefficient & ec_gfc, const Coefficient & j)
+ElectrolyteConcentration::Update()
 {
   // Mass coefficient.
   Vector mass_vec({/* NE */ EPS_N /* length scaling */ * (LNE / NNE * NX),
@@ -17,14 +17,14 @@ ElectrolyteConcentration::Update(const GridFunctionCoefficient & ec_gfc, const C
                      /* PE */ (1 - TPLUS) * AP /* length scaling */ * (LPE / NPE * NX)});
 
   PWConstCoefficient source_part(source_vec);
-  ProductCoefficient source(source_part, const_cast<Coefficient &>(j));
+  ProductCoefficient source(source_part, j);
 
   // Diffusion coefficient.
   Vector D_scale_vec({/* NE */ BNE /* length scaling */ / (LNE / NNE * NX),
                       /* SEP */ BSEP /* length scaling */ / (LSEP / NSEP * NX),
                       /* PE */ BPE /* length scaling */ / (LPE / NPE * NX)});
 
-  TransformedCoefficient D_coeff(&const_cast<GridFunctionCoefficient &>(ec_gfc), DE);
+  TransformedCoefficient D_coeff(&ec_gfc, DE);
   PWConstCoefficient D_scale_coeff(D_scale_vec);
   ProductCoefficient D(D_scale_coeff, D_coeff);
 

diff --git a/equations/ElectrolyteConcentration.hpp b/equations/ElectrolyteConcentration.hpp
@@ -7,7 +7,11 @@ using namespace mfem;
 
 class ElectrolyteConcentration : public Equation
 {
+protected:
+  Coefficient & j;
+  GridFunctionCoefficient & ec_gfc;
+
 public:
-  using Equation::Equation;
-  virtual void Update(const GridFunctionCoefficient & ec_gfc, const Coefficient & j) override;
+  ElectrolyteConcentration(ParFiniteElementSpace & f, Coefficient & j, GridFunctionCoefficient & ec_gfc) : Equation(f), j(j), ec_gfc(ec_gfc) {};
+  virtual void Update();
 };
diff --git a/equations/ElectrolytePotential.cpp b/equations/ElectrolytePotential.cpp
@@ -2,30 +2,25 @@
 #include "equations/ElectrolytePotential.hpp"
 
 void
-ElectrolytePotential::Update(const GridFunctionCoefficient & ec_gfc, const Coefficient & j)
+ElectrolytePotential::Update()
 {
   // Source term.
   Vector source_vec({/* NE */ AN /* length scaling */ * (LNE / NNE * NX),
                      /* SEP */ 0.,
                      /* PE */ AP /* length scaling */ * (LPE / NPE * NX)});
 
   PWConstCoefficient source_part(source_vec);
-  ProductCoefficient source(source_part, const_cast<Coefficient &>(j));
+  ProductCoefficient source(source_part, j);
 
   Vector b_vec({/* NE */ BNE /* length scaling */ / (LNE / NNE * NX),
                 /* SEP */ BSEP /* length scaling */ / (LSEP / NSEP * NX),
                 /* PE */ BPE /* length scaling */ / (LPE / NPE * NX)});
 
   PWConstCoefficient b_part(b_vec);
 
-  TransformedCoefficient kappa(&const_cast<GridFunctionCoefficient &>(ec_gfc), Kappa);
+  TransformedCoefficient kappa(&ec_gfc, Kappa);
   ProductCoefficient kappa_eff(b_part, kappa);
 
-  GradientGridFunctionCoefficient grad_ec(ec_gfc.GetGridFunction());
-  RatioCoefficient ec_inv(2 * T * (1 - TPLUS), const_cast<GridFunctionCoefficient &>(ec_gfc));
-  ScalarVectorProductCoefficient grad_ln_ec(ec_inv, grad_ec);
-  ScalarVectorProductCoefficient grad_ln_ec_kappad(kappa_eff, grad_ln_ec);
-
   if (!K)
   {
     K = new ParBilinearForm(&fespace);
@@ -37,7 +32,7 @@ ElectrolytePotential::Update(const GridFunctionCoefficient & ec_gfc, const Coeff
   if (!Qc)
   {
     Qc = new ParLinearForm(&fespace);
-    Qc->AddDomainIntegrator(new DomainLFGradIntegrator(grad_ln_ec_kappad));
+    Qc->AddDomainIntegrator(new DomainLFGradIntegrator(gl));
     Qc->Assemble();
     Qc->ParallelAssemble(bc);
   }

diff --git a/equations/ElectrolytePotential.hpp b/equations/ElectrolytePotential.hpp
@@ -2,21 +2,25 @@
 
 #include "mfem.hpp"
 #include "equations/Equation.hpp"
+#include "coefficients/EPGradientLogCoefficient.hpp"
 
 using namespace mfem;
 
 class ElectrolytePotential : public Equation
 {
 protected:
+  Coefficient & j;
+  GridFunctionCoefficient & ec_gfc;
+  EPGradientLogCoefficient gl;
   ParLinearForm * Qc = nullptr;
   Vector bc;
 
 public:
-  ElectrolytePotential(ParFiniteElementSpace & f) : Equation(f), bc(f.GetTrueVSize())
+  ElectrolytePotential(ParFiniteElementSpace & f, Coefficient & j, GridFunctionCoefficient & ec_gfc) : Equation(f), j(j), ec_gfc(ec_gfc), gl(ec_gfc), bc(f.GetTrueVSize())
   {
     f.GetEssentialTrueDofs(Array<int>({1, 0}), ess_tdof_list);
   }
-  virtual void Update(const GridFunctionCoefficient & ec_gfc, const Coefficient & j) override;
+  virtual void Update();
   void Reset()
   {
     delete K;

diff --git a/equations/Equation.hpp b/equations/Equation.hpp
@@ -31,9 +31,6 @@ class Equation
   const HypreParMatrix & GetK() const { return Kmat; };
   const Vector & GetZ() const { return b; };
 
-  virtual void Update(const Coefficient & j) {}
-  virtual void Update(const GridFunctionCoefficient & u, const Coefficient & j) {}
-
   virtual ~Equation()
   {
     delete M;

diff --git a/equations/SolidConcentration.cpp b/equations/SolidConcentration.cpp
@@ -1,7 +1,7 @@
 #include "equations/SolidConcentration.hpp"
 
 void
-SolidConcentration::Update(const Coefficient & j)
+SolidConcentration::Update(const real_t & j)
 {
   MFEM_ASSERT(particle_region == NE || particle_region == PE, "Particle not in electrode!");
 
@@ -11,7 +11,7 @@ SolidConcentration::Update(const Coefficient & j)
 
   FunctionCoefficient r2([](const Vector & r) { return r(0) * r(0); });
   ProductCoefficient dr2(D / R / R, r2);
-  ProductCoefficient jjr2(const_cast<Coefficient &>(j), r2);
+  ProductCoefficient jjr2(j, r2);
   ProductCoefficient jr2(-1. / R / t_scale, jjr2);
 
   if (!M)

diff --git a/equations/SolidConcentration.hpp b/equations/SolidConcentration.hpp
@@ -37,7 +37,7 @@ class SolidConcentration : public Equation
   {
   }
 
-  virtual void Update(const Coefficient & j) override;
+  virtual void Update(const real_t & j);
   real_t SurfaceConcentration(const BlockVector & x);
   int GetParticleRank() { return particle_rank; }
   int GetSurfaceRank() { return surface_rank; }

diff --git a/equations/SolidPotential.cpp b/equations/SolidPotential.cpp
@@ -2,15 +2,15 @@
 #include "equations/SolidPotential.hpp"
 
 void
-SolidPotential::Update(const Coefficient & j)
+SolidPotential::Update()
 {
   // Source term.
   Vector source_vec({/* NE */ -AN /* length scaling */ * (LNE / NNE * NX),
                      /* SEP */ 0.,
                      /* PE */ -AP /* length scaling */ * (LPE / NPE * NX)});
 
   PWConstCoefficient source_part(source_vec);
-  ProductCoefficient source(source_part, const_cast<Coefficient &>(j));
+  ProductCoefficient source(source_part, j);
 
   // Effective conductivity (does not account for electrode filler).
   Vector sigma_vec({/* NE */ (1 - EPS_N) * SIGN /* length scaling */ / (LNE / NNE * NX),

diff --git a/equations/SolidPotential.hpp b/equations/SolidPotential.hpp
@@ -7,10 +7,13 @@ using namespace mfem;
 
 class SolidPotential : public Equation
 {
+protected:
+  Coefficient &j;
+
 public:
-  SolidPotential(ParFiniteElementSpace & f) : Equation(f)
+  SolidPotential(ParFiniteElementSpace & f, Coefficient & j) : Equation(f), j(j)
   {
     f.GetEssentialTrueDofs(Array<int>({1, 1}), ess_tdof_list);
   }
-  virtual void Update(const Coefficient & j) override;
+  virtual void Update();
 };
diff --git a/makefile b/makefile
@@ -24,7 +24,8 @@ CF_INC_FILES = ExchangeCurrentCoefficient.hpp \
                ReactionCurrentCoefficient.hpp \
                ElectrodeReactionCurrentCoefficient.hpp \
                OpenCircuitPotentialCoefficient.hpp \
-               OverPotentialCoefficient.hpp
+               OverPotentialCoefficient.hpp \
+               EPGradientLogCoefficient.hpp
 CL_INC_FILES = LGM50.hpp
 PR_SRC_FILES = settings.cpp
 PR_INC_FILES = settings.hpp \