wip 1D bug moments-based integration.

mfem · May 9, 2024 · 959b07c · 959b07c
1 parent 31bd94e
commit 959b07c
Show file tree

Hide file tree

Showing 2 changed files with 77 additions and 15 deletions.
diff --git a/fem/intrules_cut.cpp b/fem/intrules_cut.cpp
@@ -175,6 +175,7 @@ void MomentFittingIntRules::ComputeFaceWeights(ElementTransformation& Tr)
          local_mesh.GetElementTransformation(0, &faceTrafo);
 
          // The 3D face integrals are computed as 2D volumetric integrals.
+         // The 2D face integrals are computed as 1D volumetric integrals.
          MomentFittingIntRules FaceRules(Order, *LvlSet, lsOrder);
          IntegrationRule FaceRule;
          FaceRules.GetVolumeIntegrationRule(faceTrafo, FaceRule);
@@ -254,6 +255,55 @@ void MomentFittingIntRules::ComputeSurfaceWeights1D(ElementTransformation& Tr)
    }
 }
 
+double bisect(ElementTransformation &Tr, Coefficient *LvlSet)
+{
+   IntegrationPoint intp;
+
+   IntegrationPoint ip0;
+   ip0.x = 0.;
+   IntegrationPoint ip1;
+   ip1.x = 1.;
+   Tr.SetIntPoint(&ip0);
+   if (LvlSet->Eval(Tr, ip0) * LvlSet->Eval(Tr, ip1) < 0.)
+   {
+      IntegrationPoint ip2;
+      ip2.x = .5;
+      while (LvlSet->Eval(Tr, ip2) > 1e-12
+             || LvlSet->Eval(Tr, ip2) < -1e-12)
+      {
+         if (LvlSet->Eval(Tr, ip0) * LvlSet->Eval(Tr, ip2) < 0.)
+         {
+            ip1.x = ip2.x;
+         }
+         else
+         {
+            ip0.x = ip2.x;
+         }
+
+         ip2.x = (ip1.x + ip0.x) / 2.;
+      }
+      intp.x = ip2.x;
+      intp.weight = 1. / Tr.Weight();
+   }
+   else if (LvlSet->Eval(Tr, ip0) > 0. && LvlSet->Eval(Tr, ip1) <= 1e-12)
+   {
+      intp.x = 1.;
+      intp.weight = 1. / Tr.Weight();
+   }
+   else if (LvlSet->Eval(Tr, ip1) > 0. && LvlSet->Eval(Tr, ip0) <= 1e-12)
+   {
+      intp.x = 0.;
+      intp.weight = 1. / Tr.Weight();
+   }
+   else
+   {
+      intp.x = .5;
+      intp.weight = 0.;
+   }
+
+   return intp.x;
+}
+
 void MomentFittingIntRules::ComputeVolumeWeights1D(ElementTransformation& Tr,
                                                    const IntegrationRule* sir)
 {
@@ -272,6 +322,7 @@ void MomentFittingIntRules::ComputeVolumeWeights1D(ElementTransformation& Tr,
       if (LvlSet->Eval(Tr, ip0) > 0.)
       {
          length = sir->IntPoint(0).x;
+         //length = bisect(Tr, LvlSet);
          for (int ip = 0; ip < ir.GetNPoints(); ip++)
          {
             IntegrationPoint &intp = ir.IntPoint(ip);
@@ -282,10 +333,12 @@ void MomentFittingIntRules::ComputeVolumeWeights1D(ElementTransformation& Tr,
       else
       {
          length = 1. - sir->IntPoint(0).x;
+         //length = 1. - bisect(Tr, LvlSet);
          for (int ip = 0; ip < ir.GetNPoints(); ip++)
          {
             IntegrationPoint &intp = ir.IntPoint(ip);
             intp.x = sir->IntPoint(ip).x + ir2.IntPoint(ip).x * length;
+            //intp.x = bisect(Tr, LvlSet) + ir2.IntPoint(ip).x * length;
             intp.weight = ir2.IntPoint(ip).weight * length;
          }
       }
@@ -1490,24 +1543,32 @@ void MomentFittingIntRules::GetVolumeIntegrationRule(ElementTransformation& Tr,
       FaceWeightsComp = 0.;
    }
 
+   //std::cout << "before surface: " << ir.GetNPoints() << std::endl;
+
    IntegrationRule SIR;
-   if (sir == NULL)
-   {
-      Order++;
-      GetSurfaceIntegrationRule(Tr, SIR);
-      Order--;
-   }
-   else if ((sir->GetOrder() - 1) != ir.GetOrder())
+   // if (Tr.GetDimension() != 1)
+   if (true)
    {
-      Order++;
-      GetSurfaceIntegrationRule(Tr, SIR);
-      Order--;
-   }
-   else
-   {
-      SIR = *sir;
+      if (sir == NULL)
+      {
+         Order++;
+         GetSurfaceIntegrationRule(Tr, SIR);
+         Order--;
+      }
+      else if ((sir->GetOrder() - 1) != ir.GetOrder())
+      {
+         Order++;
+         GetSurfaceIntegrationRule(Tr, SIR);
+         Order--;
+      }
+      else
+      {
+         SIR = *sir;
+      }
    }
 
+   //std::cout << "here: " << ir.GetNPoints() << std::endl;
+
    if (Tr.GetDimension() == 1)
    {
       ComputeVolumeWeights1D(Tr, &SIR);
@@ -1521,6 +1582,7 @@ void MomentFittingIntRules::GetVolumeIntegrationRule(ElementTransformation& Tr,
       ComputeVolumeWeights3D(Tr, &SIR);
    }
 
+   //std::cout << "end: " << ir.GetNPoints() << std::endl;
    result.SetSize(ir.GetNPoints());
    for (int ip = 0; ip < ir.GetNPoints(); ip++)
    {

diff --git a/fem/intrules_cut.hpp b/fem/intrules_cut.hpp
@@ -204,7 +204,7 @@ class MomentFittingIntRules : public CutIntegrationRules
     @param [in] sir corresponding IntegrationRule on surface
     */
    void ComputeVolumeWeights1D(ElementTransformation& Tr,
-                               const IntegrationRule* sir);
+                               const IntegrationRule *sir);
 
    /**
     @brief Compute 2D quadrature weights