diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7385d711b..7abf33366 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -287,6 +287,7 @@ if(MPM_BUILD_TESTING)
     ${mpm_SOURCE_DIR}/tests/test_main.cc
     ${mpm_SOURCE_DIR}/tests/io/write_mesh_particles.cc
     ${mpm_SOURCE_DIR}/tests/io/write_mesh_particles_unitcell.cc
+    ${mpm_SOURCE_DIR}/tests/solvers/mpm_explicit_levelset_test.cc
     ${mpm_SOURCE_DIR}/tests/solvers/mpm_explicit_constraint_absorbing_test.cc
     ${mpm_SOURCE_DIR}/tests/solvers/mpm_explicit_constraint_acceleration_test.cc
     ${mpm_SOURCE_DIR}/tests/solvers/mpm_explicit_constraint_friction_test.cc
diff --git a/include/contacts/contact.h b/include/contacts/contact.h
index 98e6b1f69..27cb3a81f 100644
--- a/include/contacts/contact.h
+++ b/include/contacts/contact.h
@@ -14,12 +14,24 @@ class Contact {
   //! Default constructor with mesh class
   Contact(const std::shared_ptr<mpm::Mesh<Tdim>>& mesh);
 
-  //! Intialize
+  //! Intialise
   virtual inline void initialise(){};
 
-  //! Compute contact forces
+  //! Initialise levelset properties
+  //! \param[in] levelset_damping Levelset damping factor
+  //! \param[in] levelset_pic Particle in cell method bool for contact velocity
+  //! \param[in] levelset_violation_corrector Violation correction factor
+  virtual inline void initialise_levelset_properties(
+      double levelset_damping, bool levelset_pic,
+      double levelset_violation_corrector){};
+
+  //! Compute contact reaction forces
   virtual inline void compute_contact_forces(){};
 
+  //! Compute contact reaction forces
+  //! \param[in] dt Analysis time step
+  virtual inline void compute_contact_forces(double dt){};
+
  protected:
   //! Mesh object
   std::shared_ptr<mpm::Mesh<Tdim>> mesh_;
diff --git a/include/contacts/contact_friction.h b/include/contacts/contact_friction.h
index 5c9bf6cec..466c8dd16 100644
--- a/include/contacts/contact_friction.h
+++ b/include/contacts/contact_friction.h
@@ -15,10 +15,10 @@ class ContactFriction : public Contact<Tdim> {
   ContactFriction(const std::shared_ptr<mpm::Mesh<Tdim>>& mesh);
 
   //! Intialize
-  virtual inline void initialise() override;
+  void initialise() override;
 
-  //! Compute contact forces
-  virtual inline void compute_contact_forces() override;
+  //! Compute contact reaction forces
+  void compute_contact_forces() override;
 
  protected:
   //! Mesh object
diff --git a/include/contacts/contact_friction.tcc b/include/contacts/contact_friction.tcc
index 6097c7829..83d6097dd 100644
--- a/include/contacts/contact_friction.tcc
+++ b/include/contacts/contact_friction.tcc
@@ -4,9 +4,9 @@ mpm::ContactFriction<Tdim>::ContactFriction(
     const std::shared_ptr<mpm::Mesh<Tdim>>& mesh)
     : mpm::Contact<Tdim>(mesh) {}
 
-//! Initialize nodal properties
+//! Initialise nodal properties
 template <unsigned Tdim>
-inline void mpm::ContactFriction<Tdim>::initialise() {
+void mpm::ContactFriction<Tdim>::initialise() {
   // Initialise nodal properties
   mesh_->initialise_nodal_properties();
 
@@ -16,9 +16,9 @@ inline void mpm::ContactFriction<Tdim>::initialise() {
                 std::placeholders::_1));
 }
 
-//! Compute contact forces
+//! Compute contact reaction forces
 template <unsigned Tdim>
-inline void mpm::ContactFriction<Tdim>::compute_contact_forces() {
+void mpm::ContactFriction<Tdim>::compute_contact_forces() {
 
   // Map multimaterial properties from particles to nodes
   mesh_->iterate_over_particles(std::bind(
@@ -49,4 +49,4 @@ inline void mpm::ContactFriction<Tdim>::compute_contact_forces() {
   mesh_->iterate_over_nodes(
       std::bind(&mpm::NodeBase<Tdim>::compute_multimaterial_normal_unit_vector,
                 std::placeholders::_1));
-}
+}
\ No newline at end of file
diff --git a/include/contacts/contact_levelset.h b/include/contacts/contact_levelset.h
new file mode 100644
index 000000000..9b0cd7d1b
--- /dev/null
+++ b/include/contacts/contact_levelset.h
@@ -0,0 +1,48 @@
+#ifndef MPM_CONTACT_LEVELSET_H_
+#define MPM_CONTACT_LEVELSET_H_
+
+#include "contact.h"
+#include "mesh_levelset.h"
+#include "particle_levelset.h"
+
+namespace mpm {
+
+//! ContactLevelset class
+//! \brief ContactLevelset base class
+//! \tparam Tdim Dimension
+template <unsigned Tdim>
+class ContactLevelset : public Contact<Tdim> {
+ public:
+  //! Default constructor with mesh class
+  ContactLevelset(const std::shared_ptr<mpm::Mesh<Tdim>>& mesh);
+
+  //! Initialise levelset properties
+  //! \param[in] levelset_damping Levelset damping factor
+  //! \param[in] levelset_pic Particle in cell method bool for contact velocity
+  //! \param[in] levelset_violation_corrector Violation correction factor
+  void initialise_levelset_properties(
+      double levelset_damping, bool levelset_pic,
+      double levelset_violation_corrector) override;
+
+  //! Compute contact reaction forces
+  //! \param[in] dt Analysis time step
+  void compute_contact_forces(double dt) override;
+
+  //! Mesh levelset object
+  std::shared_ptr<mpm::MeshLevelset<Tdim>> mesh_levelset_;
+
+ protected:
+  //! levelset_damping_
+  double levelset_damping_;
+  //! levelset_pic_
+  bool levelset_pic_;
+  //! levelset_violation_corrector_
+  double levelset_violation_corrector_;
+  //! Mesh object
+  using mpm::Contact<Tdim>::mesh_;
+};  // ContactLevelset class
+}  // namespace mpm
+
+#include "contact_levelset.tcc"
+
+#endif  // MPM_CONTACT_LEVELSET_H_
diff --git a/include/contacts/contact_levelset.tcc b/include/contacts/contact_levelset.tcc
new file mode 100644
index 000000000..a519acdf2
--- /dev/null
+++ b/include/contacts/contact_levelset.tcc
@@ -0,0 +1,25 @@
+//! Constructor of contact with mesh
+template <unsigned Tdim>
+mpm::ContactLevelset<Tdim>::ContactLevelset(
+    const std::shared_ptr<mpm::Mesh<Tdim>>& mesh)
+    : mpm::Contact<Tdim>(mesh) {
+  // Assign mesh
+  mesh_ = mesh;
+}
+
+//! Initialise levelset properties
+template <unsigned Tdim>
+void mpm::ContactLevelset<Tdim>::initialise_levelset_properties(
+    double levelset_damping, bool levelset_pic,
+    double levelset_violation_corrector) {
+  mpm::ParticleLevelset<Tdim>::set_levelset_properties(
+      levelset_damping, levelset_pic, levelset_violation_corrector);
+}
+
+//! Compute contact reaction forces
+template <unsigned Tdim>
+void mpm::ContactLevelset<Tdim>::compute_contact_forces(double dt) {
+  mesh_->iterate_over_particles(
+      std::bind(&mpm::ParticleBase<Tdim>::levelset_contact_force,
+                std::placeholders::_1, dt));
+}
\ No newline at end of file
diff --git a/include/io/io_mesh.h b/include/io/io_mesh.h
index d7e6da99a..d7b9648e9 100644
--- a/include/io/io_mesh.h
+++ b/include/io/io_mesh.h
@@ -64,11 +64,11 @@ class IOMesh {
   virtual std::vector<Eigen::Matrix<double, 6, 1>> read_particles_stresses(
       const std::string& particles_stresses) = 0;
 
-  //! Read particle scalar properties
+  //! Read scalar properties for particles or nodes
   //! \param[in] scalar_file file name with particle scalar properties
-  //! \retval Vector of particles scalar properties
-  virtual std::vector<std::tuple<mpm::Index, double>>
-      read_particles_scalar_properties(const std::string& scalar_file) = 0;
+  //! \retval Vector of scalar properties for particles or nodes
+  virtual std::vector<std::tuple<mpm::Index, double>> read_scalar_properties(
+      const std::string& scalar_file) = 0;
 
   //! Read pressure constraints file
   //! \param[in] pressure_constraints_files file name with pressure constraints
@@ -123,6 +123,11 @@ class IOMesh {
       read_adhesion_constraints(
           const std::string& adhesion_constraints_file) = 0;
 
+  //! Read levelset file
+  //! \param[in] levelset_input_file file name with levelset values
+  virtual std::vector<std::tuple<mpm::Index, double, double, double, double>>
+      read_levelset_input(const std::string& levelset_input_file) = 0;
+
   //! Read forces file
   //! \param[in] forces_file file name with nodal concentrated force
   virtual std::vector<std::tuple<mpm::Index, unsigned, double>> read_forces(
diff --git a/include/io/io_mesh_ascii.h b/include/io/io_mesh_ascii.h
index 9ce82ad82..cf398d167 100644
--- a/include/io/io_mesh_ascii.h
+++ b/include/io/io_mesh_ascii.h
@@ -51,10 +51,10 @@ class IOMeshAscii : public IOMesh<Tdim> {
   std::vector<Eigen::Matrix<double, 6, 1>> read_particles_stresses(
       const std::string& particles_stresses) override;
 
-  //! Read particle scalar properties
-  //! \param[in] scalar_file file name with particle scalar properties
-  //! \retval Vector of particles scalar properties
-  std::vector<std::tuple<mpm::Index, double>> read_particles_scalar_properties(
+  //! Read scalar properties for particles or nodes
+  //! \param[in] scalar_file file name with particle or node scalar properties
+  //! \retval Vector of scalar properties for particles or nodes
+  std::vector<std::tuple<mpm::Index, double>> read_scalar_properties(
       const std::string& scalar_file) override;
 
   //! Read pressure constraints file
@@ -110,6 +110,11 @@ class IOMeshAscii : public IOMesh<Tdim> {
       read_adhesion_constraints(
           const std::string& adhesion_constraints_file) override;
 
+  //! Read levelset file
+  //! \param[in] levelset_input_file file name with levelset values
+  std::vector<std::tuple<mpm::Index, double, double, double, double>>
+      read_levelset_input(const std::string& levelset_input_file) override;
+
   //! Read traction file
   //! \param[in] forces_file file name with nodal concentrated force
   std::vector<std::tuple<mpm::Index, unsigned, double>> read_forces(
diff --git a/include/io/io_mesh_ascii.tcc b/include/io/io_mesh_ascii.tcc
index 303c632ed..4f1af17cc 100644
--- a/include/io/io_mesh_ascii.tcc
+++ b/include/io/io_mesh_ascii.tcc
@@ -254,13 +254,13 @@ std::vector<Eigen::Matrix<double, 6, 1>>
   return stresses;
 }
 
-//! Return particles scalar properties
+//! Return scalar properties for particles or nodes
 template <unsigned Tdim>
 std::vector<std::tuple<mpm::Index, double>>
-    mpm::IOMeshAscii<Tdim>::read_particles_scalar_properties(
+    mpm::IOMeshAscii<Tdim>::read_scalar_properties(
         const std::string& scalar_file) {
 
-  // Particles scalar properties
+  // Scalar properties for particles or nodes
   std::vector<std::tuple<mpm::Index, double>> scalar_properties;
 
   // input file stream
@@ -293,7 +293,7 @@ std::vector<std::tuple<mpm::Index, double>>
     }
     file.close();
   } catch (std::exception& exception) {
-    console_->error("Read particle {} #{}: {}\n", __FILE__, __LINE__,
+    console_->error("Read particle/node {} #{}: {}\n", __FILE__, __LINE__,
                     exception.what());
     file.close();
   }
@@ -641,7 +641,7 @@ std::vector<std::tuple<mpm::Index, unsigned, double>>
   return constraints;
 }
 
-//! Return friction constraints of particles
+//! Return friction constraints
 template <unsigned Tdim>
 std::vector<std::tuple<mpm::Index, unsigned, int, double>>
     mpm::IOMeshAscii<Tdim>::read_friction_constraints(
@@ -748,6 +748,61 @@ std::vector<std::tuple<mpm::Index, unsigned, int, double, double, int>>
   return constraints;
 }
 
+//! Return nodal levelset information
+template <unsigned Tdim>
+std::vector<std::tuple<mpm::Index, double, double, double, double>>
+    mpm::IOMeshAscii<Tdim>::read_levelset_input(
+        const std::string& levelset_input_file) {
+
+  // Nodal levelset information
+  std::vector<std::tuple<mpm::Index, double, double, double, double>>
+      levelset_inputs;
+  levelset_inputs.clear();
+
+  // input file stream
+  std::fstream file;
+  file.open(levelset_input_file.c_str(), std::ios::in);
+
+  try {
+    if (file.is_open() && file.good()) {
+      // Line
+      std::string line;
+      while (std::getline(file, line)) {
+        boost::algorithm::trim(line);
+        std::istringstream istream(line);
+        // ignore comment lines (# or !) or blank lines
+        if ((line.find('#') == std::string::npos) &&
+            (line.find('!') == std::string::npos) && (line != "")) {
+          // ID
+          mpm::Index id;
+          // Levelset value
+          double levelset;
+          // Friction
+          double levelset_mu;
+          // Adhesion coefficient
+          double levelset_alpha;
+          // Barrier stiffness
+          double barrier_stiffness;
+          while (istream.good()) {
+            // Read stream
+            istream >> id >> levelset >> levelset_mu >> levelset_alpha >>
+                barrier_stiffness;
+            levelset_inputs.emplace_back(std::make_tuple(
+                id, levelset, levelset_mu, levelset_alpha, barrier_stiffness));
+          }
+        }
+      }
+    } else {
+      throw std::runtime_error("File not open or not good!");
+    }
+    file.close();
+  } catch (std::exception& exception) {
+    console_->error("Read levelset inputs: {}", exception.what());
+    file.close();
+  }
+  return levelset_inputs;
+}
+
 //! Return particles force
 template <unsigned Tdim>
 std::vector<std::tuple<mpm::Index, unsigned, double>>
diff --git a/include/mesh/mesh.h b/include/mesh/mesh.h
index 59c315d24..2c350c009 100644
--- a/include/mesh/mesh.h
+++ b/include/mesh/mesh.h
@@ -164,7 +164,7 @@ class Mesh {
   //! \param[in] cells Node ids of cells
   //! \param[in] check_duplicates Parameter to check duplicates
   //! \retval status Create cells status
-  bool create_cells(mpm::Index gnid,
+  bool create_cells(mpm::Index gcid,
                     const std::shared_ptr<mpm::Element<Tdim>>& element,
                     const std::vector<std::vector<mpm::Index>>& cells,
                     bool check_duplicates = true);
@@ -263,7 +263,7 @@ class Mesh {
   //! Locate particles in a cell
   //! Iterate over all cells in a mesh to find the cell in which particles
   //! are located.
-  //! \retval particles Particles which cannot be located in the mesh
+  //! \retval particles which cannot be located in the mesh
   std::vector<std::shared_ptr<mpm::ParticleBase<Tdim>>> locate_particles_mesh();
 
   //! Iterate over particles
@@ -475,7 +475,7 @@ class Mesh {
   //! Create map of vector of particles in sets
   //! \param[in] map of particles ids in sets
   //! \param[in] check_duplicates Parameter to check duplicates
-  //! \retval status Status of  create particle sets
+  //! \retval status Status of create particle sets
   bool create_particle_sets(
       const tsl::robin_map<mpm::Index, std::vector<mpm::Index>>& particle_sets,
       bool check_duplicates);
@@ -518,11 +518,33 @@ class Mesh {
   void inject_particles(double current_time);
 
   // Create the nodal properties' map
-  void create_nodal_properties();
+  virtual void create_nodal_properties();
 
   // Initialise the nodal properties' map
   void initialise_nodal_properties();
 
+  /**
+   * \defgroup Levelset Functions
+   */
+  /**@{*/
+
+  //! Assign nodal levelset values
+  //! \ingroup Levelset
+  //! \param[in] levelset Levelset value at the particle
+  //! \param[in] levelset_mu Levelset friction coefficient
+  //! \param[in] levelset_alpha Levelset adhesion coefficient
+  //! \param[in] barrier_stiffness Barrier stiffness
+  virtual bool assign_nodal_levelset_values(
+      const std::vector<std::tuple<mpm::Index, double, double, double, double>>&
+          levelset_input_file) {
+    throw std::runtime_error(
+        "Calling the base class function (assign_nodal_levelset_values) in "
+        "Mesh:: illegal operation!");
+    return false;
+  };
+
+  /**@}*/
+
   /**
    * \defgroup MultiPhase Functions dealing with multi-phase MPM
    */
@@ -673,10 +695,11 @@ class Mesh {
                            const Json& generator, unsigned pset_id);
 
   // Locate a particle in mesh cells
+  //! \param[in] particle of interest
   bool locate_particle_cells(
       const std::shared_ptr<mpm::ParticleBase<Tdim>>& particle);
 
- private:
+ protected:
   //! mesh id
   unsigned id_{std::numeric_limits<unsigned>::max()};
   //! Isoparametric mesh
diff --git a/include/mesh/mesh_levelset.h b/include/mesh/mesh_levelset.h
new file mode 100644
index 000000000..d7fe40cde
--- /dev/null
+++ b/include/mesh/mesh_levelset.h
@@ -0,0 +1,76 @@
+#ifndef MPM_MESH_LEVELSET_H_
+#define MPM_MESH_LEVELSET_H_
+
+#include "logger.h"
+#include "mesh.h"
+
+namespace mpm {
+
+//! Levelset subclass
+//! \brief subclass that stores the information about levelset mesh
+//! \tparam Tdim Dimension
+template <unsigned Tdim>
+class MeshLevelset : public Mesh<Tdim> {
+
+ public:
+  //! Define a vector of size dimension
+  using VectorDim = Eigen::Matrix<double, Tdim, 1>;
+
+  // Construct a mesh with a global unique id
+  //! \param[in] id Global mesh id
+  //! \param[in] isoparametric Mesh is isoparametric
+  MeshLevelset(unsigned id, bool isoparametric)
+      : mpm::Mesh<Tdim>(id, isoparametric),
+        id_{id},
+        isoparametric_{isoparametric} {
+    // Check if the dimension is between 1 & 3
+    static_assert((Tdim >= 1 && Tdim <= 3), "Invalid global dimension");
+    //! Logger
+    std::string logger = "meshlevelset::" + std::to_string(id);
+    console_ = std::make_unique<spdlog::logger>(logger, mpm::stdout_sink);
+  }
+
+  //! Default destructor
+  ~MeshLevelset() = default;
+
+  //! Delete copy constructor
+  MeshLevelset(const MeshLevelset<Tdim>&) = delete;
+
+  //! Delete assignement operator
+  MeshLevelset& operator=(const MeshLevelset<Tdim>&) = delete;
+
+  //! Assign mesh levelset values to nodes
+  //! \param[in] levelset Levelset value at the particle
+  //! \param[in] levelset_mu Levelset friction coefficient
+  //! \param[in] levelset_alpha Levelset adhesion coefficient
+  //! \param[in] barrier_stiffness Barrier stiffness
+  bool assign_nodal_levelset_values(
+      const std::vector<std::tuple<mpm::Index, double, double, double, double>>&
+          levelset_input_file) override;
+
+  // Create the nodal properties' map
+  void create_nodal_properties() override;
+
+ private:
+  //! mesh id
+  unsigned id_{std::numeric_limits<unsigned>::max()};
+  //! Isoparametric mesh
+  bool isoparametric_{true};
+  //! Logger
+  std::unique_ptr<spdlog::logger> console_;
+  //! Vector of nodes
+  using mpm::Mesh<Tdim>::nodes_;
+  //! Map of nodes for fast retrieval
+  using mpm::Mesh<Tdim>::map_nodes_;
+  //! Materials
+  using mpm::Mesh<Tdim>::materials_;
+  //! Nodal property pool
+  using mpm::Mesh<Tdim>::nodal_properties_;
+
+};  // MeshLevelset class
+
+}  // namespace mpm
+
+#include "mesh_levelset.tcc"
+
+#endif  // MPM_MESH_LEVELSET_H_
\ No newline at end of file
diff --git a/include/mesh/mesh_levelset.tcc b/include/mesh/mesh_levelset.tcc
new file mode 100644
index 000000000..405674bbe
--- /dev/null
+++ b/include/mesh/mesh_levelset.tcc
@@ -0,0 +1,87 @@
+//! Assign mesh levelset values to nodes
+template <unsigned Tdim>
+bool mpm::MeshLevelset<Tdim>::assign_nodal_levelset_values(
+    const std::vector<std::tuple<mpm::Index, double, double, double, double>>&
+        levelset_input_file) {
+  bool status = true;
+  try {
+    if (!nodes_.size())
+      throw std::runtime_error(
+          "No nodes have been assigned in mesh, cannot assign levelset values");
+
+    for (const auto& levelset_inputs : levelset_input_file) {
+      // Node id
+      mpm::Index nid = std::get<0>(levelset_inputs);
+      // Levelset
+      double levelset = std::get<1>(levelset_inputs);
+      // Levelset friction coefficient
+      double levelset_mu = std::get<2>(levelset_inputs);
+      // Levelset adhesion coefficient
+      double levelset_alpha = std::get<3>(levelset_inputs);
+      // Barrier stiffness
+      double barrier_stiffness = std::get<4>(levelset_inputs);
+
+      if (map_nodes_.find(nid) != map_nodes_.end())
+        status = map_nodes_[nid]->assign_levelset(
+            levelset, levelset_mu, levelset_alpha, barrier_stiffness);
+
+      if (!status)
+        throw std::runtime_error("Cannot assign invalid nodal levelset values");
+    }
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+    status = false;
+  }
+  return status;
+}
+
+// Create the nodal properties' map
+template <unsigned Tdim>
+void mpm::MeshLevelset<Tdim>::create_nodal_properties() {
+  // Initialise the shared pointer to nodal properties
+  nodal_properties_ = std::make_shared<mpm::NodalProperties>();
+
+  // Check if nodes_ and materials_is empty (if empty throw runtime error)
+  if (nodes_.size() != 0 && materials_.size() != 0) {
+    // Compute number of rows in nodal properties for vector entities
+    const unsigned nrows = nodes_.size() * Tdim;
+    // Create pool data for each property in the nodal properties struct
+    // object. Properties must be named in the plural form
+    nodal_properties_->create_property("masses", nodes_.size(),
+                                       materials_.size());
+    nodal_properties_->create_property("momenta", nrows, materials_.size());
+    nodal_properties_->create_property("change_in_momenta", nrows,
+                                       materials_.size());
+    nodal_properties_->create_property("displacements", nrows,
+                                       materials_.size());
+    nodal_properties_->create_property("separation_vectors", nrows,
+                                       materials_.size());
+    nodal_properties_->create_property("domain_gradients", nrows,
+                                       materials_.size());
+    nodal_properties_->create_property("normal_unit_vectors", nrows,
+                                       materials_.size());
+    nodal_properties_->create_property("wave_velocities", nrows,
+                                       materials_.size());
+    nodal_properties_->create_property("density", nodes_.size(),
+                                       materials_.size());
+    // levelset properties
+    nodal_properties_->create_property("levelsets", nodes_.size(),
+                                       materials_.size());
+    nodal_properties_->create_property("levelset_mus", nodes_.size(),
+                                       materials_.size());
+    nodal_properties_->create_property("levelset_alphas", nodes_.size(),
+                                       materials_.size());
+    nodal_properties_->create_property("barrier_stiffnesses", nodes_.size(),
+                                       materials_.size());
+    nodal_properties_->create_property("levelset_mp_radii", nodes_.size(),
+                                       materials_.size());
+
+    // Iterate over all nodes to initialise the property handle in each node
+    // and assign its node id as the prop id in the nodal property data pool
+    for (auto nitr = nodes_.cbegin(); nitr != nodes_.cend(); ++nitr)
+      (*nitr)->initialise_property_handle((*nitr)->id(), nodal_properties_);
+  } else {
+    throw std::runtime_error(
+        "Number of nodes or number of materials is zero (levelset)");
+  }
+}
\ No newline at end of file
diff --git a/include/nodes/node.h b/include/nodes/node.h
index aa182e563..bc9769c33 100644
--- a/include/nodes/node.h
+++ b/include/nodes/node.h
@@ -332,11 +332,11 @@ class Node : public NodeBase<Tdim> {
    */
   /**@{*/
   //! Initialise nodal properties for implicit solver
-  //! \ingroup Impolicit
+  //! \ingroup Implicit
   void initialise_implicit() noexcept override;
 
   //! Initialise nodal forces
-  //! \ingroup Impolicit
+  //! \ingroup Implicit
   void initialise_force() noexcept override;
 
   //! Update inertia at the nodes
diff --git a/include/nodes/node_base.h b/include/nodes/node_base.h
index af4763e16..bd6fbc9c5 100644
--- a/include/nodes/node_base.h
+++ b/include/nodes/node_base.h
@@ -549,6 +549,59 @@ class NodeBase {
 
   /**@}*/
 
+  /**
+   * \defgroup Levelset Functions
+   */
+  /**@{*/
+
+  // Assign levelset values to nodes
+  //! \param[in] levelset Levelset value at the particle
+  //! \param[in] levelset_mu Levelset friction coefficient
+  //! \param[in] levelset_alpha Levelset adhesion coefficient
+  //! \param[in] barrier_stiffness Barrier stiffness
+  virtual bool assign_levelset(double levelset, double levelset_mu,
+                               double levelset_alpha,
+                               double barrier_stiffness) {
+    throw std::runtime_error(
+        "Calling the base class function (assign_levelset) in NodeBase:: "
+        "illegal operation!");
+    return false;
+  };
+
+  //! Return levelset value
+  virtual double levelset() const {
+    throw std::runtime_error(
+        "Calling the base class function (levelset) in NodeBase:: illegal "
+        "operation!");
+    return 0.;
+  }
+
+  //! Return levelset friction coefficient
+  virtual double levelset_mu() const {
+    throw std::runtime_error(
+        "Calling the base class function (levelset_mu) in NodeBase:: illegal "
+        "operation!");
+    return 0.;
+  }
+
+  //! Return levelset adhesion coefficient
+  virtual double levelset_alpha() const {
+    throw std::runtime_error(
+        "Calling the base class function (levelset_alpha) in NodeBase:: "
+        "illegal operation!");
+    return 0.;
+  }
+
+  //! Return barrier stiffness
+  virtual double barrier_stiffness() const {
+    throw std::runtime_error(
+        "Calling the base class function (barrier_stiffness) in NodeBase:: "
+        "illegal operation!");
+    return 0.;
+  }
+
+  /**@}*/
+
 };  // NodeBase class
 }  // namespace mpm
 
diff --git a/include/nodes/node_levelset.h b/include/nodes/node_levelset.h
new file mode 100644
index 000000000..13dae9bee
--- /dev/null
+++ b/include/nodes/node_levelset.h
@@ -0,0 +1,76 @@
+#ifndef MPM_NODE_LEVELSET_H_
+#define MPM_NODE_LEVELSET_H_
+
+#include "node.h"
+#include "node_base.h"
+
+namespace mpm {
+
+//! Levelset subclass
+//! \brief subclass that stores the information about levelset node
+//! \tparam Tdim Dimension
+template <unsigned Tdim, unsigned Tdof, unsigned Tnphases>
+class NodeLevelset : public Node<Tdim, Tdof, Tnphases> {
+
+ public:
+  //! Define a vector of size dimension
+  using VectorDim = Eigen::Matrix<double, Tdim, 1>;
+
+  //! Constructor with id, coordinates and dof
+  //! \param[in] id Node id
+  //! \param[in] coord coordinates of the node
+  NodeLevelset(Index id, const VectorDim& coord)
+      : mpm::Node<Tdim, Tdof, Tnphases>(id, coord) {
+    console_ =
+        std::make_unique<spdlog::logger>("NodeLevelset", mpm::stdout_sink);
+  };
+
+  //! Virtual destructor
+  ~NodeLevelset() override{};
+
+  //! Delete copy constructor
+  NodeLevelset(const NodeLevelset<Tdim, Tdof, Tnphases>&) = delete;
+
+  //! Delete assignement operator
+  NodeLevelset& operator=(const NodeLevelset<Tdim, Tdof, Tnphases>&) = delete;
+
+  // Assign levelset values to nodes
+  //! \param[in] levelset Levelset value at the particle
+  //! \param[in] levelset_mu Levelset friction coefficient
+  //! \param[in] levelset_alpha Levelset adhesion coefficient
+  //! \param[in] barrier_stiffness Barrier stiffness
+  bool assign_levelset(double levelset, double levelset_mu,
+                       double levelset_alpha,
+                       double barrier_stiffness) override;
+
+  //! Return levelset value
+  inline double levelset() const override { return levelset_; }
+
+  //! Return levelset friction coefficient
+  inline double levelset_mu() const override { return levelset_mu_; }
+
+  //! Return levelset adhesion coefficient
+  inline double levelset_alpha() const override { return levelset_alpha_; }
+
+  //! Return barrier stiffness
+  inline double barrier_stiffness() const override {
+    return barrier_stiffness_;
+  }
+
+ private:
+  //! Logger
+  std::unique_ptr<spdlog::logger> console_;
+  //! Levelset value
+  double levelset_;
+  //! Levelset friction coefficient
+  double levelset_mu_;
+  //! Levelset adhesion coefficient
+  double levelset_alpha_;
+  //! Barrier stiffness
+  double barrier_stiffness_;
+};  // NodeLevelset class
+}  // namespace mpm
+
+#include "node_levelset.tcc"
+
+#endif  // MPM_NODE_LEVELSET_H_
\ No newline at end of file
diff --git a/include/nodes/node_levelset.tcc b/include/nodes/node_levelset.tcc
new file mode 100644
index 000000000..90923b396
--- /dev/null
+++ b/include/nodes/node_levelset.tcc
@@ -0,0 +1,24 @@
+// Assign levelset values to the nodes
+template <unsigned Tdim, unsigned Tdof, unsigned Tnphases>
+bool mpm::NodeLevelset<Tdim, Tdof, Tnphases>::assign_levelset(
+    double levelset, double levelset_mu, double levelset_alpha,
+    double barrier_stiffness) {
+  bool status = true;
+  try {
+    if ((levelset_mu < 0.) || (levelset_alpha < 0.))
+      throw std::runtime_error("Levelset mu and alpha cannot be negative");
+    if ((barrier_stiffness <= 0.))
+      throw std::runtime_error("Barrier stiffness must be greater than zero");
+
+    // Set variables
+    this->levelset_ = levelset;
+    this->levelset_mu_ = levelset_mu;
+    this->levelset_alpha_ = levelset_alpha;
+    this->barrier_stiffness_ = barrier_stiffness;
+
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+    status = false;
+  }
+  return status;
+}
\ No newline at end of file
diff --git a/include/particles/anti_locking/particle_bbar.h b/include/particles/anti_locking/particle_bbar.h
index a2eb5cd1b..ec7877e2d 100644
--- a/include/particles/anti_locking/particle_bbar.h
+++ b/include/particles/anti_locking/particle_bbar.h
@@ -61,7 +61,6 @@ class ParticleBbar : public mpm::Particle<Tdim> {
 
  protected:
   //! Compute strain rate
-  //! \ingroup Implicit
   //! \param[in] dn_dx The spatial gradient of shape function
   //! \param[in] phase Index to indicate phase
   //! \retval strain rate at particle inside a cell
diff --git a/include/particles/anti_locking/particle_levelset_bbar.h b/include/particles/anti_locking/particle_levelset_bbar.h
new file mode 100644
index 000000000..f3af3bd24
--- /dev/null
+++ b/include/particles/anti_locking/particle_levelset_bbar.h
@@ -0,0 +1,92 @@
+#ifndef MPM_PARTICLE_LEVELSET_BBAR_H_
+#define MPM_PARTICLE_LEVELSET_BBAR_H_
+
+#include <array>
+#include <limits>
+#include <memory>
+#include <string>
+#include <vector>
+
+#include "logger.h"
+#include "particle.h"
+#include "particle_levelset.h"
+
+namespace mpm {
+
+//! ParticleLevelsetBbar class
+//! \brief Class that introduce levelset B-Bar formulation to the standard
+//! single-phase particle \tparam Tdim Dimension
+template <unsigned Tdim>
+class ParticleLevelsetBbar : public mpm::ParticleLevelset<Tdim> {
+ public:
+  //! Define a vector of size dimension
+  using VectorDim = Eigen::Matrix<double, Tdim, 1>;
+
+  //! Construct a levelset particle using B-bar method with id and coordinates
+  //! \param[in] id Particle id
+  //! \param[in] coord Coordinates of the particles
+  ParticleLevelsetBbar(Index id, const VectorDim& coord);
+
+  //! Construct a levelset particle using B-bar method with id, coordinates and
+  //! status \param[in] id Particle id \param[in] coord coordinates of the
+  //! particle \param[in] status Particle status (active / inactive)
+  ParticleLevelsetBbar(Index id, const VectorDim& coord, bool status);
+
+  //! Destructor
+  ~ParticleLevelsetBbar() override{};
+
+  //! Delete copy constructor
+  ParticleLevelsetBbar(const ParticleLevelsetBbar<Tdim>&) = delete;
+
+  //! Delete assignment operator
+  ParticleLevelsetBbar& operator=(const ParticleLevelsetBbar<Tdim>&) = delete;
+
+  //! Map internal force
+  inline void map_internal_force() noexcept override;
+
+  //! Type of particle
+  std::string type() const override {
+    return (Tdim == 2) ? "P2DLSBBAR" : "P3DLSBBAR";
+  }
+
+ protected:
+  //! Compute strain rate
+  //! \param[in] dn_dx The spatial gradient of shape function
+  //! \param[in] phase Index to indicate phase
+  //! \retval strain rate at particle inside a cell
+  inline Eigen::Matrix<double, 6, 1> compute_strain_rate(
+      const Eigen::MatrixXd& dn_dx, unsigned phase) noexcept override;
+
+ protected:
+  //! Nodes
+  using ParticleBase<Tdim>::nodes_;
+  //! Volume
+  using Particle<Tdim>::volume_;
+  //! Stresses
+  using Particle<Tdim>::stress_;
+  //! Strains
+  using Particle<Tdim>::strain_;
+  //! dvolumetric strain
+  using Particle<Tdim>::dvolumetric_strain_;
+  //! Strain rate
+  using Particle<Tdim>::strain_rate_;
+  //! dstrains
+  using Particle<Tdim>::dstrain_;
+  //! Velocity
+  using Particle<Tdim>::velocity_;
+  //! Displacement
+  using Particle<Tdim>::displacement_;
+  //! dN/dX
+  using Particle<Tdim>::dn_dx_;
+  //! dN/dX at cell centroid
+  using Particle<Tdim>::dn_dx_centroid_;
+
+  //! Logger
+  std::unique_ptr<spdlog::logger> console_;
+
+};  // ParticleLevelsetBbar class
+}  // namespace mpm
+
+#include "particle_levelset_bbar.tcc"
+
+#endif  // MPM_PARTICLE_LEVELSET_BBAR_H__
diff --git a/include/particles/anti_locking/particle_levelset_bbar.tcc b/include/particles/anti_locking/particle_levelset_bbar.tcc
new file mode 100644
index 000000000..da021d57b
--- /dev/null
+++ b/include/particles/anti_locking/particle_levelset_bbar.tcc
@@ -0,0 +1,160 @@
+//! Construct a particle with id and coordinates
+template <unsigned Tdim>
+mpm::ParticleLevelsetBbar<Tdim>::ParticleLevelsetBbar(Index id,
+                                                      const VectorDim& coord)
+    : mpm::ParticleLevelset<Tdim>(id, coord) {
+  // Logger
+  std::string logger =
+      "particle_bbar" + std::to_string(Tdim) + "d::" + std::to_string(id);
+  console_ = std::make_unique<spdlog::logger>(logger, mpm::stdout_sink);
+}
+
+//! Construct a particle with id, coordinates and status
+template <unsigned Tdim>
+mpm::ParticleLevelsetBbar<Tdim>::ParticleLevelsetBbar(Index id,
+                                                      const VectorDim& coord,
+                                                      bool status)
+    : mpm::ParticleLevelset<Tdim>(id, coord, status) {
+  //! Logger
+  std::string logger = "particle_levelset_bbar" + std::to_string(Tdim) +
+                       "d::" + std::to_string(id);
+  console_ = std::make_unique<spdlog::logger>(logger, mpm::stdout_sink);
+}
+
+// Compute strain rate of the particle
+template <>
+inline Eigen::Matrix<double, 6, 1>
+    mpm::ParticleLevelsetBbar<1>::compute_strain_rate(
+        const Eigen::MatrixXd& dn_dx, unsigned phase) noexcept {
+  // Define strain rate
+  Eigen::Matrix<double, 6, 1> strain_rate = Eigen::Matrix<double, 6, 1>::Zero();
+
+  for (unsigned i = 0; i < this->nodes_.size(); ++i) {
+    Eigen::Matrix<double, 1, 1> vel = nodes_[i]->velocity(phase);
+    strain_rate[0] += dn_dx(i, 0) * vel[0];
+  }
+
+  if (std::fabs(strain_rate(0)) < 1.E-15) strain_rate[0] = 0.;
+  return strain_rate;
+}
+
+// Compute strain rate of the particle
+template <>
+inline Eigen::Matrix<double, 6, 1>
+    mpm::ParticleLevelsetBbar<2>::compute_strain_rate(
+        const Eigen::MatrixXd& dn_dx, unsigned phase) noexcept {
+  // Define strain rate
+  Eigen::Matrix<double, 6, 1> strain_rate = Eigen::Matrix<double, 6, 1>::Zero();
+
+  for (unsigned i = 0; i < this->nodes_.size(); ++i) {
+    Eigen::Matrix<double, 2, 1> vel = nodes_[i]->velocity(phase);
+    // clang-format off
+    strain_rate[0] += (dn_dx(i, 0) + (dn_dx_centroid_(i, 0) - dn_dx(i, 0)) / 2.) * vel[0] +
+                      (dn_dx_centroid_(i, 1) - dn_dx(i, 1)) / 2. * vel[1];
+    strain_rate[1] += (dn_dx_centroid_(i, 0) - dn_dx(i, 0)) / 2. * vel[0] +
+                      (dn_dx(i, 1) + (dn_dx_centroid_(i, 1) - dn_dx(i, 1)) / 2.) * vel[1];
+    strain_rate[3] += dn_dx(i, 1) * vel[0] + dn_dx(i, 0) * vel[1];
+    // clang-format on
+  }
+
+  if (std::fabs(strain_rate[0]) < 1.E-15) strain_rate[0] = 0.;
+  if (std::fabs(strain_rate[1]) < 1.E-15) strain_rate[1] = 0.;
+  if (std::fabs(strain_rate[3]) < 1.E-15) strain_rate[3] = 0.;
+  return strain_rate;
+}
+
+// Compute strain rate of the particle
+template <>
+inline Eigen::Matrix<double, 6, 1>
+    mpm::ParticleLevelsetBbar<3>::compute_strain_rate(
+        const Eigen::MatrixXd& dn_dx, unsigned phase) noexcept {
+  // Define strain rate
+  Eigen::Matrix<double, 6, 1> strain_rate = Eigen::Matrix<double, 6, 1>::Zero();
+
+  for (unsigned i = 0; i < this->nodes_.size(); ++i) {
+    Eigen::Matrix<double, 3, 1> vel = nodes_[i]->velocity(phase);
+    // clang-format off
+    strain_rate[0] += (dn_dx(i, 0) + (dn_dx_centroid_(i, 0) - dn_dx(i, 0)) / 3.) * vel[0] +
+                      (dn_dx_centroid_(i, 1) - dn_dx(i, 1)) / 3. * vel[1] +
+                      (dn_dx_centroid_(i, 2) - dn_dx(i, 2)) / 3. * vel[2];
+    strain_rate[1] += (dn_dx_centroid_(i, 0) - dn_dx(i, 0)) / 3. * vel[0] +
+                      (dn_dx(i, 1) + (dn_dx_centroid_(i, 1) - dn_dx(i, 1)) / 3.) * vel[1] +
+                      (dn_dx_centroid_(i, 2) - dn_dx(i, 2)) / 3. * vel[2];
+    strain_rate[2] += (dn_dx_centroid_(i, 0) - dn_dx(i, 0)) / 3. * vel[0] +
+                      (dn_dx_centroid_(i, 1) - dn_dx(i, 1)) / 3. * vel[1] +
+                      (dn_dx(i, 2) + (dn_dx_centroid_(i, 2) - dn_dx(i, 2)) / 3.) * vel[2];
+    strain_rate[3] += dn_dx(i, 1) * vel[0] + dn_dx(i, 0) * vel[1];
+    strain_rate[4] += dn_dx(i, 2) * vel[1] + dn_dx(i, 1) * vel[2];
+    strain_rate[5] += dn_dx(i, 2) * vel[0] + dn_dx(i, 0) * vel[2];
+    // clang-format on
+  }
+
+  for (unsigned i = 0; i < strain_rate.size(); ++i)
+    if (std::fabs(strain_rate[i]) < 1.E-15) strain_rate[i] = 0.;
+  return strain_rate;
+}
+
+//! Map internal force
+template <>
+inline void mpm::ParticleLevelsetBbar<1>::map_internal_force() noexcept {
+  // Compute nodal internal forces
+  for (unsigned i = 0; i < nodes_.size(); ++i) {
+    // Compute force: -pstress * volume
+    Eigen::Matrix<double, 1, 1> force;
+    force[0] = -1. * dn_dx_(i, 0) * volume_ * stress_[0];
+
+    nodes_[i]->update_internal_force(true, mpm::ParticlePhase::Solid, force);
+  }
+}
+
+//! Map internal force
+template <>
+inline void mpm::ParticleLevelsetBbar<2>::map_internal_force() noexcept {
+  // Compute nodal internal forces
+  for (unsigned i = 0; i < nodes_.size(); ++i) {
+    // Compute force: -pstress * volume
+    Eigen::Matrix<double, 2, 1> force;
+    // clang-format off
+    force[0] = (dn_dx_(i, 0) + (dn_dx_centroid_(i, 0) - dn_dx_(i, 0)) / 2.) * stress_[0] +
+               (dn_dx_centroid_(i, 0) - dn_dx_(i, 0)) / 2. * stress_[1] +
+               dn_dx_(i, 1) * stress_[3];
+    force[1] = (dn_dx_centroid_(i, 1) - dn_dx_(i, 1)) / 2. * stress_[0] +
+               (dn_dx_(i, 1) + (dn_dx_centroid_(i, 1) - dn_dx_(i, 1)) / 2.) * stress_[1] +
+               dn_dx_(i, 0) * stress_[3];
+    // clang-format on
+
+    force *= -1. * this->volume_;
+
+    nodes_[i]->update_internal_force(true, mpm::ParticlePhase::Solid, force);
+  }
+}
+
+//! Map internal force
+template <>
+inline void mpm::ParticleLevelsetBbar<3>::map_internal_force() noexcept {
+  // Compute nodal internal forces
+  for (unsigned i = 0; i < nodes_.size(); ++i) {
+    // Compute force: -pstress * volume
+    Eigen::Matrix<double, 3, 1> force;
+    // clang-format off
+    force[0] = (dn_dx_(i, 0) + (dn_dx_centroid_(i, 0) - dn_dx_(i, 0)) / 3.) * stress_[0] +
+               (dn_dx_centroid_(i, 0) - dn_dx_(i, 0)) / 3. * stress_[1] +
+               (dn_dx_centroid_(i, 0) - dn_dx_(i, 0)) / 3. * stress_[2] +
+               dn_dx_(i, 1) * stress_[3] + dn_dx_(i, 2) * stress_[5];
+
+    force[1] = (dn_dx_centroid_(i, 1) - dn_dx_(i, 1)) / 3. * stress_[0] +
+               (dn_dx_(i, 1) + (dn_dx_centroid_(i, 1) - dn_dx_(i, 1)) / 3.) * stress_[1] +
+               (dn_dx_centroid_(i, 1) - dn_dx_(i, 1)) / 3. * stress_[2] +
+               dn_dx_(i, 0) * stress_[3] + dn_dx_(i, 2) * stress_[4];
+
+    force[2] = (dn_dx_centroid_(i, 2) - dn_dx_(i, 2)) / 3. * stress_[0] +
+               (dn_dx_centroid_(i, 2) - dn_dx_(i, 2)) / 3. * stress_[1] +
+               (dn_dx_(i, 2) + (dn_dx_centroid_(i, 2) - dn_dx_(i, 2)) / 3.) * stress_[2] +
+               dn_dx_(i, 1) * stress_[4] + dn_dx_(i, 0) * stress_[5];
+    // clang-format on
+
+    force *= -1. * this->volume_;
+
+    nodes_[i]->update_internal_force(true, mpm::ParticlePhase::Solid, force);
+  }
+}
\ No newline at end of file
diff --git a/include/particles/particle.h b/include/particles/particle.h
index 53c5920e2..4849dd8f5 100644
--- a/include/particles/particle.h
+++ b/include/particles/particle.h
@@ -630,7 +630,7 @@ class Particle : public ParticleBase<Tdim> {
   //! Volume
   double volume_{0.};
   //! Size of particle
-  Eigen::Matrix<double, Tdim, 1> size_;
+  double size_;
   //! Size of particle in natural coordinates
   Eigen::Matrix<double, Tdim, 1> natural_size_;
   //! Stresses
diff --git a/include/particles/particle.tcc b/include/particles/particle.tcc
index cabb0aed2..bd1b8d9ea 100644
--- a/include/particles/particle.tcc
+++ b/include/particles/particle.tcc
@@ -40,18 +40,17 @@ bool mpm::Particle<Tdim>::initialise_particle(PODParticle& particle) {
   this->mass_ = particle.mass;
   // Volume
   this->volume_ = particle.volume;
-  // Compute size of particle in each direction
-  const double length = std::pow(this->volume_, static_cast<double>(1. / Tdim));
-  // Set particle size as length on each side
-  this->size_.fill(length);
 
   // Mass Density
   this->mass_density_ = particle.mass / particle.volume;
   // Set local size of particle
-  Eigen::Vector3d psize;
-  psize << particle.nsize_x, particle.nsize_y, particle.nsize_z;
+  this->size_ = particle.size;
+  Eigen::Vector3d pnsize;
+  pnsize << particle.nsize_x, particle.nsize_y, particle.nsize_z;
   // Initialise particle size
-  for (unsigned i = 0; i < Tdim; ++i) this->natural_size_(i) = psize(i);
+  for (unsigned i = 0; i < Tdim; ++i) {
+    this->natural_size_(i) = pnsize(i);
+  }
 
   // Coordinates
   Eigen::Vector3d coordinates;
@@ -200,8 +199,9 @@ std::shared_ptr<void> mpm::Particle<Tdim>::pod() const {
   // Particle local size
   Eigen::Vector3d nsize;
   nsize.setZero();
-  Eigen::VectorXd size = this->natural_size();
-  for (unsigned j = 0; j < Tdim; ++j) nsize[j] = size[j];
+  for (unsigned j = 0; j < Tdim; ++j) {
+    nsize[j] = natural_size_[j];
+  }
 
   Eigen::Matrix<double, 6, 1> stress = this->stress_;
 
@@ -234,6 +234,8 @@ std::shared_ptr<void> mpm::Particle<Tdim>::pod() const {
   particle_data->displacement_y = displacement[1];
   particle_data->displacement_z = displacement[2];
 
+  particle_data->size = this->size_;
+
   particle_data->nsize_x = nsize[0];
   particle_data->nsize_y = nsize[1];
   particle_data->nsize_z = nsize[2];
@@ -313,7 +315,7 @@ void mpm::Particle<Tdim>::initialise() {
   mass_ = 0.;
   natural_size_.setZero();
   set_traction_ = false;
-  size_.setZero();
+  size_ = 0.0;
   strain_rate_.setZero();
   strain_.setZero();
   previous_stress_.setZero();
@@ -561,11 +563,8 @@ bool mpm::Particle<Tdim>::assign_volume(double volume) {
       throw std::runtime_error("Particle volume cannot be negative");
 
     this->volume_ = volume;
-    // Compute size of particle in each direction
-    const double length =
-        std::pow(this->volume_, static_cast<double>(1. / Tdim));
-    // Set particle size as length on each side
-    this->size_.fill(length);
+    // Calculate initial particle size from volume
+    this->size_ = std::pow(this->volume_, static_cast<double>(1. / Tdim));
 
     if (cell_ != nullptr) {
       // Get element ptr of a cell
@@ -967,8 +966,11 @@ bool mpm::Particle<Tdim>::assign_traction(unsigned direction, double traction) {
       throw std::runtime_error(
           "Particle traction property: volume / direction is invalid");
     }
+    // Updated size
+    const double updated_size =
+        std::pow(this->volume_, static_cast<double>(1. / Tdim));
     // Assign traction
-    traction_(direction) = traction * this->volume_ / this->size_(direction);
+    traction_(direction) = traction * this->volume_ / updated_size;
     status = true;
     this->set_traction_ = true;
   } catch (std::exception& exception) {
@@ -1285,8 +1287,8 @@ int mpm::Particle<Tdim>::compute_pack_size() const {
   MPI_Pack_size(3 * 1, MPI_DOUBLE, MPI_COMM_WORLD, &partial_size);
   total_size += partial_size;
 
-  // Coordinates, displacement, natural size, velocity, acceleration
-  MPI_Pack_size(5 * Tdim, MPI_DOUBLE, MPI_COMM_WORLD, &partial_size);
+  // Coordinates, displacement, size, natural size, velocity, acceleration
+  MPI_Pack_size(6 * Tdim, MPI_DOUBLE, MPI_COMM_WORLD, &partial_size);
   total_size += partial_size;
   // Stress & strain
   MPI_Pack_size(6 * 2, MPI_DOUBLE, MPI_COMM_WORLD, &partial_size);
@@ -1376,6 +1378,9 @@ std::vector<uint8_t> mpm::Particle<Tdim>::serialize() {
   // Displacement
   MPI_Pack(displacement_.data(), Tdim, MPI_DOUBLE, data_ptr, data.size(),
            &position, MPI_COMM_WORLD);
+  // Size
+  MPI_Pack(&size_, 1, MPI_DOUBLE, data_ptr, data.size(), &position,
+           MPI_COMM_WORLD);
   // Natural size
   MPI_Pack(natural_size_.data(), Tdim, MPI_DOUBLE, data_ptr, data.size(),
            &position, MPI_COMM_WORLD);
@@ -1480,6 +1485,9 @@ void mpm::Particle<Tdim>::deserialize(
   // Displacement
   MPI_Unpack(data_ptr, data.size(), &position, displacement_.data(), Tdim,
              MPI_DOUBLE, MPI_COMM_WORLD);
+  // Size
+  MPI_Unpack(data_ptr, data.size(), &position, &size_, 1, MPI_DOUBLE,
+             MPI_COMM_WORLD);
   // Natural size
   MPI_Unpack(data_ptr, data.size(), &position, natural_size_.data(), Tdim,
              MPI_DOUBLE, MPI_COMM_WORLD);
diff --git a/include/particles/particle_base.h b/include/particles/particle_base.h
index c941a5578..1f574a7b9 100644
--- a/include/particles/particle_base.h
+++ b/include/particles/particle_base.h
@@ -456,6 +456,31 @@ class ParticleBase {
   //! \ingroup AdvancedMapping
   virtual Eigen::MatrixXd mapping_matrix() const = 0;
 
+  //! Levelset functions--------------------------------------------------------
+  //! Update contact reaction force due to levelset interface
+  //! \param[in] dt Analysis time step
+  virtual void levelset_contact_force(double dt) {
+    throw std::runtime_error(
+        "Calling the base class function (levelset_contact_force) "
+        "in ParticleBase:: illegal operation!");
+  };
+
+  //! Return levelset value
+  virtual double levelset() const {
+    throw std::runtime_error(
+        "Calling the base class function (levelset) in "
+        "ParticleBase:: illegal operation!");
+    return 0;
+  }
+
+  //! Return levelset contact reaction force
+  virtual VectorDim reaction_force() const {
+    throw std::runtime_error(
+        "Calling the base class function (reaction_force) in "
+        "ParticleBase:: illegal operation!");
+    return VectorDim::Zero();
+  }
+
   //! Navier-Stokes functions----------------------------------
   //! Assigning beta parameter to particle
   //! \param[in] parameter parameter determining type of projection
diff --git a/include/particles/particle_levelset.h b/include/particles/particle_levelset.h
new file mode 100644
index 000000000..e43b0c76e
--- /dev/null
+++ b/include/particles/particle_levelset.h
@@ -0,0 +1,135 @@
+#ifndef MPM_PARTICLE_LEVELSET_H_
+#define MPM_PARTICLE_LEVELSET_H_
+
+#include "logger.h"
+#include "math_utility.h"
+#include "particle.h"
+
+#include <cmath>
+
+namespace mpm {
+
+//! Levelset subclass
+//! \brief subclass that stores the information about levelset particle
+//! \tparam Tdim Dimension
+template <unsigned Tdim>
+class ParticleLevelset : public Particle<Tdim> {
+
+ public:
+  //! Define a vector of size dimension
+  using VectorDim = Eigen::Matrix<double, Tdim, 1>;
+
+  //! Setter to initialize levelset static properties
+  //! \param[in] levelset_damping Levelset damping factor
+  //! \param[in] levelset_pic Particle in cell method bool for contact velocity
+  //! \param[in] levelset_violation_corrector Violation correction factor
+  static void set_levelset_properties(double levelset_damping,
+                                      bool levelset_pic,
+                                      double levelset_violation_corrector) {
+    levelset_damping_ = levelset_damping;
+    levelset_pic_ = levelset_pic;
+    levelset_violation_corrector_ = levelset_violation_corrector;
+  }
+
+  //! Construct a levelset particle with id and coordinates
+  //! \param[in] id Particle id
+  //! \param[in] coord Coordinates of the particles
+  ParticleLevelset(Index id, const VectorDim& coord);
+
+  //! Construct a levelset particle with id, coordinates and status
+  //! \param[in] id Particle id
+  //! \param[in] coord coordinates of the particle
+  //! \param[in] status Particle status (active / inactive)
+  ParticleLevelset(Index id, const VectorDim& coord, bool status);
+
+  //! Destructor
+  ~ParticleLevelset() override{};
+
+  //! Delete copy constructor
+  ParticleLevelset(const ParticleLevelset<Tdim>&) = delete;
+
+  //! Delete assignment operator
+  ParticleLevelset& operator=(const ParticleLevelset<Tdim>&) = delete;
+
+  //! Type of particle
+  std::string type() const override { return (Tdim == 2) ? "P2DLS" : "P3DLS"; }
+
+  //! Initialise particle levelset
+  void initialise() override;
+
+  //! Update contact reaction force due to levelset interface
+  //! \param[in] dt Analysis time step
+  void levelset_contact_force(double dt) override;
+
+  //! Return levelset value
+  double levelset() const { return levelset_; }
+
+  //! Return levelset contact reaction force
+  VectorDim reaction_force() const { return reaction_force_; }
+
+ protected:
+  //! Map levelset to particle
+  void map_levelset_to_particle() noexcept;
+
+  //! Check if particle in contact with levelset interface
+  //! \param[in] init_radius Particle initial radius
+  bool is_levelset_contact(double init_radius);
+
+  //! Compute levelset contact reaction force at particle
+  //! \param[in] dt Analysis time step
+  //! \param[in] init_radius Particle initial radius
+  void compute_particle_contact_force(double dt, double init_radius) noexcept;
+
+  //! Map levelset contact reaction force to nodes
+  void map_contact_force_to_nodes() noexcept;
+
+ protected:
+  //! Static levelset variables
+  static double levelset_damping_;
+  static bool levelset_pic_;
+  static double levelset_violation_corrector_;
+
+  //! Logger
+  std::unique_ptr<spdlog::logger> console_;
+  //! levelset value
+  double levelset_{0.};
+  //! levelset friction coefficient
+  double levelset_mu_{0.};
+  //! levelset adhesion coefficient
+  double levelset_alpha_{0.};
+  //! barrier stiffness
+  double barrier_stiffness_{0.};
+  //! levelset gradient
+  VectorDim levelset_gradient_{VectorDim::Zero()};
+  //! contact velocity
+  VectorDim contact_vel_{VectorDim::Zero()};
+  //! levelset normal
+  VectorDim levelset_normal_{VectorDim::Zero()};
+  //! levelset tangent
+  VectorDim levelset_tangent_{VectorDim::Zero()};
+  //! reaction force
+  VectorDim reaction_force_{VectorDim::Zero()};
+  //! Nodes
+  using Particle<Tdim>::nodes_;
+  //! Cell
+  using Particle<Tdim>::cell_;
+  //! Shape functions
+  using Particle<Tdim>::shapefn_;
+  //! dN/dX
+  using Particle<Tdim>::dn_dx_;
+  //! Velocity
+  using Particle<Tdim>::velocity_;
+  //! Volume
+  using Particle<Tdim>::volume_;
+  //! Mass
+  using Particle<Tdim>::mass_;
+  //! Size of particle
+  using Particle<Tdim>::size_;
+  //! particleBase id
+  using Particle<Tdim>::id_;
+};  // Particle_Levelset class
+}  // namespace mpm
+
+#include "particle_levelset.tcc"
+
+#endif  // MPM_PARTICLE_LEVELSET_H__
\ No newline at end of file
diff --git a/include/particles/particle_levelset.tcc b/include/particles/particle_levelset.tcc
new file mode 100644
index 000000000..8d1fb3a33
--- /dev/null
+++ b/include/particles/particle_levelset.tcc
@@ -0,0 +1,190 @@
+//! Static levelset variables
+template <unsigned Tdim>
+double mpm::ParticleLevelset<Tdim>::levelset_damping_;
+template <unsigned Tdim>
+bool mpm::ParticleLevelset<Tdim>::levelset_pic_;
+template <unsigned Tdim>
+double mpm::ParticleLevelset<Tdim>::levelset_violation_corrector_;
+
+//! Construct a particle with id and coordinates
+template <unsigned Tdim>
+mpm::ParticleLevelset<Tdim>::ParticleLevelset(Index id, const VectorDim& coord)
+    : mpm::Particle<Tdim>(id, coord) {
+  this->initialise();
+  // Clear cell ptr
+  cell_ = nullptr;
+  // Nodes
+  nodes_.clear();
+  // Set material containers
+  this->initialise_material(1);
+  // Logger
+  std::string logger =
+      "particlelevelset" + std::to_string(Tdim) + "d::" + std::to_string(id);
+  console_ = std::make_unique<spdlog::logger>(logger, mpm::stdout_sink);
+}
+
+//! Construct a particle with id, coordinates and status
+template <unsigned Tdim>
+mpm::ParticleLevelset<Tdim>::ParticleLevelset(Index id, const VectorDim& coord,
+                                              bool status)
+    : mpm::Particle<Tdim>(id, coord, status) {
+  this->initialise();
+  // Clear cell ptr
+  cell_ = nullptr;
+  // Nodes
+  nodes_.clear();
+  // Set material containers
+  this->initialise_material(1);
+  //! Logger
+  std::string logger =
+      "particlelevelset" + std::to_string(Tdim) + "d::" + std::to_string(id);
+  console_ = std::make_unique<spdlog::logger>(logger, mpm::stdout_sink);
+}
+
+// Initialise particle levelset properties
+template <unsigned Tdim>
+void mpm::ParticleLevelset<Tdim>::initialise() {
+  mpm::Particle<Tdim>::initialise();
+
+  // Initialize scalar and vector data properties
+  this->scalar_properties_["levelset"] = [this]() { return levelset(); };
+  this->vector_properties_["levelset_reaction_forces"] = [this]() {
+    return reaction_force();
+  };
+}
+
+//! Update contact reaction force due to levelset interface
+template <unsigned Tdim>
+void mpm::ParticleLevelset<Tdim>::levelset_contact_force(double dt) {
+
+  // Calculate initial radius from initial size
+  double init_vol = std::pow(this->size_, Tdim);
+  double init_radius = 0.0;  // constant radial influence
+  if (Tdim == 2) init_radius = std::sqrt(init_vol / M_PI);  // unit cylinder
+  if (Tdim == 3) init_radius = std::cbrt(init_vol * 0.75 / M_PI);  // sphere
+
+  // Map levelset to particle
+  map_levelset_to_particle();
+
+  // Check if particle in contact with levelset interface
+  if (is_levelset_contact(init_radius)) {
+
+    // Compute levelset contact reaction force at particle
+    compute_particle_contact_force(dt, init_radius);
+
+    // Map levelset contact reaction force to nodes
+    map_contact_force_to_nodes();
+  }
+}
+
+//! Map levelset to particle
+template <unsigned Tdim>
+void mpm::ParticleLevelset<Tdim>::map_levelset_to_particle() noexcept {
+  // Reset mapped levelset values at particle, per step
+  levelset_ = 0.;
+  levelset_mu_ = 0.;
+  levelset_alpha_ = 0.;
+  barrier_stiffness_ = 0.;
+  levelset_gradient_ = VectorDim::Zero();
+  contact_vel_ = VectorDim::Zero();
+
+  // Reset levelset vtk data properties
+  reaction_force_ = VectorDim::Zero();
+
+  // Map levelset to particle
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    levelset_ += shapefn_[i] * nodes_[i]->levelset();
+  }
+}
+
+//! Check if particle in contact with levelset interface
+template <unsigned Tdim>
+bool mpm::ParticleLevelset<Tdim>::is_levelset_contact(double init_radius) {
+  // Check particle levelset minimum value
+  if (levelset_ < std::numeric_limits<double>::epsilon()) {
+    console_->warn("Levelset particle {} violates interface", id_);
+    levelset_ = levelset_violation_corrector_ * init_radius;
+  }
+
+  if ((levelset_ < init_radius) && (levelset_ > 0.))
+    return true;
+  else
+    return false;
+}
+
+//! Compute levelset contact reaction force at particle
+template <unsigned Tdim>
+void mpm::ParticleLevelset<Tdim>::compute_particle_contact_force(
+    double dt, double init_radius) noexcept {
+
+  // Global contact velocity update scheme
+  if (!levelset_pic_) contact_vel_ = velocity_;
+
+  // Map other levelset values to particle
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    levelset_gradient_ += dn_dx_.row(i).transpose() * nodes_[i]->levelset();
+    levelset_mu_ += shapefn_[i] * nodes_[i]->levelset_mu();
+    levelset_alpha_ += shapefn_[i] * nodes_[i]->levelset_alpha();
+    barrier_stiffness_ += shapefn_[i] * nodes_[i]->barrier_stiffness();
+
+    // PIC contact velocity update scheme (map contact velocity from nodes)
+    if (levelset_pic_)
+      contact_vel_ +=
+          shapefn_[i] * nodes_[i]->velocity(mpm::ParticlePhase::Solid);
+  }
+
+  // Compute normals
+  levelset_normal_ = levelset_gradient_.normalized();
+
+  // Calculate normal reaction force magnitude
+  double reaction_normal_mag =
+      barrier_stiffness_ * (levelset_ - init_radius) *
+      (2. * log(levelset_ / init_radius) - (init_radius / levelset_) + 1.);
+
+  // Calculate normal reaction force vector
+  VectorDim reaction_force_normal = reaction_normal_mag * levelset_normal_;
+
+  // Calculate levelset tangential unit vector
+  double vel_n = contact_vel_.dot(levelset_normal_);
+  VectorDim tangent_calc = contact_vel_ - (vel_n * levelset_normal_);
+  if (tangent_calc.norm() > std::numeric_limits<double>::epsilon())
+    levelset_tangent_ = tangent_calc.normalized();
+
+  // Calculate friction tangential reaction force magnitude
+  double tangent_friction = levelset_mu_ * reaction_normal_mag;
+
+  // Calculate adhesion tangential reaction force magnitude
+  double contact_area = 0.0;  // changing rectangular influence
+  if (Tdim == 2) contact_area = std::sqrt(volume_);  // unit hexahedron
+  if (Tdim == 3) contact_area = std::cbrt(volume_);  // cube
+  double tangent_adhesion = levelset_alpha_ * contact_area;
+
+  // Calculate tangential reaction force magntiude
+  double reaction_tangent_mag = tangent_friction + tangent_adhesion;
+
+  // Calculate tangential contact force magnitude
+  double contact_tangent_mag =
+      (mass_ * contact_vel_ / dt).dot(levelset_tangent_);
+
+  // Reaction force must not exceed cancellation of contact tangential force
+  reaction_tangent_mag = std::min(reaction_tangent_mag, contact_tangent_mag);
+
+  // Calculate tangential reaction force vector
+  VectorDim reaction_force_tangent = -levelset_tangent_ * reaction_tangent_mag;
+
+  // Calculate total reaction force vector
+  reaction_force_ = reaction_force_normal + reaction_force_tangent;
+
+  // Damp reaction force if particle moving away from interface
+  if ((contact_vel_.dot(levelset_normal_)) >= 0.)
+    reaction_force_ = (1. - levelset_damping_) * reaction_force_;
+}
+
+//! Map levelset contact reaction force to nodes
+template <unsigned Tdim>
+void mpm::ParticleLevelset<Tdim>::map_contact_force_to_nodes() noexcept {
+  for (unsigned i = 0; i < nodes_.size(); ++i) {
+    nodes_[i]->update_external_force(true, mpm::ParticlePhase::Solid,
+                                     (shapefn_[i] * reaction_force_));
+  }
+}
\ No newline at end of file
diff --git a/include/particles/particle_twophase.tcc b/include/particles/particle_twophase.tcc
index 7370520f8..2561b7db4 100644
--- a/include/particles/particle_twophase.tcc
+++ b/include/particles/particle_twophase.tcc
@@ -57,8 +57,9 @@ std::shared_ptr<void> mpm::TwoPhaseParticle<Tdim>::pod() const {
   // Particle local size
   Eigen::Vector3d nsize;
   nsize.setZero();
-  Eigen::VectorXd size = this->natural_size();
-  for (unsigned j = 0; j < Tdim; ++j) nsize[j] = size[j];
+  for (unsigned j = 0; j < Tdim; ++j) {
+    nsize[j] = natural_size_[j];
+  }
 
   Eigen::Matrix<double, 6, 1> stress = this->stress_;
 
@@ -91,6 +92,8 @@ std::shared_ptr<void> mpm::TwoPhaseParticle<Tdim>::pod() const {
   particle_data->displacement_y = displacement[1];
   particle_data->displacement_z = displacement[2];
 
+  particle_data->size = this->size_;
+
   particle_data->nsize_x = nsize[0];
   particle_data->nsize_y = nsize[1];
   particle_data->nsize_z = nsize[2];
@@ -1009,16 +1012,19 @@ bool mpm::TwoPhaseParticle<Tdim>::assign_traction(unsigned direction,
       throw std::runtime_error(
           "Particle traction property: volume / direction is invalid");
     }
+    // Updated size
+    const double updated_size =
+        std::pow(this->volume_, static_cast<double>(1. / Tdim));
     // Assign mixture traction
     if (direction < Tdim) {
       this->set_traction_ = true;
-      traction_(direction) = traction * this->volume_ / this->size_(direction);
+      traction_(direction) = traction * this->volume_ / updated_size;
     }
     // Assign liquid traction
     else {
       this->set_liquid_traction_ = true;
       liquid_traction_(direction - Tdim) =
-          traction * this->volume_ / this->size_(direction - Tdim);
+          traction * this->volume_ / updated_size;
     }
     status = true;
   } catch (std::exception& exception) {
@@ -1118,6 +1124,9 @@ std::vector<uint8_t> mpm::TwoPhaseParticle<Tdim>::serialize() {
   // Displacement
   MPI_Pack(displacement_.data(), Tdim, MPI_DOUBLE, data_ptr, data.size(),
            &position, MPI_COMM_WORLD);
+  // Size
+  MPI_Pack(&size_, 1, MPI_DOUBLE, data_ptr, data.size(), &position,
+           MPI_COMM_WORLD);
   // Natural size
   MPI_Pack(natural_size_.data(), Tdim, MPI_DOUBLE, data_ptr, data.size(),
            &position, MPI_COMM_WORLD);
@@ -1267,6 +1276,9 @@ void mpm::TwoPhaseParticle<Tdim>::deserialize(
   // Displacement
   MPI_Unpack(data_ptr, data.size(), &position, displacement_.data(), Tdim,
              MPI_DOUBLE, MPI_COMM_WORLD);
+  // Size
+  MPI_Unpack(data_ptr, data.size(), &position, &size_, 1, MPI_DOUBLE,
+             MPI_COMM_WORLD);
   // Natural size
   MPI_Unpack(data_ptr, data.size(), &position, natural_size_.data(), Tdim,
              MPI_DOUBLE, MPI_COMM_WORLD);
diff --git a/include/particles/pod_particles/pod_particle.h b/include/particles/pod_particles/pod_particle.h
index 392d0cd3a..a230b5b7b 100644
--- a/include/particles/pod_particles/pod_particle.h
+++ b/include/particles/pod_particles/pod_particle.h
@@ -22,7 +22,8 @@ typedef struct PODParticle {
   double coord_x, coord_y, coord_z;
   // Displacement
   double displacement_x, displacement_y, displacement_z;
-  // Natural particle size
+  // Initial particle size and Natural particle size
+  double size;
   double nsize_x, nsize_y, nsize_z;
   // Velocity
   double velocity_x, velocity_y, velocity_z;
@@ -59,7 +60,7 @@ typedef struct PODParticle {
 
 namespace pod {
 namespace particle {
-const hsize_t NFIELDS = 74;
+const hsize_t NFIELDS = 75;
 
 const size_t dst_size = sizeof(PODParticle);
 
diff --git a/include/particles/pod_particles/pod_particle_twophase.h b/include/particles/pod_particles/pod_particle_twophase.h
index 19f0bcb55..108515323 100644
--- a/include/particles/pod_particles/pod_particle_twophase.h
+++ b/include/particles/pod_particles/pod_particle_twophase.h
@@ -27,7 +27,7 @@ typedef struct PODParticleTwoPhase : PODParticle {
 
 namespace pod {
 namespace particletwophase {
-const hsize_t NFIELDS = 87;
+const hsize_t NFIELDS = 88;
 
 const size_t dst_size = sizeof(PODParticleTwoPhase);
 
diff --git a/include/solvers/mpm_base.h b/include/solvers/mpm_base.h
index 67d26e9fe..1f79da159 100644
--- a/include/solvers/mpm_base.h
+++ b/include/solvers/mpm_base.h
@@ -24,6 +24,7 @@
 #include "constraints.h"
 #include "contact.h"
 #include "contact_friction.h"
+#include "contact_levelset.h"
 #include "mpm.h"
 #include "mpm_scheme.h"
 #include "mpm_scheme_musl.h"
@@ -221,6 +222,24 @@ class MPMBase : public MPM {
   //! Initialise particle types
   void initialise_particle_types();
 
+  /**
+   * \defgroup Interface Functions (includes multimaterial and levelset)
+   */
+  /**@{*/
+  //! \ingroup Interface
+  //! Return if interface and levelset are active
+  //! \retval levelset status of mesh
+  bool is_levelset();
+
+  //! \ingroup Interface
+  //! Nodal levelset inputs
+  //! \param[in] mesh_prop Mesh properties
+  //! \param[in] mesh_io Mesh IO handle
+  void interface_inputs(const Json& mesh_prop,
+                        const std::shared_ptr<mpm::IOMesh<Tdim>>& mesh_io);
+
+  /**@}*/
+
   /**
    * \defgroup Implicit Functions dealing with implicit MPM
    */
@@ -259,8 +278,6 @@ class MPMBase : public MPM {
   std::string stress_update_{"usf"};
   //! Stress update scheme
   std::shared_ptr<mpm::MPMScheme<Tdim>> mpm_scheme_{nullptr};
-  //! Interface scheme
-  std::shared_ptr<mpm::Contact<Tdim>> contact_{nullptr};
   //! Velocity update method
   mpm::VelocityUpdate velocity_update_{mpm::VelocityUpdate::FLIP};
   //! FLIP-PIC blending ratio
@@ -294,6 +311,24 @@ class MPMBase : public MPM {
   //! Absorbing Boundary Variables
   bool absorbing_boundary_{false};
 
+  /**
+   * \defgroup Interface Variables (includes multimaterial and levelset)
+   * @{
+   */
+  //! Interface scheme
+  std::shared_ptr<mpm::Contact<Tdim>> contact_{nullptr};
+  //! Interface bool
+  bool interface_{false};
+  //! Interface type
+  std::string interface_type_{"none"};
+  //! Levelset damping factor
+  double levelset_damping_{0.05};
+  //! Levelset PIC contact velocity
+  bool levelset_pic_{true};
+  //! Levelset violation correction factor
+  double levelset_violation_corrector_{0.01};
+  /**@}*/
+
   /**
    * \defgroup Nonlocal Variables for nonlocal MPM
    * @{
diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc
index 8d1fa8ff5..87fb906f2 100644
--- a/include/solvers/mpm_base.tcc
+++ b/include/solvers/mpm_base.tcc
@@ -13,12 +13,17 @@ mpm::MPMBase<Tdim>::MPMBase(const std::shared_ptr<IO>& io) : mpm::MPM(io) {
   // Set mesh as isoparametric
   bool isoparametric = is_isoparametric();
 
-  mesh_ = std::make_shared<mpm::Mesh<Tdim>>(id, isoparametric);
+  // Construct mesh, use levelset mesh if levelset active
+  if (is_levelset()) {
+    mesh_ = std::make_shared<mpm::MeshLevelset<Tdim>>(id, isoparametric);
+  } else {
+    mesh_ = std::make_shared<mpm::Mesh<Tdim>>(id, isoparametric);
+  }
 
   // Create constraints
   constraints_ = std::make_shared<mpm::Constraints<Tdim>>(mesh_);
 
-  // Empty all materials
+  // Empty all materSials
   materials_.clear();
 
   // Variable list
@@ -29,10 +34,12 @@ mpm::MPMBase<Tdim>::MPMBase(const std::shared_ptr<IO>& io) : mpm::MPM(io) {
       {"mass", VariableType::Scalar},
       {"volume", VariableType::Scalar},
       {"mass_density", VariableType::Scalar},
+      {"levelset", VariableType::Scalar},
       // Vector variables
       {"displacements", VariableType::Vector},
       {"velocities", VariableType::Vector},
       {"normals", VariableType::Vector},
+      {"levelset_reaction_forces", VariableType::Vector},
       // Tensor variables
       {"strains", VariableType::Tensor},
       {"stresses", VariableType::Tensor}};
@@ -291,6 +298,9 @@ void mpm::MPMBase<Tdim>::initialise_mesh() {
   // Read and assign absorbing constraintes
   this->nodal_absorbing_constraints(mesh_props, mesh_io);
 
+  // Read and assign interface (includes multimaterial and levelset)
+  this->interface_inputs(mesh_props, mesh_io);
+
   // Initialise cell
   auto cells_begin = std::chrono::steady_clock::now();
   // Shape function name
@@ -768,6 +778,23 @@ bool mpm::MPMBase<Tdim>::is_isoparametric() {
   return isoparametric;
 }
 
+//! Return if interface and levelset are active
+template <unsigned Tdim>
+bool mpm::MPMBase<Tdim>::is_levelset() {
+  bool levelset_active = true;
+
+  const auto mesh_props = io_->json_object("mesh");
+  if (mesh_props.find("interface") != mesh_props.end()) {
+    this->interface_ = true;
+    this->interface_type_ =
+        mesh_props["interface"]["interface_type"].template get<std::string>();
+    if (interface_type_ != "levelset") levelset_active = false;
+  } else
+    levelset_active = false;
+
+  return levelset_active;
+}
+
 //! Initialise loads
 template <unsigned Tdim>
 void mpm::MPMBase<Tdim>::initialise_loads() {
@@ -1255,6 +1282,109 @@ void mpm::MPMBase<Tdim>::nodal_adhesional_constraints(
   }
 }
 
+// Interface inputs (includes multimaterial and levelset)
+template <unsigned Tdim>
+void mpm::MPMBase<Tdim>::interface_inputs(
+    const Json& mesh_props, const std::shared_ptr<mpm::IOMesh<Tdim>>& mesh_io) {
+  try {
+    // Read and assign interface inputs
+    if (mesh_props.find("interface") != mesh_props.end()) {
+      this->interface_ = true;
+      // Get interface type
+      this->interface_type_ =
+          mesh_props["interface"]["interface_type"].template get<std::string>();
+      if (interface_type_ == "levelset") {
+        // Check if levelset inputs are specified in a file
+        if (mesh_props["interface"].find("location") !=
+            mesh_props["interface"].end()) {
+          // Retrieve the file name
+          std::string levelset_input_file =
+              mesh_props["interface"]["location"].template get<std::string>();
+          // Retrieve levelset info from file
+          bool levelset_inputs =
+              mesh_->assign_nodal_levelset_values(mesh_io->read_levelset_input(
+                  io_->file_name(levelset_input_file)));
+          if (!levelset_inputs) {
+            throw std::runtime_error(
+                "Levelset interface is undefined; Levelset inputs are not "
+                "properly assigned");
+          }
+        } else {
+          throw std::runtime_error(
+              "Levelset interface is undefined; Levelset location is not "
+              "specified");
+        }
+        // Check if levelset damping factor is specified
+        if (mesh_props["interface"].find("damping") !=
+            mesh_props["interface"].end()) {
+          // Retrieve levelset damping factor
+          levelset_damping_ =
+              mesh_props["interface"]["damping"].template get<double>();
+          if ((levelset_damping_ < 0.) || (levelset_damping_ > 1.)) {
+            levelset_damping_ = 0.05;
+            throw std::runtime_error(
+                "Levelset damping factor is not properly specified, using "
+                "0.05 as default");
+          }
+        } else {
+          throw std::runtime_error(
+              "Levelset damping is not specified, using 0.05 as default");
+        }
+        // Check if levelset contact velocity update scheme is specified
+        if (mesh_props["interface"].find("velocity_update") !=
+            mesh_props["interface"].end()) {
+          // Retrieve levelset damping factor
+          std::string levelset_velocity_update_ =
+              mesh_props["interface"]["velocity_update"]
+                  .template get<std::string>();
+          if (levelset_velocity_update_ == "global")
+            levelset_pic_ = false;
+          else if (levelset_velocity_update_ != "pic") {
+            throw std::runtime_error(
+                "Levelset contact velocity update is not properly specified, "
+                " using \"pic\" as default");
+          }
+        } else {
+          throw std::runtime_error(
+              "Levelset contact velocity update is not specified, "
+              " using \"pic\" as default");
+        }
+        // Check if levelset violation corrector is specified
+        if (mesh_props["interface"].find("violation_corrector") !=
+            mesh_props["interface"].end()) {
+          // Retrieve levelset violation corrector
+          levelset_violation_corrector_ =
+              mesh_props["interface"]["violation_corrector"]
+                  .template get<double>();
+          if ((levelset_violation_corrector_ < 0.) ||
+              (levelset_violation_corrector_ > 1.)) {
+            levelset_violation_corrector_ = 0.01;
+            throw std::runtime_error(
+                "Levelset violation corrector is not properly specified, using "
+                "0.01 as default");
+          }
+        } else {
+          throw std::runtime_error(
+              "Levelset violation corrector is not specified, using 0.01 as "
+              "default");
+        }
+      } else if (interface_type_ == "multimaterial") {
+        throw std::runtime_error(
+            "Interfaces are undefined; Interface type \"multimaterial\" not "
+            "supported");
+      } else {
+        throw std::runtime_error(
+            "Interfaces are undefined; Interface type not properly specified");
+      }
+    } else {
+      throw std::runtime_error(
+          "Interfaces are undefined; Interfaces JSON data not found");
+    }
+  } catch (std::exception& exception) {
+    console_->warn("#{}: {}", __LINE__, exception.what());
+  }
+}
+
 // Nodal pressure constraints
 template <unsigned Tdim>
 void mpm::MPMBase<Tdim>::nodal_pressure_constraints(
@@ -1570,7 +1700,7 @@ void mpm::MPMBase<Tdim>::particles_pore_pressures(
         if (!io_->file_name(fparticles_pore_pressures).empty()) {
           // Read and assign particles pore pressures
           if (!mesh_->assign_particles_pore_pressures(
-                  particle_io->read_particles_scalar_properties(
+                  particle_io->read_scalar_properties(
                       io_->file_name(fparticles_pore_pressures))))
             throw std::runtime_error(
                 "Particles pore pressures are not properly assigned");
diff --git a/include/solvers/mpm_explicit.h b/include/solvers/mpm_explicit.h
index c08de3e1e..e4eb494fb 100644
--- a/include/solvers/mpm_explicit.h
+++ b/include/solvers/mpm_explicit.h
@@ -78,12 +78,25 @@ class MPMExplicit : public MPMBase<Tdim> {
   //! Absorbing Boundary
   using mpm::MPMBase<Tdim>::absorbing_boundary_;
 
+  /**
+   * \defgroup Interface Variables (includes multimaterial and levelset)
+   * @{
+   */
+  //! Interface boolean
+  using mpm::MPMBase<Tdim>::interface_;
+  //! Interface type
+  using mpm::MPMBase<Tdim>::interface_type_;
+  //! Levelset damping factor
+  using mpm::MPMBase<Tdim>::levelset_damping_;
+  //! Levelset PIC contact velocity
+  using mpm::MPMBase<Tdim>::levelset_pic_;
+  //! Levelset violation correction factor
+  using mpm::MPMBase<Tdim>::levelset_violation_corrector_;
+  /**@}*/
+
  private:
   //! Pressure smoothing
   bool pressure_smoothing_{false};
-  //! Interface
-  bool interface_{false};
-
 };  // MPMExplicit class
 }  // namespace mpm
 
diff --git a/include/solvers/mpm_explicit.tcc b/include/solvers/mpm_explicit.tcc
index 57317f040..c7b5ef5b9 100644
--- a/include/solvers/mpm_explicit.tcc
+++ b/include/solvers/mpm_explicit.tcc
@@ -12,9 +12,13 @@ mpm::MPMExplicit<Tdim>::MPMExplicit(const std::shared_ptr<IO>& io)
   else
     mpm_scheme_ = std::make_shared<mpm::MPMSchemeUSF<Tdim>>(mesh_, dt_);
 
-  //! Interface scheme
   if (this->interface_)
-    contact_ = std::make_shared<mpm::ContactFriction<Tdim>>(mesh_);
+    if (this->interface_type_ == "multimaterial")
+      contact_ = std::make_shared<mpm::ContactFriction<Tdim>>(mesh_);
+    else if (this->interface_type_ == "levelset")
+      contact_ = std::make_shared<mpm::ContactLevelset<Tdim>>(mesh_);
+    else  // default is "none"
+      contact_ = std::make_shared<mpm::Contact<Tdim>>(mesh_);
   else
     contact_ = std::make_shared<mpm::Contact<Tdim>>(mesh_);
 }
@@ -55,9 +59,6 @@ bool mpm::MPMExplicit<Tdim>::solve() {
   // Pressure smoothing
   pressure_smoothing_ = io_->analysis_bool("pressure_smoothing");
 
-  // Interface
-  interface_ = io_->analysis_bool("interface");
-
   // Initialise material
   this->initialise_materials();
 
@@ -101,6 +102,12 @@ bool mpm::MPMExplicit<Tdim>::solve() {
   // Create nodal properties
   if (interface_ or absorbing_boundary_) mesh_->create_nodal_properties();
 
+  // Initialise levelset properties
+  if (this->interface_type_ == "levelset")
+    contact_->initialise_levelset_properties(
+        this->levelset_damping_, this->levelset_pic_,
+        this->levelset_violation_corrector_);
+
   // Initialise loading conditions
   this->initialise_loads();
 
@@ -127,14 +134,15 @@ bool mpm::MPMExplicit<Tdim>::solve() {
     // Initialise nodes, cells and shape functions
     mpm_scheme_->initialise();
 
-    // Initialise nodal properties and append material ids to node
-    contact_->initialise();
+    // Initialise multimaterial nodal properties
+    contact_->initialise();  // multimaterial interface
 
     // Mass momentum and compute velocity at nodes
     mpm_scheme_->compute_nodal_kinematics(velocity_update_, phase);
 
-    // Map material properties to nodes
-    contact_->compute_contact_forces();
+    // Contact reaction forces at nodes
+    contact_->compute_contact_forces(dt_);  // levelset interface
+    contact_->compute_contact_forces();     // multimaterial interface
 
     // Update stress first
     mpm_scheme_->precompute_stress_strain(phase, pressure_smoothing_);
diff --git a/src/node.cc b/src/node.cc
index ea93ebb86..dd0e82b98 100644
--- a/src/node.cc
+++ b/src/node.cc
@@ -1,6 +1,7 @@
 #include "node.h"
 #include "factory.h"
 #include "node_base.h"
+#include "node_levelset.h"
 
 // Node2D (2 DoF, 1 Phase)
 static Register<mpm::NodeBase<2>, mpm::Node<2, 2, 1>, mpm::Index,
@@ -20,4 +21,14 @@ static Register<mpm::NodeBase<2>, mpm::Node<2, 2, 2>, mpm::Index,
 // Node3D (3 DoF, 2 Phase)
 static Register<mpm::NodeBase<3>, mpm::Node<3, 3, 2>, mpm::Index,
                 const Eigen::Matrix<double, 3, 1>&>
-    node3d2phase("N3D2P");
\ No newline at end of file
+    node3d2phase("N3D2P");
+
+// Node2D (2 DoF, 1 Phase)
+static Register<mpm::NodeBase<2>, mpm::NodeLevelset<2, 2, 1>, mpm::Index,
+                const Eigen::Matrix<double, 2, 1>&>
+    node2dlevelset("N2DLS");
+
+// Node3D (3 DoF, 1 Phase)
+static Register<mpm::NodeBase<3>, mpm::NodeLevelset<3, 3, 1>, mpm::Index,
+                const Eigen::Matrix<double, 3, 1>&>
+    node3dlevelset("N3DLS");
\ No newline at end of file
diff --git a/src/particle.cc b/src/particle.cc
index 3b2878735..e0e689e10 100644
--- a/src/particle.cc
+++ b/src/particle.cc
@@ -4,18 +4,22 @@
 #include "particle_bbar.h"
 #include "particle_finite_strain.h"
 #include "particle_fluid.h"
+#include "particle_levelset.h"
+#include "particle_levelset_bbar.h"
 #include "particle_twophase.h"
 
 namespace mpm {
 // ParticleType
 std::map<std::string, int> ParticleType = {
-    {"P2D", 0},       {"P3D", 1},       {"P2DFLUID", 2}, {"P3DFLUID", 3},
-    {"P2D2PHASE", 4}, {"P3D2PHASE", 5}, {"P2DBBAR", 6},  {"P3DBBAR", 7},
-    {"P2DFS", 8},     {"P3DFS", 9}};
+    {"P2D", 0},        {"P3D", 1},       {"P2DFLUID", 2}, {"P3DFLUID", 3},
+    {"P2D2PHASE", 4},  {"P3D2PHASE", 5}, {"P2DBBAR", 6},  {"P3DBBAR", 7},
+    {"P2DFS", 8},      {"P3DFS", 9},     {"P2DLS", 10},   {"P3DLS", 11},
+    {"P2DLSBBAR", 12}, {"P3DLSBBAR", 13}};
 std::map<int, std::string> ParticleTypeName = {
-    {0, "P2D"},       {1, "P3D"},       {2, "P2DFLUID"}, {3, "P3DFLUID"},
-    {4, "P2D2PHASE"}, {5, "P3D2PHASE"}, {6, "P2DBBAR"},  {7, "P3DBBAR"},
-    {8, "P2DFS"},     {9, "P3DFS"}};
+    {0, "P2D"},        {1, "P3D"},       {2, "P2DFLUID"}, {3, "P3DFLUID"},
+    {4, "P2D2PHASE"},  {5, "P3D2PHASE"}, {6, "P2DBBAR"},  {7, "P3DBBAR"},
+    {8, "P2DFS"},      {9, "P3DFS"},     {10, "P2DLS"},   {11, "P3DLS"},
+    {12, "P2DLSBBAR"}, {13, "P3DLSBBAR"}};
 std::map<std::string, std::string> ParticlePODTypeName = {
     {"P2D", "particles"},
     {"P3D", "particles"},
@@ -26,7 +30,11 @@ std::map<std::string, std::string> ParticlePODTypeName = {
     {"P2DBBAR", "particles"},
     {"P3DBBAR", "particles"},
     {"P2DFS", "particles"},
-    {"P3DFS", "particles"}};
+    {"P3DFS", "particles"},
+    {"P2DLS", "particles"},
+    {"P3DLS", "particles"},
+    {"P2DLSBBAR", "particles"},
+    {"P3DLSBBAR", "particles"}};
 }  // namespace mpm
 
 // Particle2D (2 Dim)
@@ -77,4 +85,24 @@ static Register<mpm::ParticleBase<2>, mpm::ParticleFiniteStrain<2>, mpm::Index,
 // Particle3D with finite strain formulation (3 Dim)
 static Register<mpm::ParticleBase<3>, mpm::ParticleFiniteStrain<3>, mpm::Index,
                 const Eigen::Matrix<double, 3, 1>&>
-    particle3dfinitestrain("P3DFS");
\ No newline at end of file
+    particle3dfinitestrain("P3DFS");
+
+// Particle2D with Levelset (2 Dim)
+static Register<mpm::ParticleBase<2>, mpm::ParticleLevelset<2>, mpm::Index,
+                const Eigen::Matrix<double, 2, 1>&>
+    particle2dlevelset("P2DLS");
+
+// Particle3D with Levelset (3 Dim)
+static Register<mpm::ParticleBase<3>, mpm::ParticleLevelset<3>, mpm::Index,
+                const Eigen::Matrix<double, 3, 1>&>
+    particle3dlevelset("P3DLS");
+
+// Particle2D with Levelset and B-bar method (2 Dim)
+static Register<mpm::ParticleBase<2>, mpm::ParticleLevelsetBbar<2>, mpm::Index,
+                const Eigen::Matrix<double, 2, 1>&>
+    particle2dlevelsetbbar("P2DLSBBAR");
+
+// Particle3D with Levelset and B-bar method (3 Dim)
+static Register<mpm::ParticleBase<3>, mpm::ParticleLevelsetBbar<3>, mpm::Index,
+                const Eigen::Matrix<double, 3, 1>&>
+    particle3dlevelsetbbar("P3DLSBBAR");
\ No newline at end of file
diff --git a/src/pod_particle.cc b/src/pod_particle.cc
index c2f46bdec..e0a4e8a2a 100644
--- a/src/pod_particle.cc
+++ b/src/pod_particle.cc
@@ -13,6 +13,7 @@ const size_t dst_offset[NFIELDS] = {
     HOFFSET(PODParticle, displacement_x),
     HOFFSET(PODParticle, displacement_y),
     HOFFSET(PODParticle, displacement_z),
+    HOFFSET(PODParticle, size),
     HOFFSET(PODParticle, nsize_x),
     HOFFSET(PODParticle, nsize_y),
     HOFFSET(PODParticle, nsize_z),
@@ -92,6 +93,7 @@ const size_t dst_sizes[NFIELDS] = {
     sizeof(particle.displacement_x),
     sizeof(particle.displacement_y),
     sizeof(particle.displacement_z),
+    sizeof(particle.size),
     sizeof(particle.nsize_x),
     sizeof(particle.nsize_y),
     sizeof(particle.nsize_z),
@@ -170,6 +172,7 @@ const char* field_names[NFIELDS] = {
     "displacement_x",
     "displacement_y",
     "displacement_z",
+    "size",
     "nsize_x",
     "nsize_y",
     "nsize_z",
@@ -248,15 +251,15 @@ const hid_t field_type[NFIELDS] = {
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
-    H5T_NATIVE_HBOOL,  H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
+    H5T_NATIVE_DOUBLE, H5T_NATIVE_HBOOL,  H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
-    H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_HBOOL,  H5T_NATIVE_LLONG,
-    H5T_NATIVE_UINT,   H5T_NATIVE_UINT,   H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
+    H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_HBOOL,
+    H5T_NATIVE_LLONG,  H5T_NATIVE_UINT,   H5T_NATIVE_UINT,   H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
-    H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE};
+    H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE};
 }  // namespace particle
 }  // namespace pod
 }  // namespace mpm
diff --git a/src/pod_particle_twophase.cc b/src/pod_particle_twophase.cc
index 9330382e1..10462953b 100644
--- a/src/pod_particle_twophase.cc
+++ b/src/pod_particle_twophase.cc
@@ -14,6 +14,7 @@ const size_t dst_offset[NFIELDS] = {
     HOFFSET(PODParticleTwoPhase, displacement_x),
     HOFFSET(PODParticleTwoPhase, displacement_y),
     HOFFSET(PODParticleTwoPhase, displacement_z),
+    HOFFSET(PODParticleTwoPhase, size),
     HOFFSET(PODParticleTwoPhase, nsize_x),
     HOFFSET(PODParticleTwoPhase, nsize_y),
     HOFFSET(PODParticleTwoPhase, nsize_z),
@@ -108,6 +109,7 @@ const size_t dst_sizes[NFIELDS] = {
     sizeof(particle.displacement_x),
     sizeof(particle.displacement_y),
     sizeof(particle.displacement_z),
+    sizeof(particle.size),
     sizeof(particle.nsize_x),
     sizeof(particle.nsize_y),
     sizeof(particle.nsize_z),
@@ -201,6 +203,7 @@ const char* field_names[NFIELDS] = {
     "displacement_x",
     "displacement_y",
     "displacement_z",
+    "size",
     "nsize_x",
     "nsize_y",
     "nsize_z",
@@ -293,18 +296,18 @@ const hid_t field_type[NFIELDS] = {
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
-    H5T_NATIVE_HBOOL,  H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
+    H5T_NATIVE_DOUBLE, H5T_NATIVE_HBOOL,  H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
-    H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_HBOOL,  H5T_NATIVE_LLONG,
-    H5T_NATIVE_UINT,   H5T_NATIVE_UINT,   H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
+    H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_HBOOL,
+    H5T_NATIVE_LLONG,  H5T_NATIVE_UINT,   H5T_NATIVE_UINT,   H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
     H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
-    H5T_NATIVE_UINT,   H5T_NATIVE_UINT,   H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
-    H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE};
+    H5T_NATIVE_DOUBLE, H5T_NATIVE_UINT,   H5T_NATIVE_UINT,   H5T_NATIVE_DOUBLE,
+    H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE};
 }  // namespace particletwophase
 }  // namespace pod
 }  // namespace mpm
diff --git a/tests/include/write_mesh_particles.h b/tests/include/write_mesh_particles.h
index d71ec3074..ea3e58d16 100644
--- a/tests/include/write_mesh_particles.h
+++ b/tests/include/write_mesh_particles.h
@@ -13,6 +13,11 @@ bool write_json_warnings(unsigned dim, bool material_sets, bool math_functions,
                          std::string math_functions_type, int grav_dim,
                          std::string file_name);
 
+// Write JSON Configuration file for levelset boundary
+bool write_json_levelset(unsigned dim, bool resume, const std::string& analysis,
+                         const std::string& file_name, bool BBar,
+                         bool levelset_pic);
+
 // Write JSON Configuration file for absorbing boundary
 bool write_json_absorbing(unsigned dim, bool resume,
                           const std::string& analysis,
diff --git a/tests/io/io_mesh_ascii_test.cc b/tests/io/io_mesh_ascii_test.cc
index 58356a53c..3c174ffeb 100644
--- a/tests/io/io_mesh_ascii_test.cc
+++ b/tests/io/io_mesh_ascii_test.cc
@@ -6,7 +6,7 @@ using Json = nlohmann::json;
 
 #include "io_mesh_ascii.h"
 
-// Check IOMeshAscii
+//! \brief Check IOMeshAscii for 2D
 TEST_CASE("IOMeshAscii is checked for 2D", "[IOMesh][IOMeshAscii][2D]") {
 
   // Dimension
@@ -191,6 +191,74 @@ TEST_CASE("IOMeshAscii is checked for 2D", "[IOMesh][IOMeshAscii][2D]") {
     }
   }
 
+  SECTION("Check levelset file") {
+    // Vector of levelset inputs
+    std::vector<std::tuple<mpm::Index, double, double, double, double>>
+        levelset_input_file;
+
+    // Inputs
+    levelset_input_file.emplace_back(
+        std::make_tuple(0, 0.00, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(1, 0.00, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(2, 0.50, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(3, 0.50, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(4, 0.00, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(5, 0.50, 0.1, 1000., 1.0E+06));
+
+    // Dump levelset inputs as a file to be read
+    std::ofstream file;
+    file.open("levelset-nodal-values-2d.txt");
+    // Write particle coordinates
+    for (const auto& levelset_inputs : levelset_input_file) {
+      file << std::get<0>(levelset_inputs) << "\t";
+      file << std::get<1>(levelset_inputs) << "\t";
+      file << std::get<2>(levelset_inputs) << "\t";
+      file << std::get<3>(levelset_inputs) << "\t";
+      file << std::get<4>(levelset_inputs) << "\t";
+
+      file << "\n";
+    }
+
+    file.close();
+  }
+
+  // Check read levelset inputs
+  SECTION("Check levelset inputs") {
+    // Create a read_mesh object
+    auto read_mesh = std::make_unique<mpm::IOMeshAscii<dim>>();
+
+    // Try to read inputs from a non-existant file
+    auto levelset_inputs =
+        read_mesh->read_levelset_input("levelset-nodal-values-missing.txt");
+    // Check number of inputs
+    REQUIRE(levelset_inputs.size() == 0);
+
+    // Check inputs
+    levelset_inputs =
+        read_mesh->read_levelset_input("levelset-nodal-values-2d.txt");
+    // Check number of particles
+    REQUIRE(levelset_inputs.size() == levelset_inputs.size());
+
+    // Check coordinates of nodes
+    for (unsigned i = 0; i < levelset_inputs.size(); ++i) {
+      REQUIRE(std::get<0>(levelset_inputs.at(i)) ==
+              Approx(std::get<0>(levelset_inputs.at(i))).epsilon(Tolerance));
+      REQUIRE(std::get<1>(levelset_inputs.at(i)) ==
+              Approx(std::get<1>(levelset_inputs.at(i))).epsilon(Tolerance));
+      REQUIRE(std::get<2>(levelset_inputs.at(i)) ==
+              Approx(std::get<2>(levelset_inputs.at(i))).epsilon(Tolerance));
+      REQUIRE(std::get<3>(levelset_inputs.at(i)) ==
+              Approx(std::get<3>(levelset_inputs.at(i))).epsilon(Tolerance));
+      REQUIRE(std::get<4>(levelset_inputs.at(i)) ==
+              Approx(std::get<4>(levelset_inputs.at(i))).epsilon(Tolerance));
+    }
+  }
+
   SECTION("Check displacement constraints file") {
     // Vector of particle coordinates
     std::vector<std::tuple<mpm::Index, unsigned, double>>
@@ -843,14 +911,14 @@ TEST_CASE("IOMeshAscii is checked for 2D", "[IOMesh][IOMeshAscii][2D]") {
       auto io_mesh = std::make_unique<mpm::IOMeshAscii<dim>>();
 
       // Try to read particles scalar properties from a non-existant file
-      auto read_particles_scalars = io_mesh->read_particles_scalar_properties(
-          "particles-scalar-missing.txt");
+      auto read_particles_scalars =
+          io_mesh->read_scalar_properties("particles-scalar-missing.txt");
       // Check number of particles scalar properties
       REQUIRE(read_particles_scalars.size() == 0);
 
       // Check particles scalar properties
       read_particles_scalars =
-          io_mesh->read_particles_scalar_properties("particles-scalars-2d.txt");
+          io_mesh->read_scalar_properties("particles-scalars-2d.txt");
 
       // Check number of particles
       REQUIRE(read_particles_scalars.size() == particles_scalars.size());
@@ -912,7 +980,7 @@ TEST_CASE("IOMeshAscii is checked for 2D", "[IOMesh][IOMeshAscii][2D]") {
   }
 }
 
-// Check IOMeshAscii
+//! \brief Check IOMeshAscii for 3D
 TEST_CASE("IOMeshAscii is checked for 3D", "[IOMesh][IOMeshAscii][3D]") {
 
   // Dimension
@@ -1139,6 +1207,86 @@ TEST_CASE("IOMeshAscii is checked for 3D", "[IOMesh][IOMeshAscii][3D]") {
     }
   }
 
+  SECTION("Check levelset file") {
+    // Vector of levelset inputs
+    std::vector<std::tuple<mpm::Index, double, double, double, double>>
+        levelset_input_file;
+
+    // Inputs
+    levelset_input_file.emplace_back(
+        std::make_tuple(0, 0.00, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(1, 0.00, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(2, 0.50, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(3, 0.50, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(4, 0.00, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(5, 0.00, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(6, 0.50, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(7, 0.50, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(8, 0.00, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(9, 0.50, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(10, 0.00, 0.1, 1000., 1.0E+06));
+    levelset_input_file.emplace_back(
+        std::make_tuple(11, 0.50, 0.1, 1000., 1.0E+06));
+
+    // Dump levelset inputs as a file to be read
+    std::ofstream file;
+    file.open("levelset-nodal-values-3d.txt");
+    // Write particle coordinates
+    for (const auto& levelset_inputs : levelset_input_file) {
+      file << std::get<0>(levelset_inputs) << "\t";
+      file << std::get<1>(levelset_inputs) << "\t";
+      file << std::get<2>(levelset_inputs) << "\t";
+      file << std::get<3>(levelset_inputs) << "\t";
+      file << std::get<4>(levelset_inputs) << "\t";
+
+      file << "\n";
+    }
+
+    file.close();
+  }
+
+  // Check read levelset inputs
+  SECTION("Check levelset inputs") {
+    // Create a read_mesh object
+    auto read_mesh = std::make_unique<mpm::IOMeshAscii<dim>>();
+
+    // Try to read inputs from a non-existant file
+    auto levelset_inputs =
+        read_mesh->read_levelset_input("levelset-nodal-values-missing.txt");
+    // Check number of inputs
+    REQUIRE(levelset_inputs.size() == 0);
+
+    // Check inputs
+    levelset_inputs =
+        read_mesh->read_levelset_input("levelset-nodal-values-3d.txt");
+    // Check number of particles
+    REQUIRE(levelset_inputs.size() == levelset_inputs.size());
+
+    // Check coordinates of nodes
+    for (unsigned i = 0; i < levelset_inputs.size(); ++i) {
+      REQUIRE(std::get<0>(levelset_inputs.at(i)) ==
+              Approx(std::get<0>(levelset_inputs.at(i))).epsilon(Tolerance));
+      REQUIRE(std::get<1>(levelset_inputs.at(i)) ==
+              Approx(std::get<1>(levelset_inputs.at(i))).epsilon(Tolerance));
+      REQUIRE(std::get<2>(levelset_inputs.at(i)) ==
+              Approx(std::get<2>(levelset_inputs.at(i))).epsilon(Tolerance));
+      REQUIRE(std::get<3>(levelset_inputs.at(i)) ==
+              Approx(std::get<3>(levelset_inputs.at(i))).epsilon(Tolerance));
+      REQUIRE(std::get<4>(levelset_inputs.at(i)) ==
+              Approx(std::get<4>(levelset_inputs.at(i))).epsilon(Tolerance));
+    }
+  }
+
   SECTION("Check displacement constraints file") {
     // Vector of particle coordinates
     std::vector<std::tuple<mpm::Index, unsigned, double>>
@@ -1792,14 +1940,14 @@ TEST_CASE("IOMeshAscii is checked for 3D", "[IOMesh][IOMeshAscii][3D]") {
       auto io_mesh = std::make_unique<mpm::IOMeshAscii<dim>>();
 
       // Try to read particles scalar properties from a non-existant file
-      auto read_particles_scalars = io_mesh->read_particles_scalar_properties(
-          "particles-scalar-missing.txt");
+      auto read_particles_scalars =
+          io_mesh->read_scalar_properties("particles-scalar-missing.txt");
       // Check number of particles scalar properties
       REQUIRE(read_particles_scalars.size() == 0);
 
       // Check particles scalar properties
       read_particles_scalars =
-          io_mesh->read_particles_scalar_properties("particles-scalars-3d.txt");
+          io_mesh->read_scalar_properties("particles-scalars-3d.txt");
 
       // Check number of particles
       REQUIRE(read_particles_scalars.size() == particles_scalars.size());
diff --git a/tests/io/write_mesh_particles.cc b/tests/io/write_mesh_particles.cc
index 023a94b67..5b35a8b95 100644
--- a/tests/io/write_mesh_particles.cc
+++ b/tests/io/write_mesh_particles.cc
@@ -201,6 +201,109 @@ bool write_json_warnings(unsigned dim, bool material_sets, bool math_functions,
   return true;
 }
 
+// Write JSON Configuration file for levelset boundary
+bool write_json_levelset(unsigned dim, bool resume, const std::string& analysis,
+                         const std::string& file_name, bool BBar,
+                         bool levelset_pic) {
+  // Make json object with input files
+  // 2D
+  std::string dimension = "2d";
+  auto particle_type = "P2DLS";
+  auto node_type = "N2DLS";
+  auto cell_type = "ED2Q4";
+  auto io_type = "Ascii2D";
+  std::string material = "LinearElastic2D";
+  std::vector<double> gravity{{0., -9.81}};
+  unsigned material_id = 0;
+  std::vector<double> xvalues{{0.0, 0.5, 1.0}};
+  std::vector<double> fxvalues{{0.0, 1.0, 1.0}};
+  if (BBar) {
+    particle_type = "P2DLSBBAR";
+  }
+
+  // 3D
+  if (dim == 3) {
+    dimension = "3d";
+    particle_type = "P3DLS";
+    node_type = "N3DLS";
+    cell_type = "ED3H8";
+    io_type = "Ascii3D";
+    material = "LinearElastic3D";
+    gravity.clear();
+    gravity = {0., -9.81, 0.};
+    if (BBar) {
+      particle_type = "P3DLSBBAR";
+    }
+  }
+
+  auto location = "levelset-nodal-values-" + dimension + ".txt";
+  auto levelset_velocity_update = "global";
+  if (levelset_pic) {
+    levelset_velocity_update = "pic";
+  }
+
+  Json json_file = {{"title", "Example JSON Input for MPM"},
+                    {"mesh",
+                     {{"mesh", "mesh-" + dimension + ".txt"},
+                      {"io_type", io_type},
+                      {"check_duplicates", true},
+                      {"isoparametric", false},
+                      {"node_type", node_type},
+                      {"cell_type", cell_type},
+                      {
+                          "interface",
+                          {{"interface_type", "levelset"},
+                           {"location", location},
+                           {"damping", 0.05},
+                           {"velocity_update", levelset_velocity_update},
+                           {"violation_corrector", 0.001}},
+                      }}},
+                    {"particles",
+                     {{{"generator",
+                        {{"type", "file"},
+                         {"material_id", material_id},
+                         {"pset_id", 0},
+                         {"io_type", io_type},
+                         {"particle_type", particle_type},
+                         {"check_duplicates", true},
+                         {"location", "particles-" + dimension + ".txt"}}}}}},
+                    {"materials",
+                     {
+                         {{"id", 0},
+                          {"type", material},
+                          {"density", 1000.},
+                          {"youngs_modulus", 1.0E+6},
+                          {"poisson_ratio", 0.0}},
+                     }},
+                    {"external_loading_conditions", {{"gravity", gravity}}},
+                    {"analysis",
+                     {{"type", analysis},
+                      {"mpm_scheme", "usf"},
+                      {"velocity_update", "flip"},
+                      {"locate_particles", true},
+                      {"dt", 0.001},
+                      {"uuid", file_name + "-" + dimension},
+                      {"nsteps", 5},
+                      {"resume",
+                       {{"resume", resume},
+                        {"uuid", file_name + "-" + dimension},
+                        {"step", 5}}},
+                      {"nload_balance_steps", 1000}}},
+                    {"post_processing",
+                     {{"path", "results/"},
+                      {"vtk", {"id", "levelset", "levelset_reaction_forces"}},
+                      {"output_steps", 5}}}};
+
+  // Dump JSON as an input file to be read
+  std::string fname = (file_name + "-" + dimension + ".json").c_str();
+  std::ofstream file;
+  file.open(fname, std::ios_base::out);
+  file << json_file.dump(2);
+  file.close();
+
+  return true;
+}
+
 // Write JSON Configuration file for absorbing boundary
 bool write_json_absorbing(unsigned dim, bool resume,
                           const std::string& analysis,
diff --git a/tests/mesh/mesh_test_2d.cc b/tests/mesh/mesh_test_2d.cc
index 9a7512462..1443059f5 100644
--- a/tests/mesh/mesh_test_2d.cc
+++ b/tests/mesh/mesh_test_2d.cc
@@ -17,11 +17,13 @@
 #include "hexahedron_element.h"
 #include "linear_function.h"
 #include "mesh.h"
+#include "mesh_levelset.h"
 #include "node.h"
+#include "node_levelset.h"
 #include "partio_writer.h"
 #include "quadrilateral_element.h"
 
-//! Check mesh class for 2D case
+//! \brief Check mesh class for 2D case
 TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") {
   // Dimension
   const unsigned Dim = 2;
@@ -1119,6 +1121,7 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") {
                           particles_stresses) == false);
             }
 
+            //! Boundary Conditions from Entity Sets
             // Test assign particles velocity constraints
             SECTION("Check assign particles velocity constraints") {
               tsl::robin_map<mpm::Index, std::vector<mpm::Index>> particle_sets;
@@ -1459,6 +1462,7 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") {
                       mfunction, -1, 0, pressure) == true);
         }
 
+        //! Boundary Conditions from Files
         // Test assign acceleration constraints to nodes
         SECTION("Check assign acceleration constraints to nodes") {
           // Vector of particle coordinates
@@ -1697,6 +1701,191 @@ TEST_CASE("Mesh is checked for 2D case", "[mesh][2D]") {
     }
   }
 
+  //! Levelset Mesh
+  //! Check create nodes and cells in a levelset mesh
+  SECTION("Check create levelset nodes and cells") {
+    // Vector of nodal coordinates
+    std::vector<Eigen::Matrix<double, Dim, 1>> coordinates;
+
+    // Nodal coordinates
+    Eigen::Matrix<double, Dim, 1> node;
+
+    // Cell 0
+    // Node 0
+    node << 0., 0.;
+    coordinates.emplace_back(node);
+    // Node 1
+    node << 0.5, 0.;
+    coordinates.emplace_back(node);
+    // Node 2
+    node << 0.5, 0.5;
+    coordinates.emplace_back(node);
+    // Node 3
+    node << 0., 0.5;
+    coordinates.emplace_back(node);
+
+    // Cell 1
+    // Node 4
+    node << 1.0, 0.;
+    coordinates.emplace_back(node);
+    // Node 5
+    node << 1.0, 0.5;
+    coordinates.emplace_back(node);
+
+    //! Create a new levelset mesh
+    unsigned meshid = 0;
+    auto mesh = std::make_shared<mpm::MeshLevelset<Dim>>(meshid, false);
+
+    SECTION("Check creation of levelset nodes") {
+      //! Node type 2DLS
+      const std::string node_type = "N2DLS";
+      // Global node index
+      mpm::Index gnid = 0;
+      mesh->create_nodes(gnid, node_type, coordinates, false);
+      // Check if mesh has added nodes
+      REQUIRE(mesh->nnodes() == coordinates.size());
+      // Try again this shouldn't add more coordinates
+      mesh->create_nodes(gnid, node_type, coordinates);
+      // Check if mesh has added nodes
+      REQUIRE(mesh->nnodes() == coordinates.size());
+      // Clear coordinates and try creating a list of nodes with an empty list
+      unsigned nnodes = coordinates.size();
+      coordinates.clear();
+      // This fails with empty list error in node creation
+      mesh->create_nodes(gnid, node_type, coordinates);
+      REQUIRE(mesh->nnodes() == nnodes);
+
+      SECTION("Check creation of levelset cells") {
+        // Cell with node ids
+        std::vector<std::vector<mpm::Index>> cells{// cell #0
+                                                   {0, 1, 2, 3},
+                                                   // cell #1
+                                                   {1, 4, 5, 2}};
+        // Assign 4-noded quadrilateral element to cell
+        std::shared_ptr<mpm::Element<Dim>> element =
+            Factory<mpm::Element<Dim>>::instance()->create("ED2Q4");
+
+        // Global cell index
+        mpm::Index gcid = 0;
+        mesh->create_cells(gcid, element, cells, false);
+        // Check if mesh has added cells
+        REQUIRE(mesh->ncells() == cells.size());
+        // Try again this shouldn't add more cells
+        mesh->create_cells(gcid, element, cells);
+        // Check if mesh has added cells
+        REQUIRE(mesh->ncells() == cells.size());
+        // Clear cells and try creating a list of empty cells
+        unsigned ncells = cells.size();
+        cells.clear();
+        // This fails with empty list error in node creation
+        gcid = 10;
+        mesh->create_cells(gcid, element, cells);
+        REQUIRE(mesh->ncells() == ncells);
+
+        // Try with invalid node ids
+        cells = {// cell #0
+                 {10, 11, 12, 13},
+                 // cell #1
+                 {11, 14, 15, 12}};
+        gcid = 20;
+        mesh->create_cells(gcid, element, cells);
+        REQUIRE(mesh->ncells() == ncells);
+
+        SECTION("Check creation of levelset particles") {
+          // Vector of particle coordinates
+          std::vector<Eigen::Matrix<double, Dim, 1>> coordinates;
+          coordinates.clear();
+
+          // Particle coordinates
+          Eigen::Matrix<double, Dim, 1> particle;
+
+          // Cell 0
+          // Particle 0
+          particle << 0.125, 0.125;
+          coordinates.emplace_back(particle);
+          // Particle 1
+          particle << 0.25, 0.125;
+          coordinates.emplace_back(particle);
+          // Particle 2
+          particle << 0.25, 0.25;
+          coordinates.emplace_back(particle);
+          // Particle 3
+          particle << 0.125, 0.25;
+          coordinates.emplace_back(particle);
+
+          // Cell 1
+          // Particle 4
+          particle << 0.675, 0.125;
+          coordinates.emplace_back(particle);
+          // Particle 5
+          particle << 0.85, 0.125;
+          coordinates.emplace_back(particle);
+          // Particle 6
+          particle << 0.85, 0.25;
+          coordinates.emplace_back(particle);
+          // Particle 7
+          particle << 0.675, 0.25;
+          coordinates.emplace_back(particle);
+
+          // Initialise material models in mesh
+          mesh->initialise_material_models(materials);
+
+          SECTION("Check addition of levelset particles to levelset mesh") {
+            //! Particle type 2DLS
+            const std::string particle_type = "P2DLS";
+            // Create particles from file
+            mesh->create_particles(particle_type, coordinates, mids, 0, false);
+            // Check if mesh has added particles
+            REQUIRE(mesh->nparticles() == coordinates.size());
+          }
+        }
+
+        //! Boundary Conditions from Files
+        // Test assign levelset inputs to nodes from file
+        SECTION("Check assign levelset inputs to nodes") {
+          // Levelset inputs
+          std::vector<std::tuple<mpm::Index, double, double, double, double>>
+              levelset_inputs;
+
+          levelset_inputs.emplace_back(
+              std::make_tuple(0, 0.05, 0.1, 1000., 1.0E+06));
+          levelset_inputs.emplace_back(
+              std::make_tuple(1, 0.05, 0.1, 1000., 1.0E+06));
+          levelset_inputs.emplace_back(
+              std::make_tuple(2, 0.05, 0.1, 1000., 1.0E+06));
+          levelset_inputs.emplace_back(
+              std::make_tuple(3, 0.05, 0.1, 1000., 1.0E+06));
+
+          REQUIRE(mesh != nullptr);
+          REQUIRE(mesh->assign_nodal_levelset_values(levelset_inputs) == true);
+
+          // Check assigned values
+          for (const auto& levelset_input : levelset_inputs) {
+            auto node_id = std::get<0>(levelset_input);
+            auto node = mesh->node(node_id);
+            REQUIRE(node->levelset() ==
+                    Approx(std::get<1>(levelset_input)).epsilon(Tolerance));
+            REQUIRE(node->levelset_mu() ==
+                    Approx(std::get<2>(levelset_input)).epsilon(Tolerance));
+            REQUIRE(node->levelset_alpha() ==
+                    Approx(std::get<3>(levelset_input)).epsilon(Tolerance));
+            REQUIRE(node->barrier_stiffness() ==
+                    Approx(std::get<4>(levelset_input)).epsilon(Tolerance));
+          }
+
+          // Failing inputs (negative mu and alpha)
+          levelset_inputs.emplace_back(
+              std::make_tuple(3, 0.05, -0.1, -1000., 1.0E+06));
+          REQUIRE(mesh->assign_nodal_levelset_values(levelset_inputs) == false);
+          // Failing inputs (negative barrier stiffness)
+          levelset_inputs.emplace_back(
+              std::make_tuple(3, 0.05, 0.1, 1000., -1.0E+06));
+          REQUIRE(mesh->assign_nodal_levelset_values(levelset_inputs) == false);
+        }
+      }
+    }
+  }
+
   //! Check if nodal properties is initialised
   SECTION("Check nodal properties initialisation") {
     // Create the different meshes
diff --git a/tests/mesh/mesh_test_3d.cc b/tests/mesh/mesh_test_3d.cc
index 14ebd0959..3d270875a 100644
--- a/tests/mesh/mesh_test_3d.cc
+++ b/tests/mesh/mesh_test_3d.cc
@@ -17,7 +17,9 @@
 #include "hexahedron_element.h"
 #include "linear_function.h"
 #include "mesh.h"
+#include "mesh_levelset.h"
 #include "node.h"
+#include "node_levelset.h"
 #include "partio_writer.h"
 #include "quadrilateral_element.h"
 
@@ -1244,6 +1246,7 @@ TEST_CASE("Mesh is checked for 3D case", "[mesh][3D]") {
                           particles_stresses) == false);
             }
 
+            //! Boundary Conditions from Entity Sets
             // Test assign particles velocity constraints
             SECTION("Check assign particles velocity constraints") {
               tsl::robin_map<mpm::Index, std::vector<mpm::Index>> particle_sets;
@@ -1557,6 +1560,7 @@ TEST_CASE("Mesh is checked for 3D case", "[mesh][3D]") {
                       set_id, adhesion_constraint) == false);
         }
 
+        //! Boundary Conditions from Files
         // Test assign acceleration constraints to nodes
         SECTION("Check assign acceleration constraints to nodes") {
           // Vector of particle coordinates
@@ -1794,6 +1798,218 @@ TEST_CASE("Mesh is checked for 3D case", "[mesh][3D]") {
     }
   }
 
+  //! Levelset Mesh
+  //! Check create nodes and cells in a levelset mesh
+  SECTION("Check create levelset nodes and cells") {
+    // Vector of nodal coordinates
+    std::vector<Eigen::Matrix<double, Dim, 1>> coordinates;
+
+    // Nodal coordinates
+    Eigen::Matrix<double, Dim, 1> node;
+
+    // Cell 0
+    // Node 0
+    node << 0., 0., 0.;
+    coordinates.emplace_back(node);
+    // Node 1
+    node << 0.5, 0., 0.;
+    coordinates.emplace_back(node);
+    // Node 2
+    node << 0.5, 0.5, 0.;
+    coordinates.emplace_back(node);
+    // Node 3
+    node << 0., 0.5, 0.;
+    coordinates.emplace_back(node);
+    // Node 4
+    node << 0., 0., 0.5;
+    coordinates.emplace_back(node);
+    // Node 5
+    node << 0.5, 0., 0.5;
+    coordinates.emplace_back(node);
+    // Node 6
+    node << 0.5, 0.5, 0.5;
+    coordinates.emplace_back(node);
+    // Node 7
+    node << 0., 0.5, 0.5;
+    coordinates.emplace_back(node);
+
+    // Cell 1
+    // Node 8
+    node << 1.0, 0., 0.;
+    coordinates.emplace_back(node);
+    // Node 9
+    node << 1.0, 0.5, 0.;
+    coordinates.emplace_back(node);
+    // Node 10
+    node << 1.0, 0., 0.5;
+    coordinates.emplace_back(node);
+    // Node 11
+    node << 1.0, 0.5, 0.5;
+    coordinates.emplace_back(node);
+
+    //! Create a new levelset mesh
+    unsigned meshid = 0;
+    auto mesh = std::make_shared<mpm::MeshLevelset<Dim>>(meshid, false);
+
+    SECTION("Check creation of levelset nodes") {
+      //! Node type 3DLS
+      const std::string node_type = "N3DLS";
+      // Global node index
+      mpm::Index gnid = 0;
+      mesh->create_nodes(gnid, node_type, coordinates);
+      // Check if mesh has added nodes
+      REQUIRE(mesh->nnodes() == coordinates.size());
+      // Try again this shouldn't add more coordinates
+      mesh->create_nodes(gnid, node_type, coordinates);
+      // Check if mesh has added nodes
+      REQUIRE(mesh->nnodes() == coordinates.size());
+      // Clear coordinates and try creating a list of nodes with an empty list
+      unsigned nnodes = coordinates.size();
+      coordinates.clear();
+      // This fails with empty list error in node creation
+      mesh->create_nodes(gnid, node_type, coordinates);
+      REQUIRE(mesh->nnodes() == nnodes);
+
+      SECTION("Check creation of levelset cells") {
+        // Cell with node ids
+        std::vector<std::vector<mpm::Index>> cells{// cell #0
+                                                   {0, 1, 2, 3, 4, 5, 6, 7},
+                                                   // cell #1
+                                                   {1, 8, 9, 2, 5, 10, 11, 6}};
+        // Assign 8-noded hexahedron element to cell
+        std::shared_ptr<mpm::Element<Dim>> element =
+            Factory<mpm::Element<Dim>>::instance()->create("ED3H8");
+
+        // Global cell index
+        mpm::Index gcid = 0;
+        mesh->create_cells(gcid, element, cells);
+        // Check if mesh has added cells
+        REQUIRE(mesh->ncells() == cells.size());
+        // Try again this shouldn't add more cells
+        mesh->create_cells(gcid, element, cells);
+        // Check if mesh has added cells
+        REQUIRE(mesh->ncells() == cells.size());
+        // Clear cells and try creating a list of empty cells
+        unsigned ncells = cells.size();
+        cells.clear();
+        // This fails with empty list error in node creation
+        gcid = 100;
+        mesh->create_cells(gcid, element, cells);
+        REQUIRE(mesh->ncells() == ncells);
+
+        // Try with invalid node ids
+        cells = {// cell #0
+                 {90, 91, 92, 93, 94, 95, 96, 97},
+                 // cell #1
+                 {71, 88, 89, 82, 85, 80, 81, 86}};
+        gcid = 200;
+        mesh->create_cells(gcid, element, cells);
+        REQUIRE(mesh->ncells() == ncells);
+
+        SECTION("Check creation of levelset particles") {
+          // Vector of particle coordinates
+          std::vector<Eigen::Matrix<double, Dim, 1>> coordinates;
+          coordinates.clear();
+
+          // Particle coordinates
+          Eigen::Matrix<double, Dim, 1> particle;
+
+          // Cell 0
+          // Particle 0
+          particle << 0.125, 0.125, 0.125;
+          coordinates.emplace_back(particle);
+          // Particle 1
+          particle << 0.25, 0.125, 0.125;
+          coordinates.emplace_back(particle);
+          // Particle 2
+          particle << 0.25, 0.25, 0.125;
+          coordinates.emplace_back(particle);
+          // Particle 3
+          particle << 0.125, 0.25, 0.125;
+          coordinates.emplace_back(particle);
+          // Particle 4
+          particle << 0.125, 0.125, 0.25;
+          coordinates.emplace_back(particle);
+          // Particle 5
+          particle << 0.25, 0.125, 0.25;
+          coordinates.emplace_back(particle);
+          // Particle 6
+          particle << 0.25, 0.25, 0.25;
+          coordinates.emplace_back(particle);
+          // Particle 7
+          particle << 0.125, 0.25, 0.25;
+          coordinates.emplace_back(particle);
+
+          // Cell 1
+          // Particle 8
+          particle << 0.675, 0.125, 0.125;
+          coordinates.emplace_back(particle);
+          // Particle 9
+          particle << 0.85, 0.125, 0.125;
+          coordinates.emplace_back(particle);
+          // Particle 10
+          particle << 0.85, 0.25, 0.125;
+          coordinates.emplace_back(particle);
+          // Particle 11
+          particle << 0.675, 0.25, 0.125;
+          coordinates.emplace_back(particle);
+          // Particle 12
+          particle << 0.675, 0.125, 0.25;
+          coordinates.emplace_back(particle);
+          // Particle 13
+          particle << 0.85, 0.125, 0.25;
+          coordinates.emplace_back(particle);
+          // Particle 14
+          particle << 0.85, 0.25, 0.25;
+          coordinates.emplace_back(particle);
+          // Particle 15
+          particle << 0.675, 0.25, 0.25;
+          coordinates.emplace_back(particle);
+
+          // Initialise material models in mesh
+          mesh->initialise_material_models(materials);
+
+          SECTION("Check addition of levelset particles to levelset mesh") {
+            //! Particle type 3DLS
+            const std::string particle_type = "P3DLS";
+            // Create particles from file
+            mesh->create_particles(particle_type, coordinates, mids, 0, false);
+            // Check if mesh has added particles
+            REQUIRE(mesh->nparticles() == coordinates.size());
+          }
+        }
+
+        //! Boundary Conditions from Files
+        // Test assign levelset inputs to nodes from file
+        SECTION("Check assign levelset inputs to nodes") {
+          // Levelset inputs
+          std::vector<std::tuple<mpm::Index, double, double, double, double>>
+              levelset_inputs;
+
+          levelset_inputs.emplace_back(
+              std::make_tuple(0, 0.05, 0.1, 1000., 1.0E+06));
+          levelset_inputs.emplace_back(
+              std::make_tuple(1, 0.05, 0.1, 1000., 1.0E+06));
+          levelset_inputs.emplace_back(
+              std::make_tuple(2, 0.05, 0.1, 1000., 1.0E+06));
+          levelset_inputs.emplace_back(
+              std::make_tuple(3, 0.05, 0.1, 1000., 1.0E+06));
+
+          REQUIRE(mesh != nullptr);
+          REQUIRE(mesh->assign_nodal_levelset_values(levelset_inputs) == true);
+          // When inputs fail (negative mu and alpha)
+          levelset_inputs.emplace_back(
+              std::make_tuple(3, 0.05, -0.1, -1000., 1.0E+06));
+          REQUIRE(mesh->assign_nodal_levelset_values(levelset_inputs) == false);
+          // When inputs fail (negative barrier stiffness)
+          levelset_inputs.emplace_back(
+              std::make_tuple(3, 0.05, 0.1, 1000., -1.0E+06));
+          REQUIRE(mesh->assign_nodal_levelset_values(levelset_inputs) == false);
+        }
+      }
+    }
+  }
+
   //! Check if nodal properties is initialised
   SECTION("Check nodal properties initialisation") {
     // Create the different meshes
diff --git a/tests/nodes/node_test.cc b/tests/nodes/node_test.cc
index 90e12b631..509ed9825 100644
--- a/tests/nodes/node_test.cc
+++ b/tests/nodes/node_test.cc
@@ -9,7 +9,7 @@
 #include "geometry.h"
 #include "node.h"
 
-// Check node class for 1D case
+//! \brief Check node class for 1D case
 TEST_CASE("Node is checked for 1D case", "[node][1D]") {
   const unsigned Dim = 1;
   const unsigned Dof = 1;
@@ -493,7 +493,7 @@ TEST_CASE("Node is checked for 1D case", "[node][1D]") {
   }
 }
 
-// \brief Check node class for 2D case
+//! \brief Check node class for 2D case
 TEST_CASE("Node is checked for 2D case", "[node][2D]") {
   const unsigned Dim = 2;
   const unsigned Dof = 2;
@@ -1309,7 +1309,7 @@ TEST_CASE("Node is checked for 2D case", "[node][2D]") {
   }
 }
 
-// \brief Check node class for 3D case
+//! \brief Check node class for 3D case
 TEST_CASE("Node is checked for 3D case", "[node][3D]") {
   const unsigned Dim = 3;
   const unsigned Dof = 3;
diff --git a/tests/nodes/node_twophase_test.cc b/tests/nodes/node_twophase_test.cc
index b9b74a182..5d4914028 100644
--- a/tests/nodes/node_twophase_test.cc
+++ b/tests/nodes/node_twophase_test.cc
@@ -764,7 +764,7 @@ TEST_CASE("Twophase Node is checked for 1D case", "[node][1D][2Phase]") {
   }
 }
 
-// \brief Check node class for 2D case
+//! \brief Check node class for 2D case
 TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") {
   const unsigned Dim = 2;
   const unsigned Dof = 2;
@@ -1755,7 +1755,7 @@ TEST_CASE("Twophase Node is checked for 2D case", "[node][2D][2Phase]") {
   }
 }
 
-// \brief Check node class for 3D case
+//! \brief Check node class for 3D case
 TEST_CASE("Twophase Node is checked for 3D case", "[node][3D][2Phase]") {
   const unsigned Dim = 3;
   const unsigned Dof = 3;
diff --git a/tests/particles/particle_bbar_test.cc b/tests/particles/particle_bbar_test.cc
index ed74b4617..e0c32ae2c 100644
--- a/tests/particles/particle_bbar_test.cc
+++ b/tests/particles/particle_bbar_test.cc
@@ -10,6 +10,7 @@
 #include "material.h"
 #include "node.h"
 #include "particle_bbar.h"
+#include "particle_levelset_bbar.h"
 #include "pod_particle.h"
 #include "quadrilateral_element.h"
 
@@ -37,7 +38,8 @@ TEST_CASE("ParticleBbar is checked for 2D case", "[particle][2D][Bbar]") {
     // Add particle
     mpm::Index id = 0;
     coords << 0.75, 0.75;
-    auto particle = std::make_shared<mpm::ParticleBbar<Dim>>(id, coords);
+    std::shared_ptr<mpm::ParticleBase<Dim>> particle =
+        std::make_shared<mpm::ParticleBbar<Dim>>(id, coords);
 
     // Particle type
     REQUIRE(particle->type() == "P2DBBAR");
@@ -310,36 +312,37 @@ TEST_CASE("ParticleBbar is checked for 2D case", "[particle][2D][Bbar]") {
   }
 }
 
-//! \brief Check particle bbar class for 3D case
-TEST_CASE("ParticleBbar is checked for 3D case", "[particle][3D][Bbar]") {
+//! \brief Check particle levelset bbar class for 2D case
+TEST_CASE("ParticleLevelsetBbar is checked for 2D case",
+          "[particle][2D][Levelset][Bbar]") {
   // Dimension
-  const unsigned Dim = 3;
-  // Dimension
-  const unsigned Dof = 6;
+  const unsigned Dim = 2;
+  // Degree of freedom
+  const unsigned Dof = 2;
   // Number of nodes per cell
-  const unsigned Nnodes = 8;
+  const unsigned Nnodes = 4;
   // Number of phases
   const unsigned Nphases = 1;
+  // Phase
+  const unsigned phase = 0;
   // Tolerance
   const double Tolerance = 1.E-7;
 
   // Coordinates
-  Eigen::Vector3d coords;
+  Eigen::Vector2d coords;
   coords.setZero();
 
-  //! Test particle, cell and node functions
-  SECTION("Test particle, cell and node functions") {
+  //! Test levelset particle, cell and node functions
+  SECTION("Test levelset particle, cell and node functions") {
     // Add particle
     mpm::Index id = 0;
-    coords << 1.5, 1.5, 1.5;
-    std::shared_ptr<mpm::ParticleBase<Dim>> particle =
-        std::make_shared<mpm::ParticleBbar<Dim>>(id, coords);
+    coords << 0.75, 0.75;
+    auto particle =
+        std::make_shared<mpm::ParticleLevelsetBbar<Dim>>(id, coords);
 
     // Particle type
-    REQUIRE(particle->type() == "P3DBBAR");
+    REQUIRE(particle->type() == "P2DLSBBAR");
 
-    // Phase
-    const unsigned phase = 0;
     // Time-step
     const double dt = 0.1;
 
@@ -348,87 +351,47 @@ TEST_CASE("ParticleBbar is checked for 3D case", "[particle][3D][Bbar]") {
     for (unsigned i = 0; i < coordinates.size(); ++i)
       REQUIRE(coordinates(i) == Approx(coords(i)).epsilon(Tolerance));
 
-    // Assign hexahedron shape function
+    // Shape function
     std::shared_ptr<mpm::Element<Dim>> element =
-        std::make_shared<mpm::HexahedronElement<Dim, 8>>();
+        std::make_shared<mpm::QuadrilateralElement<Dim, 4>>();
 
     // Create cell
     auto cell = std::make_shared<mpm::Cell<Dim>>(10, Nnodes, element);
-    // Add nodes
-    coords << 0, 0, 0;
+    // Add nodes to cell
+    coords << 0.5, 0.5;
     std::shared_ptr<mpm::NodeBase<Dim>> node0 =
         std::make_shared<mpm::Node<Dim, Dof, Nphases>>(0, coords);
 
-    coords << 2, 0, 0;
+    coords << 1.5, 0.5;
     std::shared_ptr<mpm::NodeBase<Dim>> node1 =
         std::make_shared<mpm::Node<Dim, Dof, Nphases>>(1, coords);
 
-    coords << 2, 2, 0;
+    coords << 1.5, 1.5;
     std::shared_ptr<mpm::NodeBase<Dim>> node2 =
         std::make_shared<mpm::Node<Dim, Dof, Nphases>>(2, coords);
 
-    coords << 0, 2, 0;
+    coords << 0.5, 1.5;
     std::shared_ptr<mpm::NodeBase<Dim>> node3 =
         std::make_shared<mpm::Node<Dim, Dof, Nphases>>(3, coords);
-
-    coords << 0, 0, 2;
-    std::shared_ptr<mpm::NodeBase<Dim>> node4 =
-        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(4, coords);
-
-    coords << 2, 0, 2;
-    std::shared_ptr<mpm::NodeBase<Dim>> node5 =
-        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(5, coords);
-
-    coords << 2, 2, 2;
-    std::shared_ptr<mpm::NodeBase<Dim>> node6 =
-        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(6, coords);
-
-    coords << 0, 2, 2;
-    std::shared_ptr<mpm::NodeBase<Dim>> node7 =
-        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(7, coords);
+    cell->add_node(0, node0);
+    cell->add_node(1, node1);
+    cell->add_node(2, node2);
+    cell->add_node(3, node3);
+    REQUIRE(cell->nnodes() == 4);
 
     std::vector<std::shared_ptr<mpm::NodeBase<Dim>>> nodes;
     nodes.emplace_back(node0);
     nodes.emplace_back(node1);
     nodes.emplace_back(node2);
     nodes.emplace_back(node3);
-    nodes.emplace_back(node4);
-    nodes.emplace_back(node5);
-    nodes.emplace_back(node6);
-    nodes.emplace_back(node7);
-
-    cell->add_node(0, node0);
-    cell->add_node(1, node1);
-    cell->add_node(2, node2);
-    cell->add_node(3, node3);
-    cell->add_node(4, node4);
-    cell->add_node(5, node5);
-    cell->add_node(6, node6);
-    cell->add_node(7, node7);
-    REQUIRE(cell->nnodes() == 8);
 
     // Initialise cell properties
     cell->initialise();
 
-    // Check if cell is initialised
-    REQUIRE(cell->is_initialised() == true);
-
     // Add cell to particle
     REQUIRE(cell->status() == false);
-    // Check compute shape functions of a particle
-    // TODO Assert: REQUIRE(particle->compute_shapefn() == false);
     // Compute reference location should throw
     REQUIRE(particle->compute_reference_location() == false);
-    // Compute updated particle location should fail
-    // TODO Assert: REQUIRE(particle->compute_updated_position(dt) == false);
-    // Compute updated particle location from nodal velocity should fail
-    // TODO Assert: REQUIRE_NOTHROW(particle->compute_updated_position(dt,
-    // true));
-
-    // Compute volume
-    // TODO Assert: REQUIRE(particle->compute_volume() == false);
-    // Update volume should fail
-    // TODO Assert: REQUIRE(particle->update_volume() == false);
 
     REQUIRE(particle->assign_cell(cell) == true);
     REQUIRE(cell->status() == true);
@@ -449,17 +412,17 @@ TEST_CASE("ParticleBbar is checked for 3D case", "[particle][3D][Bbar]") {
     // Compute volume
     REQUIRE_NOTHROW(particle->compute_volume());
     // Check volume
-    REQUIRE(particle->volume() == Approx(8.0).epsilon(Tolerance));
+    REQUIRE(particle->volume() == Approx(1.0).epsilon(Tolerance));
 
     // Check reference location
-    coords << 0.5, 0.5, 0.5;
+    coords << -0.5, -0.5;
     REQUIRE(particle->compute_reference_location() == true);
     auto ref_coordinates = particle->reference_location();
     for (unsigned i = 0; i < ref_coordinates.size(); ++i)
       REQUIRE(ref_coordinates(i) == Approx(coords(i)).epsilon(Tolerance));
 
     // Assign material
-    unsigned mid = 0;
+    unsigned mid = 1;
     // Initialise material
     Json jmaterial;
     jmaterial["density"] = 1000.;
@@ -468,19 +431,13 @@ TEST_CASE("ParticleBbar is checked for 3D case", "[particle][3D][Bbar]") {
 
     auto material =
         Factory<mpm::Material<Dim>, unsigned, const Json&>::instance()->create(
-            "LinearElastic3D", std::move(mid), jmaterial);
-
-    // Check compute mass before material and volume
-    // TODO Assert: REQUIRE(particle->compute_mass() == false);
-
-    // Test compute stress before material assignment
-    // TODO Assert: REQUIRE(particle->compute_stress() == false);
+            "LinearElastic2D", std::move(mid), jmaterial);
 
     // Assign material properties
     REQUIRE(particle->assign_material(material) == true);
 
-    // Check material id from particle
-    REQUIRE(particle->material_id() == 0);
+    // Check material id
+    REQUIRE(particle->material_id() == 1);
 
     // Compute volume
     REQUIRE_NOTHROW(particle->compute_volume());
@@ -488,14 +445,665 @@ TEST_CASE("ParticleBbar is checked for 3D case", "[particle][3D][Bbar]") {
     // Compute mass
     REQUIRE_NOTHROW(particle->compute_mass());
     // Mass
-    REQUIRE(particle->mass() == Approx(8000.).epsilon(Tolerance));
+    REQUIRE(particle->mass() == Approx(1000.).epsilon(Tolerance));
+
+    // Map particle mass to nodes
+    particle->assign_mass(std::numeric_limits<double>::max());
+
+    // Assign mass to nodes
+    REQUIRE(particle->compute_reference_location() == true);
+    REQUIRE_NOTHROW(particle->compute_shapefn());
+
+    // Check velocity
+    Eigen::VectorXd velocity;
+    velocity.resize(Dim);
+    for (unsigned i = 0; i < velocity.size(); ++i) velocity(i) = i;
+    REQUIRE(particle->assign_velocity(velocity) == true);
+    for (unsigned i = 0; i < velocity.size(); ++i)
+      REQUIRE(particle->velocity()(i) == Approx(i).epsilon(Tolerance));
+
+    REQUIRE_NOTHROW(particle->compute_mass());
+    REQUIRE_NOTHROW(particle->map_mass_momentum_to_nodes());
+
+    REQUIRE(particle->compute_pressure_smoothing() == false);
+
+    // Values of nodal mass
+    std::array<double, 4> nodal_mass{562.5, 187.5, 62.5, 187.5};
+    // Check nodal mass
+    for (unsigned i = 0; i < nodes.size(); ++i)
+      REQUIRE(nodes.at(i)->mass(phase) ==
+              Approx(nodal_mass.at(i)).epsilon(Tolerance));
+
+    // Compute nodal velocity
+    for (const auto& node : nodes) node->compute_velocity();
+
+    // Values of nodal momentum
+    Eigen::Matrix<double, 4, 2> nodal_momentum;
+    // clang-format off
+    nodal_momentum << 0., 562.5,
+                      0., 187.5,
+                      0., 62.5,
+                      0., 187.5;
+    // clang-format on
+    // Check nodal momentum
+    for (unsigned i = 0; i < nodal_momentum.rows(); ++i)
+      for (unsigned j = 0; j < nodal_momentum.cols(); ++j)
+        REQUIRE(nodes.at(i)->momentum(phase)(j) ==
+                Approx(nodal_momentum(i, j)).epsilon(Tolerance));
+
+    // Values of nodal velocity
+    Eigen::Matrix<double, 4, 2> nodal_velocity;
+    // clang-format off
+    nodal_velocity << 0., 1.,
+                      0., 1.,
+                      0., 1.,
+                      0., 1.;
+    // clang-format on
+    // Check nodal velocity
+    for (unsigned i = 0; i < nodal_velocity.rows(); ++i)
+      for (unsigned j = 0; j < nodal_velocity.cols(); ++j)
+        REQUIRE(nodes.at(i)->velocity(phase)(j) ==
+                Approx(nodal_velocity(i, j)).epsilon(Tolerance));
+
+    // Set momentum to get non-zero strain
+    // clang-format off
+    nodal_momentum << 0., 562.5 * 1.,
+                      0., 187.5 * 2.,
+                      0.,  62.5 * 3.,
+                      0., 187.5 * 4.;
+    // clang-format on
+    for (unsigned i = 0; i < nodes.size(); ++i)
+      REQUIRE_NOTHROW(
+          nodes.at(i)->update_momentum(false, phase, nodal_momentum.row(i)));
+
+    // nodal velocity
+    // clang-format off
+    nodal_velocity << 0., 1.,
+                      0., 2.,
+                      0., 3.,
+                      0., 4.;
+    // clang-format on
+    // Compute nodal velocity
+    for (const auto& node : nodes) node->compute_velocity();
+    // Check nodal velocity
+    for (unsigned i = 0; i < nodal_velocity.rows(); ++i)
+      for (unsigned j = 0; j < nodal_velocity.cols(); ++j)
+        REQUIRE(nodes.at(i)->velocity(phase)(j) ==
+                Approx(nodal_velocity(i, j)).epsilon(Tolerance));
+
+    // Check pressure
+    REQUIRE(std::isnan(particle->pressure()) == true);
+
+    // Compute strain
+    particle->compute_strain(dt);
+    // Strain
+    Eigen::Matrix<double, 6, 1> strain;
+    strain << -0.025, 0.225, 0., 0.050, 0., 0.;
+    // Check strains
+    for (unsigned i = 0; i < strain.rows(); ++i)
+      REQUIRE(particle->strain()(i) == Approx(strain(i)).epsilon(Tolerance));
+
+    // Check updated pressure
+    REQUIRE(std::isnan(particle->pressure()) == true);
+
+    // Update volume strain rate
+    REQUIRE(particle->volume() == Approx(1.0).epsilon(Tolerance));
+    particle->compute_strain(dt);
+    REQUIRE_NOTHROW(particle->update_volume());
+    REQUIRE(particle->volume() == Approx(1.2).epsilon(Tolerance));
+
+    // Compute stress
+    REQUIRE_NOTHROW(particle->compute_stress(dt));
+
+    Eigen::Matrix<double, 6, 1> stress;
+    // clang-format off
+    stress <<  961538.4615384613,
+              2884615.3846153845,
+              1153846.1538461535,
+                192307.692307692,
+                    0.0000000000,
+                    0.0000000000;
+    // clang-format on
+    // Check stress
+    for (unsigned i = 0; i < stress.rows(); ++i)
+      REQUIRE(particle->stress()(i) == Approx(stress(i)).epsilon(Tolerance));
+
+    // Internal force
+    Eigen::Matrix<double, 4, 2> internal_force;
+    // clang-format off
+    internal_force <<  384615.3846153844,  1826923.0769230768,
+                      -192307.6923076921,  1057692.3076923077,
+                      -769230.7692307691,          -1250000.0,
+                       576923.0769230768, -1634615.3846153845;
+    // clang-format on
+
+    // Map particle internal force
+    particle->assign_volume(1.0);
+    particle->map_internal_force();
+
+    // Check nodal internal force
+    for (unsigned i = 0; i < internal_force.rows(); ++i)
+      for (unsigned j = 0; j < internal_force.cols(); ++j)
+        REQUIRE(nodes[i]->internal_force(phase)[j] ==
+                Approx(internal_force(i, j)).epsilon(Tolerance));
+  }
+}
+
+//! \brief Check particle bbar class for 3D case
+TEST_CASE("ParticleBbar is checked for 3D case", "[particle][3D][Bbar]") {
+  // Dimension
+  const unsigned Dim = 3;
+  // Dimension
+  const unsigned Dof = 6;
+  // Number of nodes per cell
+  const unsigned Nnodes = 8;
+  // Number of phases
+  const unsigned Nphases = 1;
+  // Tolerance
+  const double Tolerance = 1.E-7;
+
+  // Coordinates
+  Eigen::Vector3d coords;
+  coords.setZero();
+
+  //! Test particle, cell and node functions
+  SECTION("Test particle, cell and node functions") {
+    // Add particle
+    mpm::Index id = 0;
+    coords << 1.5, 1.5, 1.5;
+    std::shared_ptr<mpm::ParticleBase<Dim>> particle =
+        std::make_shared<mpm::ParticleBbar<Dim>>(id, coords);
+
+    // Particle type
+    REQUIRE(particle->type() == "P3DBBAR");
+
+    // Phase
+    const unsigned phase = 0;
+    // Time-step
+    const double dt = 0.1;
+
+    // Check particle coordinates
+    auto coordinates = particle->coordinates();
+    for (unsigned i = 0; i < coordinates.size(); ++i)
+      REQUIRE(coordinates(i) == Approx(coords(i)).epsilon(Tolerance));
+
+    // Assign hexahedron shape function
+    std::shared_ptr<mpm::Element<Dim>> element =
+        std::make_shared<mpm::HexahedronElement<Dim, 8>>();
+
+    // Create cell
+    auto cell = std::make_shared<mpm::Cell<Dim>>(10, Nnodes, element);
+    // Add nodes
+    coords << 0, 0, 0;
+    std::shared_ptr<mpm::NodeBase<Dim>> node0 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(0, coords);
+
+    coords << 2, 0, 0;
+    std::shared_ptr<mpm::NodeBase<Dim>> node1 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(1, coords);
+
+    coords << 2, 2, 0;
+    std::shared_ptr<mpm::NodeBase<Dim>> node2 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(2, coords);
+
+    coords << 0, 2, 0;
+    std::shared_ptr<mpm::NodeBase<Dim>> node3 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(3, coords);
+
+    coords << 0, 0, 2;
+    std::shared_ptr<mpm::NodeBase<Dim>> node4 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(4, coords);
+
+    coords << 2, 0, 2;
+    std::shared_ptr<mpm::NodeBase<Dim>> node5 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(5, coords);
+
+    coords << 2, 2, 2;
+    std::shared_ptr<mpm::NodeBase<Dim>> node6 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(6, coords);
+
+    coords << 0, 2, 2;
+    std::shared_ptr<mpm::NodeBase<Dim>> node7 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(7, coords);
+
+    std::vector<std::shared_ptr<mpm::NodeBase<Dim>>> nodes;
+    nodes.emplace_back(node0);
+    nodes.emplace_back(node1);
+    nodes.emplace_back(node2);
+    nodes.emplace_back(node3);
+    nodes.emplace_back(node4);
+    nodes.emplace_back(node5);
+    nodes.emplace_back(node6);
+    nodes.emplace_back(node7);
+
+    cell->add_node(0, node0);
+    cell->add_node(1, node1);
+    cell->add_node(2, node2);
+    cell->add_node(3, node3);
+    cell->add_node(4, node4);
+    cell->add_node(5, node5);
+    cell->add_node(6, node6);
+    cell->add_node(7, node7);
+    REQUIRE(cell->nnodes() == 8);
+
+    // Initialise cell properties
+    cell->initialise();
+
+    // Check if cell is initialised
+    REQUIRE(cell->is_initialised() == true);
+
+    // Add cell to particle
+    REQUIRE(cell->status() == false);
+    // Check compute shape functions of a particle
+    // TODO Assert: REQUIRE(particle->compute_shapefn() == false);
+    // Compute reference location should throw
+    REQUIRE(particle->compute_reference_location() == false);
+    // Compute updated particle location should fail
+    // TODO Assert: REQUIRE(particle->compute_updated_position(dt) == false);
+    // Compute updated particle location from nodal velocity should fail
+    // TODO Assert: REQUIRE_NOTHROW(particle->compute_updated_position(dt,
+    // true));
+
+    // Compute volume
+    // TODO Assert: REQUIRE(particle->compute_volume() == false);
+    // Update volume should fail
+    // TODO Assert: REQUIRE(particle->update_volume() == false);
+
+    REQUIRE(particle->assign_cell(cell) == true);
+    REQUIRE(cell->status() == true);
+    REQUIRE(particle->cell_id() == 10);
+
+    // Check if cell is initialised
+    REQUIRE(cell->is_initialised() == true);
+
+    // Check compute shape functions of a particle
+    REQUIRE_NOTHROW(particle->compute_shapefn());
+
+    // Assign volume
+    REQUIRE(particle->assign_volume(0.0) == false);
+    REQUIRE(particle->assign_volume(-5.0) == false);
+    REQUIRE(particle->assign_volume(2.0) == true);
+    // Check volume
+    REQUIRE(particle->volume() == Approx(2.0).epsilon(Tolerance));
+    // Compute volume
+    REQUIRE_NOTHROW(particle->compute_volume());
+    // Check volume
+    REQUIRE(particle->volume() == Approx(8.0).epsilon(Tolerance));
+
+    // Check reference location
+    coords << 0.5, 0.5, 0.5;
+    REQUIRE(particle->compute_reference_location() == true);
+    auto ref_coordinates = particle->reference_location();
+    for (unsigned i = 0; i < ref_coordinates.size(); ++i)
+      REQUIRE(ref_coordinates(i) == Approx(coords(i)).epsilon(Tolerance));
+
+    // Assign material
+    unsigned mid = 0;
+    // Initialise material
+    Json jmaterial;
+    jmaterial["density"] = 1000.;
+    jmaterial["youngs_modulus"] = 1.0E+7;
+    jmaterial["poisson_ratio"] = 0.3;
+
+    auto material =
+        Factory<mpm::Material<Dim>, unsigned, const Json&>::instance()->create(
+            "LinearElastic3D", std::move(mid), jmaterial);
+
+    // Check compute mass before material and volume
+    // TODO Assert: REQUIRE(particle->compute_mass() == false);
+
+    // Test compute stress before material assignment
+    // TODO Assert: REQUIRE(particle->compute_stress() == false);
+
+    // Assign material properties
+    REQUIRE(particle->assign_material(material) == true);
+
+    // Check material id from particle
+    REQUIRE(particle->material_id() == 0);
+
+    // Compute volume
+    REQUIRE_NOTHROW(particle->compute_volume());
+
+    // Compute mass
+    REQUIRE_NOTHROW(particle->compute_mass());
+    // Mass
+    REQUIRE(particle->mass() == Approx(8000.).epsilon(Tolerance));
+
+    // Map particle mass to nodes
+    particle->assign_mass(std::numeric_limits<double>::max());
+    // TODO Assert: REQUIRE(particle->map_mass_momentum_to_nodes() == false);
+
+    // Map particle pressure to nodes
+    // TODO Assert: REQUIRE(particle->map_pressure_to_nodes() == false);
+
+    // Assign mass to nodes
+    REQUIRE(particle->compute_reference_location() == true);
+    REQUIRE_NOTHROW(particle->compute_shapefn());
+
+    // Check velocity
+    Eigen::VectorXd velocity;
+    velocity.resize(Dim);
+    for (unsigned i = 0; i < velocity.size(); ++i) velocity(i) = i;
+    REQUIRE(particle->assign_velocity(velocity) == true);
+    for (unsigned i = 0; i < velocity.size(); ++i)
+      REQUIRE(particle->velocity()(i) == Approx(i).epsilon(Tolerance));
+
+    REQUIRE_NOTHROW(particle->compute_mass());
+    REQUIRE_NOTHROW(particle->map_mass_momentum_to_nodes());
+
+    // TODO Assert: REQUIRE(particle->map_pressure_to_nodes() == false);
+    REQUIRE(particle->compute_pressure_smoothing() == false);
+
+    // Values of nodal mass
+    std::array<double, 8> nodal_mass{125., 375.,  1125., 375.,
+                                     375., 1125., 3375., 1125.};
+    // Check nodal mass
+    for (unsigned i = 0; i < nodes.size(); ++i)
+      REQUIRE(nodes.at(i)->mass(phase) ==
+              Approx(nodal_mass.at(i)).epsilon(Tolerance));
+
+    // Compute nodal velocity
+    for (const auto& node : nodes) node->compute_velocity();
+
+    // Values of nodal momentum
+    Eigen::Matrix<double, 8, 3> nodal_momentum;
+    // clang-format off
+    nodal_momentum << 0.,  125.,  250.,
+                      0.,  375.,  750.,
+                      0., 1125., 2250.,
+                      0.,  375.,  750.,
+                      0.,  375.,  750.,
+                      0., 1125., 2250.,
+                      0., 3375., 6750.,
+                      0., 1125., 2250.;
+    // clang-format on
+
+    // Check nodal momentum
+    for (unsigned i = 0; i < nodal_momentum.rows(); ++i)
+      for (unsigned j = 0; j < nodal_momentum.cols(); ++j)
+        REQUIRE(nodes.at(i)->momentum(phase)[j] ==
+                Approx(nodal_momentum(i, j)).epsilon(Tolerance));
+
+    // Values of nodal velocity
+    Eigen::Matrix<double, 8, 3> nodal_velocity;
+    // clang-format off
+    nodal_velocity << 0., 1., 2.,
+                      0., 1., 2.,
+                      0., 1., 2.,
+                      0., 1., 2.,
+                      0., 1., 2.,
+                      0., 1., 2.,
+                      0., 1., 2.,
+                      0., 1., 2.;
+    // clang-format on
+    // Check nodal velocity
+    for (unsigned i = 0; i < nodal_velocity.rows(); ++i)
+      for (unsigned j = 0; j < nodal_velocity.cols(); ++j)
+        REQUIRE(nodes.at(i)->velocity(phase)(j) ==
+                Approx(nodal_velocity(i, j)).epsilon(Tolerance));
+
+    // Set momentum to get non-zero strain
+    // clang-format off
+    nodal_momentum << 0.,  125. * 1.,  250. * 1.,
+                      0.,  375. * 2.,  750. * 2.,
+                      0., 1125. * 3., 2250. * 3.,
+                      0.,  375. * 4.,  750. * 4.,
+                      0.,  375. * 5.,  750. * 5.,
+                      0., 1125. * 6., 2250. * 6.,
+                      0., 3375. * 7., 6750. * 7.,
+                      0., 1125. * 8., 2250. * 8.;
+    // clang-format on
+    for (unsigned i = 0; i < nodes.size(); ++i)
+      REQUIRE_NOTHROW(
+          nodes.at(i)->update_momentum(false, phase, nodal_momentum.row(i)));
+
+    // nodal velocity
+    // clang-format off
+    nodal_velocity << 0., 1.,  2.,
+                      0., 2.,  4.,
+                      0., 3.,  6.,
+                      0., 4.,  8.,
+                      0., 5., 10.,
+                      0., 6., 12.,
+                      0., 7., 14.,
+                      0., 8., 16.;
+    // clang-format on
+    // Compute nodal velocity
+    for (const auto& node : nodes) node->compute_velocity();
+    // Check nodal velocity
+    for (unsigned i = 0; i < nodal_velocity.rows(); ++i)
+      for (unsigned j = 0; j < nodal_velocity.cols(); ++j)
+        REQUIRE(nodes.at(i)->velocity(phase)(j) ==
+                Approx(nodal_velocity(i, j)).epsilon(Tolerance));
+
+    // Check pressure
+    REQUIRE(std::isnan(particle->pressure()) == true);
+
+    // Compute strain
+    particle->compute_strain(dt);
+    // Strain
+    Eigen::Matrix<double, 6, 1> strain;
+    strain << 0.0083333333, 0.0833333333, 0.4083333333, -0.02500, 0.35000,
+        -0.05000;
+
+    // Check strains
+    for (unsigned i = 0; i < strain.rows(); ++i)
+      REQUIRE(particle->strain()(i) == Approx(strain(i)).epsilon(Tolerance));
+
+    // Check updated pressure
+    REQUIRE(std::isnan(particle->pressure()) == true);
+
+    // Update volume strain rate
+    REQUIRE(particle->volume() == Approx(8.0).epsilon(Tolerance));
+    particle->compute_strain(dt);
+    REQUIRE_NOTHROW(particle->update_volume());
+    REQUIRE(particle->volume() == Approx(12.0).epsilon(Tolerance));
+
+    // Compute stress
+    REQUIRE_NOTHROW(particle->compute_stress(dt));
+    Eigen::Matrix<double, 6, 1> stress;
+    // clang-format off
+    stress << 2948717.9487179490,
+              3525641.0256410260,
+              6025641.0256410269,
+               -96153.8461538463,
+              1346153.8461538465,
+              -192307.6923076927;
+    // clang-format on
+    // Check stress
+    for (unsigned i = 0; i < stress.rows(); ++i)
+      REQUIRE(particle->stress()(i) == Approx(stress(i)).epsilon(Tolerance));
+
+    // Internal force
+    Eigen::Matrix<double, 8, 3> internal_force;
+    // clang-format off
+    internal_force <<  3790064.1025641034,   4318910.256410257,   4919871.7948717959,
+                      -4078525.6410256424,   4719551.282051282,   6618589.7435897458,
+                      -3613782.0512820524,  -584935.8974358966,   7483974.3589743618,
+                       3133012.8205128210, -3068910.2564102570,   5080128.2051282059,
+                       3229166.6666666670,  3277243.5897435895,  -3766025.6410256415,
+                      -3325320.5128205139,  1786858.9743589731,  -2387820.5128205139,
+                       -777243.5897435879, -5536858.9743589740, -10945512.8205128238,
+                       1642628.2051282048, -4911858.9743589750,  -7003205.1282051317;
+    // clang-format on
+
+    // Map particle internal force
+    particle->assign_volume(8.0);
+    particle->map_internal_force();
+
+    // Check nodal internal force
+    for (unsigned i = 0; i < internal_force.rows(); ++i)
+      for (unsigned j = 0; j < internal_force.cols(); ++j)
+        REQUIRE(nodes[i]->internal_force(phase)[j] ==
+                Approx(internal_force(i, j)).epsilon(Tolerance));
+  }
+}
+
+//! \brief Check particle levelset bbar class for 3D case
+TEST_CASE("ParticleLevelsetBbar is checked for 3D case",
+          "[particle][3D][Levelset][Bbar]") {
+  // Dimension
+  const unsigned Dim = 3;
+  // Dimension
+  const unsigned Dof = 6;
+  // Number of nodes per cell
+  const unsigned Nnodes = 8;
+  // Number of phases
+  const unsigned Nphases = 1;
+  // Tolerance
+  const double Tolerance = 1.E-7;
+
+  // Coordinates
+  Eigen::Vector3d coords;
+  coords.setZero();
+
+  //! Test levelset particle, cell and node functions
+  SECTION("Test levelset particle, cell and node functions") {
+    // Add particle
+    mpm::Index id = 0;
+    coords << 1.5, 1.5, 1.5;
+    std::shared_ptr<mpm::ParticleBase<Dim>> particle =
+        std::make_shared<mpm::ParticleLevelsetBbar<Dim>>(id, coords);
+
+    // Particle type
+    REQUIRE(particle->type() == "P3DLSBBAR");
+
+    // Phase
+    const unsigned phase = 0;
+    // Time-step
+    const double dt = 0.1;
+
+    // Check particle coordinates
+    auto coordinates = particle->coordinates();
+    for (unsigned i = 0; i < coordinates.size(); ++i)
+      REQUIRE(coordinates(i) == Approx(coords(i)).epsilon(Tolerance));
+
+    // Assign hexahedron shape function
+    std::shared_ptr<mpm::Element<Dim>> element =
+        std::make_shared<mpm::HexahedronElement<Dim, 8>>();
+
+    // Create cell
+    auto cell = std::make_shared<mpm::Cell<Dim>>(10, Nnodes, element);
+    // Add nodes
+    coords << 0, 0, 0;
+    std::shared_ptr<mpm::NodeBase<Dim>> node0 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(0, coords);
+
+    coords << 2, 0, 0;
+    std::shared_ptr<mpm::NodeBase<Dim>> node1 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(1, coords);
+
+    coords << 2, 2, 0;
+    std::shared_ptr<mpm::NodeBase<Dim>> node2 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(2, coords);
+
+    coords << 0, 2, 0;
+    std::shared_ptr<mpm::NodeBase<Dim>> node3 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(3, coords);
+
+    coords << 0, 0, 2;
+    std::shared_ptr<mpm::NodeBase<Dim>> node4 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(4, coords);
+
+    coords << 2, 0, 2;
+    std::shared_ptr<mpm::NodeBase<Dim>> node5 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(5, coords);
+
+    coords << 2, 2, 2;
+    std::shared_ptr<mpm::NodeBase<Dim>> node6 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(6, coords);
+
+    coords << 0, 2, 2;
+    std::shared_ptr<mpm::NodeBase<Dim>> node7 =
+        std::make_shared<mpm::Node<Dim, Dof, Nphases>>(7, coords);
+
+    std::vector<std::shared_ptr<mpm::NodeBase<Dim>>> nodes;
+    nodes.emplace_back(node0);
+    nodes.emplace_back(node1);
+    nodes.emplace_back(node2);
+    nodes.emplace_back(node3);
+    nodes.emplace_back(node4);
+    nodes.emplace_back(node5);
+    nodes.emplace_back(node6);
+    nodes.emplace_back(node7);
+
+    cell->add_node(0, node0);
+    cell->add_node(1, node1);
+    cell->add_node(2, node2);
+    cell->add_node(3, node3);
+    cell->add_node(4, node4);
+    cell->add_node(5, node5);
+    cell->add_node(6, node6);
+    cell->add_node(7, node7);
+    REQUIRE(cell->nnodes() == 8);
+
+    // Initialise cell properties
+    cell->initialise();
+
+    // Check if cell is initialised
+    REQUIRE(cell->is_initialised() == true);
+
+    // Add cell to particle
+    REQUIRE(cell->status() == false);
+    // Compute reference location should throw
+    REQUIRE(particle->compute_reference_location() == false);
+
+    REQUIRE(particle->assign_cell(cell) == true);
+    REQUIRE(cell->status() == true);
+    REQUIRE(particle->cell_id() == 10);
+
+    // Check if cell is initialised
+    REQUIRE(cell->is_initialised() == true);
+
+    // Check compute shape functions of a particle
+    REQUIRE_NOTHROW(particle->compute_shapefn());
+
+    // Assign volume
+    REQUIRE(particle->assign_volume(0.0) == false);
+    REQUIRE(particle->assign_volume(-5.0) == false);
+    REQUIRE(particle->assign_volume(2.0) == true);
+    // Check volume
+    REQUIRE(particle->volume() == Approx(2.0).epsilon(Tolerance));
+    // Compute volume
+    REQUIRE_NOTHROW(particle->compute_volume());
+    // Check volume
+    REQUIRE(particle->volume() == Approx(8.0).epsilon(Tolerance));
+
+    // Check reference location
+    coords << 0.5, 0.5, 0.5;
+    REQUIRE(particle->compute_reference_location() == true);
+    auto ref_coordinates = particle->reference_location();
+    for (unsigned i = 0; i < ref_coordinates.size(); ++i)
+      REQUIRE(ref_coordinates(i) == Approx(coords(i)).epsilon(Tolerance));
+
+    // Assign material
+    unsigned mid = 0;
+    // Initialise material
+    Json jmaterial;
+    jmaterial["density"] = 1000.;
+    jmaterial["youngs_modulus"] = 1.0E+7;
+    jmaterial["poisson_ratio"] = 0.3;
+
+    auto material =
+        Factory<mpm::Material<Dim>, unsigned, const Json&>::instance()->create(
+            "LinearElastic3D", std::move(mid), jmaterial);
+
+    // Assign material properties
+    REQUIRE(particle->assign_material(material) == true);
+
+    // Check material id from particle
+    REQUIRE(particle->material_id() == 0);
+
+    // Compute volume
+    REQUIRE_NOTHROW(particle->compute_volume());
+
+    // Compute mass
+    REQUIRE_NOTHROW(particle->compute_mass());
+    // Mass
+    REQUIRE(particle->mass() == Approx(8000.).epsilon(Tolerance));
 
     // Map particle mass to nodes
     particle->assign_mass(std::numeric_limits<double>::max());
-    // TODO Assert: REQUIRE(particle->map_mass_momentum_to_nodes() == false);
-
-    // Map particle pressure to nodes
-    // TODO Assert: REQUIRE(particle->map_pressure_to_nodes() == false);
 
     // Assign mass to nodes
     REQUIRE(particle->compute_reference_location() == true);
@@ -512,7 +1120,6 @@ TEST_CASE("ParticleBbar is checked for 3D case", "[particle][3D][Bbar]") {
     REQUIRE_NOTHROW(particle->compute_mass());
     REQUIRE_NOTHROW(particle->map_mass_momentum_to_nodes());
 
-    // TODO Assert: REQUIRE(particle->map_pressure_to_nodes() == false);
     REQUIRE(particle->compute_pressure_smoothing() == false);
 
     // Values of nodal mass
diff --git a/tests/particles/particle_test.cc b/tests/particles/particle_test.cc
index 9b7c914f3..18406ded9 100644
--- a/tests/particles/particle_test.cc
+++ b/tests/particles/particle_test.cc
@@ -10,7 +10,9 @@
 #include "material.h"
 #include "math_utility.h"
 #include "node.h"
+#include "node_levelset.h"
 #include "particle.h"
+#include "particle_levelset.h"
 #include "pod_particle.h"
 #include "quadrilateral_element.h"
 
@@ -189,6 +191,34 @@ TEST_CASE("Particle is checked for 1D case", "[particle][1D]") {
       REQUIRE(std::isnan(particle->tensor_data("invalid")(i)) == true);
   }
 
+  //! Test levelset particles scalar, vector and tensor data
+  SECTION("Check levelset particle scalar, vector, and tensor data") {
+    mpm::Index id = 0;
+    const double Tolerance = 1.E-7;
+    bool status = true;
+
+    //! Create levelset particle
+    std::shared_ptr<mpm::ParticleBase<Dim>> base_particle =
+        std::make_shared<mpm::ParticleLevelset<Dim>>(id, coords, status);
+    auto particle =
+        std::dynamic_pointer_cast<mpm::ParticleLevelset<Dim>>(base_particle);
+
+    particle->initialise();
+
+    // Check scalar data: levelset
+    REQUIRE(particle->scalar_data("levelset") ==
+            Approx(0.0).epsilon(Tolerance));
+
+    // Check vector data: levelset_reaction_forces
+    Eigen::VectorXd levelset_reaction;
+    levelset_reaction.resize(Dim);
+
+    REQUIRE(particle->vector_data("levelset_reaction_forces").size() == Dim);
+    for (unsigned i = 0; i < levelset_reaction.size(); ++i)
+      REQUIRE(particle->vector_data("levelset_reaction_forces")(i) ==
+              Approx(0.).epsilon(Tolerance));
+  }
+
   //! Test particles velocity constraints
   SECTION("Particle with velocity constraints") {
     mpm::Index id = 0;
@@ -867,6 +897,122 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") {
     REQUIRE(particle->velocity()(1) == Approx(-12.5).epsilon(Tolerance));
   }
 
+  //! Test particles levelset
+  SECTION("Particle with levelset") {
+    const double Tolerance = 1.E-7;
+    bool status = true;
+
+    //! Levelset static variables
+    mpm::ParticleLevelset<Dim>::set_levelset_properties(0.05, false, 0.01);
+
+    // Particle ID and Coordinates
+    mpm::Index id1 = 0;
+    mpm::Index id2 = 1;
+    Eigen::Vector2d coords1;
+    Eigen::Vector2d coords2;
+    coords1 << 0.5, 0.2;
+    coords2 << 0.5, 0.4;
+
+    //! Create levelset particle
+    std::shared_ptr<mpm::ParticleLevelset<Dim>> particle1 =
+        std::make_shared<mpm::ParticleLevelset<Dim>>(id1, coords1, status);
+    std::shared_ptr<mpm::ParticleLevelset<Dim>> particle2 =
+        std::make_shared<mpm::ParticleLevelset<Dim>>(id2, coords2, status);
+
+    // Initialise particle
+    particle1->initialise();
+    particle2->initialise();
+
+    // Assign particle properties
+    Eigen::Vector2d velocity;
+    velocity << 1., 0.1;
+    REQUIRE(particle1->assign_velocity(velocity) == true);
+    REQUIRE(particle1->assign_volume(0.25) == true);
+    particle1->assign_mass(2000.);
+
+    // Check particle type
+    REQUIRE(particle1->type() == "P2DLS");
+    REQUIRE(particle2->type() == "P2DLS");
+
+    // Check initial levelset value
+    REQUIRE(particle1->levelset() == Approx(0.0).epsilon(Tolerance));
+    REQUIRE(particle2->levelset() == Approx(0.0).epsilon(Tolerance));
+
+    // Check initial reaction force
+    auto reaction_force1 = particle1->reaction_force();
+    auto reaction_force2 = particle2->reaction_force();
+    REQUIRE(reaction_force1.size() == Dim);
+    REQUIRE(reaction_force2.size() == Dim);
+    for (unsigned i = 0; i < Dim; ++i) {
+      REQUIRE(reaction_force1(i) == Approx(0.0).epsilon(Tolerance));
+      REQUIRE(reaction_force2(i) == Approx(0.0).epsilon(Tolerance));
+    }
+
+    // Create cell
+    const unsigned Nnodes = 4;
+    auto element = std::make_shared<mpm::QuadrilateralElement<Dim, Nnodes>>();
+    auto cell = std::make_shared<mpm::Cell<Dim>>(10, Nnodes, element);
+
+    // Add nodes to the cell
+    Eigen::Vector2d node_coords;
+    std::vector<std::shared_ptr<mpm::NodeBase<Dim>>> nodes;
+    node_coords << 0.0, 0.0;
+    nodes.emplace_back(
+        std::make_shared<mpm::NodeLevelset<Dim, Dim, 1>>(0, node_coords));
+    node_coords << 1.0, 0.0;
+    nodes.emplace_back(
+        std::make_shared<mpm::NodeLevelset<Dim, Dim, 1>>(1, node_coords));
+    node_coords << 1.0, 1.0;
+    nodes.emplace_back(
+        std::make_shared<mpm::NodeLevelset<Dim, Dim, 1>>(2, node_coords));
+    node_coords << 0.0, 1.0;
+    nodes.emplace_back(
+        std::make_shared<mpm::NodeLevelset<Dim, Dim, 1>>(3, node_coords));
+
+    for (unsigned i = 0; i < nodes.size(); ++i) cell->add_node(i, nodes[i]);
+
+    // Initialise cell
+    cell->initialise();
+
+    // Assign the cell to the particle
+    REQUIRE(particle1->assign_cell(cell) == true);
+    REQUIRE(particle2->assign_cell(cell) == true);
+
+    // Assign levelset values to nodes
+    REQUIRE(nodes[0]->assign_levelset(0.0, 0.1, 1000, 1E+06) == true);
+    REQUIRE(nodes[1]->assign_levelset(0.0, 0.1, 1000, 1E+06) == true);
+    REQUIRE(nodes[2]->assign_levelset(1.0, 0.1, 1000, 1E+06) == true);
+    REQUIRE(nodes[3]->assign_levelset(1.0, 0.1, 1000, 1E+06) == true);
+
+    // Compute shape functions and map levelset to particle
+    particle1->compute_shapefn();
+    particle2->compute_shapefn();
+    REQUIRE_NOTHROW(particle1->levelset_contact_force(0.1));
+    REQUIRE_NOTHROW(particle2->levelset_contact_force(0.1));
+
+    // Check mapped levelset value
+    REQUIRE(particle1->levelset() == Approx(0.20).epsilon(Tolerance));
+    REQUIRE(particle2->levelset() == Approx(0.40).epsilon(Tolerance));
+
+    // Check reaction force for particle not in contact
+    reaction_force2 = particle2->reaction_force();
+    for (unsigned i = 0; i < Dim; ++i) {
+      REQUIRE(reaction_force2(i) == Approx(0.0).epsilon(Tolerance));
+    }
+
+    // Check reaction force after contact reaction force computation
+    reaction_force1 = particle1->reaction_force();
+    REQUIRE(reaction_force1(0) == Approx(-9040.8465928880).epsilon(Tolerance));
+    REQUIRE(reaction_force1(1) == Approx(85658.4659288802).epsilon(Tolerance));
+
+    // Check mapped nodal contact reaction force
+    for (const auto& node : nodes) {
+      auto external_force = node->external_force(0);
+      REQUIRE(external_force(0) != Approx(0.0).epsilon(Tolerance));
+      REQUIRE(external_force(1) != Approx(0.0).epsilon(Tolerance));
+    }
+  }
+
   //! Test particle, cell and node functions
   SECTION("Test particle, cell and node functions") {
     // Add particle
@@ -2326,6 +2472,140 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") {
     REQUIRE(particle->velocity()(2) == Approx(14.5).epsilon(Tolerance));
   }
 
+  //! Test particles levelset
+  SECTION("Particle with levelset") {
+    const double Tolerance = 1.E-7;
+    bool status = true;
+
+    //! Levelset static variables
+    mpm::ParticleLevelset<Dim>::set_levelset_properties(0.05, false, 0.01);
+
+    // Particle ID and Coordinates
+    mpm::Index id1 = 0;
+    mpm::Index id2 = 1;
+    Eigen::Vector3d coords1;
+    Eigen::Vector3d coords2;
+    coords1 << 0.5, 0.2, 0.5;
+    coords2 << 0.5, 0.4, 0.5;
+
+    //! Create levelset particle
+    std::shared_ptr<mpm::ParticleLevelset<Dim>> particle1 =
+        std::make_shared<mpm::ParticleLevelset<Dim>>(id1, coords1, status);
+    std::shared_ptr<mpm::ParticleLevelset<Dim>> particle2 =
+        std::make_shared<mpm::ParticleLevelset<Dim>>(id2, coords2, status);
+
+    // Initialise particle
+    particle1->initialise();
+    particle2->initialise();
+
+    // Assign particle properties
+    Eigen::Vector3d velocity;
+    velocity << 1., 0.1, 1.;
+    REQUIRE(particle1->assign_velocity(velocity) == true);
+    REQUIRE(particle1->assign_volume(0.4) == true);
+    particle1->assign_mass(2000.);
+
+    // Check particle type
+    REQUIRE(particle1->type() == "P3DLS");
+    REQUIRE(particle2->type() == "P3DLS");
+
+    // Check initial levelset value
+    REQUIRE(particle1->levelset() == Approx(0.0).epsilon(Tolerance));
+    REQUIRE(particle2->levelset() == Approx(0.0).epsilon(Tolerance));
+
+    // Check initial reaction force
+    auto reaction_force1 = particle1->reaction_force();
+    auto reaction_force2 = particle2->reaction_force();
+    REQUIRE(reaction_force1.size() == Dim);
+    REQUIRE(reaction_force2.size() == Dim);
+    for (unsigned i = 0; i < Dim; ++i) {
+      REQUIRE(reaction_force1(i) == Approx(0.0).epsilon(Tolerance));
+      REQUIRE(reaction_force2(i) == Approx(0.0).epsilon(Tolerance));
+    }
+
+    // Create cell
+    const unsigned Nnodes = 8;
+    auto element = std::make_shared<mpm::HexahedronElement<Dim, Nnodes>>();
+    auto cell = std::make_shared<mpm::Cell<Dim>>(10, Nnodes, element);
+
+    // Add nodes to the cell
+    Eigen::Vector3d node_coords;
+    std::vector<std::shared_ptr<mpm::NodeBase<Dim>>> nodes;
+    node_coords << 0.0, 0.0, 0.0;
+    nodes.emplace_back(
+        std::make_shared<mpm::NodeLevelset<Dim, Dim, 1>>(0, node_coords));
+    node_coords << 1.0, 0.0, 0.0;
+    nodes.emplace_back(
+        std::make_shared<mpm::NodeLevelset<Dim, Dim, 1>>(1, node_coords));
+    node_coords << 1.0, 1.0, 0.0;
+    nodes.emplace_back(
+        std::make_shared<mpm::NodeLevelset<Dim, Dim, 1>>(2, node_coords));
+    node_coords << 0.0, 1.0, 0.0;
+    nodes.emplace_back(
+        std::make_shared<mpm::NodeLevelset<Dim, Dim, 1>>(3, node_coords));
+    node_coords << 0.0, 0.0, 1.0;
+    nodes.emplace_back(
+        std::make_shared<mpm::NodeLevelset<Dim, Dim, 1>>(4, node_coords));
+    node_coords << 1.0, 0.0, 1.0;
+    nodes.emplace_back(
+        std::make_shared<mpm::NodeLevelset<Dim, Dim, 1>>(5, node_coords));
+    node_coords << 1.0, 1.0, 1.0;
+    nodes.emplace_back(
+        std::make_shared<mpm::NodeLevelset<Dim, Dim, 1>>(6, node_coords));
+    node_coords << 0.0, 1.0, 1.0;
+    nodes.emplace_back(
+        std::make_shared<mpm::NodeLevelset<Dim, Dim, 1>>(7, node_coords));
+
+    for (unsigned i = 0; i < nodes.size(); ++i) cell->add_node(i, nodes[i]);
+
+    // Initialise cell
+    cell->initialise();
+
+    // Assign the cell to the particle
+    REQUIRE(particle1->assign_cell(cell) == true);
+    REQUIRE(particle2->assign_cell(cell) == true);
+
+    // Assign levelset values to nodes
+    REQUIRE(nodes[0]->assign_levelset(0.0, 0.1, 1000, 1E+05) == true);
+    REQUIRE(nodes[1]->assign_levelset(0.0, 0.1, 1000, 1E+05) == true);
+    REQUIRE(nodes[2]->assign_levelset(1.0, 0.1, 1000, 1E+05) == true);
+    REQUIRE(nodes[3]->assign_levelset(1.0, 0.1, 1000, 1E+05) == true);
+    REQUIRE(nodes[4]->assign_levelset(0.0, 0.1, 1000, 1E+05) == true);
+    REQUIRE(nodes[5]->assign_levelset(0.0, 0.1, 1000, 1E+05) == true);
+    REQUIRE(nodes[6]->assign_levelset(1.0, 0.1, 1000, 1E+05) == true);
+    REQUIRE(nodes[7]->assign_levelset(1.0, 0.1, 1000, 1E+05) == true);
+
+    // Compute shape functions and map levelset to particle
+    particle1->compute_shapefn();
+    particle2->compute_shapefn();
+    REQUIRE_NOTHROW(particle1->levelset_contact_force(0.1));
+    REQUIRE_NOTHROW(particle2->levelset_contact_force(0.1));
+
+    // Check mapped levelset value
+    REQUIRE(particle1->levelset() == Approx(0.20).epsilon(Tolerance));
+    REQUIRE(particle2->levelset() == Approx(0.40).epsilon(Tolerance));
+
+    // Check reaction force for particle not in contact
+    reaction_force2 = particle2->reaction_force();
+    for (unsigned i = 0; i < Dim; ++i) {
+      REQUIRE(reaction_force2(i) == Approx(0.0).epsilon(Tolerance));
+    }
+
+    // Check reaction force after contact reaction force computation
+    reaction_force1 = particle1->reaction_force();
+    REQUIRE(reaction_force1(0) == Approx(-5569.4624469601).epsilon(Tolerance));
+    REQUIRE(reaction_force1(1) == Approx(71764.4334287699).epsilon(Tolerance));
+    REQUIRE(reaction_force1(2) == Approx(-5569.4624469601).epsilon(Tolerance));
+
+    // Check mapped nodal contact reaction force
+    for (const auto& node : nodes) {
+      auto external_force = node->external_force(0);
+      REQUIRE(external_force(0) != Approx(0.0).epsilon(Tolerance));
+      REQUIRE(external_force(1) != Approx(0.0).epsilon(Tolerance));
+      REQUIRE(external_force(2) != Approx(0.0).epsilon(Tolerance));
+    }
+  }
+
   //! Test particle, cell and node functions
   SECTION("Test particle, cell and node functions") {
     // Add particle
diff --git a/tests/solvers/mpm_explicit_levelset_test.cc b/tests/solvers/mpm_explicit_levelset_test.cc
new file mode 100644
index 000000000..fc63f8584
--- /dev/null
+++ b/tests/solvers/mpm_explicit_levelset_test.cc
@@ -0,0 +1,139 @@
+#include "catch.hpp"
+
+//! Alias for JSON
+#include "json.hpp"
+using Json = nlohmann::json;
+
+#include "contact_levelset.h"
+#include "mpm_explicit.h"
+#include "particle_levelset.h"
+#include "write_mesh_particles.h"
+
+//! \brief Check MPM Explicit Levelset Interface (2D)
+TEST_CASE("MPM 2D Explicit Levelset Interface is checked",
+          "[MPM][2D][Explicit][Levelset]") {
+  // Dimension
+  const unsigned Dim = 2;
+  // Tolerance
+  const double Tolerance = 1.E-7;
+
+  // Write JSON file
+  const std::string fname = "mpm-levelset-test";
+  const std::string analysis = "MPMExplicit2D";
+
+  // Write Mesh with levelset boundaries
+  REQUIRE(mpm_test::write_mesh_2d() == true);
+
+  // Write Particles
+  REQUIRE(mpm_test::write_particles_2d() == true);
+
+  // Assign argc and argv to input arguments of MPM
+  int argc = 5;
+  // clang-format off
+  char* argv[] = {(char*)"./mpm",
+                  (char*)"-f",  (char*)"./",
+                  (char*)"-i",  (char*)"mpm-levelset-test-2d.json"};
+  // clang-format on
+
+  SECTION("Check initialisation") {
+    // Write custom JSON for levelset test with interface_type: "levelset"
+    REQUIRE(mpm_test::write_json_levelset(Dim, false, analysis, fname, false,
+                                          false) == true);
+    REQUIRE(mpm_test::write_json_levelset(Dim, false, analysis, fname, true,
+                                          true) == true);
+
+    // Check JSON information
+    std::ifstream json_file("mpm-levelset-test-2d.json");
+    REQUIRE(json_file.is_open());
+
+    Json json_data;
+    json_file >> json_data;
+    json_file.close();
+
+    REQUIRE(json_data.contains("mesh"));
+    REQUIRE(json_data["mesh"].contains("interface"));
+    const auto& interface = json_data["mesh"]["interface"];
+    REQUIRE(interface.contains("interface_type"));
+    REQUIRE(interface["interface_type"] == "levelset");
+    REQUIRE(interface.contains("location"));
+    REQUIRE(interface["location"] == "levelset-nodal-values-2d.txt");
+    REQUIRE(interface.contains("damping"));
+    REQUIRE(interface["damping"] == Approx(0.05).epsilon(Tolerance));
+    REQUIRE(interface.contains("velocity_update"));
+    REQUIRE(interface["velocity_update"] == "pic");
+    REQUIRE(interface.contains("violation_corrector"));
+    REQUIRE(interface["violation_corrector"] ==
+            Approx(0.001).epsilon(Tolerance));
+
+    // Create an IO object and run explicit MPM
+    auto io = std::make_unique<mpm::IO>(argc, argv);
+    auto mpm = std::make_unique<mpm::MPMExplicit<Dim>>(std::move(io));
+
+    // Initialise and check that the solver runs with levelset boundary
+    REQUIRE_NOTHROW(mpm->initialise_materials());
+    REQUIRE_NOTHROW(mpm->initialise_mesh());
+    REQUIRE_NOTHROW(mpm->initialise_particles());
+    REQUIRE_NOTHROW(mpm->initialise_loads());
+  }
+
+  SECTION("Check solver") {
+    // Create an IO object and run explicit MPM
+    auto io = std::make_unique<mpm::IO>(argc, argv);
+    auto mpm = std::make_unique<mpm::MPMExplicit<Dim>>(std::move(io));
+
+    // Check that the solver runs without throwing
+    REQUIRE(mpm->solve() == true);
+  }
+}
+
+//! \brief Check MPM Explicit Levelset Interface (3D)
+TEST_CASE("MPM 3D Explicit Levelset Interface is checked",
+          "[MPM][3D][Explicit][Levelset]") {
+  // Dimension
+  const unsigned Dim = 3;
+
+  // Write JSON file
+  const std::string fname = "mpm-levelset-test";
+  const std::string analysis = "MPMExplicit3D";
+
+  // Write Mesh with levelset boundaries
+  REQUIRE(mpm_test::write_mesh_3d() == true);
+
+  // Write Particles
+  REQUIRE(mpm_test::write_particles_3d() == true);
+
+  // Assign argc and argv to input arguments of MPM
+  int argc = 5;
+  // clang-format off
+  char* argv[] = {(char*)"./mpm",
+            (char*)"-f",  (char*)"./",
+            (char*)"-i",  (char*)"mpm-levelset-test-3d.json"};
+  // clang-format on
+
+  SECTION("Check initialisation and solver") {
+    // Write custom JSON for levelset test with interface_type: "levelset"
+    REQUIRE(mpm_test::write_json_levelset(Dim, false, analysis, fname, false,
+                                          false) == true);
+    REQUIRE(mpm_test::write_json_levelset(Dim, false, analysis, fname, true,
+                                          true) == true);
+
+    // Create an IO object and run explicit MPM
+    auto io = std::make_unique<mpm::IO>(argc, argv);
+    auto mpm = std::make_unique<mpm::MPMExplicit<Dim>>(std::move(io));
+
+    // Initialise and check that the solver runs with levelset boundary
+    REQUIRE_NOTHROW(mpm->initialise_materials());
+    REQUIRE_NOTHROW(mpm->initialise_mesh());
+    REQUIRE_NOTHROW(mpm->initialise_particles());
+    REQUIRE_NOTHROW(mpm->initialise_loads());
+  }
+
+  SECTION("Check solver") {
+    // Create an IO object and run explicit MPM
+    auto io = std::make_unique<mpm::IO>(argc, argv);
+    auto mpm = std::make_unique<mpm::MPMExplicit<Dim>>(std::move(io));
+
+    // Check that the solver runs without throwing
+    REQUIRE(mpm->solve() == true);
+  }
+}