DGLaxFriedrichsFluxIntegrator interface gradient verified.

LLNL · May 7, 2024 · e1d6dee · e1d6dee
1 parent 4537d56
commit e1d6dee
Show file tree

Hide file tree

Showing 4 changed files with 401 additions and 0 deletions.
diff --git a/include/interfaceinteg.hpp b/include/interfaceinteg.hpp
@@ -304,6 +304,20 @@ class DGLaxFriedrichsFluxIntegrator : public InterfaceNonlinearFormIntegrator
                            FaceElementTransformations &T, const IntegrationPoint &ip,
                            const Vector &x, DenseMatrix &jac) override;
 
+   void AssembleInterfaceVector(const FiniteElement &el1,
+                                 const FiniteElement &el2,
+                                 FaceElementTransformations &Tr1,
+                                 FaceElementTransformations &Tr2,
+                                 const Vector &elfun1, const Vector &elfun2,
+                                 Vector &elvect1, Vector &elvect2) override;
+
+   void AssembleInterfaceGrad(const FiniteElement &el1,
+                              const FiniteElement &el2,
+                              FaceElementTransformations &Tr1,
+                              FaceElementTransformations &Tr2,
+                              const Vector &elfun1, const Vector &elfun2,
+                              Array2D<DenseMatrix *> &elmats) override;
+
 private:
    void AppendPrecomputeFaceCoeffs(const FiniteElementSpace *fes, 
                                     FaceElementTransformations *T,

diff --git a/src/interfaceinteg.cpp b/src/interfaceinteg.cpp
@@ -2058,4 +2058,242 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleGrad_Fast(
    }
 }
 
+void DGLaxFriedrichsFluxIntegrator::AssembleInterfaceVector(
+   const FiniteElement &el1, const FiniteElement &el2,
+   FaceElementTransformations &Tr1, FaceElementTransformations &Tr2,
+   const Vector &elfun1, const Vector &elfun2,
+   Vector &elvect1, Vector &elvect2)
+{
+   dim = el1.GetDim();
+   ndofs1 = el1.GetDof();
+   ndofs2 = el2.GetDof();
+   nvdofs = dim * (ndofs1 + ndofs2);
+   elvect1.SetSize(dim * ndofs1);
+   elvect2.SetSize(dim * ndofs2);
+
+   nor.SetSize(dim);
+   flux.SetSize(dim);
+
+   udof1.UseExternalData(elfun1.GetData(), ndofs1, dim);
+   elv1.UseExternalData(elvect1.GetData(), ndofs1, dim);
+   shape1.SetSize(ndofs1);
+   u1.SetSize(dim);
+
+   udof2.UseExternalData(elfun2.GetData(), ndofs2, dim);
+   elv2.UseExternalData(elvect2.GetData(), ndofs2, dim);
+   shape2.SetSize(ndofs2);
+   u2.SetSize(dim);
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      // a simple choice for the integration order; is this OK?
+      const int order = (int)(ceil(1.5 * (2 * max(el1.GetOrder(), el2.GetOrder()) - 1)));
+      ir = &IntRules.Get(Tr1.GetGeometryType(), order);
+   }
+
+   elvect1 = 0.0; elvect2 = 0.0;
+   for (int pind = 0; pind < ir->GetNPoints(); ++pind)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(pind);
+
+      // Set the integration point in the face and the neighboring elements
+      Tr1.SetAllIntPoints(&ip);
+      Tr2.SetAllIntPoints(&ip);
+
+      // Access the neighboring elements' integration points
+      // Note: only Element1 of each Tr is valid.
+      const IntegrationPoint &eip1 = Tr1.GetElement1IntPoint();
+      const IntegrationPoint &eip2 = Tr2.GetElement1IntPoint();
+
+      el1.CalcShape(eip1, shape1);
+      udof1.MultTranspose(shape1, u1);
+
+      if (dim == 1)
+      {
+         nor(0) = 2*eip1.x - 1.0;
+      }
+      else
+      {
+         CalcOrtho(Tr1.Jacobian(), nor);
+      }
+
+      el2.CalcShape(eip2, shape2);
+      udof2.MultTranspose(shape2, u2);
+
+      w = 0.5 * ip.weight * Tr1.Weight();
+      if (Q) { w *= Q->Eval(*Tr1.Elem1, eip1); }
+
+      nor *= w;
+
+      un1 = nor * u1;
+      un2 = (ndofs2) ? nor * u2 : 0.0;
+
+      flux.Set(un1, u1);
+      flux.Add(un2, u2);
+
+      un = 2.0 * max(abs(un1), abs(un2));
+      flux.Add(un, u1);
+      if (ndofs2)
+         flux.Add(-un, u2);
+
+      AddMultVWt(shape1, flux, elv1);
+      if (ndofs2)
+         AddMult_a_VWt(-1.0, shape2, flux, elv2);
+   }
+}
+
+void DGLaxFriedrichsFluxIntegrator::AssembleInterfaceGrad(
+   const FiniteElement &el1, const FiniteElement &el2,
+   FaceElementTransformations &Tr1, FaceElementTransformations &Tr2,
+   const Vector &elfun1, const Vector &elfun2, Array2D<DenseMatrix *> &elmats)
+{
+   assert(elmats.NumRows() == 2);
+   assert(elmats.NumCols() == 2);
+   for (int i = 0; i < 2; i++)
+      for (int j = 0; j < 2; j++) assert(elmats(i, j));
+
+   dim = el1.GetDim();
+   ndofs1 = el1.GetDof();
+   ndofs2 = el2.GetDof();
+
+   elmats(0, 0)->SetSize(dim * ndofs1);
+   elmats(0, 1)->SetSize(dim * ndofs1, dim * ndofs2);
+   elmats(1, 0)->SetSize(dim * ndofs2, dim * ndofs1);
+   elmats(1, 1)->SetSize(dim * ndofs2);
+
+   elmat_comp11.SetSize(ndofs1);
+   elmat_comp12.SetSize(ndofs1, ndofs2);
+   elmat_comp21.SetSize(ndofs2, ndofs1);
+   elmat_comp22.SetSize(ndofs2);
+
+   nor.SetSize(dim);
+   flux.SetSize(dim);
+
+   udof1.UseExternalData(elfun1.GetData(), ndofs1, dim);
+   shape1.SetSize(ndofs1);
+   u1.SetSize(dim);  
+
+   udof2.UseExternalData(elfun2.GetData(), ndofs2, dim);
+   shape2.SetSize(ndofs2);
+   u2.SetSize(dim);
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      // a simple choice for the integration order; is this OK?
+      const int order = (int)(ceil(1.5 * (2 * max(el1.GetOrder(), ndofs2 ? el2.GetOrder() : 0) - 1)));
+      ir = &IntRules.Get(Tr1.GetGeometryType(), order);
+   } 
+
+   for (int i = 0; i < 2; i++)
+      for (int j = 0; j < 2; j++) *(elmats(i, j)) = 0.0;
+
+   for (int pind = 0; pind < ir->GetNPoints(); ++pind)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(pind);
+
+      // Set the integration point in the face and the neighboring elements
+      Tr1.SetAllIntPoints(&ip);
+      Tr2.SetAllIntPoints(&ip);
+
+      // Access the neighboring elements' integration points
+      // Note: eip2 will only contain valid data if Elem2 exists
+      const IntegrationPoint &eip1 = Tr1.GetElement1IntPoint();
+      const IntegrationPoint &eip2 = Tr2.GetElement1IntPoint();
+
+      el1.CalcShape(eip1, shape1);
+      udof1.MultTranspose(shape1, u1);
+
+      if (dim == 1)
+      {
+         nor(0) = 2*eip1.x - 1.0;
+      }
+      else
+      {
+         CalcOrtho(Tr1.Jacobian(), nor);
+      }
+
+      el2.CalcShape(eip2, shape2);
+      udof2.MultTranspose(shape2, u2);
+
+      w = ip.weight * Tr1.Weight();
+      if (Q) { w *= Q->Eval(*Tr1.Elem1, eip1); }
+
+      nor *= w;
+      // Just for the average operator gradient.
+      nor *= 0.5;
+
+      un1 = nor * u1;
+      un2 = (ndofs2) ? nor * u2 : 0.0;
+
+      MultVVt(shape1, elmat_comp11);
+      MultVWt(shape1, shape2, elmat_comp12);
+      elmat_comp21.Transpose(elmat_comp12);
+      MultVVt(shape2, elmat_comp22);
+
+      for (int di = 0; di < dim; di++)
+      {
+         elmats(0, 0)->AddMatrix(un1, elmat_comp11, di * ndofs1, di * ndofs1);
+         elmats(1, 0)->AddMatrix(-un1, elmat_comp21, di * ndofs2, di * ndofs1);
+         elmats(0, 1)->AddMatrix(un2, elmat_comp12, di * ndofs1, di * ndofs2);
+         elmats(1, 1)->AddMatrix(-un2, elmat_comp22, di * ndofs2, di * ndofs2);
+
+         for (int dj = 0; dj < dim; dj++)
+         {
+            elmats(0, 0)->AddMatrix(u1(di) * nor(dj), elmat_comp11, di * ndofs1, dj * ndofs1);
+            elmats(1, 0)->AddMatrix(-u1(di) * nor(dj), elmat_comp21, di * ndofs2, dj * ndofs1);
+            elmats(0, 1)->AddMatrix(u2(di) * nor(dj), elmat_comp12, di * ndofs1, dj * ndofs2);
+            elmats(1, 1)->AddMatrix(-u2(di) * nor(dj), elmat_comp22, di * ndofs2, dj * ndofs2);
+         }
+      }
+
+      // Recover 1/2 factor on normal vector.
+      if (ndofs2)
+         nor *= 2.0;
+
+      // un1 as maximum absolute eigenvalue, un2 as the sign.
+      un = max( abs(un1), abs(un2) );
+      if (ndofs2) un *= 2.0;
+
+      double sgn = 0.0;
+      bool u1_lg_u2 = (abs(un1) >= abs(un2));
+      if (u1_lg_u2)
+      {
+         sgn = (un1 >= 0.0) ? 1.0 : -1.0;
+      }
+      else
+      {
+         sgn = (un2 >= 0.0) ? 1.0 : -1.0;
+      }
+
+      // [ u ] = u1 - u2
+      if (ndofs2)
+         u1.Add(-1.0, u2);
+
+      for (int di = 0; di < dim; di++)
+      {
+         elmats(0, 0)->AddMatrix(un, elmat_comp11, di * ndofs1, di * ndofs1);
+         elmats(1, 0)->AddMatrix(-un, elmat_comp21, di * ndofs2, di * ndofs1);
+         elmats(0, 1)->AddMatrix(-un, elmat_comp12, di * ndofs1, di * ndofs2);
+         elmats(1, 1)->AddMatrix(un, elmat_comp22, di * ndofs2, di * ndofs2);
+
+         // remember u1 in this loop is in fact [ u ] = u1 - u2.
+         for (int dj = 0; dj < dim; dj++)
+         {
+            if (u1_lg_u2)
+            {
+               elmats(0, 0)->AddMatrix(sgn * u1(di) * nor(dj), elmat_comp11, di * ndofs1, dj * ndofs1);
+               elmats(1, 0)->AddMatrix(-sgn * u1(di) * nor(dj), elmat_comp21, di * ndofs2, dj * ndofs1);
+            }
+            else
+            {
+               elmats(0, 1)->AddMatrix(sgn * u1(di) * nor(dj), elmat_comp12, di * ndofs1, dj * ndofs2);
+               elmats(1, 1)->AddMatrix(-sgn * u1(di) * nor(dj), elmat_comp22, di * ndofs2, dj * ndofs2);
+            }
+         }
+      }
+   }
+}
+
 }
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
@@ -65,6 +65,8 @@ add_executable(test_block_smoother test_block_smoother.cpp $<TARGET_OBJECTS:scal
 
 add_executable(nonlinear_integ_grad nonlinear_integ_grad.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 
+add_executable(interfaceinteg_grad interfaceinteg_grad.cpp $<TARGET_OBJECTS:scaleupROMObj>)
+
 add_executable(test_linalg_utils test_linalg_utils.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 
 add_executable(test_rom_nonlinearform test_rom_nonlinearform.cpp $<TARGET_OBJECTS:scaleupROMObj>)