Ghost-penalization

miniapps/shifted/CMakeLists.txt

@@ -14,13 +14,16 @@ if (MFEM_USE_MPI)
         sbm_solver.cpp
         marking.cpp
         extrapolator.cpp
-	integ_algoim.cpp)
+        integ_algoim.cpp
+        cut_integrators.cpp)
+
   list(APPEND DIST_COMMON_HEADERS
         sbm_solver.hpp
         sbm_aux.hpp
         marking.hpp
         extrapolator.hpp
-	integ_algoim.hpp)
+        integ_algoim.hpp
+        cut_integrators.hpp)
 
   convert_filenames_to_full_paths(DIST_COMMON_SOURCES)
   convert_filenames_to_full_paths(DIST_COMMON_HEADERS)
@@ -49,6 +52,11 @@ if (MFEM_USE_MPI)
     ${DIST_COMMON_FILES}
     LIBRARIES mfem mfem-common)
 
+  add_mfem_miniapp(cut_integral
+    MAIN cut_integral.cpp
+    ${DIST_COMMON_FILES}
+    LIBRARIES mfem mfem-common)
+
   if (MFEM_ENABLE_TESTING)
     add_test(NAME shifted_distance_np${MFEM_MPI_NP}
       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}

miniapps/shifted/cut_integral.cpp

@@ -0,0 +1,170 @@
+#include "mfem.hpp"
+#include "marking.hpp"
+#include "cut_integrators.hpp"
+
+using namespace mfem;
+
+//Level set function for sphere in 3D and circle in 2D
+double sphere_ls(const Vector &x)
+{
+   double r2= x*x;
+   return -sqrt(r2)+1.0;//the radius is 1.0
+}
+
+
+int main(int argc, char *argv[])
+{
+    // Initialize MPI and HYPRE.
+    Mpi::Init(argc, argv);
+    int myrank = Mpi::WorldRank();
+    Hypre::Init();
+
+    // Parse command-line options.
+    const char *mesh_file = "../../data/star-q3.mesh";
+    int order = 2;
+    bool visualization = true;
+    int ser_ref_levels = 0;
+    int par_ref_levels = 0;
+
+    OptionsParser args(argc, argv);
+    args.AddOption(&mesh_file, "-m", "--mesh",
+                   "Mesh file to use.");
+    args.AddOption(&order, "-o", "--order",
+                   "Finite element order (polynomial degree) or -1 for"
+                   " isoparametric space.");
+    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                   "--no-visualization",
+                   "Enable or disable GLVis visualization.");
+    args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                   "Number of times to refine the mesh uniformly in serial.");
+    args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
+                   "Number of times to refine the mesh uniformly in parallel.");
+    args.Parse();
+    if (!args.Good())
+    {
+       if (myrank == 0)
+       {
+          args.PrintUsage(std::cout);
+       }
+       return 1;
+    }
+    if (myrank == 0) { args.PrintOptions(std::cout); }
+
+
+
+    // Enable hardware devices such as GPUs, and programming models such as CUDA,
+    // OCCA, RAJA and OpenMP based on command line options.
+    Device device("cpu");
+    if (myrank == 0) { device.Print(); }
+
+    // Refine the mesh.
+    Mesh mesh(mesh_file, 1, 1);
+    int dim = mesh.Dimension();
+    for (int lev = 0; lev < ser_ref_levels; lev++) { mesh.UniformRefinement(); }
+    if (myrank == 0)
+    {
+       std::cout << "Number of elements: " << mesh.GetNE() << std::endl;
+    }
+
+    Coefficient *ls_coeff = nullptr;
+    ls_coeff=new FunctionCoefficient(sphere_ls);
+
+    // MPI distribution.
+    ParMesh pmesh(MPI_COMM_WORLD, mesh);
+    for (int lev = 0; lev < par_ref_levels; lev++) { pmesh.UniformRefinement(); }
+    mesh.Clear();
+
+    // Define a finite element space on the mesh. Here we use continuous Lagrange
+    // finite elements of the specified order. If order < 1, we fix it to 1.
+    if (order < 1) { order = 1; }
+    H1_FECollection fec(order, dim);
+    ParFiniteElementSpace sfespace(&pmesh, &fec, 1); //LS fe space
+    ParGridFunction lsf(&sfespace);
+    lsf.ProjectCoefficient(*ls_coeff);
+
+    ParFiniteElementSpace pfespace(&pmesh, &fec, dim);
+    pmesh.SetNodalFESpace(&pfespace);
+    ParGridFunction nod(&pfespace);
+    pmesh.SetNodalGridFunction(&nod);
+    ParGridFunction ggf(&pfespace);
+
+    ParNonlinearForm* nf=new ParNonlinearForm(&pfespace);
+
+    ElementMarker elmark(pmesh,true,true,3);
+    elmark.SetLevelSetFunction(lsf);
+    Array<int> el_markers;
+    elmark.MarkElements(el_markers);
+    CutIntegrationRules* cint=new CutIntegrationRules(3,lsf,el_markers);
+
+    CutVolLagrangeIntegrator* intgr=new CutVolLagrangeIntegrator(3);
+    intgr->SetCutIntegration(cint);
+
+    nf->AddDomainIntegrator(intgr);
+
+    double vol=nf->GetEnergy(nod.GetTrueVector());
+    if(myrank==0){ std::cout<<"vol="<<vol<<" err="<<vol-M_PI<<std::endl;}
+
+    Vector grd; grd.SetSize(pfespace.TrueVSize()); grd=0.0;
+    nf->Mult(nod.GetTrueVector(),grd); ggf.SetFromTrueDofs(grd);
+
+    //check the gradients by FD
+    {
+        Vector vtmp;
+        Vector prtv;
+        Vector tmpv;
+
+        prtv.SetSize(pfespace.TrueVSize());
+        tmpv.SetSize(pfespace.TrueVSize());
+        vtmp.SetSize(pfespace.TrueVSize());
+        nod.GetTrueDofs(vtmp);
+
+
+        prtv.Randomize();
+
+        double nd=mfem::InnerProduct(pmesh.GetComm(), prtv,prtv);
+        double td=mfem::InnerProduct(pmesh.GetComm(),prtv,grd);
+
+        td=td/nd;
+        double lsc=1.0;
+        double lqoi;
+
+        for(int l=0;l<10;l++){
+            lsc/=10.0;
+            prtv/=10.0;
+            add(prtv,vtmp,tmpv);
+
+            //calculate new integration rule
+            nod.SetFromTrueDofs(tmpv);
+            pmesh.DeleteGeometricFactors();
+            lqoi=nf->GetEnergy(tmpv);
+            double ld=(lqoi-vol)/lsc;
+            if(myrank==0){
+                std::cout << "dx=" << lsc <<" FD approximation=" << ld/nd
+                                     << " true gradient=" << td<<" r="<< ld/(td*nd)
+                                     << " err=" << std::fabs(ld/nd-td) << std::endl;
+            }
+        }
+
+    }
+
+    delete cint;
+    delete nf;
+    delete ls_coeff;
+
+    //dump out the data for ParaView
+    if(visualization)
+    {
+        mfem::ParaViewDataCollection paraview_dc("TopOpt", &pmesh);
+        paraview_dc.SetPrefixPath("ParaView");
+        paraview_dc.SetLevelsOfDetail(order);
+        paraview_dc.SetDataFormat(mfem::VTKFormat::BINARY);
+        paraview_dc.SetHighOrderOutput(true);
+        paraview_dc.SetCycle(0);
+        paraview_dc.SetTime(0.0);
+        paraview_dc.RegisterField("lsf",&lsf);
+        paraview_dc.RegisterField("grd",&ggf);
+        paraview_dc.RegisterField("nod",&nod);
+        paraview_dc.Save();
+    }
+
+}

miniapps/shifted/cut_integrators.cpp

@@ -0,0 +1,236 @@
+#include "cut_integrators.hpp"
+
+namespace mfem {
+
+void CutVolLagrangeIntegrator::AssembleElementVector(const FiniteElement &el,
+                           ElementTransformation &tr,
+                           const Vector &elfun, Vector &elvect)
+{
+    elvect.SetSize(elfun.Size());
+    elvect=0.0;
+
+    if(cint==nullptr){return;}
+
+    int elmark=cint->GetElementMarker(tr.ElementNo);
+
+    //only cut elements will have gradients
+    if(elmark!=ElementMarker::SBElementType::OUTSIDE){
+        int ndof = el.GetDof();
+        int ndim = tr.GetSpaceDim();
+
+        DenseMatrix bmat(ndof,ndim); //gradients of the shape functions in isoparametric space
+        DenseMatrix pmat(ndof,ndim);
+        Vector sh(ndof);
+        double f=1.0;
+
+        DenseMatrix gp; gp.Reset(elvect.GetData(),ndof,ndim);
+
+
+        const IntegrationRule* ir;
+        if(elmark==ElementMarker::SBElementType::CUT){
+            ir=cint->GetVolIntegrationRule(tr.ElementNo);
+        }else{
+            ir=&IntRules.Get(el.GetGeomType(),int_order);
+        }
+
+        for (int j = 0; j < ir->GetNPoints(); j++)
+        {
+            const IntegrationPoint &ip = ir->IntPoint(j);
+            tr.SetIntPoint(&ip);
+            el.CalcDShape(ip,bmat);
+            Mult(bmat,tr.AdjugateJacobian(),pmat);
+            if(coeff!=nullptr){f = coeff->Eval(tr,ip);}
+            gp.Add(ip.weight*f,pmat);
+
+        }
+    }
+}
+
+
+/// Compute the local energy, i.e., the volume
+double CutVolLagrangeIntegrator::GetElementEnergy(const FiniteElement &el,
+                            ElementTransformation &tr,
+                            const Vector &elfun)
+{
+    double rez=0.0;
+    if(cint==nullptr){return rez;}
+
+    int elmark=cint->GetElementMarker(tr.ElementNo);
+    if(elmark!=ElementMarker::SBElementType::OUTSIDE){
+
+        const IntegrationRule* ir;
+        if(elmark==ElementMarker::SBElementType::CUT){
+            ir=cint->GetVolIntegrationRule(tr.ElementNo);
+        }else{
+            ir=&IntRules.Get(el.GetGeomType(),int_order);
+        }
+
+        double f=1.0;
+        double w;
+        for (int j = 0; j < ir->GetNPoints(); j++)
+        {
+            const IntegrationPoint &ip = ir->IntPoint(j);
+            tr.SetIntPoint(&ip);
+            w = ip.weight ;
+            if(coeff!=nullptr){f = coeff->Eval(tr,ip);}
+
+            rez=rez+w*f*tr.Weight();
+
+            //std::cout<<"Tr.Weight="<<tr.Weight()<<std::endl;
+        }
+
+    }
+    return rez;
+}
+
+
+void VolGhostPenaltyIntegrator::AssembleFaceMatrix(const FiniteElement &fe1,
+                                                const FiniteElement &fe2,
+                                                FaceElementTransformations &Tr,
+                                                DenseMatrix &elmat)
+{
+    const int ndim=Tr.GetSpaceDim();
+    const int vdim=fe1.GetVDim();
+    int elem2 = Tr.Elem2No;
+    if(elem2<0){
+        elmat.SetSize(fe1.GetDof()*vdim);
+        elmat=0.0;
+        return;
+    }
+
+    const int ndof1 = fe1.GetDof();
+    const int ndof2 = fe2.GetDof();
+    const int ndofe = ndof1+ndof2;
+
+    elmat.SetSize(ndofe*vdim);
+    elmat=0.0;
+
+    int order=std::max(fe1.GetOrder(), fe2.GetOrder());
+
+    int ndofg;
+    if(ndim==1){ndofg=order+1;}
+    else if(ndim==2){ ndofg=(order+1)*(order+2)/2;}
+    else if(ndim==3){ ndofg=(order+1)*(order+2)*(order+3)/6;}
+
+    Vector sh1(ndof1);
+    Vector sh2(ndof2);
+    Vector shg(ndofg);
+
+    Vector xx(ndim);
+
+    DenseMatrix Mge(ndofg,ndofe); Mge=0.0;
+    DenseMatrix Mgg(ndofg,ndofg); Mgg=0.0;
+    DenseMatrix Mee(ndofe,ndofe); Mee=0.0;
+
+    ElementTransformation &Tr1 = Tr.GetElement1Transformation();
+    ElementTransformation &Tr2 = Tr.GetElement2Transformation();
+
+    //elements' volumes
+    double vol1=0.0;
+    double vol2=0.0;
+
+    const IntegrationRule* ir;
+
+    //element 1
+    double w;
+    ir=&IntRules.Get(Tr1.GetGeometryType(), 2*order+2);
+
+    for(int ii=0;ii<ir->GetNPoints();ii++){
+        const IntegrationPoint &ip = ir->IntPoint(ii);
+        Tr1.SetIntPoint(&ip);
+        Tr1.Transform(ip,xx);
+        fe1.CalcPhysShape(Tr1,sh1);
+        Shape(xx,ndim,order,shg);
+
+        w = Tr1.Weight();
+        w = ip.weight * w;
+        vol1=vol1+w;
+        for(int i=0;i<ndofg;i++){
+            for(int j=0;j<i;j++){
+                Mgg(i,j)=Mgg(i,j)+shg(i)*shg(j)*w;
+                Mgg(j,i)=Mgg(j,i)+shg(i)*shg(j)*w;
+            }
+            Mgg(i,i)=Mgg(i,i)+shg(i)*shg(i)*w;
+        }
+
+        for(int i=0;i<ndof1;i++){
+            for(int j=0;j<i;j++){
+                Mee(i,j)=Mee(i,j)+sh1(i)*sh1(j)*w;
+                Mee(j,i)=Mee(j,i)+sh1(i)*sh1(j)*w;
+            }
+            Mee(i,i)=Mee(i,i)+sh1(i)*sh1(i)*w;
+        }
+
+        for(int i=0;i<ndof1;i++){
+            for(int j=0;j<ndofg;j++){
+                Mge(j,i)=Mge(j,i)+shg(j)*sh1(i)*w;
+            }}
+    }
+
+    //element 2
+    ir=&IntRules.Get(Tr2.GetGeometryType(), 2*order+2);
+    for(int ii=0;ii<ir->GetNPoints();ii++){
+        const IntegrationPoint &ip = ir->IntPoint(ii);
+        Tr2.SetIntPoint(&ip);
+        Tr2.Transform(ip,xx);
+
+        fe2.CalcPhysShape(Tr2,sh2);
+        Shape(xx,ndim,order,shg);
+
+        w = Tr2.Weight();
+        w = ip.weight * w;
+        vol2=vol2+w;
+        for(int i=0;i<ndofg;i++){
+            for(int j=0;j<i;j++){
+                Mgg(i,j)=Mgg(i,j)+shg(i)*shg(j)*w;
+                Mgg(j,i)=Mgg(j,i)+shg(i)*shg(j)*w;
+            }
+            Mgg(i,i)=Mgg(i,i)+shg(i)*shg(i)*w;
+        }
+
+        for(int i=0;i<ndof2;i++){
+            for(int j=0;j<i;j++){
+                Mee(ndof1+i,ndof1+j)=Mee(ndof1+i,ndof1+j)+sh2(i)*sh2(j)*w;
+                Mee(ndof1+j,ndof1+i)=Mee(ndof1+j,ndof1+i)+sh2(i)*sh2(j)*w;
+            }
+            Mee(ndof1+i,ndof1+i)=Mee(ndof1+i,ndof1+i)+sh2(i)*sh2(i)*w;
+        }
+
+        for(int i=0;i<ndof2;i++){
+            for(int j=0;j<ndofg;j++){
+                Mge(j,ndof1+i)=Mge(j,ndof1+i)+shg(j)*sh2(i)*w;
+            }}
+    }
+
+
+    DenseMatrixInverse Mii(Mgg);
+    DenseMatrix Mre(ndofg,ndofe);
+    DenseMatrix Mff(ndofe,ndofe);
+    Mii.Mult(Mge,Mre);
+    MultAtB(Mge,Mre,Mff);
+
+
+    double hh;
+    if(ndim==1){
+        hh=(vol1+vol2)/2.0;
+    }else if(ndim==2){
+        hh=std::sqrt((vol1+vol2)/2.0);
+    }else{
+        hh=std::cbrt((vol1+vol2)/2.0);
+    }
+
+    double tv;
+
+    hh=penal/(hh*hh);
+
+    for(int i=0;i<ndofe;i++){
+        for(int j=0;j<ndofe;j++){
+            tv=hh*(Mee(i,j)+Mee(j,i)-Mff(i,j)-Mff(j,i))/(2.0);
+            for(int d=0;d<vdim;d++){
+                elmat(i+d*ndofe,j+d*ndofe)=tv;
+            }
+        }
+    }
+}
+
+}

miniapps/shifted/cut_integrators.hpp

@@ -0,0 +1,153 @@
+// Copyright (c) 2010-2023, 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.
+
+#ifndef MFEM_CUT_INTEGRATORS_HPP
+#define MFEM_CUT_INTEGRATORS_HPP
+
+#include "mfem.hpp"
+#include "marking.hpp"
+
+namespace mfem
+{
+
+/// Cut integrator for evaluating a volume and its gradients
+/// with respect to displacements of the mesh nodes
+class CutVolLagrangeIntegrator: public NonlinearFormIntegrator
+{
+
+public:
+    CutVolLagrangeIntegrator(int io=2)
+    {
+        coeff=nullptr;
+        cint=nullptr;
+        int_order=io;
+    }
+
+    virtual
+    ~CutVolLagrangeIntegrator()
+    {
+
+    }
+
+    void SetCutIntegration(CutIntegrationRules* cint_)
+    {
+        cint=cint_;
+    }
+
+    void SetCoefficient(Coefficient& cf){
+        coeff=&cf;
+    }
+
+    void SetIntOrder(int io=2)
+    {
+        int_order=io;
+    }
+
+    /// Perform the local action of the NonlinearFormIntegrator -
+    /// evaluates the gradients with respect to displacements of
+    /// the mesh nodes
+    virtual void AssembleElementVector(const FiniteElement &el,
+                                       ElementTransformation &Tr,
+                                       const Vector &elfun, Vector &elvect);
+
+
+    /// Compute the local energy, i.e., the volume
+    virtual double GetElementEnergy(const FiniteElement &el,
+                                        ElementTransformation &Tr,
+                                        const Vector &elfun);
+
+private:
+    CutIntegrationRules* cint; //cut integration rulle
+    Coefficient* coeff;
+    int int_order;
+};
+
+/// Ghost-penalty volume integrator
+class VolGhostPenaltyIntegrator:public BilinearFormIntegrator
+{
+public:
+    VolGhostPenaltyIntegrator(double penal_=1.0):penal(penal_)
+    {
+
+    }
+
+    virtual
+    ~VolGhostPenaltyIntegrator()
+    {
+
+    }
+
+    virtual
+    void AssembleFaceMatrix(const FiniteElement &fe1,
+                            const FiniteElement &fe2,
+                            FaceElementTransformations &Tr,
+                            DenseMatrix &elmat);
+
+
+
+protected:
+    double penal;
+
+    void Shape1D(double x, int order, Vector& sh){
+        sh.SetSize(order+1);
+        sh[0]=1.0;
+        for(int i=0;i<order;i++){sh[i+1]=sh[i]*x;}
+    }
+
+    void Shape2D(double x, double y, int order, Vector& sh)
+    {
+        Vector shx(order+1); Shape1D(x,order,shx);
+        Vector shy(order+1); Shape1D(y,order,shy);
+        sh.SetSize((order+1)*(order+2)/2);
+        int k=0;
+        for(int i=0;i<order+1;i++){
+        for(int j=0;j<order+1;j++){
+            if((i+j)<(order+1)){
+                sh[k]=shx[i]*shy[j];
+                k=k+1;
+            }
+        }}
+    }
+
+    void Shape3D(double x, double y, double z, int order, Vector& sh)
+    {
+        Vector shx(order+1); Shape1D(x,order,shx);
+        Vector shy(order+1); Shape1D(y,order,shy);
+        Vector shz(order+1); Shape1D(z,order,shz);
+        sh.SetSize((order+1)*(order+2)*(order+3)/6);
+        int p=0;
+        for(int i=0;i<order+1;i++){
+        for(int j=0;j<order+1;j++){
+        for(int k=0;k<order+1;k++){
+            if((i+j+k)<(order+1)){
+                sh[p]=shx[i]*shy[j]*shz[k];
+                p=p+1;
+            }
+        }}}
+    }
+
+    void Shape(Vector& xx, int dim, int order, Vector& sh)
+    {
+        if(dim==1){
+            Shape1D(xx[0],order,sh);
+        }else
+        if(dim==2){
+            Shape2D(xx[0],xx[1],order,sh);
+        }else
+        if(dim==3){
+            Shape3D(xx[0],xx[1],xx[2],order,sh);
+        }
+    }
+
+};
+
+}
+#endif

miniapps/shifted/integ_algoim.cpp

@@ -36,6 +36,139 @@ AlgoimIntegrationRule::AlgoimIntegrationRule(int o, const FiniteElement &el,
    T.Mult(lsfun,lsvec);
 }
 
+AlgoimIntegrationRule::AlgoimIntegrationRule(int o,ElementTransformation &trans,
+                                             Coefficient& lsfun, int lso)
+{
+    int_order=o;
+    vir=nullptr;
+    sir=nullptr;
+
+    if(trans.GetGeometryType()==Geometry::Type::SQUARE)
+    {
+        if(lso==-1){
+            pe=new H1Pos_QuadrilateralElement(trans.Order());
+        }else{
+            pe=new H1Pos_QuadrilateralElement(lso);
+        }
+    }
+    else if (trans.GetGeometryType()==Geometry::Type::CUBE)
+    {
+        if(lso==-1){
+            pe=new H1Pos_HexahedronElement(trans.Order());
+        }else{
+            pe=new H1Pos_HexahedronElement(lso);
+        }
+    }
+    else
+    {
+        MFEM_ABORT("Currently MFEM + Algoim supports only quads and hexes.");
+    }
+
+    //evaluate the level-set function
+    DenseMatrix M; M.SetSize(pe->GetDof()); M=0.0;
+    Vector rhs; rhs.SetSize(pe->GetDof()); rhs=0.0;
+    Vector shape; shape.SetSize(pe->GetDof()); shape=0.0;
+
+    const IntegrationRule* ir= nullptr;
+    int io=2*pe->GetOrder()+trans.OrderJ();
+    ir=&IntRules.Get(trans.GetGeometryType(),io);
+
+    lsvec.SetSize(pe->GetDof());
+
+    // form the mass matrix and the rhs
+    double w;
+    for(int i=0;i<ir->GetNPoints();i++){
+        const IntegrationPoint &ip = ir->IntPoint(i);
+        trans.SetIntPoint(&ip);
+
+        w=trans.Weight();
+        w = ip.weight * w;
+
+        pe->CalcPhysShape(trans,shape);
+
+        for(int ii=0;ii<pe->GetDof();ii++){
+                rhs[ii]=rhs[ii]+w*shape[ii]*lsfun.Eval(trans,ip);
+                M(ii,ii)=M(ii,ii)+w*shape[ii]*shape[ii];
+            for(int jj=0;jj<ii;jj++){
+                M(ii,jj)=M(ii,jj)+w*shape[ii]*shape[jj];
+        }}
+    }
+
+    for(int ii=0;ii<pe->GetDof();ii++){
+        for(int jj=0;jj<ii;jj++){
+            M(jj,ii)=M(ii,jj);
+        }
+    }
+
+    DenseMatrixInverse im(&M,true);
+    im.Factor();
+    im.Mult(rhs,lsvec);
+}
+
+AlgoimIntegrationRule::   AlgoimIntegrationRule(int o, ElementTransformation &trans,
+                                                GridFunction& lsfun)
+{
+
+    int lso= lsfun.FESpace()->GetElementOrder(trans.ElementNo);
+
+    int_order=o;
+    vir=nullptr;
+    sir=nullptr;
+
+    if(trans.GetGeometryType()==Geometry::Type::SQUARE)
+    {
+        pe=new H1Pos_QuadrilateralElement(lso);
+    }
+    else if (trans.GetGeometryType()==Geometry::Type::CUBE)
+    {
+        pe=new H1Pos_HexahedronElement(lso);
+    }
+    else
+    {
+        MFEM_ABORT("Currently MFEM + Algoim supports only quads and hexes.");
+    }
+
+    //evaluate the level-set function
+    DenseMatrix M; M.SetSize(pe->GetDof()); M=0.0;
+    Vector rhs; rhs.SetSize(pe->GetDof()); rhs=0.0;
+    Vector shape; shape.SetSize(pe->GetDof()); shape=0.0;
+
+    const IntegrationRule* ir= nullptr;
+    int io=2*pe->GetOrder()+trans.OrderJ();
+    ir=&IntRules.Get(trans.GetGeometryType(),io);
+
+    lsvec.SetSize(pe->GetDof());
+
+    // form the mass matrix and the rhs
+    double w;
+    for(int i=0;i<ir->GetNPoints();i++){
+        const IntegrationPoint &ip = ir->IntPoint(i);
+        trans.SetIntPoint(&ip);
+
+        w=trans.Weight();
+        w = ip.weight * w;
+
+        pe->CalcPhysShape(trans,shape);
+
+        for(int ii=0;ii<pe->GetDof();ii++){
+                rhs[ii]=rhs[ii]+w*shape[ii]*lsfun.GetValue(trans,ip);
+                M(ii,ii)=M(ii,ii)+w*shape[ii]*shape[ii];
+            for(int jj=0;jj<ii;jj++){
+                M(ii,jj)=M(ii,jj)+w*shape[ii]*shape[jj];
+        }}
+    }
+
+    for(int ii=0;ii<pe->GetDof();ii++){
+        for(int jj=0;jj<ii;jj++){
+            M(jj,ii)=M(ii,jj);
+        }
+    }
+
+    DenseMatrixInverse im(&M,true);
+    im.Factor();
+    im.Mult(rhs,lsvec);
+}
+
 
 const IntegrationRule* AlgoimIntegrationRule::GetVolumeIntegrationRule()
 {
@@ -113,6 +246,63 @@ const IntegrationRule* AlgoimIntegrationRule::GetSurfaceIntegrationRule()
    return sir;
 }
 
+
+double AlgoimIntegrationRule::SurfaceWeight(ElementTransformation &trans,
+                                            const IntegrationPoint& ip)
+{
+    DenseMatrix bmat; //gradients of the shape functions in isoparametric space
+    DenseMatrix pmat; //gradients of the shape functions in physical space
+    Vector inormal; //normal to the level set in isoparametric space
+    Vector tnormal; //normal to the level set in physical space
+
+    bmat.SetSize(pe->GetDof(),pe->GetDim());
+    pmat.SetSize(pe->GetDof(),pe->GetDim());
+    inormal.SetSize(pe->GetDim());
+    tnormal.SetSize(pe->GetDim());
+
+    pe->CalcDShape(ip,bmat);
+    Mult(bmat, trans.AdjugateJacobian(), pmat);
+    //compute the normal to the LS in isoparametric space
+    bmat.MultTranspose(lsvec,inormal);
+    //compute the normal to the LS in physical space
+    pmat.MultTranspose(lsvec,tnormal);
+
+    return tnormal.Norml2() / inormal.Norml2();
+}
+
+const Array<double>* AlgoimIntegrationRule::GetSurfaceWeights(
+                                            ElementTransformation &trans)
+{
+    const IntegrationRule* lsir=GetSurfaceIntegrationRule();
+
+    sweight.SetSize(lsir->GetNPoints());
+
+    DenseMatrix bmat; //gradients of the shape functions in isoparametric space
+    DenseMatrix pmat; //gradients of the shape functions in physical space
+    Vector inormal; //normal to the level set in isoparametric space
+    Vector tnormal; //normal to the level set in physical space
+
+    bmat.SetSize(pe->GetDof(),pe->GetDim());
+    pmat.SetSize(pe->GetDof(),pe->GetDim());
+    inormal.SetSize(pe->GetDim());
+    tnormal.SetSize(pe->GetDim());
+
+    for(int i=0;i<lsir->GetNPoints();i++){
+        const IntegrationPoint &ip = lsir->IntPoint(i);
+        trans.SetIntPoint(&ip);
+        pe->CalcDShape(ip,bmat);
+        Mult(bmat, trans.AdjugateJacobian(), pmat);
+        //compute the normal to the LS in isoparametric space
+        bmat.MultTranspose(lsvec,inormal);
+        //compute the normal to the LS in physical space
+        pmat.MultTranspose(lsvec,tnormal);
+        sweight[i]=tnormal.Norml2() / inormal.Norml2();
+    }
+
+    return &sweight;
+}
+
+
 }
 
 #endif

miniapps/shifted/integ_algoim.hpp

@@ -117,6 +117,21 @@ public:
                          ElementTransformation &trans, const Vector &lsfun);
 
 
+   /// Construct Algoim object using  a specified integration order o,
+   /// an element transformation, level set function supplied in the form
+   /// of a coefficient and approximation order of the level-set. If the
+   /// approximation order is set to -1, the order will be set to be equal
+   /// to the element order from the element transformation.
+   AlgoimIntegrationRule(int o, ElementTransformation &trans,
+                         Coefficient& lsfun, int lso=-1);
+
+   /// Construct Algoim object using  a specified integration order o,
+   /// an element transformation, level set function supplied in the form
+   /// of a grid function.
+   AlgoimIntegrationRule(int o, ElementTransformation &trans,
+                         GridFunction& lsfun);
+
+
    /// Destructor of the Algoim object
    ~AlgoimIntegrationRule()
    {
@@ -133,6 +148,18 @@ public:
    /// level-set function.
    const IntegrationRule* GetSurfaceIntegrationRule();
 
+   /// Returns the surface weight for surface integration at
+   /// integration point ip. The total weight should be computed
+   /// as the product of the integration point weight ip.weight
+   /// and the result returned by SurfaceWeight.
+   double SurfaceWeight(ElementTransformation &trans,
+                        const IntegrationPoint& ip);
+
+
+   /// Returns the surface weights for the current surface integration
+   /// rule. The weights are recomputed for every call to the method.
+   const Array<double>* GetSurfaceWeights(ElementTransformation &trans);
+
 
 private:
 
@@ -269,11 +296,16 @@ private:
    };
 
 
+   Array<double> sweight; // Surface weight evaluated for
+                          // every surface integration point.
+
    IntegrationRule* sir; // Surface integration rule. Owned.
    IntegrationRule* vir; // Volumetric integration rule. Owned.
    PositiveTensorFiniteElement *pe;
    Vector lsvec; //level-set in Bernstein bases
    int int_order; //integration order
+
+
 };
 
 }

miniapps/shifted/lsf_integral.cpp

@@ -182,26 +182,14 @@ int main(int argc, char *argv[])
 
 #ifdef MFEM_USE_ALGOIM
    double area=0.0;
-   DenseMatrix bmat; //gradients of the shape functions in isoparametric space
-   DenseMatrix pmat; //gradients of the shape functions in physical space
-   Vector inormal; //normal to the level set in isoparametric space
-   Vector tnormal; //normal to the level set in physical space
-   Vector lsfun; //level set function restricted to an element
-   DofTransformation *doftrans;
-   Array<int> vdofs;
    for (int i=0; i<fespace.GetNE(); i++)
    {
-      const FiniteElement* el=fespace.GetFE(i);
-
       //get the element transformation
       trans = fespace.GetElementTransformation(i);
 
-      //extract the element vector from the level-set
-      doftrans = fespace.GetElementVDofs(i,vdofs);
-      x.GetSubVector(vdofs, lsfun);
-
       //construct Algoim integration object
-      AlgoimIntegrationRule air(aorder, *el, *trans, lsfun);
+      AlgoimIntegrationRule air(aorder, *trans, *ls_coeff, 3);
+
 
       //compute the volume contribution from the element
       ir = air.GetVolumeIntegrationRule();
@@ -212,25 +200,12 @@ int main(int argc, char *argv[])
          vol += ip.weight * trans->Weight();
       }
 
-      //compute the perimeter/area contribution from the element
-      bmat.SetSize(el->GetDof(),el->GetDim());
-      pmat.SetSize(el->GetDof(),el->GetDim());
-      inormal.SetSize(el->GetDim());
-      tnormal.SetSize(el->GetDim());
-
       ir = air.GetSurfaceIntegrationRule();
+      const Array<double>& sweights= *air.GetSurfaceWeights(*trans);
       for (int j = 0; j < ir->GetNPoints(); j++)
       {
          const IntegrationPoint &ip = ir->IntPoint(j);
-         trans->SetIntPoint(&ip);
-
-         el->CalcDShape(ip,bmat);
-         Mult(bmat, trans->AdjugateJacobian(), pmat);
-         //compute the normal to the LS in isoparametric space
-         bmat.MultTranspose(lsfun,inormal);
-         //compute the normal to the LS in physical space
-         pmat.MultTranspose(lsfun,tnormal);
-         area += ip.weight * tnormal.Norml2() / inormal.Norml2();
+         area += ip.weight * sweights[j];
       }
    }
 
@@ -245,7 +220,6 @@ int main(int argc, char *argv[])
       std::cout<<"Algoim Area="<<area<<std::endl;
    }
 #endif
-
    //Perform standard MFEM integration
    vol=0.0;
    for (int i=0; i<fespace.GetNE(); i++)

miniapps/shifted/marking.cpp

@@ -10,6 +10,7 @@
 // CONTRIBUTING.md for details.
 
 #include "marking.hpp"
+#include "integ_algoim.hpp"
 
 namespace mfem
 {
@@ -351,4 +352,354 @@ void ShiftedFaceMarker::ListShiftedFaceDofs2(const Array<int> &elem_marker,
    }
 }
 
+void ElementMarker::SetLevelSetFunction(Coefficient &ls_fun)
+{
+    FiniteElementCollection* fec=new H1_FECollection(h1_order,pmesh->Dimension());
+    ParFiniteElementSpace* pfes_sltn=new ParFiniteElementSpace(pmesh,fec);
+
+    Vector vals;
+    Array<int> vdofs;
+
+    if(use_cut_marks==false){
+    if(include_cut_elements){
+        elgf=(double)(SBElementType::INSIDE);
+        for(int e=0;e<pmesh->GetNE();e++){
+            const IntegrationRule &ir = pfes_sltn->GetFE(e)->GetNodes();
+            {
+                int n = ir.GetNPoints();
+                vals.SetSize(n);
+                ElementTransformation *Tr = pfes_sltn->GetElementTransformation(e);
+                for(int k=0;k<n;k++){
+                    Tr->SetIntPoint(&ir.IntPoint(k));
+                    vals[k]=ls_fun.Eval(*Tr,ir.IntPoint(k));
+                }
+            }
+
+            int countp = 0;
+            int countn = 0;
+            for (int j = 0; j < ir.GetNPoints(); j++){
+                if (vals(j)>0.0) { countp++; }
+                else{countn++;}
+            }
+            if (countn == ir.GetNPoints()) // completely outside
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::OUTSIDE;
+            }
+        }
+    }else{//DEFAULT - do not include cuts
+        elgf=(double)(SBElementType::OUTSIDE);
+        for(int e=0;e<pmesh->GetNE();e++){
+            const IntegrationRule &ir = pfes_sltn->GetFE(e)->GetNodes();
+            {
+                int n = ir.GetNPoints();
+                vals.SetSize(n);
+                ElementTransformation *Tr = pfes_sltn->GetElementTransformation(e);
+                for(int k=0;k<n;k++){
+                    Tr->SetIntPoint(&ir.IntPoint(k));
+                    vals[k]=ls_fun.Eval(*Tr,ir.IntPoint(k));
+                }
+            }
+
+            int countp = 0;
+            int countn = 0;
+            for (int j = 0; j < ir.GetNPoints(); j++){
+                if (vals(j)>0.0) { countp++; }
+                else{countn++;}
+            }
+            if (countp == ir.GetNPoints()) // completely inside
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::INSIDE;
+            }
+        }
+    }}else{// use CUT mark
+        elgf=(double)(SBElementType::INSIDE);
+        for(int e=0;e<pmesh->GetNE();e++){
+            const IntegrationRule &ir = pfes_sltn->GetFE(e)->GetNodes();
+            {
+                int n = ir.GetNPoints();
+                vals.SetSize(n);
+                ElementTransformation *Tr = pfes_sltn->GetElementTransformation(e);
+                for(int k=0;k<n;k++){
+                    Tr->SetIntPoint(&ir.IntPoint(k));
+                    vals[k]=ls_fun.Eval(*Tr,ir.IntPoint(k));
+                }
+            }
+
+            int countp = 0;
+            int countn = 0;
+            for (int j = 0; j < ir.GetNPoints(); j++){
+                if (vals(j)>0) {countp++;}
+                else {countn++;}
+            }
+            if (countn == ir.GetNPoints()) // completely outside
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::OUTSIDE;
+            }else
+            if ((countp>0)&&(countn>0))
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::CUT;
+            }
+        }
+    }
+
+    elgf.ExchangeFaceNbrData();
+
+    delete pfes_sltn;
+    delete fec;
+
+}
+
+void ElementMarker::SetLevelSetFunction(const ParGridFunction &ls_fun)
+{
+    ParFiniteElementSpace* pfes_sltn=ls_fun.ParFESpace();
+    Vector vals;
+    Array<int> vdofs;
+
+    if(use_cut_marks==false){
+    if(include_cut_elements){
+        elgf=(double)(SBElementType::INSIDE);
+        for(int e=0;e<pmesh->GetNE();e++){
+            const IntegrationRule &ir = pfes_sltn->GetFE(e)->GetNodes();
+            ls_fun.GetValues(e, ir, vals);
+            int countn = 0;
+            for (int j = 0; j < ir.GetNPoints(); j++){
+                if (vals(j)>0.0) {}
+                else{ countn++; }
+            }
+            if (countn == ir.GetNPoints()) // completely outside
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::OUTSIDE;
+            }
+        }
+    }else{//DEFAULT - do not include cuts
+        elgf=(double)(SBElementType::OUTSIDE);
+        for(int e=0;e<pmesh->GetNE();e++){
+            const IntegrationRule &ir = pfes_sltn->GetFE(e)->GetNodes();
+            ls_fun.GetValues(e, ir, vals);
+            int countp = 0;
+            for (int j = 0; j < ir.GetNPoints(); j++){
+                if (vals(j)>0.0) { countp++; }
+            }
+            if (countp == ir.GetNPoints()) // completely inside
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::INSIDE;
+            }
+        }
+    }}else{//use CUT marks
+        elgf=(double)(SBElementType::INSIDE);
+        for(int e=0;e<pmesh->GetNE();e++){
+            const IntegrationRule &ir = pfes_sltn->GetFE(e)->GetNodes();
+            ls_fun.GetValues(e, ir, vals);
+            int countp = 0;
+            int countn = 0;
+            for (int j = 0; j < ir.GetNPoints(); j++){
+                if (vals(j)>0) {countp++;}
+                else {countn++;}
+            }
+            if (countn == ir.GetNPoints()) // completely outside
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::OUTSIDE;
+            }else
+            if ((countp>0)&&(countn>0))
+            {
+                elfes->GetElementVDofs(e,vdofs);
+                elgf[vdofs[0]] = SBElementType::CUT;
+            }
+        }
+    }
+    if(pmesh->GetNRanks()>0){
+    elgf.ExchangeFaceNbrData();}
+
+}
+
+void ElementMarker::MarkElements(Array<int> &elem_marker)
+{
+    elem_marker.SetSize(pmesh->GetNE());
+    for(int e=0;e<pmesh->GetNE();e++)
+    {
+        ElementTransformation* tr=elfes->GetElementTransformation(e);
+        IntegrationPoint ip; ip.Init(0);
+        elem_marker[e] = elgf.GetValue(*tr, ip);
+    }
+}
+
+void ElementMarker::MarkGhostPenaltyFaces(Array<int> &face_marker)
+{
+    face_marker.SetSize(pmesh->GetNumFaces());
+    face_marker=SBFaceType::UNDEFINED;
+    IntegrationPoint ip; ip.Init(0);
+
+    for(int f=0;f<pmesh->GetNumFaces();f++){
+        auto *ft = pmesh->GetFaceElementTransformations(f, 3);
+        if (ft->Elem2No < 0) { continue; } //do not mark boundary faces
+        const int attr1 = elgf.GetValue(*ft->Elem1,ip);
+        const int attr2 = elgf.GetValue(*ft->Elem2,ip);
+        if((attr1==SBElementType::CUT)&&(attr2!=SBElementType::OUTSIDE))
+        {
+            face_marker[f]=SBFaceType::GHOSTP;
+        }else
+        if((attr1!=SBElementType::OUTSIDE)&&(attr2==SBElementType::CUT))
+        {
+            face_marker[f]=SBFaceType::GHOSTP;
+        }
+    }
+
+    elgf.ExchangeFaceNbrData();
+
+    for (int f = 0; f < pmesh->GetNSharedFaces(); f++)
+    {
+        auto *ftr = pmesh->GetSharedFaceTransformations(f, true);
+        const int attr1 = elgf.GetValue(*ftr->Elem1, ip);
+        const int attr2 = elgf.GetValue(*ftr->Elem2, ip);
+        int faceno = pmesh->GetSharedFace(f);
+        if((attr1==SBElementType::CUT)&&(attr2!=SBElementType::OUTSIDE))
+        {
+            face_marker[faceno]=SBFaceType::GHOSTP;
+        }else
+        if((attr1!=SBElementType::OUTSIDE)&&(attr2==SBElementType::CUT))
+        {
+            face_marker[faceno]=SBFaceType::GHOSTP;
+        }
+    }
+}
+
+void ElementMarker::MarkFaces(Array<int> &face_marker)
+{
+    face_marker.SetSize(pmesh->GetNumFaces());
+    face_marker=SBFaceType::UNDEFINED;
+    IntegrationPoint ip; ip.Init(0);
+
+    if(include_cut_elements==true){
+        for(int f=0;f<pmesh->GetNumFaces();f++){
+            auto *ft = pmesh->GetFaceElementTransformations(f, 3);
+            if (ft->Elem2No < 0) { continue; } //do not mark boundary faces
+            const int attr1 = elgf.GetValue(*ft->Elem1,ip);
+            const int attr2 = elgf.GetValue(*ft->Elem2,ip);
+            if((attr1==SBElementType::OUTSIDE)||(attr2==SBElementType::OUTSIDE)){
+                if(attr1!=attr2){
+                    face_marker[f]=SBFaceType::SURROGATE;
+                }
+            }
+        }
+    }else{
+        for(int f=0;f<pmesh->GetNumFaces();f++){
+            auto *ft = pmesh->GetFaceElementTransformations(f, 3);
+            if (ft->Elem2No < 0) { continue; } //do not mark boundary faces
+            const int attr1 = elgf.GetValue(*ft->Elem1,ip);
+            const int attr2 = elgf.GetValue(*ft->Elem2,ip);
+            if((attr1==SBElementType::INSIDE)||(attr2==SBElementType::INSIDE)){
+                if(attr1!=attr2){
+                    face_marker[f]=SBFaceType::SURROGATE;
+                }
+            }
+        }
+    }
+
+    elgf.ExchangeFaceNbrData();
+    if(include_cut_elements==true){
+        for (int f = 0; f < pmesh->GetNSharedFaces(); f++)
+        {
+            auto *ftr = pmesh->GetSharedFaceTransformations(f, true);
+            const int attr1 = elgf.GetValue(*ftr->Elem1, ip);
+            const int attr2 = elgf.GetValue(*ftr->Elem2, ip);
+            int faceno = pmesh->GetSharedFace(f);
+            if((attr1==SBElementType::OUTSIDE)||(attr2==SBElementType::OUTSIDE)){
+                if(attr1!=attr2){
+                    face_marker[faceno]=SBFaceType::SURROGATE;
+                }
+            }
+        }
+    }else{
+        for (int f = 0; f < pmesh->GetNSharedFaces(); f++)
+        {
+            auto *ftr = pmesh->GetSharedFaceTransformations(f, true);
+            const int attr1 = elgf.GetValue(*ftr->Elem1, ip);
+            const int attr2 = elgf.GetValue(*ftr->Elem2, ip);
+            int faceno = pmesh->GetSharedFace(f);
+            if((attr1==SBElementType::INSIDE)||(attr2==SBElementType::INSIDE)){
+                if(attr1!=attr2){
+                    face_marker[faceno]=SBFaceType::SURROGATE;
+                }
+            }
+        }
+
+    }
+}
+
+void ElementMarker::ListEssentialTDofs(const Array<int> &elem_marker,
+                                       ParFiniteElementSpace &lfes,
+                                       Array<int> &ess_tdof_list) const
+{
+    Array<int> dofs;
+
+    mfem::Vector vvdof; vvdof.SetSize(lfes.GetVSize()); vvdof=0.0;
+
+    for(int i=0;i<lfes.GetNE();i++)
+    {
+        if(elem_marker[i]==SBElementType::INSIDE){
+            lfes.GetElementVDofs(i,dofs);
+            for(int j=0;j<dofs.Size();j++){
+                vvdof[dofs[j]]=1.0;
+            }
+        }
+
+        if(include_cut_elements==true){
+            if(elem_marker[i]==SBElementType::CUT){
+                lfes.GetElementVDofs(i,dofs);
+                for(int j=0;j<dofs.Size();j++){
+                    vvdof[dofs[j]]=1.0;
+                }
+            }
+        }
+    }
+
+    Array<int> tdof_mark; tdof_mark.SetSize(lfes.GetTrueVSize());
+    Vector vtdof; vtdof.SetSize(lfes.GetTrueVSize()); vtdof=0.0;
+    lfes.GetProlongationMatrix()->MultTranspose(vvdof,vtdof);
+    for(int i=0;i<vtdof.Size();i++){
+        if(vtdof[i]<1.0){tdof_mark[i]=1;}
+        else{tdof_mark[i]=0;}
+    }
+
+    lfes.MarkerToList(tdof_mark, ess_tdof_list);
+}
+
+#ifdef MFEM_USE_ALGOIM
+
+CutIntegrationRules::CutIntegrationRules(int int_order,
+                                         GridFunction& lsf,
+                                         Array<int>& elm_markers):markers(elm_markers)
+{
+
+     FiniteElementSpace* fespace=lsf.FESpace();
+     ElementTransformation *Tr;
+
+     air.SetSize(fespace->GetNE()); air=nullptr;
+
+     for (int i=0; i<fespace->GetNE(); i++){
+     if(elm_markers[i]==ElementMarker::SBElementType::CUT){
+         Tr=fespace->GetElementTransformation(i);
+         air[i]=new AlgoimIntegrationRule(int_order,*Tr,lsf);
+     }}
+
+}
+
+CutIntegrationRules::~CutIntegrationRules()
+{
+    for(int i=0;i<air.Size();i++)
+    {
+        delete air[i];
+    }
+}
+
+#endif
+
+
 }

miniapps/shifted/marking.hpp

@@ -13,6 +13,7 @@
 #define MFEM_MARKING_HPP
 
 #include "mfem.hpp"
+#include "integ_algoim.hpp"
 
 namespace mfem
 {
@@ -77,6 +78,135 @@ public:
                            Array<int> &ess_shift_bdr) const;
 };
 
+
+
+// Marking operations for elements, faces, dofs, etc, related to cut, shifted
+// boundary and interface methods.
+class ElementMarker{
+public:
+    enum SBElementType {INSIDE = 0, OUTSIDE = 1, CUT = 2};
+
+    enum SBFaceType {UNDEFINED = 0, SURROGATE = 1, GHOSTP = 2};
+
+    ///Defines element marker class with options to include the cut elements
+    /// (include_cut=true) or to mark the cut elements as SBElementType::CUT.
+    /// If use_cut=false the marking will use only INSIDE/OUTSIDE marks.
+    /// The last integer argument determines the order of the surrogate H1 field
+    /// for checking if an element is cut by a zero level set of an implicit
+    /// material distribution.
+    ElementMarker(ParMesh& mesh, bool include_cut=false,
+                  bool use_cut=false, int h1_order_=2)
+    {
+        pmesh=&mesh;
+        const int dim=pmesh->SpaceDimension();
+        elfec=new L2_FECollection(0,dim);
+        elfes=new ParFiniteElementSpace(pmesh,elfec,1);
+        elgf.SetSpace(elfes);
+        include_cut_elements=include_cut;
+        use_cut_marks=use_cut;
+
+        h1_order=h1_order_;
+    }
+
+    /// Destructor of the ElementMarker class
+    ~ElementMarker()
+    {
+        delete elfes;
+        delete elfec;
+    }
+
+    /// Mark elements according to the specified level-set
+    /// function.
+    void SetLevelSetFunction(const ParGridFunction& ls_fun);
+
+    /// Mark the elements according to the specified coefficient.
+    void SetLevelSetFunction(Coefficient& ls_fun);
+
+    /// Returns the marking of all the elements
+    /// in the mesh using the @a SBElementType
+    void MarkElements(Array<int> &elem_marker);
+
+    /// Returns the marking of all faces in the
+    /// mesh using  the @a SBFaceType
+    void MarkFaces(Array<int> &face_marker);
+
+    /// Returns the marking of all faces in the
+    /// mesh using  the @a SBFaceType.
+    /// The marks of all cut and faces between
+    /// cut and inside elements are set to GHOSTP
+    void MarkGhostPenaltyFaces(Array<int> &face_marker);
+
+    /// Lists all inactive dofs, i.e.,
+    ///  all dofs in the outside region.
+    void ListEssentialTDofs(const Array<int> &elem_marker,
+                            ParFiniteElementSpace &lfes,
+                            Array<int> &ess_tdof_list) const;
+private:
+    ParMesh* pmesh;
+    FiniteElementCollection* elfec;
+    ParFiniteElementSpace* elfes;
+    ParGridFunction elgf;
+
+    bool include_cut_elements;
+    bool use_cut_marks;
+
+    int h1_order; //order of the H1 FE space for level set functions defined by coefficient
+};
+
+#ifdef MFEM_USE_ALGOIM
+
+/// The class generates volumetric and surface
+/// body fitted integration rules for all zero level-set cut elements.
+class CutIntegrationRules
+{
+public:
+
+    /// Contructs the set of integration rules for a given grid function and a set
+    ///  of markers.
+    CutIntegrationRules(int int_order, GridFunction& lsf, Array<int>& elm_markers);
+
+    /// Frees the integration rules
+    ~CutIntegrationRules();
+
+
+    /// Returns the volumetric integration rule for element el.
+    const IntegrationRule* GetVolIntegrationRule(int el){
+        if(air[el]==nullptr){return nullptr;}
+        else{
+            return air[el]->GetVolumeIntegrationRule();
+        }
+    }
+
+    /// Returns the marker associated with element el.
+    ElementMarker::SBElementType GetElementMarker(int el){
+        return ElementMarker::SBElementType(markers[el]);
+    }
+
+    /// Returns surface integration rule for element el.
+    const IntegrationRule* GetSurfIntegrationRule(int el)
+    {
+        if(air[el]==nullptr){return nullptr;}
+        else{
+            return air[el]->GetSurfaceIntegrationRule();
+        }
+    }
+
+    /// Returns surface integration weights for element el
+    const Array<double>* GetSurfaceWeights(ElementTransformation& trans)
+    {
+        if(air[trans.ElementNo]==nullptr){return nullptr;}
+        else{
+            return air[trans.ElementNo]->GetSurfaceWeights(trans);
+        }
+    }
+
+private:
+    Array<AlgoimIntegrationRule*> air; //holds all Algoim integration rules
+    Array<int>& markers; //holds a reference to the markers
+};
+#endif
+
+
 } // namespace mfem
 
 #endif