Merge pull request #105 from awslabs/sjg/coeff-ip-opt

Remove some unneeded `SetIntPoint` calls for `Coefficient` evaluation
awslabs · Oct 12, 2023 · d55abfa · d55abfa
2 parents 35bad95 + 7aea5c7
commit d55abfa
Show file tree

Hide file tree

Showing 2 changed files with 17 additions and 22 deletions.
diff --git a/palace/fem/coefficient.hpp b/palace/fem/coefficient.hpp
@@ -196,8 +196,15 @@ class BdrGridFunctionCoefficient
   {
   }
 
-  static void GetNormal(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip,
-                        mfem::Vector &normal);
+  // Return normal vector to the boundary element at an integration point (it is assumed
+  // that the element transformation has already been configured at the integration point of
+  // interest).
+  static void GetNormal(mfem::ElementTransformation &T, mfem::Vector &normal)
+  {
+    normal.SetSize(T.GetSpaceDim());
+    mfem::CalcOrtho(T.Jacobian(), normal);
+    normal /= normal.Norml2();
+  }
 };
 
 // Computes surface current J_s = n x H on boundaries from B as a vector grid function
@@ -242,7 +249,7 @@ class BdrCurrentVectorCoefficient : public mfem::VectorCoefficient,
     }
 
     // Orient with normal pointing into el1.
-    GetNormal(T, ip, nor);
+    GetNormal(T, nor);
     V.SetSize(vdim);
     if (C1 * nor < 0.0)
     {
@@ -296,7 +303,7 @@ class BdrChargeCoefficient : public mfem::Coefficient, public BdrGridFunctionCoe
     }
 
     // Orient with normal pointing into el1.
-    GetNormal(T, ip, nor);
+    GetNormal(T, nor);
     return (C1 * nor < 0.0) ? -(VU * nor) : VU * nor;
   }
 };
@@ -338,7 +345,7 @@ class BdrFluxCoefficient : public mfem::Coefficient, public BdrGridFunctionCoeff
     }
 
     // Orient sign with the global direction.
-    GetNormal(T, ip, nor);
+    GetNormal(T, nor);
     return (dir * nor < 0.0) ? -(V * nor) : V * nor;
   }
 };
@@ -431,7 +438,7 @@ inline double DielectricInterfaceCoefficient<DielectricInterfaceType::MA>::Eval(
 {
   // Get single-sided solution and neighboring element attribute.
   Initialize(T, ip, V);
-  GetNormal(T, ip, nor);
+  GetNormal(T, nor);
 
   // Metal-air interface: 0.5 * t / ε_MA * |E_n|² .
   double Vn = V * nor;
@@ -444,7 +451,7 @@ inline double DielectricInterfaceCoefficient<DielectricInterfaceType::MS>::Eval(
 {
   // Get single-sided solution and neighboring element attribute.
   int attr = Initialize(T, ip, V);
-  GetNormal(T, ip, nor);
+  GetNormal(T, nor);
 
   // Metal-substrate interface: 0.5 * t * (ε_S)² / ε_MS * |E_n|² .
   const double Vn = V * nor;
@@ -458,7 +465,7 @@ inline double DielectricInterfaceCoefficient<DielectricInterfaceType::SA>::Eval(
 {
   // Get single-sided solution and neighboring element attribute.
   Initialize(T, ip, V);
-  GetNormal(T, ip, nor);
+  GetNormal(T, nor);
 
   // Substrate-air interface: 0.5 * t * (ε_SA * |E_t|² + 1 / ε_MS * |E_n|²) .
   double Vn = V * nor;
@@ -509,7 +516,6 @@ class EnergyDensityCoefficient : public mfem::Coefficient, public BdrGridFunctio
   {
     if (T.ElementType == mfem::ElementTransformation::ELEMENT)
     {
-      T.SetIntPoint(&ip);
       return GetLocalEnergyDensity(T, ip, T.Attribute);
     }
     if (T.ElementType == mfem::ElementTransformation::BDR_ELEMENT)
@@ -673,7 +679,7 @@ class NormalProjectedCoefficient : public mfem::Coefficient
   double Eval(mfem::ElementTransformation &T, const mfem::IntegrationPoint &ip) override
   {
     c->Eval(K, T, ip);
-    BdrGridFunctionCoefficient::GetNormal(T, ip, nor);
+    BdrGridFunctionCoefficient::GetNormal(T, nor);
     return K.InnerProduct(nor, nor);
   }
 };
@@ -1030,17 +1036,6 @@ inline void BdrGridFunctionCoefficient::GetElementTransformations(
   }
 }
 
-inline void BdrGridFunctionCoefficient::GetNormal(mfem::ElementTransformation &T,
-                                                  const mfem::IntegrationPoint &ip,
-                                                  mfem::Vector &normal)
-{
-  // Return normal vector to the boundary element at the provided integration point.
-  normal.SetSize(T.GetSpaceDim());
-  T.SetIntPoint(&ip);
-  mfem::CalcOrtho(T.Jacobian(), normal);
-  normal /= normal.Norml2();
-}
-
 }  // namespace palace
 
 #endif  // PALACE_FEM_COEFFICIENT_HPP
diff --git a/palace/models/waveportoperator.cpp b/palace/models/waveportoperator.cpp
@@ -430,6 +430,7 @@ class BdrSubmeshHVectorCoefficient : public mfem::VectorCoefficient
       else
       {
         submesh_T = submesh.GetElementTransformation(it->second);
+        submesh_T->SetIntPoint(&ip);
       }
 
       int i, o, iel1, iel2;
@@ -456,7 +457,6 @@ class BdrSubmeshHVectorCoefficient : public mfem::VectorCoefficient
 
     // Compute Re/Im{-1/i (ikₙ Eₜ + ∇ₜ Eₙ)}.
     mfem::Vector U;
-    submesh_T->SetIntPoint(&ip);
     if constexpr (RealPart)
     {
       Et.real().GetVectorValue(*submesh_T, ip, U);