integ_algoim.cpp

Go to the documentation of this file.

// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
 
#include "integ_algoim.hpp"
 
#ifdef MFEM_USE_ALGOIM
 
namespace mfem
{
 
AlgoimIntegrationRule::AlgoimIntegrationRule(int o, const FiniteElement &el,
                                             ElementTransformation &trans,
                                             const Vector &lsfun)
{
   int_order=o;
   vir=nullptr;
   sir=nullptr;
 
   if (el.GetGeomType()==Geometry::Type::SQUARE)
   {
      pe=new H1Pos_QuadrilateralElement(el.GetOrder());
   }
   else if (el.GetGeomType()==Geometry::Type::CUBE)
   {
      pe=new H1Pos_HexahedronElement(el.GetOrder());
   }
   else
   {
      MFEM_ABORT("Currently MFEM + Algoim supports only quads and hexes.");
   }
 
   // change the basis of the level-set function
   // from Lagrangian to Bernstein (positive)
   lsvec.SetSize(pe->GetDof());
   DenseMatrix T(pe->GetDof());
   pe->Project(el,trans,T);
   T.Mult(lsfun,lsvec);
}
 
const IntegrationRule* AlgoimIntegrationRule::GetVolumeIntegrationRule()
{
   if (vir!=nullptr) {return vir;}
 
   const int dim=pe->GetDim();
   int np1d=int_order/2+1;
   if (dim==2)
   {
      LevelSet2D ls(pe,lsvec);
      auto q = Algoim::quadGen<2>(ls,Algoim::BoundingBox<real_t,2>(0.0,1.0),
                                  -1, -1, np1d);
 
      vir=new IntegrationRule(q.nodes.size());
      vir->SetOrder(int_order);
      for (size_t i=0; i<q.nodes.size(); i++)
      {
         IntegrationPoint& ip=vir->IntPoint(i);
         ip.Set2w(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].w);
      }
   }
   else
   {
      LevelSet3D ls(pe,lsvec);
      auto q = Algoim::quadGen<3>(ls,Algoim::BoundingBox<real_t,3>(0.0,1.0),
                                  -1, -1, np1d);
 
      vir=new IntegrationRule(q.nodes.size());
      vir->SetOrder(int_order);
      for (size_t i=0; i<q.nodes.size(); i++)
      {
         IntegrationPoint& ip=vir->IntPoint(i);
         ip.Set(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].x(2),q.nodes[i].w);
      }
   }
 
   return vir;
}
 
const IntegrationRule* AlgoimIntegrationRule::GetSurfaceIntegrationRule()
{
   if (sir!=nullptr) {return sir;}
 
   int np1d=int_order/2+1;
   const int dim=pe->GetDim();
   if (dim==2)
   {
      LevelSet2D ls(pe,lsvec);
      auto q = Algoim::quadGen<2>(ls,Algoim::BoundingBox<real_t,2>(0.0,1.0),
                                  2, -1, np1d);
 
      sir=new IntegrationRule(q.nodes.size());
      sir->SetOrder(int_order);
      for (size_t i=0; i<q.nodes.size(); i++)
      {
         IntegrationPoint& ip=sir->IntPoint(i);
         ip.Set2w(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].w);
      }
   }
   else
   {
      LevelSet3D ls(pe,lsvec);
      auto q = Algoim::quadGen<3>(ls,Algoim::BoundingBox<real_t,3>(0.0,1.0),
                                  3, -1, np1d);
 
      sir=new IntegrationRule(q.nodes.size());
      sir->SetOrder(int_order);
      for (size_t i=0; i<q.nodes.size(); i++)
      {
         IntegrationPoint& ip=sir->IntPoint(i);
         ip.Set(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].x(2),q.nodes[i].w);
      }
   }
 
   return sir;
}
 
}
 
#endif