3D test diffusion

fem/intrules_cut.cpp

@@ -140,6 +140,32 @@ void AlgoimIntegrationRules::GetSurfaceWeights(ElementTransformation &Tr,
 
 }
 
+void AlgoimIntegrationRules::GetSurfaceNormal(ElementTransformation &Tr,
+                                               const IntegrationRule &sir,
+                                               DenseMatrix &normal)
+{
+   GenerateLSVector(Tr,LvlSet);
+   DenseMatrix pmat; // gradients of the shape functions in physcial space 
+   Vector tnormal; // normal to the level set in physcial space 
+   pmat.SetSize(pe->GetDof(),pe->GetDim());
+   tnormal.SetSize(pe->GetDim());
+   normal.SetSize(sir.GetNPoints(),pe->GetDim());
+
+   for (int j = 0; j < sir.GetNPoints(); j++)
+   {
+      const IntegrationPoint &ip = sir.IntPoint(j);
+      Tr.SetIntPoint(&ip);
+      pe->CalcPhysDShape(Tr,pmat);
+      // compute the normal to the LS in physcial space 
+      pmat.MultTranspose(lsvec,tnormal);
+      for (int k = 0; k<pe->GetDim(); k++)
+      {
+         normal(j,k) = -tnormal(k)/tnormal.Norml2();
+      }
+   }
+
+}
+
 void AlgoimIntegrationRules::GenerateLSVector(ElementTransformation &Tr,
                                               Coefficient* lvlset)
 {

fem/intrules_cut.hpp

@@ -317,6 +317,20 @@ public:
                           const IntegrationRule &sir,
                           Vector &weights) override;
 
+   /**
+   @brief Compute normal vectors for surface integration.
+
+
+   @param [in]  Tr      Specifies the IntegrationRule's associated element.
+   @param [in]  sir     IntegrationRule defining the IntegrationPoints
+   @param [out] normal  Matrix where each row contain the normal for the corrsponding 
+   integration point.
+   */
+
+   void GetSurfaceNormal(ElementTransformation &Tr,
+                          const IntegrationRule &sir,
+                          DenseMatrix &normal);
+
 private:
 
    /// projects the lvlset coefficient onto the lsvec,

miniapps/CMakeLists.txt

@@ -37,3 +37,4 @@ add_subdirectory(tribol)
 add_subdirectory(hooke)
 add_subdirectory(dpg)
 add_subdirectory(hdiv-linear-solver)
+add_subdirectory(cutfem)

miniapps/cutfem/CMakeLists.txt

@@ -0,0 +1,93 @@
+# 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.
+
+list(APPEND SEQCUT_COMMON_SOURCES
+	cut_marking.cpp
+	my_integrators.cpp)
+
+list(APPEND SEQCUT_COMMON_HEADERS
+	cut_marking.hpp
+	my_integrators.hpp)
+
+convert_filenames_to_full_paths(SEQCUT_COMMON_SOURCES)
+convert_filenames_to_full_paths(SEQCUT_COMMON_HEADERS)
+
+set(SEQCUT_COMMON_FILES
+	EXTRA_SOURCES ${SEQCUT_COMMON_SOURCES}
+	EXTRA_HEADERS ${SEQCUT_COMMON_HEADERS})
+
+add_mfem_miniapp(seqdiff
+	MAIN diffusion.cpp
+	${SEQCUT_COMMON_FILES}
+	LIBRARIES mfem)
+
+add_mfem_miniapp(seqdiffn
+	MAIN diffusion_nitsche.cpp
+	${SEQCUT_COMMON_FILES}
+        LIBRARIES mfem)
+
+add_mfem_miniapp(seqstokes
+	MAIN stokes.cpp
+	${SEQCUT_COMMON_FILES}
+        LIBRARIES mfem)
+
+add_mfem_miniapp(seq_mark_test
+        MAIN seq_mark_test.cpp
+        ${SEQCUT_COMMON_FILES}
+        LIBRARIES mfem)
+
+add_mfem_miniapp(seq_diffusion_cut
+        MAIN diffusion_cut.cpp
+        ${SEQCUT_COMMON_FILES}
+        LIBRARIES mfem)
+
+add_mfem_miniapp(seq_stokes_cut
+	MAIN stokes_cut.cpp
+	${SEQCUT_COMMON_FILES}
+	LIBRARIES mfem)
+
+add_mfem_miniapp(seq_vector
+	MAIN vector_diffusion.cpp
+	${SEQCUT_COMMON_FILES}
+	LIBRARIES mfem)
+
+add_mfem_miniapp(seq_vector_cut
+	MAIN vector_diffusion_cut.cpp
+	${SEQCUT_COMMON_FILES}
+	LIBRARIES mfem)
+
+if(MFEM_USE_MPI)
+
+list(APPEND PARCUT_COMMON_SOURCES)
+list(APPEND PARCUT_COMMON_HEADERS)
+
+convert_filenames_to_full_paths(PARCUT_COMMON_SOURCES)
+convert_filenames_to_full_paths(PARCUT_COMMON_HEADERS)
+
+set(PARCUT_COMMON_FILES
+	EXTRA_SOURCES ${PARCUT_COMMON_SOURCES} ${SEQCUT_COMMON_SOURCES}
+	EXTRA_HEADERS ${PARCUT_COMMON_HEADERS} ${SEQCUT_COMMON_HEADERS})
+
+# message(STATUS "PARCUT_COMMON_FILES: ${PARCUT_COMMON_FILES}")
+# message(STATUS "SEQCUT_COMMON_FILES: ${SEQCUT_COMMON_FILES}")
+
+add_mfem_miniapp(par_diffusion_cut
+        MAIN par_diffusion_cut.cpp
+        ${PARCUT_COMMON_FILES}
+        LIBRARIES mfem)
+
+add_mfem_miniapp(par_diffusion_cut_3D
+	MAIN par_diffusion_cut_3D.cpp
+        ${PARCUT_COMMON_FILES}
+        LIBRARIES mfem)
+
+
+endif ()

miniapps/cutfem/cut_marking.cpp

@@ -0,0 +1,614 @@
+#include "cut_marking.hpp"
+
+namespace mfem {
+
+void ElementMarker::SetLevelSetFunction(Coefficient &ls_fun)
+{
+    FiniteElementCollection* fec=new H1_FECollection(h1_order,smesh->Dimension());
+    FiniteElementSpace* pfes_sltn=new FiniteElementSpace(smesh,fec);
+
+    Vector vals;
+    Array<int> vdofs;
+
+    if(use_cut_marks==false){
+        if(include_cut_elements){
+            elgf=(double)(ElementType::INSIDE);
+            for(int e=0;e<smesh->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]] = ElementType::OUTSIDE;
+                }
+            }
+        }else{//DEFAULT - do not include cuts
+            elgf=(double)(ElementType::OUTSIDE);
+            for(int e=0;e<smesh->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]] = ElementType::INSIDE;
+                }
+            }
+        }
+    }else{// use CUT mark
+        elgf=(double)(ElementType::INSIDE);
+        for(int e=0;e<smesh->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]] = ElementType::OUTSIDE;
+            }else
+                if ((countp>0)&&(countn>0))
+                {
+                    elfes->GetElementVDofs(e,vdofs);
+                    elgf[vdofs[0]] = ElementType::CUT;
+                }
+        }
+    }
+
+    delete pfes_sltn;
+    delete fec;
+}
+
+void ElementMarker::SetLevelSetFunction(const GridFunction &ls_fun)
+{
+    const FiniteElementSpace* pfes_sltn=ls_fun.FESpace();
+    Vector vals;
+    Array<int> vdofs;
+
+    if(use_cut_marks==false){
+        if(include_cut_elements){
+            elgf=(double)(ElementType::INSIDE);
+            for(int e=0;e<smesh->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]] = ElementType::OUTSIDE;
+                }
+            }
+        }else{//DEFAULT - do not include cuts
+            elgf=(double)(ElementType::OUTSIDE);
+            for(int e=0;e<smesh->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]] = ElementType::INSIDE;
+                }
+            }
+        }
+    }else{//use CUT marks
+        elgf=(double)(ElementType::INSIDE);
+        for(int e=0;e<smesh->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]] = ElementType::OUTSIDE;
+            }else
+                if ((countp>0)&&(countn>0))
+                {
+                    elfes->GetElementVDofs(e,vdofs);
+                    elgf[vdofs[0]] = ElementType::CUT;
+                }
+        }
+    }
+}
+
+void ElementMarker::MarkElements(Array<int> &elem_marker)
+{
+    elem_marker.SetSize(smesh->GetNE());
+    for(int e=0;e<smesh->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(smesh->GetNumFaces());
+    face_marker=FaceType::UNDEFINED;
+    IntegrationPoint ip; ip.Init(0);
+
+    for(int f=0;f<smesh->GetNumFaces();f++){
+        auto *ft = smesh->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==ElementType::CUT)&&(attr2!=ElementType::OUTSIDE))
+        {
+            face_marker[f]=FaceType::GHOSTP;
+        }else
+        if((attr1!=ElementType::OUTSIDE)&&(attr2==ElementType::CUT))
+        {
+            face_marker[f]=FaceType::GHOSTP;
+        }
+    }
+}
+
+void ElementMarker::MarkFaces(Array<int> &face_marker)
+{
+    face_marker.SetSize(smesh->GetNumFaces());
+    face_marker=FaceType::UNDEFINED;
+    IntegrationPoint ip; ip.Init(0);
+
+    if(include_cut_elements==true){
+        for(int f=0;f<smesh->GetNumFaces();f++){
+            auto *ft = smesh->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==ElementType::OUTSIDE)||(attr2==ElementType::OUTSIDE)){
+                if(attr1!=attr2){
+                    face_marker[f]=FaceType::SURROGATE;
+                }
+            }
+        }
+    }else{
+        for(int f=0;f<smesh->GetNumFaces();f++){
+            auto *ft = smesh->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==ElementType::INSIDE)||(attr2==ElementType::INSIDE)){
+                if(attr1!=attr2){
+                    face_marker[f]=FaceType::SURROGATE;
+                }
+            }
+        }
+    }
+}
+
+void ElementMarker::ListEssentialTDofs(const Array<int> &elem_marker,
+                                       FiniteElementSpace &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]==ElementType::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]==ElementType::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;
+    if(lfes.GetProlongationMatrix()!=nullptr){
+        vtdof.SetSize(lfes.GetTrueVSize()); vtdof=0.0;
+        lfes.GetProlongationMatrix()->MultTranspose(vvdof,vtdof);
+    }else{
+        vtdof=vvdof;
+    }
+
+    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_MPI
+
+void ParElementMarker::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)(ElementMarker::ElementType::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]] = ElementMarker::ElementType::OUTSIDE;
+                }
+            }
+        }else{//DEFAULT - do not include cuts
+            elgf=(double)(ElementMarker::ElementType::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]] = ElementMarker::ElementType::INSIDE;
+                }
+            }
+        }
+    }else{// use CUT mark
+        elgf=(double)(ElementMarker::ElementType::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]] = ElementMarker::ElementType::OUTSIDE;
+            }else
+                if ((countp>0)&&(countn>0))
+                {
+                    elfes->GetElementVDofs(e,vdofs);
+                    elgf[vdofs[0]] = ElementMarker::ElementType::CUT;
+                }
+        }
+    }
+
+    if(pmesh->GetNRanks()>0){
+        elgf.ExchangeFaceNbrData();
+    }
+
+    delete pfes_sltn;
+    delete fec;
+}
+
+
+void ParElementMarker::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)(ElementMarker::ElementType::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]] = ElementMarker::ElementType::OUTSIDE;
+                }
+            }
+        }else{//DEFAULT - do not include cuts
+            elgf=(double)(ElementMarker::ElementType::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]] = ElementMarker::ElementType::INSIDE;
+                }
+            }
+        }}else{//use CUT marks
+        elgf=(double)(ElementMarker::ElementType::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]] = ElementMarker::ElementType::OUTSIDE;
+            }else
+                if ((countp>0)&&(countn>0))
+                {
+                    elfes->GetElementVDofs(e,vdofs);
+                    elgf[vdofs[0]] = ElementMarker::ElementType::CUT;
+                }
+        }
+    }
+
+    if(pmesh->GetNRanks()>0){
+        elgf.ExchangeFaceNbrData();}
+}
+
+void ParElementMarker::MarkElements(Array<int> &elem_marker)
+{
+    elem_marker.SetSize(pmesh->GetNE());
+
+    IntegrationPoint ip; ip.Init(0);
+
+    for(int e=0;e<pmesh->GetNE();e++)
+    {
+        ElementTransformation* tr=elfes->GetElementTransformation(e);
+        elem_marker[e] = elgf.GetValue(*tr, ip);
+    }
+}
+
+void ParElementMarker::MarkGhostPenaltyFaces(Array<int> &face_marker)
+{
+    face_marker.SetSize(pmesh->GetNumFaces());
+    face_marker=ElementMarker::FaceType::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==ElementMarker::ElementType::CUT)&&(attr2!=ElementMarker::ElementType::OUTSIDE))
+        {
+            face_marker[f]=ElementMarker::FaceType::GHOSTP;
+        }else
+        if((attr1!=ElementMarker::ElementType::OUTSIDE)&&(attr2==ElementMarker::ElementType::CUT))
+        {
+            face_marker[f]=ElementMarker::FaceType::GHOSTP;
+        }
+    }
+
+    if(pmesh->GetNRanks()>0){
+        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==ElementMarker::ElementType::CUT)&&(attr2!=ElementMarker::ElementType::OUTSIDE))
+        {
+            face_marker[faceno]=ElementMarker::FaceType::GHOSTP;
+        }else
+        if((attr1!=ElementMarker::ElementType::OUTSIDE)&&(attr2==ElementMarker::ElementType::CUT))
+        {
+            face_marker[faceno]=ElementMarker::FaceType::GHOSTP;
+        }
+    }
+}
+
+void ParElementMarker::MarkFaces(Array<int> &face_marker)
+{
+    face_marker.SetSize(pmesh->GetNumFaces());
+    face_marker=ElementMarker::FaceType::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==ElementMarker::ElementType::OUTSIDE)||(attr2==ElementMarker::ElementType::OUTSIDE)){
+                if(attr1!=attr2){
+                    face_marker[f]=ElementMarker::FaceType::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==ElementMarker::ElementType::INSIDE)||(attr2==ElementMarker::ElementType::INSIDE)){
+                if(attr1!=attr2){
+                    face_marker[f]=ElementMarker::FaceType::SURROGATE;
+                }
+            }
+        }
+    }
+
+    if(pmesh->GetNRanks()>0){
+        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==ElementMarker::ElementType::OUTSIDE)||(attr2==ElementMarker::ElementType::OUTSIDE)){
+                if(attr1!=attr2){
+                    face_marker[faceno]=ElementMarker::FaceType::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==ElementMarker::ElementType::INSIDE)||(attr2==ElementMarker::ElementType::INSIDE)){
+                if(attr1!=attr2){
+                    face_marker[faceno]=ElementMarker::FaceType::SURROGATE;
+                }
+            }
+        }
+
+    }
+}
+
+void ParElementMarker::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]==ElementMarker::ElementType::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]==ElementMarker::ElementType::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);
+}
+
+
+#endif
+
+}

miniapps/cutfem/cut_marking.hpp

@@ -0,0 +1,185 @@
+// 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.
+
+#ifndef MFEM_MARKING_HPP
+#define MFEM_MARKING_HPP
+
+#include "mfem.hpp"
+
+
+
+namespace mfem{
+
+/// Marking operations for elements, faces, dofs, etc, related to CutFEM.
+class ElementMarker{
+public:
+    enum ElementType {INSIDE = 0, OUTSIDE = 1, CUT = 2};
+
+    enum FaceType {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(Mesh& mesh, bool include_cut=false,
+                  bool use_cut=false, int h1_order_=2)
+    {
+        const int dim=mesh.SpaceDimension();
+        elfec=new L2_FECollection(0,dim);
+        smesh=&mesh;
+
+        elfes=new FiniteElementSpace(smesh,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 GridFunction& 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 ElementType
+    void MarkElements(Array<int> &elem_marker);
+
+    /// Returns the marking of all faces in the
+    /// mesh using  the @a FaceType
+    void MarkFaces(Array<int> &face_marker);
+
+    /// Returns the marking of all faces in the
+    /// mesh using  the @a FaceType.
+    /// 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,
+                            FiniteElementSpace &lfes,
+                            Array<int> &ess_tdof_list) const;
+
+protected:
+    FiniteElementCollection* elfec;
+    Mesh* smesh;
+
+    FiniteElementSpace* elfes;
+    GridFunction 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_MPI
+
+/// Marking operations for elements, faces, dofs, etc,
+/// related to parallel implementations of CutFEM.
+class ParElementMarker
+{
+public:
+    ///Defines parallel element marker class with options to include the cut
+    ///  elements (include_cut=true) or to mark the cut elements as
+    /// ElementType::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.
+    ParElementMarker(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
+    ~ParElementMarker()
+    {
+        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 ElementMarker::ElementType
+    void MarkElements(Array<int> &elem_marker);
+
+    /// Returns the marking of all faces in the
+    /// mesh using  the @a ElementMarker::FaceType
+    void MarkFaces(Array<int> &face_marker);
+
+    /// Returns the marking of all faces in the
+    /// mesh using  the @a ElementMarker::FaceType.
+    /// 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;
+
+    // order of the H1 FE space for level set functions
+    // defined by coefficient
+    int h1_order;
+};
+
+
+
+#endif
+
+
+
+
+}
+
+
+
+
+#endif

miniapps/cutfem/diffusion.cpp

@@ -0,0 +1,137 @@
+//                                MFEM Example 0 with some small changes 
+//
+// 
+//// Description: This example code demonstrates the most basic usage of MFEM to
+//              define a simple finite element discretization of the Laplace
+//              problem -Delta u = 1 with Dirichlet boundary conditions.
+
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "my_integrators.hpp"
+
+
+using namespace std;
+using namespace mfem;
+
+
+real_t f_rhs(const Vector &x);
+real_t u_ex(const Vector &x);
+void u_grad_exact(const Vector &x, Vector &u);
+real_t koeff(const Vector &x);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command line options.
+   string mesh_file = "../../data/star.mesh";
+   int order = 3;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.ParseCheck();
+
+   // 2. Read the mesh from the given mesh file, and refine once uniformly.
+   Mesh mesh(mesh_file);
+   mesh.UniformRefinement();
+	mesh.UniformRefinement();
+
+
+   // 3. Define a finite element space on the mesh. Here we use H1 continuous
+   //    high-order Lagrange finite elements of the given order.
+   H1_FECollection fec(order, mesh.Dimension());
+   FiniteElementSpace fespace(&mesh, &fec);
+   cout << "Number of unknowns: " << fespace.GetTrueVSize() << endl;
+
+   // 4. Extract the list of all the boundary DOFs. These will be marked as
+   //    Dirichlet in order to enforce zero boundary conditions.
+   Array<int> boundary_dofs;
+   fespace.GetBoundaryTrueDofs(boundary_dofs);
+
+   // 5. Define the solution x as a finite element grid function in fespace. Set
+   //    the initial guess to the function defining boundary conditions.
+   GridFunction x(&fespace);
+   FunctionCoefficient bc (u_ex);
+   x.ProjectCoefficient(bc);
+
+   // 6. Set up the linear form b(.) corresponding to the right-hand side.
+   FunctionCoefficient f (f_rhs);
+   LinearForm b(&fespace);
+   b.AddDomainIntegrator(new DomainLFIntegrator(f));
+   b.Assemble();
+
+   // 7. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
+   BilinearForm a(&fespace);
+   FunctionCoefficient coeff(koeff);
+
+   a.AddDomainIntegrator(new MyDiffusionIntegrator(coeff));
+   a.Assemble();
+
+   // 8. Form the linear system A X = B. This includes eliminating boundary
+   //    conditions, applying AMR constraints, and other transformations.
+   SparseMatrix A;
+   Vector B, X;
+   a.FormLinearSystem(boundary_dofs, x, b, A, X, B);
+
+   // 9. Solve the system using PCG with symmetric Gauss-Seidel preconditioner.
+   GSSmoother M(A);
+   PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
+
+   // 10. Recover the solution x as a grid function and save to file
+   a.RecoverFEMSolution(X, b, x);
+
+
+   // compute errors:     
+   VectorFunctionCoefficient u_grad(mesh.Dimension(), u_grad_exact);
+
+   cout << "\n|| u_h - u ||_{L^2} = " << x.ComputeL2Error(bc) << '\n' << endl;
+   cout << "\n|| grad u_h - grad u ||_{L^2} = " << x.ComputeH1Error(&bc,&u_grad) << '\n' << endl;
+
+
+   // save solution with paraview 
+   ParaViewDataCollection paraview_dc("diffusion", &mesh);
+   paraview_dc.SetPrefixPath("ParaView");
+   paraview_dc.SetLevelsOfDetail(order);
+   paraview_dc.SetCycle(0);
+   paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   paraview_dc.SetHighOrderOutput(true);
+   paraview_dc.SetTime(0.0); // set the time
+   paraview_dc.RegisterField("velocity",&x);
+   paraview_dc.Save();
+
+   return 0;
+}
+
+real_t f_rhs(const Vector &x)
+{
+   return sin(x(0)) * sin( x(1)) + (2*x(0)+4)*cos(x(0))*sin(x(1));
+}
+
+real_t u_ex(const Vector &x)
+{
+    return cos(x(0))*sin(x(1));
+}
+
+
+void u_grad_exact(const Vector &x, Vector &u)
+{
+    u(0) =  - sin(x(0)) * sin( x(1));
+    u(1) =  cos(x(0)) * cos( x(1));
+}
+
+real_t koeff(const Vector &x)
+{
+   return x(0) + 2; 
+}
+
+// Q = 1
+// real_t f_rhs(const Vector &x)
+// {
+//    return 2*cos(x(0))*sin(x(1));
+// }
+
+// real_t u_ex(const Vector &x)
+// {
+//     return cos(x(0))*sin(x(1));
+// }

miniapps/cutfem/diffusion_cut.cpp

@@ -0,0 +1,211 @@
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "my_integrators.hpp"
+
+
+using namespace std;
+using namespace mfem;
+
+
+real_t f_rhs(const Vector &x);
+real_t u_ex(const Vector &x);
+real_t g_neumann(const Vector &x);
+void u_grad_exact(const Vector &x, Vector &u);
+real_t circle_func(const Vector &x);
+real_t ellipsoide_func(const Vector &x);
+
+
+
+int main(int argc, char *argv[])
+{
+
+   int n = 50;
+   Mesh mesh = Mesh::MakeCartesian2D( n*2, n , mfem::Element::Type::QUADRILATERAL, true, 1, 0.5);
+
+   int order = 2;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.ParseCheck();
+
+
+   double h_min, h_max, kappa_min, kappa_max;
+   mesh.GetCharacteristics(h_min, h_max, kappa_min, kappa_max);
+
+   H1_FECollection fec(order, mesh.Dimension());
+   FiniteElementSpace fespace(&mesh, &fec);
+   cout << "Number of unknowns: " << fespace.GetTrueVSize() << endl;
+
+   ConstantCoefficient one(1.0);
+
+   GridFunction x(&fespace);
+   FunctionCoefficient bc (u_ex);
+   FunctionCoefficient f (f_rhs);
+   x.ProjectCoefficient(bc);
+   FunctionCoefficient neumann(g_neumann);
+
+
+   // level set function 
+   GridFunction cgf(&fespace);
+   FunctionCoefficient circle(ellipsoide_func);
+   cgf.ProjectCoefficient(circle);
+
+   // mark elements and outside DOFs
+   Array<int> boundary_dofs;
+   fespace.GetBoundaryTrueDofs(boundary_dofs);
+   Array<int> outside_dofs;
+   Array<int> marks;
+   Array<int> face_marks;
+   {
+       Array<int> outside_dofs;
+       ElementMarker* elmark=new ElementMarker(mesh,true,true); 
+       elmark->SetLevelSetFunction(cgf);
+       elmark->MarkElements(marks);
+       elmark->MarkGhostPenaltyFaces(face_marks);
+       elmark->ListEssentialTDofs(marks,fespace,outside_dofs);
+       delete elmark;
+   }
+   outside_dofs.Append(boundary_dofs);
+   outside_dofs.Sort();
+   outside_dofs.Unique();
+
+   int otherorder = 2;
+   int aorder = 2; // Algoim integration points
+   AlgoimIntegrationRules* air=new AlgoimIntegrationRules(aorder,circle,otherorder);
+   real_t gp = 0.1/(h_min*h_min);
+
+   // 6. Set up the linear form b(.) corresponding to the right-hand side.
+
+   LinearForm b(&fespace);
+
+   b.AddDomainIntegrator(new CutDomainLFIntegrator(f,&marks,air));
+   b.AddDomainIntegrator(new CutUnfittedBoundaryLFIntegrator(neumann,&marks,air));
+   b.Assemble();
+
+   // 7. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
+   BilinearForm a(&fespace);
+   a.AddDomainIntegrator(new CutDiffusionIntegrator(one,&marks,air));
+   a.AddInteriorFaceIntegrator(new CutGhostPenaltyIntegrator(gp,&face_marks));
+   a.Assemble();
+
+
+   SparseMatrix A;
+   Vector B, X;
+   a.FormLinearSystem(boundary_dofs, x, b, A, X, B);
+
+   // 9. Solve the system using PCG with symmetric Gauss-Seidel preconditioner.
+   GSSmoother M(A);
+   PCG(A,M, B, X, 2, 2000, 1e-16, 0.0);
+
+   // 10. Recover the solution x as a grid function and save to file
+   a.RecoverFEMSolution(X, b, x);
+
+
+   //compute the error
+   {
+       NonlinearForm* nf=new NonlinearForm(&fespace);
+       nf->AddDomainIntegrator(new CutScalarErrorIntegrator(bc,&marks,air));
+
+       cout << "\n|| u_h - u ||_{L^2} = " << nf->GetEnergy(x.GetTrueVector())<< std::endl;
+
+       delete nf;
+   }
+
+   // to visualize level set and markings
+   L2_FECollection* l2fec= new L2_FECollection(0,mesh.Dimension());
+   FiniteElementSpace* l2fes= new FiniteElementSpace(&mesh,l2fec,1);
+   GridFunction mgf(l2fes);
+   for(int i=0;i<marks.Size();i++){
+      mgf[i]=marks[i];
+   }
+   GridFunction exact_sol(&fespace);
+   exact_sol.ProjectCoefficient(bc);
+
+
+   GridFunction error(&fespace);
+   error = x;
+   error -= exact_sol;
+
+    // GLVis 
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock(vishost, visport);
+   sol_sock.precision(8);
+   sol_sock << "solution\n" << mesh << x << flush;
+   //    // save solution with paraview 
+   ParaViewDataCollection paraview_dc("diffusion_cut", &mesh);
+   paraview_dc.SetPrefixPath("ParaView");
+   paraview_dc.SetLevelsOfDetail(order);
+   paraview_dc.SetCycle(0);
+   paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   paraview_dc.SetHighOrderOutput(true);
+   paraview_dc.SetTime(0.0); // set the time
+   paraview_dc.RegisterField("solution",&x);
+   paraview_dc.RegisterField("marks", &mgf);
+   paraview_dc.RegisterField("level_set",&cgf);
+   paraview_dc.RegisterField("exact_sol",&exact_sol);
+   paraview_dc.RegisterField("error",&error);
+   paraview_dc.Save();
+
+
+   delete l2fec;
+   delete l2fes;
+   delete air;
+
+   cout << "h:" << h_min<<endl;
+   return 0;
+}
+
+
+real_t f_rhs(const Vector &x)
+{
+
+    return  2*M_PI*M_PI*(sin(M_PI*x(0))*cos((M_PI*x(1))));
+}
+
+
+real_t u_ex(const Vector &x)
+{
+       real_t x0 = 0.5;
+    real_t y0 = 0.5;
+   return sin(M_PI*x(0))*cos((M_PI*x(1)));
+}
+
+void u_grad_exact(const Vector &x, Vector &u)
+{
+    u(0) =M_PI*cos(M_PI*x(0)) *cos(M_PI* x(1));
+    u(1) =  -M_PI*sin(M_PI*x(0)) * sin(M_PI* x(1));
+}
+
+
+real_t g_neumann(const Vector &x)
+{
+    real_t a = 1.5;
+    real_t b = 0.5;
+    real_t x0 = 0.5;
+    real_t y0 = 0.25;
+    real_t xx = x(0)-x0;
+    real_t y = x(1)-y0;
+    real_t normalize = sqrt((xx*xx)/(a*a*a*a) + y*y/(b*b*b*b));
+   return M_PI*cos(M_PI*x(0)) *cos(M_PI* x(1))*xx/(a*a*normalize) -M_PI*sin(M_PI*x(0)) * sin(M_PI* x(1))*y/(b*b*normalize);
+}
+
+real_t circle_func(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0.5;
+    real_t r = 0.4;
+    return -(x(0)-x0)*(x(0)-x0) - (x(1)-y0)*(x(1)-y0) + r*r;  
+}
+
+real_t ellipsoide_func(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0.25;
+    real_t r = 0.35;
+    real_t xx = x(0)-x0;
+    real_t y = x(1)-y0;
+    return -(xx)*(xx)/(1.5*1.5) - (y)*(y)/(0.5*0.5)+ r*r; // + 0.25*cos(atan2(x(1)-y0,x(0)-x0))*cos(atan2(x(1)-y0,x(0)-x0));  
+}

miniapps/cutfem/diffusion_nitsche.cpp

@@ -0,0 +1,118 @@
+//       Solving the  Laplace   problem -Delta u = 1 with 
+//       Dirichlet boundary conditions weakly inforced. 
+//
+
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "my_integrators.hpp"
+
+
+using namespace std;
+using namespace mfem;
+
+
+real_t f_rhs(const Vector &x);
+real_t u_ex(const Vector &x);
+void u_grad_exact(const Vector &x, Vector &u);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command line options.
+   string mesh_file = "../../data/star.mesh";
+   int order = 2;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.ParseCheck();
+
+   // 2. Read the mesh from the given mesh file, and refine once uniformly.
+   Mesh mesh(mesh_file);
+   mesh.UniformRefinement();
+
+   // 3. Define a finite element space on the mesh. 
+   H1_FECollection fec(order, mesh.Dimension());
+   FiniteElementSpace fespace(&mesh, &fec);
+   cout << "Number of unknowns: " << fespace.GetTrueVSize() << endl;
+
+   // 4. Array of boundary dofs, empty but needed in FormLinearSystem
+   Array<int> boundary_dofs;
+
+   // 5. Define the solution x as a finite element grid function in fespace. 
+   GridFunction x(&fespace);
+
+   // 6. Set up the linear form b(.) corresponding to the right-hand side.
+   LinearForm b(&fespace);
+
+   FunctionCoefficient f (f_rhs);
+   ConstantCoefficient one(1.0);
+   FunctionCoefficient bc (u_ex);
+   real_t sigma = -1.0; // IP 
+   real_t lambda = 10.0; // Nitsche penalty parameter
+
+   b.AddDomainIntegrator(new DomainLFIntegrator(f));
+   // rhs Nitsche terms 
+   b.AddBdrFaceIntegrator(new DGDirichletLFIntegrator(bc, one, sigma, lambda));
+   b.Assemble();
+
+   // 7. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
+   BilinearForm a(&fespace);
+
+   a.AddDomainIntegrator(new MyDiffusionIntegrator);
+   // Nitsche terms 
+    a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, lambda)); 
+    // a.AddBdrFaceIntegrator(new MyNitscheBilinIntegrator(one, lambda)); //not working
+   a.Assemble();
+
+   // 8. Form the linear system A X = B. 
+   SparseMatrix A;
+   Vector B, X;
+   a.FormLinearSystem(boundary_dofs, x, b, A, X, B);
+
+   // 9. Solve the system using PCG with symmetric Gauss-Seidel preconditioner.
+   GSSmoother M(A);
+   PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
+
+   // 10. Recover the solution x as a grid function 
+   a.RecoverFEMSolution(X, b, x);
+
+
+   VectorFunctionCoefficient u_grad(mesh.Dimension(), u_grad_exact);
+
+   cout << "\n|| u_h - u ||_{L^2} = " << x.ComputeL2Error(bc) << '\n' << endl;
+   cout << "\n|| grad u_h - grad u ||_{L^2} = " << x.ComputeH1Error(&bc,&u_grad) << '\n' << endl;
+
+
+   ParaViewDataCollection paraview_dc("diffusionnit", &mesh);
+   paraview_dc.SetPrefixPath("ParaView");
+   paraview_dc.SetLevelsOfDetail(order);
+   paraview_dc.SetCycle(0);
+   paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   paraview_dc.SetHighOrderOutput(true);
+   paraview_dc.SetTime(0.0); // set the time
+   paraview_dc.RegisterField("velocity",&x);
+   paraview_dc.Save();
+
+   return 0;
+}
+
+real_t f_rhs(const Vector &x)
+{
+   return 2*cos(x(0))*sin(x(1));
+}
+
+real_t u_ex(const Vector &x)
+{
+    return cos(x(0))*sin(x(1));
+}
+
+
+void u_grad_exact(const Vector &x, Vector &u)
+{
+    u(0) =  - sin(x(0)) * sin( x(1));
+    u(1) =  cos(x(0)) * cos( x(1));
+}
+
+

miniapps/cutfem/makefile

@@ -0,0 +1,79 @@
+# 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.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+MFEM_INSTALL_DIR ?= ../../mfem
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/cutfem/,)
+CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
+   $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
+
+# Include defaults.mk to get XLINKER
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+include $(DEFAULTS_MK)
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+cutfem_COMMON_SRC = my_integrators.cpp 
+
+cutfem_COMMON_OBJ = $(cutfem_COMMON_SRC:.cpp=.o)
+
+SEQ_MINIAPPS = diffusion
+PAR_MINIAPPS = parheat
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS = $(SEQ_MINIAPPS)
+else
+   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+%: %.o $(cutfem_COMMON_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $^ -o $@ $(MFEM_LIBS)
+
+%.o: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
+
+all: $(MINIAPPS)
+
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+# Testing: Parallel vs. serial runs
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+TEST_NAME := cutfem miniapp
+%-test-par: %
+	@$(call mfem-test,$<, $(RUN_MPI), $(TEST_NAME))
+%-test-seq: %
+	@$(call mfem-test,$<,, $(TEST_NAME))
+
+# Testing: "test" target and mfem-test* variables are defined in config/test.mk
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -rf SeqHeat* ParHeat*

miniapps/cutfem/par_diffusion_cut.cpp

@@ -0,0 +1,398 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+#include "my_integrators.hpp"
+
+using namespace std;
+using namespace mfem;
+
+real_t f_rhs(const Vector &x);
+real_t u_ex(const Vector &x);
+real_t bcf(const Vector &x);
+real_t g_neumann(const Vector &x);
+real_t circle_func(const Vector &x);
+real_t ellipsoide_func(const Vector &x);
+real_t sphere_func(const Vector &x);
+real_t new_func(const Vector &xx);
+real_t new_func3d(const Vector &xx);
+real_t new_func2d(const Vector &xx);
+real_t koeff(const Vector &x);
+
+
+// solves the diffusion problem Delta u = f, with either Dirichlet conditions weakly imposed, 
+// or neumann + Dirichlet conditions 
+int main(int argc, char *argv[])
+{
+    // 1. Initialize MPI and HYPRE.
+    Mpi::Init();
+    int num_procs = Mpi::WorldSize();
+    int myid = Mpi::WorldRank();
+    Hypre::Init();
+
+    // 2. Parse command-line options.
+    const char *mesh_file = "../data/star.mesh";
+    int order = 1;
+    bool static_cond = false;
+    bool pa = false;
+    bool fa = false;
+    const char *device_config = "cpu";
+    bool visualization = false;
+    bool visualization_paraview = true;
+    bool algebraic_ceed = false;
+    int ser_ref_levels = 1;
+    int aorder = 8; // Algoim integration points
+    real_t g = 1;
+
+    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(&static_cond, "-sc", "--static-condensation", "-no-sc", 
+                   "--no-static-condensation", "Enable static condensation.");  //these three not valid options now
+    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                   "--no-partial-assembly", "Enable Partial Assembly.");
+    args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
+                   "--no-full-assembly", "Enable Full Assembly.");
+    args.AddOption(&device_config, "-d", "--device",
+                   "Device configuration string, see Device::Configure().");
+    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                   "--no-visualization",
+                   "Enable or disable GLVis visualization.");
+    args.AddOption(&visualization_paraview, "-vispv", "--visualizationpv", "-no-vispv",
+                   "--no-visualizationpv",
+                   "Enable or disable ParaView visualization.");
+    args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+    args.AddOption(&aorder, "-ao", "--aorder",
+                  "Set order for alogim integration");
+    args.AddOption(&g, "-g", "--ghost penalty constant",
+                  "Ghost penalty constant");
+    args.Parse();
+    if (!args.Good())
+    {
+        if (myid == 0)
+        {
+            args.PrintUsage(cout);
+        }
+        return 1;
+    }
+    if (myid == 0)
+    {
+        args.PrintOptions(cout);
+    }
+    //Mesh mesh(mesh_file, 1, 1);
+
+    Mesh mesh = Mesh::MakeCartesian2D(2, 2, mfem::Element::Type::QUADRILATERAL, true, 1, 1);
+    // Mesh mesh = Mesh::MakeCartesian3D( 4, 4 ,4, mfem::Element::Type::HEXAHEDRON, 2.2, 2.2,2.2 );
+
+
+    int dim = mesh.Dimension();
+    {
+        for (int l = 0; l < ser_ref_levels; l++)
+        {
+            mesh.UniformRefinement();
+        }
+    }
+
+    ParMesh pmesh(MPI_COMM_WORLD, mesh);
+    mesh.Clear();
+    {
+        int par_ref_levels = 1;
+        for (int l = 0; l < par_ref_levels; l++)
+        {
+            pmesh.UniformRefinement();
+        }
+    }
+
+    std::cout<<"id="<<myid<<" "<<pmesh.GetNE()<<" "<<pmesh.GetNumFaces()
+            <<" "<<pmesh.GetNSharedFaces()<<std::endl; std::cout.flush();
+
+    double h_min, h_max, kappa_min, kappa_max;
+    pmesh.GetCharacteristics(h_min, h_max, kappa_min, kappa_max);
+
+    FiniteElementCollection *fec;
+    fec = new H1_FECollection(order, dim);
+    ParFiniteElementSpace fespace(&pmesh, fec);
+    HYPRE_BigInt size = fespace.GlobalTrueVSize();
+    if (myid == 0)
+    {
+        cout << "Number of finite element unknowns: " << size << endl;
+    }
+
+    ConstantCoefficient one(1.0);
+    // FunctionCoefficient coeff(koeff);
+    ParGridFunction x(&fespace);
+    FunctionCoefficient bc (u_ex);
+    FunctionCoefficient f (f_rhs);
+    x.ProjectCoefficient(bc);
+    FunctionCoefficient neumann(g_neumann);
+
+    // level set function
+    ParGridFunction lsgf(&fespace);
+    FunctionCoefficient level_set(circle_func);
+    lsgf.ProjectCoefficient(level_set);
+
+    // mark elements and outside DOFs
+    Array<int> boundary_dofs;
+    // fespace.GetBoundaryTrueDofs(boundary_dofs);
+    Array<int> outside_dofs;
+    Array<int> marks;
+    Array<int> face_marks;
+    {
+        ParElementMarker* elmark=new ParElementMarker(pmesh,true,true);
+        elmark->SetLevelSetFunction(lsgf);
+        elmark->MarkElements(marks);
+        elmark->MarkGhostPenaltyFaces(face_marks);
+        elmark->ListEssentialTDofs(marks,fespace,outside_dofs);
+        delete elmark;
+    }
+    outside_dofs.Append(boundary_dofs);
+    outside_dofs.Sort();
+    outside_dofs.Unique();
+
+
+    std::cout<<"myid="<<myid<<" marks_size="<<marks.Size()<<std::endl;
+    std::cout.flush();
+
+    int otherorder = 2;
+    AlgoimIntegrationRules* air=new AlgoimIntegrationRules(aorder,level_set,otherorder);
+    real_t gp = g/(h_min*h_min);
+    real_t lambda = 10/h_min;
+
+    ParLinearForm b(&fespace);
+
+    b.AddDomainIntegrator(new CutDomainLFIntegrator(f,&marks,air));
+    // b.AddDomainIntegrator(new CutUnfittedBoundaryLFIntegrator(neumann,&marks,air)); //when neumann condition 
+    b.AddDomainIntegrator(new  CutUnfittedNitscheLFIntegrator(bc,one,lambda,&marks,air));
+
+    b.Assemble();
+
+    ParBilinearForm a(&fespace);
+    a.AddDomainIntegrator(new CutDiffusionIntegrator(one,&marks,air,false));
+    a.AddInteriorFaceIntegrator(new CutGhostPenaltyIntegrator(gp,&face_marks));
+    a.AddDomainIntegrator(new CutNitscheIntegrator(one,lambda,&marks,air));
+
+    a.Assemble();
+
+    OperatorPtr A;
+    Vector B, X;
+    a.FormLinearSystem(outside_dofs, x, b, A, X, B);
+
+    Solver *prec = nullptr;
+    prec = new HypreBoomerAMG;
+    CGSolver cg(MPI_COMM_WORLD);
+    cg.SetRelTol(1e-20);
+    cg.SetMaxIter(1500);
+    cg.SetPrintLevel(1);
+    if (prec) { cg.SetPreconditioner(*prec); }
+    cg.SetOperator(*A);
+    cg.Mult(B, X);
+    delete prec;
+
+    a.RecoverFEMSolution(X, b, x);
+
+    FunctionCoefficient uex (u_ex);
+   ParGridFunction exact_sol(&fespace);
+   exact_sol.ProjectCoefficient(uex);
+
+    // to visualize level set and markings
+    L2_FECollection* l2fec= new L2_FECollection(0,pmesh.Dimension());
+    ParFiniteElementSpace* l2fes= new ParFiniteElementSpace(&pmesh,l2fec,1);
+    ParGridFunction mgf(l2fes); mgf=0.0;
+    ParGridFunction par(l2fes); par=0.0;
+    for(int i=0;i<marks.Size();i++){
+        mgf[i]=marks[i];
+        par[i]=myid;
+    }
+
+
+   {
+       ParNonlinearForm* nf=new ParNonlinearForm(&fespace);
+       nf->AddDomainIntegrator(new CutScalarErrorIntegrator(uex,&marks,air));
+        real_t error_squared = nf->GetEnergy(x.GetTrueVector());
+        if (myid==0)
+        {
+       cout << "\n|| u_h - u ||_{L^2} = " << sqrt(error_squared)<< std::endl;
+        cout << "h: " << h_min<< std::endl;
+        }
+       delete nf;
+   }
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << pmesh << x << flush;
+   }
+
+
+
+    if (visualization_paraview)
+    {
+
+        ParGridFunction error(&fespace);
+        error = x;
+        error -= exact_sol;
+
+        ParaViewDataCollection paraview_dc("diffusion_cut", &pmesh);
+        paraview_dc.SetPrefixPath("ParaView");
+        paraview_dc.SetLevelsOfDetail(order);
+        paraview_dc.SetCycle(0);
+        paraview_dc.SetDataFormat(VTKFormat::BINARY);
+        paraview_dc.SetHighOrderOutput(true);
+        paraview_dc.SetTime(0.0); // set the time
+        paraview_dc.RegisterField("u",&x);
+        paraview_dc.RegisterField("marks", &mgf);
+        paraview_dc.RegisterField("parts", &par);
+        paraview_dc.RegisterField("level_set",&lsgf);
+        paraview_dc.RegisterField("error",&error);
+        paraview_dc.RegisterField("u_ex",&exact_sol);
+        paraview_dc.Save();
+    }
+
+    delete l2fes;
+    delete l2fec;
+    delete air;
+    delete fec;
+
+    return 0;
+
+}
+
+real_t f_rhs(const Vector &x)
+{
+
+    return  2*M_PI*M_PI*(sin(M_PI*x(0))*cos((M_PI*x(1))));
+}
+
+
+
+real_t u_ex(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0.5;
+    return sin(M_PI*x(0))*cos((M_PI*x(1)));
+}
+
+real_t bcf(const Vector &x)
+{
+    return cos(M_PI*x(0))*cos((M_PI*x(1)))+ sin(M_PI*x(0))*sin((M_PI*x(1)));;//sin(M_PI*x(0))*cos((M_PI*x(1)));
+}
+
+
+real_t g_neumann(const Vector &xx)
+{    
+    real_t x = xx(0);
+    real_t y = xx(1);
+
+    real_t dldx =- M_PI*(4*sin(x*2* M_PI)*cos(x*2* M_PI) +cos(x*M_PI/2)/2 );
+    real_t dldy = 1 ;
+
+    real_t dudx = M_PI*(-sin(M_PI*xx(0)) * cos(M_PI* xx(1)) + cos(M_PI*xx(0)) * sin(M_PI* xx(1)));
+    real_t dudy = M_PI*( - cos(M_PI*xx(0)) * sin(M_PI* xx(1)) + sin(M_PI*xx(0)) * cos(M_PI* xx(1)));
+
+    real_t normalize = sqrt(dldx*dldx + dldy*dldy);
+   return dudx *dldx/normalize + dudy *dldy/normalize;
+}
+
+
+
+real_t g_neumann3d(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0;
+    real_t z0 = 0.5;
+    real_t xx = x(0)-x0;
+    real_t y = x(1)-y0;
+    real_t z = x(2)-z0;
+    real_t normalize = sqrt(xx*xx + y*y + z*z);
+    return M_PI*cos(M_PI*x(0)) *cos(M_PI* x(1))*cos((x(2)))*xx/(normalize) -M_PI*sin(M_PI*x(0)) * sin(M_PI* x(1))*cos((x(2)))*y/(normalize)- sin(M_PI*x(0)) * cos(M_PI* x(1))*sin((x(2)))*z/(normalize);
+}
+
+
+real_t circle_func(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0.5;
+    real_t r = 0.4;
+    return -(x(0)-x0)*(x(0)-x0) - (x(1)-y0)*(x(1)-y0) + r*r;
+}
+
+real_t ellipsoide_func(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0.25;
+    real_t r = 0.351;
+    real_t xx = x(0)-x0;
+    real_t y = x(1)-y0;
+    return -(xx)*(xx)/(1.5*1.5) - (y)*(y)/(0.5*0.5)+ r*r; // + 0.25*cos(atan2(x(1)-y0,x(0)-x0))*cos(atan2(x(1)-y0,x(0)-x0));
+}
+
+
+real_t sphere_func(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0.5;
+    real_t z0 = 0.5;
+    real_t r = 0.4;
+    return -(x(0)-x0)*(x(0)-x0) - (x(1)-y0)*(x(1)-y0)  - (x(2)-z0)*(x(2)-z0) + r*r;
+}
+
+
+real_t new_func(const Vector &xx)
+{
+
+    real_t x = xx(0);
+    real_t y = xx(1);
+
+    return sin(x*2* M_PI)*sin(x*2* M_PI)+sin(x*M_PI/2)-y;
+
+}
+
+real_t new_func3d(const Vector &xx)
+{
+
+    real_t x = xx(0)-1;
+    real_t y = xx(1)-1;
+    real_t z = xx(2)-1;
+    real_t r = 0.5;
+    real_t r0 = 3.5;
+
+    return -(sqrt(x*x + y*y + z*z) - r + r/r0*cos(5*atan2(y,x))*cos(M_PI*z));
+
+}
+
+
+real_t new_func2d(const Vector &xx)
+{
+
+    real_t x = xx(0)-1.1;
+    real_t y = xx(1)-1.1;
+    real_t r = 0.6;
+    real_t r0 = 0.2;
+// - r0*cos(atan2(y,x)))
+    return -(sqrt(x*x + y*y) - r- r0*cos(5*atan2(y,x)));
+
+}
+
+// real_t f_rhs(const Vector &x) //koeff
+// {
+//    return sin(x(0)) * sin( x(1)) + (2*x(0)+4)*cos(x(0))*sin(x(1));
+// }
+
+
+real_t koeff(const Vector &x)
+{
+   return x(0) + 2; 
+}
+
+// real_t u_ex(const Vector &x)
+// {
+//     return cos(x(0))*sin(x(1));
+// }

miniapps/cutfem/par_diffusion_cut_3D.cpp

@@ -0,0 +1,399 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+#include "my_integrators.hpp"
+
+using namespace std;
+using namespace mfem;
+
+real_t f_rhs(const Vector &x);
+real_t u_ex(const Vector &x);
+real_t bcf(const Vector &x);
+real_t g_neumann(const Vector &x);
+real_t g_neumann3d(const Vector &x);
+real_t circle_func(const Vector &x);
+real_t ellipsoide_func(const Vector &x);
+real_t sphere_func(const Vector &x);
+real_t new_func(const Vector &xx);
+real_t new_func3d(const Vector &xx);
+real_t new_func2d(const Vector &xx);
+real_t koeff(const Vector &x);
+
+
+// solves the diffusion problem Delta u = f, with either Dirichlet conditions weakly imposed, 
+// or neumann + Dirichlet conditions 
+int main(int argc, char *argv[])
+{
+    // 1. Initialize MPI and HYPRE.
+    Mpi::Init();
+    int num_procs = Mpi::WorldSize();
+    int myid = Mpi::WorldRank();
+    Hypre::Init();
+
+    // 2. Parse command-line options.
+    const char *mesh_file = "../data/star.mesh";
+    int order = 1;
+    bool static_cond = false;
+    bool pa = false;
+    bool fa = false;
+    const char *device_config = "cpu";
+    bool visualization = false;
+    bool visualization_paraview = true;
+    bool algebraic_ceed = false;
+    int ser_ref_levels = 1;
+    int aorder = 2; // Algoim integration points
+    real_t g = 1;
+
+    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(&static_cond, "-sc", "--static-condensation", "-no-sc", 
+                   "--no-static-condensation", "Enable static condensation.");  //these three not valid options now
+    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                   "--no-partial-assembly", "Enable Partial Assembly.");
+    args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
+                   "--no-full-assembly", "Enable Full Assembly.");
+    args.AddOption(&device_config, "-d", "--device",
+                   "Device configuration string, see Device::Configure().");
+    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                   "--no-visualization",
+                   "Enable or disable GLVis visualization.");
+    args.AddOption(&visualization_paraview, "-vispv", "--visualizationpv", "-no-vispv",
+                   "--no-visualizationpv",
+                   "Enable or disable ParaView visualization.");
+    args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+    args.AddOption(&aorder, "-ao", "--aorder",
+                  "Set order for alogim integration");
+    args.AddOption(&g, "-g", "--ghost penalty constant",
+                  "Ghost penalty constant");
+    args.Parse();
+    if (!args.Good())
+    {
+        if (myid == 0)
+        {
+            args.PrintUsage(cout);
+        }
+        return 1;
+    }
+    if (myid == 0)
+    {
+        args.PrintOptions(cout);
+    }
+    //Mesh mesh(mesh_file, 1, 1);
+
+    //Mesh mesh = Mesh::MakeCartesian2D(2, 2, mfem::Element::Type::QUADRILATERAL, true, 1, 1);
+    Mesh mesh = Mesh::MakeCartesian3D( 4, 4 ,4, mfem::Element::Type::HEXAHEDRON, 2.2, 2.2,2.2 );
+
+
+    int dim = mesh.Dimension();
+    {
+        for (int l = 0; l < ser_ref_levels; l++)
+        {
+            mesh.UniformRefinement();
+        }
+    }
+
+    ParMesh pmesh(MPI_COMM_WORLD, mesh);
+    mesh.Clear();
+    {
+        int par_ref_levels = 1;
+        for (int l = 0; l < par_ref_levels; l++)
+        {
+            pmesh.UniformRefinement();
+        }
+    }
+
+    std::cout<<"id="<<myid<<" "<<pmesh.GetNE()<<" "<<pmesh.GetNumFaces()
+            <<" "<<pmesh.GetNSharedFaces()<<std::endl; std::cout.flush();
+
+    double h_min, h_max, kappa_min, kappa_max;
+    pmesh.GetCharacteristics(h_min, h_max, kappa_min, kappa_max);
+
+    FiniteElementCollection *fec;
+    fec = new H1_FECollection(order, dim);
+    ParFiniteElementSpace fespace(&pmesh, fec);
+    HYPRE_BigInt size = fespace.GlobalTrueVSize();
+    if (myid == 0)
+    {
+        cout << "Number of finite element unknowns: " << size << endl;
+    }
+
+    ConstantCoefficient one(1.0);
+    // FunctionCoefficient coeff(koeff);
+    ParGridFunction x(&fespace);
+    FunctionCoefficient bc (u_ex);
+    FunctionCoefficient f (f_rhs);
+    x.ProjectCoefficient(bc);
+    FunctionCoefficient neumann(g_neumann3d);
+
+    // level set function
+    ParGridFunction lsgf(&fespace);
+    FunctionCoefficient level_set(sphere_func);
+    lsgf.ProjectCoefficient(level_set);
+
+    // mark elements and outside DOFs
+    Array<int> boundary_dofs;
+    // fespace.GetBoundaryTrueDofs(boundary_dofs);
+    Array<int> outside_dofs;
+    Array<int> marks;
+    Array<int> face_marks;
+    {
+        ParElementMarker* elmark=new ParElementMarker(pmesh,true,true);
+        elmark->SetLevelSetFunction(lsgf);
+        elmark->MarkElements(marks);
+        elmark->MarkGhostPenaltyFaces(face_marks);
+        elmark->ListEssentialTDofs(marks,fespace,outside_dofs);
+        delete elmark;
+    }
+    outside_dofs.Append(boundary_dofs);
+    outside_dofs.Sort();
+    outside_dofs.Unique();
+
+
+    std::cout<<"myid="<<myid<<" marks_size="<<marks.Size()<<std::endl;
+    std::cout.flush();
+
+    int otherorder = 2;
+    AlgoimIntegrationRules* air=new AlgoimIntegrationRules(aorder,level_set,otherorder);
+    real_t gp = g/(h_min*h_min);
+    real_t lambda = 10/h_min;
+
+    ParLinearForm b(&fespace);
+
+    b.AddDomainIntegrator(new CutDomainLFIntegrator(f,&marks,air));
+    // b.AddDomainIntegrator(new CutUnfittedBoundaryLFIntegrator(neumann,&marks,air)); //when neumann condition 
+    b.AddDomainIntegrator(new  CutUnfittedNitscheLFIntegrator(bc,one,lambda,&marks,air));
+
+    b.Assemble();
+
+    ParBilinearForm a(&fespace);
+    a.AddDomainIntegrator(new CutDiffusionIntegrator(one,&marks,air,false));
+    a.AddInteriorFaceIntegrator(new CutGhostPenaltyIntegrator(gp,&face_marks));
+    a.AddDomainIntegrator(new CutNitscheIntegrator(one,lambda,&marks,air));
+
+    a.Assemble();
+
+    OperatorPtr A;
+    Vector B, X;
+    a.FormLinearSystem(outside_dofs, x, b, A, X, B);
+
+    Solver *prec = nullptr;
+    prec = new HypreBoomerAMG;
+    CGSolver cg(MPI_COMM_WORLD);
+    cg.SetRelTol(1e-20);
+    cg.SetMaxIter(1500);
+    cg.SetPrintLevel(1);
+    if (prec) { cg.SetPreconditioner(*prec); }
+    cg.SetOperator(*A);
+    cg.Mult(B, X);
+    delete prec;
+
+    a.RecoverFEMSolution(X, b, x);
+
+    FunctionCoefficient uex (u_ex);
+   ParGridFunction exact_sol(&fespace);
+   exact_sol.ProjectCoefficient(uex);
+
+    // to visualize level set and markings
+    L2_FECollection* l2fec= new L2_FECollection(0,pmesh.Dimension());
+    ParFiniteElementSpace* l2fes= new ParFiniteElementSpace(&pmesh,l2fec,1);
+    ParGridFunction mgf(l2fes); mgf=0.0;
+    ParGridFunction par(l2fes); par=0.0;
+    for(int i=0;i<marks.Size();i++){
+        mgf[i]=marks[i];
+        par[i]=myid;
+    }
+
+
+   {
+       ParNonlinearForm* nf=new ParNonlinearForm(&fespace);
+       nf->AddDomainIntegrator(new CutScalarErrorIntegrator(uex,&marks,air));
+        real_t error_squared = nf->GetEnergy(x.GetTrueVector());
+        if (myid==0)
+        {
+       cout << "\n|| u_h - u ||_{L^2} = " << sqrt(error_squared)<< std::endl;
+        cout << "h: " << h_min<< std::endl;
+        }
+       delete nf;
+   }
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << pmesh << x << flush;
+   }
+
+
+
+    if (visualization_paraview)
+    {
+
+        ParGridFunction error(&fespace);
+        error = x;
+        error -= exact_sol;
+
+        ParaViewDataCollection paraview_dc("diffusion_cut", &pmesh);
+        paraview_dc.SetPrefixPath("ParaView");
+        paraview_dc.SetLevelsOfDetail(order);
+        paraview_dc.SetCycle(0);
+        paraview_dc.SetDataFormat(VTKFormat::BINARY);
+        paraview_dc.SetHighOrderOutput(true);
+        paraview_dc.SetTime(0.0); // set the time
+        paraview_dc.RegisterField("u",&x);
+        paraview_dc.RegisterField("marks", &mgf);
+        paraview_dc.RegisterField("parts", &par);
+        paraview_dc.RegisterField("level_set",&lsgf);
+        paraview_dc.RegisterField("error",&error);
+        paraview_dc.RegisterField("u_ex",&exact_sol);
+        paraview_dc.Save();
+    }
+
+    delete l2fes;
+    delete l2fec;
+    delete air;
+    delete fec;
+
+    return 0;
+
+}
+
+real_t f_rhs(const Vector &x)
+{
+
+    return  2*M_PI*M_PI*(sin(M_PI*x(0))*cos((M_PI*x(1))));
+}
+
+
+
+real_t u_ex(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0.5;
+    return sin(M_PI*x(0))*cos((M_PI*x(1)));
+}
+
+real_t bcf(const Vector &x)
+{
+    return cos(M_PI*x(0))*cos((M_PI*x(1)))+ sin(M_PI*x(0))*sin((M_PI*x(1)));;//sin(M_PI*x(0))*cos((M_PI*x(1)));
+}
+
+
+real_t g_neumann(const Vector &xx)
+{    
+    real_t x = xx(0);
+    real_t y = xx(1);
+
+    real_t dldx =- M_PI*(4*sin(x*2* M_PI)*cos(x*2* M_PI) +cos(x*M_PI/2)/2 );
+    real_t dldy = 1 ;
+
+    real_t dudx = M_PI*(-sin(M_PI*xx(0)) * cos(M_PI* xx(1)) + cos(M_PI*xx(0)) * sin(M_PI* xx(1)));
+    real_t dudy = M_PI*( - cos(M_PI*xx(0)) * sin(M_PI* xx(1)) + sin(M_PI*xx(0)) * cos(M_PI* xx(1)));
+
+    real_t normalize = sqrt(dldx*dldx + dldy*dldy);
+   return dudx *dldx/normalize + dudy *dldy/normalize;
+}
+
+
+
+real_t g_neumann3d(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0;
+    real_t z0 = 0.5;
+    real_t xx = x(0)-x0;
+    real_t y = x(1)-y0;
+    real_t z = x(2)-z0;
+    real_t normalize = sqrt(xx*xx + y*y + z*z);
+    return M_PI*cos(M_PI*x(0)) *cos(M_PI* x(1))*cos((x(2)))*xx/(normalize) -M_PI*sin(M_PI*x(0)) * sin(M_PI* x(1))*cos((x(2)))*y/(normalize)- sin(M_PI*x(0)) * cos(M_PI* x(1))*sin((x(2)))*z/(normalize);
+}
+
+
+real_t circle_func(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0.5;
+    real_t r = 0.4;
+    return -(x(0)-x0)*(x(0)-x0) - (x(1)-y0)*(x(1)-y0) + r*r;
+}
+
+real_t ellipsoide_func(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0.25;
+    real_t r = 0.351;
+    real_t xx = x(0)-x0;
+    real_t y = x(1)-y0;
+    return -(xx)*(xx)/(1.5*1.5) - (y)*(y)/(0.5*0.5)+ r*r; // + 0.25*cos(atan2(x(1)-y0,x(0)-x0))*cos(atan2(x(1)-y0,x(0)-x0));
+}
+
+
+real_t sphere_func(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0.5;
+    real_t z0 = 0.5;
+    real_t r = 0.4;
+    return -(x(0)-x0)*(x(0)-x0) - (x(1)-y0)*(x(1)-y0)  - (x(2)-z0)*(x(2)-z0) + r*r;
+}
+
+
+real_t new_func(const Vector &xx)
+{
+
+    real_t x = xx(0);
+    real_t y = xx(1);
+
+    return sin(x*2* M_PI)*sin(x*2* M_PI)+sin(x*M_PI/2)-y;
+
+}
+
+real_t new_func3d(const Vector &xx)
+{
+
+    real_t x = xx(0)-1;
+    real_t y = xx(1)-1;
+    real_t z = xx(2)-1;
+    real_t r = 0.5;
+    real_t r0 = 3.5;
+
+    return -(sqrt(x*x + y*y + z*z) - r + r/r0*cos(5*atan2(y,x))*cos(M_PI*z));
+
+}
+
+
+real_t new_func2d(const Vector &xx)
+{
+
+    real_t x = xx(0)-1.1;
+    real_t y = xx(1)-1.1;
+    real_t r = 0.6;
+    real_t r0 = 0.2;
+// - r0*cos(atan2(y,x)))
+    return -(sqrt(x*x + y*y) - r- r0*cos(5*atan2(y,x)));
+
+}
+
+// real_t f_rhs(const Vector &x) //koeff
+// {
+//    return sin(x(0)) * sin( x(1)) + (2*x(0)+4)*cos(x(0))*sin(x(1));
+// }
+
+
+real_t koeff(const Vector &x)
+{
+   return x(0) + 2; 
+}
+
+// real_t u_ex(const Vector &x)
+// {
+//     return cos(x(0))*sin(x(1));
+// }

miniapps/cutfem/seq_mark_test.cpp

@@ -0,0 +1,125 @@
+
+
+#include "mfem.hpp"
+#include <iostream>
+#include "cut_marking.hpp"
+
+using namespace mfem;
+using namespace std;
+
+
+
+class GyroidCoeff:public Coefficient
+{
+public:
+    GyroidCoeff(double cell_size=1.0){
+        ll=cell_size;
+    }
+
+    virtual
+        double Eval(ElementTransformation &T,
+             const IntegrationPoint &ip)
+    {
+        //evaluate the true coordinate of the ip
+        Vector xx; xx.SetSize(T.GetDimension());
+        T.Transform(ip,xx);
+        double x = xx[0]*ll;
+        double y = xx[1]*ll;
+        double z = (xx.Size()==3) ? xx[2]*ll : 0.0;
+
+        double r=std::sin(x)*std::cos(y) +
+                   std::sin(y)*std::cos(z) +
+                   std::sin(z)*std::cos(x) ;
+
+        return r;
+    }
+
+private:
+    double ll;
+};
+
+
+class BinaryGyroidCoeff:public GyroidCoeff
+{
+public:
+    BinaryGyroidCoeff(double cell_size=1.0):GyroidCoeff(cell_size)
+    {
+
+    }
+
+    virtual
+    double Eval(ElementTransformation &T,
+             const IntegrationPoint &ip)
+    {
+        double r=GyroidCoeff::Eval(T,ip);
+        if(r>0.0){return 1.0;}
+        return -1.0;
+    }
+
+};
+
+
+int main(int argc, char *argv[])
+{
+    // 1. Parse command line options.
+    string mesh_file = "../../data/star.mesh";
+    int order = 3;
+    int rs_levels = 2;
+
+
+    OptionsParser args(argc, argv);
+    args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+    args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+    args.AddOption(&rs_levels, "-rs", "--refine-serial",
+                   "Number of times to refine the mesh uniformly in serial.");
+    args.ParseCheck();
+
+    Mesh mesh(mesh_file);
+    for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
+
+
+    H1_FECollection fec(order, mesh.Dimension());
+    FiniteElementSpace fespace(&mesh, &fec);
+    cout << "Number of unknowns: " << fespace.GetTrueVSize() << endl;
+
+    GridFunction cgf(&fespace);
+    // project the Gyroid coefficient onto the grid function
+    GyroidCoeff gco(2.0*M_PI);
+    cgf.ProjectCoefficient(gco);
+
+    ElementMarker* elmark=new ElementMarker(mesh,false,true);
+    elmark->SetLevelSetFunction(cgf);
+
+    Array<int> marks;
+    elmark->MarkElements(marks);
+    Array<int> ghost_penalty_marks;
+    elmark->MarkGhostPenaltyFaces(ghost_penalty_marks);
+
+
+    //Create L2 field for marking
+    L2_FECollection* l2fec=new L2_FECollection(0,mesh.Dimension());
+    FiniteElementSpace* l2fes=new FiniteElementSpace(&mesh,l2fec,1);
+    GridFunction mgf(l2fes);
+    for(int i=0;i<marks.Size();i++){
+        mgf[i]=marks[i];
+    }
+
+    delete elmark;
+
+
+    // ParaView output.
+    ParaViewDataCollection dacol("ParaViewMarking", &mesh);
+    dacol.SetLevelsOfDetail(order);
+    dacol.SetHighOrderOutput(true);
+    dacol.RegisterField("marks", &mgf);
+    dacol.RegisterField("gyroid",&cgf);
+    dacol.SetTime(1.0);
+    dacol.SetCycle(1);
+    dacol.Save();
+
+    delete l2fes;
+    delete l2fec;
+
+    return 0;
+
+}

miniapps/cutfem/stokes.cpp

@@ -0,0 +1,229 @@
+//         first test on Stokes
+
+
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+
+
+using namespace std;
+using namespace mfem;
+
+
+void f_rhs(const Vector &x,Vector &y);
+void u_ex(const Vector &x,Vector &y);
+real_t p_ex(const Vector &x);
+
+
+
+
+int main(int argc, char *argv[])
+{
+   StopWatch chrono;
+
+   // 1. Parse command line options.
+   string mesh_file = "../../data/ref-square.mesh";
+   int order = 1;
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.ParseCheck();
+
+   Mesh mesh(mesh_file);
+   mesh.UniformRefinement();
+	mesh.UniformRefinement();
+	mesh.UniformRefinement();
+   mesh.UniformRefinement();
+   // mesh.UniformRefinement(); //TODO: how to do this without repeating
+
+   H1_FECollection fec(order+1, mesh.Dimension());
+   FiniteElementSpace Vspace(&mesh, &fec,mesh.Dimension(),Ordering::byVDIM);
+   H1_FECollection fec2(order, mesh.Dimension());
+   FiniteElementSpace Qspace(&mesh, &fec2);
+   cout << "Dimension: " << mesh.Dimension() << endl;
+
+ 
+   int dim = mesh.Dimension();
+
+   const char *device_config = "cpu";
+   Device device(device_config);
+   device.Print();
+
+   Array<int> block_offsets(3); // number of variables + 1
+   block_offsets[0] = 0;
+   block_offsets[1] = Vspace.GetVSize();
+   block_offsets[2] = Qspace.GetVSize();
+   block_offsets.PartialSum();
+
+
+   //    (u,p) and the linear forms (fform, gform).
+   MemoryType mt = device.GetMemoryType();
+   BlockVector x(block_offsets, mt), rhs(block_offsets, mt);
+   // 4. Extract the list of all the boundary DOFs. 
+   Array<int> boundary_dofs;
+   Vspace.GetBoundaryTrueDofs(boundary_dofs);
+
+   GridFunction v(&Vspace);
+   VectorFunctionCoefficient v_coeff(dim, u_ex);
+   v.ProjectCoefficient(v_coeff);
+   Array<int> ess_bdr(mesh.bdr_attributes.Max());
+   ess_bdr = 1;
+
+   // 6. Set up the linear form b(.) corresponding to the right-hand side.
+   VectorFunctionCoefficient f (dim,f_rhs);
+   // LinearForm b(&fespace);
+
+   LinearForm *fform(new LinearForm);
+   fform->Update(&Vspace, rhs.GetBlock(0), 0);
+   fform->AddDomainIntegrator(new VectorDomainLFIntegrator(f));
+   fform->Assemble();
+   fform->SyncAliasMemory(rhs);
+
+   // if rhs on second block: 
+   // LinearForm *gform(new LinearForm);
+   // gform->Update(&Qspace, rhs.GetBlock(1), 0);
+   // gform->AddDomainIntegrator(new DomainLFIntegrator(zero));
+   // gform->Assemble();
+   // gform->SyncAliasMemory(rhs);
+
+
+   // Bilinear forms: 
+   BilinearForm *mVarf(new BilinearForm(&Vspace));
+   MixedBilinearForm *bVarf(new MixedBilinearForm(&Vspace, &Qspace));
+
+   ConstantCoefficient lambda(0.0);
+   ConstantCoefficient mu(1.0);
+
+   mVarf->AddDomainIntegrator(new VectorDiffusionIntegrator);
+   mVarf->Assemble();
+   mVarf->EliminateEssentialBC(ess_bdr, v, rhs.GetBlock(0));
+   mVarf->Finalize(); 
+
+   bVarf->AddDomainIntegrator(new VectorDivergenceIntegrator);
+   bVarf->Assemble();
+   bVarf->EliminateTrialDofs(ess_bdr, v, rhs.GetBlock(1));
+   bVarf->Finalize();
+
+
+   BlockOperator StokesOp(block_offsets);
+
+   TransposeOperator *Bt = NULL;
+   SparseMatrix &Ms(mVarf->SpMat());
+   SparseMatrix &Bs(bVarf->SpMat());
+   // Bs *= -1.;
+   Bt = new TransposeOperator(&Bs);
+
+   StokesOp.SetBlock(0,0, &Ms);
+   StokesOp.SetBlock(0,1, Bt);
+   StokesOp.SetBlock(1,0, &Bs);
+
+
+   // for preconditioner 
+   BilinearForm *cVarf(new BilinearForm(&Qspace));
+   ConstantCoefficient one(1.0);
+   cVarf->AddDomainIntegrator(new MassIntegrator(one));
+   cVarf->Assemble();
+
+
+   BlockDiagonalPreconditioner stokesPrec(block_offsets);
+   Solver *solM, *solI;
+
+
+   SparseMatrix &M(mVarf->SpMat());
+   SparseMatrix &I(cVarf->SpMat());
+
+
+   solM = new DSmoother(M);
+   solI = new GSSmoother(I);
+
+   solM->iterative_mode = false;
+   solI->iterative_mode = false;
+
+   stokesPrec.SetDiagonalBlock(0, solM);
+   stokesPrec.SetDiagonalBlock(1, solI);
+
+   // 11. Solve the linear system with MINRES.
+   int maxIter(1000);
+   real_t rtol(1.e-6);
+   real_t atol(1.e-10);
+
+
+   chrono.Clear();
+   chrono.Start();
+   MINRESSolver solver;
+   solver.SetAbsTol(atol);
+   solver.SetRelTol(rtol);
+   solver.SetMaxIter(maxIter);
+   solver.SetOperator(StokesOp);
+   solver.SetPreconditioner(stokesPrec);
+   solver.SetPrintLevel(1);
+   x = 0.0;
+   solver.Mult(rhs, x);
+   if (device.IsEnabled()) { x.HostRead(); }
+   chrono.Stop();
+
+
+   if (solver.GetConverged())
+   {
+      std::cout << "MINRES converged in " << solver.GetNumIterations()
+                << " iterations with a residual norm of "
+                << solver.GetFinalNorm() << ".\n";
+   }
+   else
+   {
+      std::cout << "MINRES did not converge in " << solver.GetNumIterations()
+                << " iterations. Residual norm is " << solver.GetFinalNorm()
+                << ".\n";
+   }
+   std::cout << "MINRES solver took " << chrono.RealTime() << "s.\n";
+
+
+   GridFunction u(&Vspace);
+   GridFunction p(&Qspace);
+   u.MakeTRef(&Vspace, x.GetBlock(0), 0);
+   p.MakeRef(&Qspace, x.GetBlock(1), 0);
+
+
+   int order_quad = max(2, 2*order+1);
+   const IntegrationRule *irs[Geometry::NumGeom];
+   for (int i=0; i < Geometry::NumGeom; ++i)
+   {
+      irs[i] = &(IntRules.Get(i, order_quad));
+   }
+   FunctionCoefficient press (p_ex);
+
+   cout << "\n|| u_h - u ||_{L^2} = " << u.ComputeL2Error(v_coeff,irs) << '\n' << endl;
+   cout << "\n|| p_h - p ||_{L^2} = " << p.ComputeL2Error(press,irs) << '\n' << endl;
+
+
+   ParaViewDataCollection paraview_dc("stokes", &mesh);
+   paraview_dc.SetPrefixPath("ParaView");
+   paraview_dc.SetLevelsOfDetail(order);
+   paraview_dc.SetCycle(0);
+   paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   paraview_dc.SetHighOrderOutput(true);
+   paraview_dc.SetTime(0.0); // set the time
+   paraview_dc.RegisterField("velocity",&u);
+   paraview_dc.RegisterField("pressure",&p);
+   paraview_dc.Save();
+
+   return 0;
+}
+
+void f_rhs(const Vector &x,Vector &y)
+{
+   y(0) = sin(x(1)) + 3*x(0) * x(0), y(1) = sin(x(0)) + 3*x(1) * x(1);
+}
+
+void u_ex(const Vector &x,Vector &y)
+{
+   y(0)= sin(x(1)),y(1)=sin(x(0));
+}
+
+real_t p_ex(const Vector &x)
+{
+   return - x(0) * x(0)*x(0)  - x(1) * x(1)*x(1) + 0.5;
+}
+

miniapps/cutfem/stokes_cut.cpp

@@ -0,0 +1,394 @@
+//         first test on Stokes
+
+
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "my_integrators.hpp"
+
+
+
+using namespace std;
+using namespace mfem;
+
+
+void f_rhs(const Vector &x,Vector &y);
+void u_ex(const Vector &x,Vector &y);
+real_t p_ex(const Vector &x);
+real_t ellipsoide_func(const Vector &x);
+void g_neumann(const Vector &x,Vector &z );
+
+
+
+int main(int argc, char *argv[])
+{
+   StopWatch chrono;
+
+   // 1. Parse command line options.
+   int n = 20;
+   int order = 1;
+   OptionsParser args(argc, argv);
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+      args.AddOption(&n, "-n", "--n", "n");
+
+   args.ParseCheck();
+   Mesh mesh = Mesh::MakeCartesian2D( n*2, n , mfem::Element::Type::QUADRILATERAL, true, 1, 0.5);
+
+   double h_min, h_max, kappa_min, kappa_max;
+   mesh.GetCharacteristics(h_min, h_max, kappa_min, kappa_max);
+
+   // mesh.UniformRefinement(); //TODO: how to do this without repeating
+
+   H1_FECollection fec(order+1, mesh.Dimension());
+   FiniteElementSpace Vspace(&mesh, &fec,mesh.Dimension(),Ordering::byVDIM);
+   H1_FECollection fec2(order, mesh.Dimension());
+   FiniteElementSpace Qspace(&mesh, &fec2);
+   cout << "Dimension: " << mesh.Dimension() << endl;
+
+   GridFunction cgf(&Qspace);
+   FunctionCoefficient circle(ellipsoide_func);
+   cgf.ProjectCoefficient(circle);
+ 
+
+
+
+   int dim = mesh.Dimension();
+
+   const char *device_config = "cpu";
+   Device device(device_config);
+   device.Print();
+
+   Array<int> block_offsets(3); // number of variables + 1
+   block_offsets[0] = 0;
+   block_offsets[1] = Vspace.GetVSize();
+   block_offsets[2] = Qspace.GetVSize();
+   block_offsets.PartialSum();
+
+
+   //    (u,p) and the linear forms (fform, gform).
+   MemoryType mt = device.GetMemoryType();
+   BlockVector x(block_offsets, mt), rhs(block_offsets, mt);
+   // 4. Extract the list of all the boundary DOFs. 
+
+   GridFunction v(&Vspace);
+   VectorFunctionCoefficient v_coeff(dim, u_ex);
+   v.ProjectCoefficient(v_coeff);
+
+
+   VectorFunctionCoefficient neumann(dim,g_neumann);
+   Array<int> boundary_dofs(mesh.bdr_attributes.Max());
+   boundary_dofs = 1;
+
+
+   Array<int> outside_dofs_v;
+   Array<int> outside_dofs_p;
+   Array<int> marks;
+   Array<int> face_marks;
+   {
+       ElementMarker* elmark=new ElementMarker(mesh,true,true); 
+       elmark->SetLevelSetFunction(cgf);
+       elmark->MarkElements(marks);
+       elmark->ListEssentialTDofs(marks,Vspace,outside_dofs_v);
+       elmark->ListEssentialTDofs(marks,Qspace,outside_dofs_p);
+      elmark->MarkGhostPenaltyFaces(face_marks);
+       delete elmark;
+   }
+   //  outside_dofs.Append(boundary_dofs);
+   //  outside_dofs.Sort();
+   //  outside_dofs.Unique();
+
+   int otherorder = 2;
+   int aorder = 8; // Algoim integration points
+   AlgoimIntegrationRules* air=new AlgoimIntegrationRules(aorder,circle,otherorder);
+   real_t gp1 = 0.1/(h_min*h_min);
+   real_t gp2 = -0.1;
+
+   // 6. Set up the linear form b(.) corresponding to the right-hand side.
+   VectorFunctionCoefficient f (dim,f_rhs);
+   // LinearForm b(&fespace);
+
+   LinearForm *fform(new LinearForm);
+   fform->Update(&Vspace, rhs.GetBlock(0), 0);
+   // fform->AddDomainIntegrator(new VectorDomainLFIntegrator(f));
+   fform->AddDomainIntegrator(new CutVectorDomainLFIntegrator(f,&marks,air));
+   fform->AddDomainIntegrator(new CutUnfittedVectorBoundaryLFIntegrator(neumann,&marks,air));
+
+   fform->Assemble();
+   fform->SyncAliasMemory(rhs);
+
+   ConstantCoefficient zero(0.0);
+
+   // if rhs on second block: 
+   // LinearForm *gform(new LinearForm);
+   // gform->Update(&Qspace, rhs.GetBlock(1), 0);
+   // gform->AddDomainIntegrator(new CutDomainLFIntegrator(zero,&marks,air));
+   // gform->Assemble();
+   // gform->SyncAliasMemory(rhs);
+   
+   ConstantCoefficient one(1.0);
+
+
+   // Bilinear forms: 
+   BilinearForm *mVarf(new BilinearForm(&Vspace));
+   MixedBilinearForm *bVarf(new MixedBilinearForm(&Vspace, &Qspace));
+   BilinearForm *cVarf(new BilinearForm(&Qspace));
+
+   mVarf->SetDiagonalPolicy(Operator::DiagonalPolicy::DIAG_KEEP);
+   cVarf->SetDiagonalPolicy(Operator::DiagonalPolicy::DIAG_KEEP);
+
+
+   ConstantCoefficient lambda(0.0);
+   ConstantCoefficient mu(1.0);
+
+   // mVarf->AddDomainIntegrator(new VectorDiffusionIntegrator);
+   mVarf->AddDomainIntegrator(new CutVectorDiffusionIntegrator(one,&marks,air));
+   mVarf->AddInteriorFaceIntegrator(new CutGhostPenaltyVectorIntegrator(gp1,&face_marks));
+   mVarf->Assemble(0);
+   mVarf->EliminateEssentialBC(boundary_dofs, v, rhs.GetBlock(0));
+   // smVarf->EliminateVDofs(outside_dofs_v);
+   mVarf->Finalize(0);
+
+
+
+   bVarf->AddDomainIntegrator(new CutVectorDivergenceIntegrator(one,&marks,air));
+   bVarf->Assemble(0);
+   bVarf->EliminateTrialDofs(boundary_dofs, v, rhs.GetBlock(1));
+   bVarf->Finalize(0);
+
+   cVarf->AddInteriorFaceIntegrator(new CutGhostPenaltyIntegrator(gp2,&face_marks));
+   cVarf->Assemble(0);
+   cVarf->Finalize(0);
+
+
+   BlockOperator StokesOp(block_offsets);
+
+   TransposeOperator *Bt = NULL;
+   SparseMatrix &Ms(mVarf->SpMat());
+   SparseMatrix &Bs(bVarf->SpMat());
+   SparseMatrix &C(cVarf->SpMat());
+
+
+   cout <<  outside_dofs_p.Size()<<": "<<  outside_dofs_v.Size()  << endl;
+
+   Vector before;
+   C.GetDiag(before);
+   cout <<"before"<<endl;
+   //before.Print();
+
+   for(int i=0;i<outside_dofs_v.Size();i++)
+   {
+      // cout<<outside_dofs_v[i]<<endl;
+      Ms.EliminateRowColDiag(outside_dofs_v[i],1.0);
+   }
+
+   //Ms.SetDiagIdentity();
+
+
+
+   for(int i=0;i<outside_dofs_p.Size();i++)
+   {
+      C.EliminateRowCol(outside_dofs_p[i]);
+   }
+
+
+   // Bs *= -1.;
+   Bt = new TransposeOperator(&Bs);
+
+   StokesOp.SetBlock(0,0, &Ms);
+   StokesOp.SetBlock(0,1, Bt);
+   StokesOp.SetBlock(1,0, &Bs);
+   StokesOp.SetBlock(1,1, &C);
+
+   // rhs.Print();
+   // C.Print();
+   // Bs.Print();
+
+   // for preconditioner 
+   BilinearForm *mpre(new BilinearForm(&Vspace));
+   BilinearForm *cpre(new BilinearForm(&Qspace));
+   cpre->AddDomainIntegrator(new MassIntegrator);
+   cpre->Assemble();
+   cpre->Finalize(); 
+
+   BlockDiagonalPreconditioner stokesPrec(block_offsets);
+   Solver *solM, *solI;
+
+   mpre->AddDomainIntegrator(new VectorDiffusionIntegrator);
+   mpre->Assemble();
+   mpre->Finalize(); 
+
+
+   SparseMatrix &M(mpre->SpMat());
+   SparseMatrix &I(cpre->SpMat());
+
+
+   solM = new GSSmoother(Ms);
+   solI = new GSSmoother(I);
+
+
+  Vector id_vec(bVarf->Height());
+   id_vec = 1.0;
+   SparseMatrix id_mat(id_vec);
+
+
+
+   stokesPrec.SetDiagonalBlock(0, solM);
+   stokesPrec.SetDiagonalBlock(1, &id_mat);
+
+   // 11. Solve the linear system with MINRES.
+   int maxIter(30000);
+   real_t rtol(1.e-15);
+   real_t atol(1.e-10);
+
+
+   chrono.Clear();
+   chrono.Start();
+   MINRESSolver solver;
+   solver.SetAbsTol(atol);
+   solver.SetRelTol(rtol);
+   solver.SetMaxIter(maxIter);
+   solver.SetOperator(StokesOp);
+   solver.SetPreconditioner(stokesPrec);
+   solver.SetPrintLevel(2);
+   x = 0.0;
+   solver.Mult(rhs, x);
+   if (device.IsEnabled()) { x.HostRead(); }
+   chrono.Stop();
+
+
+   if (solver.GetConverged())
+   {
+      std::cout << "MINRES converged in " << solver.GetNumIterations()
+                << " iterations with a residual norm of "
+                << solver.GetFinalNorm() << ".\n";
+   }
+   else
+   {
+      std::cout << "MINRES did not converge in " << solver.GetNumIterations()
+                << " iterations. Residual norm is " << solver.GetFinalNorm()
+                << ".\n";
+   }
+   std::cout << "MINRES solver took " << chrono.RealTime() << "s.\n";
+
+
+   GridFunction u(&Vspace);
+   GridFunction p(&Qspace);
+   u.MakeTRef(&Vspace, x.GetBlock(0), 0);
+   p.MakeRef(&Qspace, x.GetBlock(1), 0);
+
+
+   int order_quad = max(5, 2*order+1);
+   const IntegrationRule *irs[Geometry::NumGeom];
+   for (int i=0; i < Geometry::NumGeom; ++i)
+   {
+      irs[i] = &(IntRules.Get(i, order_quad));
+   }
+   FunctionCoefficient press (p_ex);
+
+      //compute the error
+   {
+       NonlinearForm* nf=new NonlinearForm(&Vspace);
+       nf->AddDomainIntegrator(new CutVectorErrorIntegrator(v_coeff,&marks,air));
+        real_t error_squared = nf->GetEnergy(u.GetTrueVector());
+       cout << "\n|| u_h - u ||_{L^2} = " << sqrt(error_squared)<< std::endl;
+
+         NonlinearForm* nf2=new NonlinearForm(&Qspace);
+       nf2->AddDomainIntegrator(new CutScalarErrorIntegrator(press,&marks,air));
+        real_t error_squared_p = nf2->GetEnergy(p.GetTrueVector());
+      cout << "\n|| p_h - p ||_{L^2} = " << sqrt(error_squared_p) << '\n' << endl;
+
+       delete nf;
+              delete nf2;
+   }
+
+   // cout << "\n|| u_h - u ||_{L^2} = " << u.ComputeL2Error(v_coeff,irs) << '\n' << endl;
+   // cout << "\n|| p_h - p ||_{L^2} = " << p.ComputeL2Error(press,irs) << '\n' << endl;
+
+   // to visualize level set and markings
+   L2_FECollection* l2fec= new L2_FECollection(0,mesh.Dimension());
+   FiniteElementSpace* l2fes= new FiniteElementSpace(&mesh,l2fec,1);
+   GridFunction mgf(l2fes);
+   for(int i=0;i<marks.Size();i++){
+      mgf[i]=marks[i];
+   }
+
+   GridFunction exact_p(&Qspace);
+   exact_p.ProjectCoefficient(press);
+
+   GridFunction error_u(&Vspace);
+   error_u = u;
+   error_u -= v;
+
+   GridFunction error_p(&Qspace);
+   error_p = p;
+   error_p -= exact_p;
+
+   ParaViewDataCollection paraview_dc("stokes_cut", &mesh);
+   paraview_dc.SetPrefixPath("ParaView");
+   paraview_dc.SetLevelsOfDetail(order);
+   paraview_dc.SetCycle(0);
+   paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   paraview_dc.SetHighOrderOutput(true);
+   paraview_dc.SetTime(0.0); // set the time
+   paraview_dc.RegisterField("velocity",&u);
+   paraview_dc.RegisterField("pressure",&p);
+   paraview_dc.RegisterField("marks", &mgf);
+   paraview_dc.RegisterField("exact_v", &v);
+      paraview_dc.RegisterField("exact_p", &exact_p);
+
+      paraview_dc.RegisterField("level_set",&cgf);
+            paraview_dc.RegisterField("error_u", &error_u);
+
+      paraview_dc.RegisterField("error_p",&error_p);
+
+   paraview_dc.Save();
+
+   return 0;
+}
+
+void f_rhs(const Vector &x,Vector &y)
+{
+   y(0) = sin(x(1)) + 3*x(0) * x(0), y(1) = sin(x(0)) + 3*x(1) * x(1);
+   // y(0) = sin(x(1)), y(1) = sin(x(0)) ;
+}
+
+void u_ex(const Vector &x,Vector &y)
+{
+   y(0)= sin(x(1)),y(1)=sin(x(0));
+}
+void g_neumann(const Vector &x,Vector &z )
+{
+   real_t a = 1.5;
+    real_t b = 0.5;
+    real_t x0 = 0.5;
+    real_t y0 = 0.25;
+    real_t xx = x(0)-x0;
+    real_t y = x(1)-y0;
+    real_t normalize = sqrt((xx*xx)/(a*a*a*a) + y*y/(b*b*b*b));
+   z(0) =cos(x(1))*y/(b*b*normalize)   + (- x(0) * x(0)*x(0)  - x(1) * x(1)*x(1) + 0.5)*xx/(a*a*normalize); 
+   z(1) =cos(x(0))*xx/(a*a*normalize)  + (-x(0) * x(0)*x(0)  - x(1) * x(1)*x(1) + 0.5)*y/(b*b*normalize);
+}
+
+real_t p_ex(const Vector &x)
+{
+   return - x(0) * x(0)*x(0)  - x(1) * x(1)*x(1) + 0.5;
+}
+
+real_t circle_func(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0.5;
+    real_t r = 0.4;
+    return -(x(0)-x0)*(x(0)-x0) - (x(1)-y0)*(x(1)-y0) + r*r;  
+}
+
+real_t ellipsoide_func(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0.25;
+    real_t r = 0.35;
+    real_t xx = x(0)-x0;
+    real_t y = x(1)-y0;
+    return -(xx)*(xx)/(1.5*1.5) - (y)*(y)/(0.5*0.5)+ r*r; // + 0.25*cos(atan2(x(1)-y0,x(0)-x0))*cos(atan2(x(1)-y0,x(0)-x0));  
+}

miniapps/cutfem/vector_diffusion.cpp

@@ -0,0 +1,134 @@
+//       Solving the vectorized version of the Laplace  problem -Delta u = 1 with 
+//       Dirichlet boundary conditions
+
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+
+
+using namespace std;
+using namespace mfem;
+
+
+void f_rhs(const Vector &x,Vector &y);
+void u_ex(const Vector &x,Vector &y);
+
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command line options.
+   string mesh_file = "../../data/star.mesh";
+   int order = 2;
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.ParseCheck();
+
+   // 2. Read the mesh from the given mesh file, and refine once uniformly.
+   Mesh mesh(mesh_file);
+   mesh.UniformRefinement();
+	mesh.UniformRefinement();
+
+   // 3. Define a finite element space on the mesh. Here we use H1 continuous
+   //    high-order Lagrange finite elements of the given order.
+   H1_FECollection fec(order, mesh.Dimension());
+   FiniteElementSpace fespace(&mesh, &fec,mesh.Dimension(),Ordering::byVDIM);
+   cout << "Dimension: " << mesh.Dimension() << endl;
+
+ 
+   int dim = mesh.Dimension();
+
+   cout << "Number of unknowns: " << fespace.GetTrueVSize() << endl;
+
+   // 4. Extract the list of all the boundary DOFs. These will be marked as
+   //    Dirichlet in order to enforce zero boundary conditions.
+   Array<int> boundary_dofs;
+   fespace.GetBoundaryTrueDofs(boundary_dofs);
+
+   // 5. Define the solution x as a finite element grid function in fespace. Set
+   //    the initial guess to zero, which also sets the boundary conditions.
+   GridFunction x(&fespace);
+   VectorFunctionCoefficient bc(dim, u_ex);
+   x.ProjectCoefficient(bc);
+
+   // 6. Set up the linear form b(.) corresponding to the right-hand side.
+   VectorFunctionCoefficient f (dim,f_rhs);
+   LinearForm b(&fespace);
+
+   ConstantCoefficient one(1.0);
+
+   b.AddDomainIntegrator(new VectorDomainLFIntegrator(f));
+   b.Assemble();
+
+   // 7. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
+   BilinearForm a(&fespace);
+
+   ConstantCoefficient lambda(0.0);
+   ConstantCoefficient mu(1.0);
+
+   a.AddDomainIntegrator(new ElasticityIntegrator(lambda,mu));
+   // a.AddDomainIntegrator(new VectorDiffusionIntegrator);
+
+   a.Assemble();
+
+   // 8. Form the linear system A X = B. This includes eliminating boundary
+   //    conditions, applying AMR constraints, and other transformations.
+   SparseMatrix A;
+   Vector B, X;
+   a.FormLinearSystem(boundary_dofs, x, b, A, X, B);
+
+   // 9. Solve the system using PCG with symmetric Gauss-Seidel preconditioner.
+   GSSmoother M(A);
+   PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
+
+   // 10. Recover the solution x as a grid function and save to file. The output
+   //     can be viewed using GLVis as follows: "glvis -m mesh.mesh -g sol.gf"
+   a.RecoverFEMSolution(X, b, x);
+   // x.Save("sol.gf");
+   // mesh.Save("mesh.mesh");
+
+
+   cout << "\n|| u_h - u ||_{L^2} = " << x.ComputeL2Error(bc) << '\n' << endl;
+   // cout << "\n|| grad u_h - grad u ||_{L^2} = " << x.ComputeH1Error(&bc,&u_grad) << '\n' << endl;
+   // todo: find out how to compute H1 for vector case
+
+   ParaViewDataCollection paraview_dc("vecdiffusion", &mesh);
+   paraview_dc.SetPrefixPath("ParaView");
+   paraview_dc.SetLevelsOfDetail(order);
+   paraview_dc.SetCycle(0);
+   paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   paraview_dc.SetHighOrderOutput(true);
+   paraview_dc.SetTime(0.0); // set the time
+   paraview_dc.RegisterField("velocity",&x);
+   paraview_dc.Save();
+
+   return 0;
+}
+// for gradient case
+// void f_rhs(const Vector &x,Vector &y)
+// {
+//    y(0) = sin(x(1)), y(1) = sin(x(0));
+// }
+
+// void u_ex(const Vector &x,Vector &y)
+// {
+//    y(0)= sin(x(1)),y(1)=sin(x(0));
+// }
+
+
+// for symmetric gradient
+void f_rhs(const Vector &x,Vector &y)
+{
+   y(0) = 3*cos(x(0))*sin(x(1)) - cos(x(0))*cos(x(1)), y(1)=3*sin(x(0))*sin(x(1))+sin(x(0))*cos(x(1));
+}
+
+void u_ex(const Vector &x,Vector &y)
+{
+   y(0)= cos(x(0))*sin(x(1)),y(1)=sin(x(0))*sin(x(1));
+}
+
+
+

miniapps/cutfem/vector_diffusion_cut.cpp

@@ -0,0 +1,234 @@
+//       Solving the vectorized version of the Laplace  problem -Delta u = 1 with 
+//       Dirichlet boundary conditions
+
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "my_integrators.hpp"
+
+
+
+using namespace std;
+using namespace mfem;
+
+
+void f_rhs(const Vector &x,Vector &y);
+void u_ex(const Vector &x,Vector &y);
+void g_neumann(const Vector &x,Vector &z);
+real_t circle_func(const Vector &x);
+real_t ellipsoide_func(const Vector &x);
+
+
+
+int main(int argc, char *argv[])
+{
+   
+
+
+   int n =30;
+
+
+   int order = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.AddOption(&n, "-n", "--n", "n");
+   args.ParseCheck();
+   Mesh mesh = Mesh::MakeCartesian2D( n, n , mfem::Element::Type::QUADRILATERAL, true, 1, 1);
+
+   int dim = mesh.Dimension();
+   double h_min, h_max, kappa_min, kappa_max;
+   mesh.GetCharacteristics(h_min, h_max, kappa_min, kappa_max);
+
+   H1_FECollection fec(order, mesh.Dimension());
+   FiniteElementSpace fespace(&mesh, &fec,mesh.Dimension(),Ordering::byNODES);
+
+   FiniteElementSpace fespacescalar(&mesh, &fec);
+
+   // cout << "Number of unknowns: " << fespace.GetTrueVSize() << endl;
+
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient zero(1.0);
+
+
+   GridFunction x(&fespace);
+   VectorFunctionCoefficient bc (dim,u_ex);
+   VectorFunctionCoefficient f (dim,f_rhs);
+   x.ProjectCoefficient(bc);
+   VectorFunctionCoefficient neumann(dim,g_neumann);
+
+
+   // level set function 
+   GridFunction cgf(&fespacescalar);
+   FunctionCoefficient circle(circle_func);
+   cgf.ProjectCoefficient(circle);
+
+   // mark elements and outside DOFs
+   Array<int> boundary_dofs;
+   fespace.GetBoundaryTrueDofs(boundary_dofs);
+   // Array<int> outside_dofs;
+   Array<int> marks;
+      Array<int> face_marks;
+   {
+      //  Array<int> outside_dofs;
+       ElementMarker* elmark=new ElementMarker(mesh,true,true); 
+       elmark->SetLevelSetFunction(cgf);
+       elmark->MarkElements(marks);
+      elmark->MarkGhostPenaltyFaces(face_marks);
+      //  elmark->ListEssentialTDofs(marks,fespace,outside_dofs);
+       delete elmark;
+   }
+   //  outside_dofs.Append(boundary_dofs);
+   //  outside_dofs.Sort();
+   //  outside_dofs.Unique();
+
+   int otherorder = 2;
+   int aorder = 4; // Algoim integration points
+   AlgoimIntegrationRules* air=new AlgoimIntegrationRules(aorder,circle,otherorder);
+   real_t gp = 10/(h_min*h_min);
+   real_t lambda = 20/h_min;
+
+   // 6. Set up the linear form b(.) corresponding to the right-hand side.
+
+   LinearForm b(&fespace);
+
+   b.AddDomainIntegrator(new CutVectorDomainLFIntegrator(f,&marks,air));
+   // b.AddDomainIntegrator(new CutUnfittedVectorBoundaryLFIntegrator(neumann,&marks,air));
+   b.AddDomainIntegrator(new  CutUnfittedNitscheSymmetricLFIntegrator(bc,one,lambda,&marks,air));
+   b.Assemble();
+
+   // 7. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
+   BilinearForm a(&fespace);
+   a.AddDomainIntegrator(new CutElasticityIntegrator(one,&marks,air));
+   // a.AddDomainIntegrator(new CutNitscheVectorIntegrator(one,lambda,&marks,air));
+      a.AddDomainIntegrator(new CutNitscheSymmetricIntegrator(one,lambda,&marks,air));
+
+   a.AddInteriorFaceIntegrator(new CutGhostPenaltyVectorIntegrator(gp,&face_marks));
+   a.Assemble();
+
+   
+   SparseMatrix A;
+   Vector B, X;
+   a.FormLinearSystem(boundary_dofs, x, b, A, X, B);
+
+   // A.Print();
+   // 9. Solve the system using PCG with symmetric Gauss-Seidel preconditioner.
+   GSSmoother M(A);
+   PCG(A,M, B, X, 2, 2000, 1e-20, 0.0);
+
+   // 10. Recover the solution x as a grid function and save to file
+   a.RecoverFEMSolution(X, b, x);
+
+
+   // compute the error
+   {
+       NonlinearForm* nf=new NonlinearForm(&fespace);
+       nf->AddDomainIntegrator(new  CutVectorErrorIntegrator(bc,&marks,air));
+
+        real_t error_squared = nf->GetEnergy(x.GetTrueVector());
+       cout << "\n|| u_h - u ||_{L^2} = " << sqrt(error_squared)<< std::endl;
+
+       delete nf;
+   }
+
+   // to visualize level set and markings
+   L2_FECollection* l2fec= new L2_FECollection(0,mesh.Dimension());
+   FiniteElementSpace* l2fes= new FiniteElementSpace(&mesh,l2fec,1);
+   GridFunction mgf(l2fes);
+   for(int i=0;i<marks.Size();i++){
+      mgf[i]=marks[i];
+   }
+   GridFunction exact_sol(&fespace);
+   exact_sol.ProjectCoefficient(bc);
+
+
+   GridFunction error(&fespace);
+   error = x;
+   error -= exact_sol;
+
+   //  // GLVis 
+   // char vishost[] = "localhost";
+   // int  visport   = 19916;
+   // socketstream sol_sock(vishost, visport);
+   // sol_sock.precision(8);
+   // sol_sock << "solution\n" << mesh << x << flush;
+   // //    // save solution with paraview 
+   ParaViewDataCollection paraview_dc("vector_cut", &mesh);
+   paraview_dc.SetPrefixPath("ParaView");
+   paraview_dc.SetLevelsOfDetail(order);
+   paraview_dc.SetCycle(0);
+   paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   paraview_dc.SetHighOrderOutput(true);
+   paraview_dc.SetTime(0.0); // set the time
+   paraview_dc.RegisterField("solution",&x);
+   paraview_dc.RegisterField("marks", &mgf);
+   paraview_dc.RegisterField("level_set",&cgf);
+   paraview_dc.RegisterField("exact_sol",&exact_sol);
+   paraview_dc.RegisterField("error",&error);
+   paraview_dc.Save();
+
+
+   delete l2fec;
+   delete l2fes;
+   delete air;
+
+   cout << "h:" << h_min<<endl;
+   return 0;  
+}
+//for gradient case
+// void f_rhs(const Vector &x,Vector &y)
+// {
+//    y(0)= sin(x(0)),y(1)=sin(x(0));//sin(x(0));
+// }
+
+// void u_ex(const Vector &x,Vector &y)
+// {
+//    y(0)= sin(x(0)),y(1)=sin(x(0));
+// }
+
+void g_neumann(const Vector &x, Vector &z)
+{
+    real_t a = 1.5;
+    real_t b = 0.5;
+    real_t x0 = 0.5;
+    real_t y0 = 0.25;
+    real_t xx = x(0)-x0;
+    real_t y = x(1)-y0;
+    real_t normalize = sqrt((xx*xx)/(a*a*a*a) + y*y/(b*b*b*b));
+   //  z(0) = cos(x(1))*y/(b*b*normalize), z(1)=cos(x(0))*xx/(a*a*normalize);
+   // z(0) =cos(x(1))*y/(b*b*normalize), z(1)=cos(x(0))*xx/(a*a*normalize);
+   z(0) =cos(x(0))*xx/(a*a*normalize), z(1)=cos(x(0))*xx/(a*a*normalize);
+}
+
+
+real_t circle_func(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0.5;
+    real_t r = 0.4;
+    return -(x(0)-x0)*(x(0)-x0) - (x(1)-y0)*(x(1)-y0) + r*r;  
+}
+
+real_t ellipsoide_func(const Vector &x)
+{
+    real_t x0 = 0.5;
+    real_t y0 = 0.25;
+    real_t r = 0.35;
+    real_t xx = x(0)-x0;
+    real_t y = x(1)-y0;
+    return -(xx)*(xx)/(1.5*1.5) - (y)*(y)/(0.5*0.5)+ r*r; // + 0.25*cos(atan2(x(1)-y0,x(0)-x0))*cos(atan2(x(1)-y0,x(0)-x0));  
+}
+// for symmetric gradient
+void f_rhs(const Vector &x,Vector &y)
+{
+   y(0) = 3*cos(x(0))*sin(x(1)) - cos(x(0))*cos(x(1)), y(1)=3*sin(x(0))*sin(x(1))+sin(x(0))*cos(x(1));
+}
+
+void u_ex(const Vector &x,Vector &y)
+{
+   y(0)= cos(x(0))*sin(x(1)),y(1)=sin(x(0))*sin(x(1));
+}
+
+
+