diff --git a/CMakeLists.txt b/CMakeLists.txt
index d1a6dfa..4639fe8 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -6,8 +6,8 @@ project(Discregrid)
 # Visual studio solution directories.
 set_property(GLOBAL PROPERTY USE_FOLDERS on)
 
-# Require C++11 compiler
-set(CMAKE_CXX_STANDARD 11)
+# Require C++20 compiler
+set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
diff --git a/discregrid/include/Discregrid/mesh/entity_iterators.hpp b/discregrid/include/Discregrid/mesh/entity_iterators.hpp
index 3c736a8..a14ea5e 100755
--- a/discregrid/include/Discregrid/mesh/entity_iterators.hpp
+++ b/discregrid/include/Discregrid/mesh/entity_iterators.hpp
@@ -2,263 +2,289 @@
 
 #include "halfedge.hpp"
 
-#include <iterator>
-#include <array>
 #include <Eigen/Core>
+#include <array>
+#include <iterator>
 
-namespace Discregrid
-{
+namespace Discregrid {
 
 class TriangleMesh;
 
 class FaceContainer;
-class FaceIterator : public
-	std::iterator<std::random_access_iterator_tag, std::array<unsigned int, 3>>
-{
+class FaceIterator {
 public:
 
-	typedef FaceIterator _Mytype;
-
-	FaceIterator() = delete;
-
-	reference operator*();
-
-	bool operator<(_Mytype const& other) const
-	{
-		return m_index < other.m_index;
-	}
-	bool operator==(_Mytype const& other) const
-	{
-		return m_index == other.m_index;
-	}
-
-	bool operator!=(_Mytype const& other) const
-	{
-		return !(*this == other);
-	}
-
-	inline _Mytype& operator++() { ++m_index; return *this; }
-	inline _Mytype& operator--() { --m_index; return *this; }
-
-	inline _Mytype operator+(_Mytype const& rhs)
-	{
-		return _Mytype(m_index + rhs.m_index, m_mesh);
-	}
-	inline difference_type operator-(_Mytype const& rhs)
-	{
-		return m_index - rhs.m_index;
-	}
-	inline _Mytype operator-(int const& rhs)
-	{
-		return _Mytype(m_index - rhs, m_mesh);
-	}
-
-	unsigned int vertex(unsigned int i) const;
-	unsigned int& vertex(unsigned int i);
+  // ITERATOR TRAITS -------------------------------
+  /// One of the @link iterator_tags tag types@endlink.
+  typedef std::random_access_iterator_tag iterator_category;
+  /// The type "pointed to" by the iterator.
+  typedef std::array<unsigned int, 3> value_type;
+  /// Distance between iterators is represented as this type.
+  typedef ptrdiff_t difference_type;
+  /// This type represents a pointer-to-value_type.
+  typedef std::array<unsigned int, 3> *pointer;
+  /// This type represents a reference-to-value_type.
+  typedef std::array<unsigned int, 3> &reference;
+  // -----------------------------------------------
+
+  typedef FaceIterator _Mytype;
+
+  FaceIterator() = delete;
+
+  reference operator*();
+
+  bool operator<(_Mytype const &other) const { return m_index < other.m_index; }
+  bool operator==(_Mytype const &other) const {
+    return m_index == other.m_index;
+  }
+
+  bool operator!=(_Mytype const &other) const { return !(*this == other); }
+
+  inline _Mytype &operator++() {
+    ++m_index;
+    return *this;
+  }
+  inline _Mytype &operator--() {
+    --m_index;
+    return *this;
+  }
+
+  inline _Mytype operator+(_Mytype const &rhs) {
+    return _Mytype(m_index + rhs.m_index, m_mesh);
+  }
+  inline difference_type operator-(_Mytype const &rhs) {
+    return m_index - rhs.m_index;
+  }
+  inline _Mytype operator-(int const &rhs) {
+    return _Mytype(m_index - rhs, m_mesh);
+  }
+
+  unsigned int vertex(unsigned int i) const;
+  unsigned int &vertex(unsigned int i);
 
 private:
+  friend class FaceContainer;
+  FaceIterator(unsigned int index, TriangleMesh *mesh)
+      : m_index(index), m_mesh(mesh) {}
 
-	friend class FaceContainer;
-	FaceIterator(unsigned int index, TriangleMesh* mesh)
-		: m_index(index), m_mesh(mesh) {}
-
-	unsigned int m_index;
-	TriangleMesh* m_mesh;
+  unsigned int m_index;
+  TriangleMesh *m_mesh;
 };
-class FaceConstIterator : public
-	std::iterator<std::random_access_iterator_tag, std::array<unsigned int, 3> const>
-{
-
+class FaceConstIterator {
 public:
 
-	typedef FaceConstIterator _Mytype;
-
-	FaceConstIterator() = delete;
-
-	reference operator*();
-
-	bool operator<(_Mytype const& other) const
-	{
-		return m_index < other.m_index;
-	}
-	bool operator==(_Mytype const& other) const
-	{
-		return m_index == other.m_index;
-	}
-
-	bool operator!=(_Mytype const& other) const
-	{
-		return !(*this == other);
-	}
-
-	inline _Mytype& operator++() { ++m_index; return *this; }
-	inline _Mytype& operator--() { --m_index; return *this; }
-
-	inline _Mytype operator+(_Mytype const& rhs) const
-	{
-		return _Mytype(m_index + rhs.m_index, m_mesh);
-	}
-	inline difference_type operator-(_Mytype const& rhs) const
-	{
-		return m_index - rhs.m_index;
-	}
-	inline _Mytype operator-(int const& rhs) const
-	{
-		return _Mytype(m_index - rhs, m_mesh);
-	}
-
-	unsigned int vertex(unsigned int i) const;
-	unsigned int& vertex(unsigned int i);
+  // ITERATOR TRAITS -------------------------------
+  /// One of the @link iterator_tags tag types@endlink.
+  typedef std::random_access_iterator_tag iterator_category;
+  /// The type "pointed to" by the iterator.
+  typedef std::array<unsigned int, 3> const value_type;
+  /// Distance between iterators is represented as this type.
+  typedef ptrdiff_t difference_type;
+  /// This type represents a pointer-to-value_type.
+  typedef std::array<unsigned int, 3> const *pointer;
+  /// This type represents a reference-to-value_type.
+  typedef std::array<unsigned int, 3> const &reference;
+  // -----------------------------------------------
+
+  typedef FaceConstIterator _Mytype;
+
+  FaceConstIterator() = delete;
+
+  reference operator*();
+
+  bool operator<(_Mytype const &other) const { return m_index < other.m_index; }
+  bool operator==(_Mytype const &other) const {
+    return m_index == other.m_index;
+  }
+
+  bool operator!=(_Mytype const &other) const { return !(*this == other); }
+
+  inline _Mytype &operator++() {
+    ++m_index;
+    return *this;
+  }
+  inline _Mytype &operator--() {
+    --m_index;
+    return *this;
+  }
+
+  inline _Mytype operator+(_Mytype const &rhs) const {
+    return _Mytype(m_index + rhs.m_index, m_mesh);
+  }
+  inline difference_type operator-(_Mytype const &rhs) const {
+    return m_index - rhs.m_index;
+  }
+  inline _Mytype operator-(int const &rhs) const {
+    return _Mytype(m_index - rhs, m_mesh);
+  }
+
+  unsigned int vertex(unsigned int i) const;
+  unsigned int &vertex(unsigned int i);
 
 private:
+  friend class FaceConstContainer;
+  FaceConstIterator(unsigned int index, TriangleMesh const *mesh)
+      : m_index(index), m_mesh(mesh) {}
 
-	friend class FaceConstContainer;
-	FaceConstIterator(unsigned int index, TriangleMesh const* mesh)
-		: m_index(index), m_mesh(mesh) {}
-
-	unsigned int m_index;
-	TriangleMesh const* m_mesh;
+  unsigned int m_index;
+  TriangleMesh const *m_mesh;
 };
 
 class IncidentFaceContainer;
-class IncidentFaceIterator : public std::iterator<std::forward_iterator_tag, Halfedge>
-{
+class IncidentFaceIterator {
 
 public:
 
-	typedef IncidentFaceIterator _Mytype;
+  // ITERATOR TRAITS -------------------------------
+  /// One of the @link iterator_tags tag types@endlink.
+  typedef std::forward_iterator_tag iterator_category;
+  /// The type "pointed to" by the iterator.
+  typedef Halfedge value_type;
+  /// Distance between iterators is represented as this type.
+  typedef ptrdiff_t difference_type;
+  /// This type represents a pointer-to-value_type.
+  typedef Halfedge *pointer;
+  /// This type represents a reference-to-value_type.
+  typedef Halfedge &reference;
+  // -----------------------------------------------
 
-	value_type operator*() { return m_h; }
-	_Mytype& operator++();
-	bool operator==(_Mytype const& other) const
-	{
-		return m_h == other.m_h;
-	}
 
-	bool operator!=(_Mytype const& other) const
-	{
-		return !(*this == other);
-	}
+  typedef IncidentFaceIterator _Mytype;
 
-private:
+  value_type operator*() { return m_h; }
+  _Mytype &operator++();
+  bool operator==(_Mytype const &other) const { return m_h == other.m_h; }
 
-	friend class IncidentFaceContainer;
-	IncidentFaceIterator(unsigned int v, TriangleMesh const* mesh);
-	IncidentFaceIterator() : m_h(), m_begin(), m_mesh(nullptr) {}
+  bool operator!=(_Mytype const &other) const { return !(*this == other); }
 
-	Halfedge m_h, m_begin;
-	TriangleMesh const* m_mesh;
-};
+private:
+  friend class IncidentFaceContainer;
+  IncidentFaceIterator(unsigned int v, TriangleMesh const *mesh);
+  IncidentFaceIterator() : m_h(), m_begin(), m_mesh(nullptr) {}
 
+  Halfedge m_h, m_begin;
+  TriangleMesh const *m_mesh;
+};
 
 class VertexContainer;
-class VertexIterator : public std::iterator<std::random_access_iterator_tag, Eigen::Vector3d>
-{
+class VertexIterator {
 
 public:
 
-	typedef VertexIterator _Mytype;
-
-	VertexIterator() = delete;
-
-	reference operator*();
-
-	bool operator<(_Mytype const& other) const
-	{
-		return m_index < other.m_index;
-	}
-	bool operator==(_Mytype const& other) const
-	{
-		return m_index == other.m_index;
-	}
-
-	bool operator!=(_Mytype const& other) const
-	{
-		return !(*this == other);
-	}
-
-	inline _Mytype& operator++() { ++m_index; return *this; }
-	inline _Mytype& operator--() { --m_index; return *this; }
-
-	inline _Mytype operator+(_Mytype const& rhs) const
-	{ 
-		return _Mytype(m_index + rhs.m_index, m_mesh); 
-	}
-	inline difference_type operator-(_Mytype const& rhs) const
-	{ 
-		return m_index - rhs.m_index; 
-	}
-	inline _Mytype operator-(int const& rhs) const
-	{ 
-		return _Mytype(m_index - rhs, m_mesh); 
-	}
-
-	unsigned int index() const;
+  // ITERATOR TRAITS -------------------------------
+  /// One of the @link iterator_tags tag types@endlink.
+  typedef std::random_access_iterator_tag iterator_category;
+  /// The type "pointed to" by the iterator.
+  typedef Eigen::Vector3d value_type;
+  /// Distance between iterators is represented as this type.
+  typedef ptrdiff_t difference_type;
+  /// This type represents a pointer-to-value_type.
+  typedef Eigen::Vector3d *pointer;
+  /// This type represents a reference-to-value_type.
+  typedef Eigen::Vector3d &reference;
+  // -----------------------------------------------
+
+
+  typedef VertexIterator _Mytype;
+
+  VertexIterator() = delete;
+
+  reference operator*();
+
+  bool operator<(_Mytype const &other) const { return m_index < other.m_index; }
+  bool operator==(_Mytype const &other) const {
+    return m_index == other.m_index;
+  }
+
+  bool operator!=(_Mytype const &other) const { return !(*this == other); }
+
+  inline _Mytype &operator++() {
+    ++m_index;
+    return *this;
+  }
+  inline _Mytype &operator--() {
+    --m_index;
+    return *this;
+  }
+
+  inline _Mytype operator+(_Mytype const &rhs) const {
+    return _Mytype(m_index + rhs.m_index, m_mesh);
+  }
+  inline difference_type operator-(_Mytype const &rhs) const {
+    return m_index - rhs.m_index;
+  }
+  inline _Mytype operator-(int const &rhs) const {
+    return _Mytype(m_index - rhs, m_mesh);
+  }
+
+  unsigned int index() const;
 
 private:
+  friend class VertexContainer;
+  VertexIterator(unsigned int index, TriangleMesh *mesh)
+      : m_index(index), m_mesh(mesh) {}
 
-	friend class VertexContainer;
-	VertexIterator(unsigned int index, TriangleMesh* mesh)
-		: m_index(index), m_mesh(mesh) {}
-
-	unsigned int m_index;
-	TriangleMesh* m_mesh;
+  unsigned int m_index;
+  TriangleMesh *m_mesh;
 };
 
-
 class VertexConstContainer;
-class VertexConstIterator : 
-	public std::iterator<std::random_access_iterator_tag, Eigen::Vector3d const>
-{
+class VertexConstIterator {
 
 public:
-
-	typedef VertexConstIterator _Mytype;
-
-	VertexConstIterator() = delete;
-
-	reference operator*();
-
-	bool operator<(_Mytype const& other) const
-	{
-		return m_index < other.m_index;
-	}
-	bool operator==(_Mytype const& other) const
-	{
-		return m_index == other.m_index;
-	}
-
-	bool operator!=(_Mytype const& other) const
-	{
-		return !(*this == other);
-	}
-
-	inline _Mytype& operator++() { ++m_index; return *this; }
-	inline _Mytype& operator--() { --m_index; return *this; }
-
-	inline _Mytype operator+(_Mytype const& rhs) const
-	{ 
-		return _Mytype(m_index + rhs.m_index, m_mesh); 
-	}
-	inline difference_type operator-(_Mytype const& rhs) const
-	{ 
-		return m_index - rhs.m_index; 
-	}
-	inline _Mytype operator-(int const& rhs) const
-	{ 
-		return _Mytype(m_index - rhs, m_mesh); 
-	}
-
-	unsigned int index() const;
+  // ITERATOR TRAITS -------------------------------
+  /// One of the @link iterator_tags tag types@endlink.
+  typedef std::random_access_iterator_tag iterator_category;
+  /// The type "pointed to" by the iterator.
+  typedef Eigen::Vector3d const value_type;
+  /// Distance between iterators is represented as this type.
+  typedef ptrdiff_t difference_type;
+  /// This type represents a pointer-to-value_type.
+  typedef Eigen::Vector3d const *pointer;
+  /// This type represents a reference-to-value_type.
+  typedef Eigen::Vector3d const &reference;
+  // -----------------------------------------------
+
+  typedef VertexConstIterator _Mytype;
+
+  VertexConstIterator() = delete;
+
+  reference operator*();
+
+  bool operator<(_Mytype const &other) const { return m_index < other.m_index; }
+  bool operator==(_Mytype const &other) const {
+    return m_index == other.m_index;
+  }
+
+  bool operator!=(_Mytype const &other) const { return !(*this == other); }
+
+  inline _Mytype &operator++() {
+    ++m_index;
+    return *this;
+  }
+  inline _Mytype &operator--() {
+    --m_index;
+    return *this;
+  }
+
+  inline _Mytype operator+(_Mytype const &rhs) const {
+    return _Mytype(m_index + rhs.m_index, m_mesh);
+  }
+  inline difference_type operator-(_Mytype const &rhs) const {
+    return m_index - rhs.m_index;
+  }
+  inline _Mytype operator-(int const &rhs) const {
+    return _Mytype(m_index - rhs, m_mesh);
+  }
+
+  unsigned int index() const;
 
 private:
+  friend class VertexConstContainer;
+  VertexConstIterator(unsigned int index, TriangleMesh const *mesh)
+      : m_index(index), m_mesh(mesh) {}
 
-	friend class VertexConstContainer;
-	VertexConstIterator(unsigned int index, TriangleMesh const* mesh)
-		: m_index(index), m_mesh(mesh) {}
-
-	unsigned int m_index;
-	TriangleMesh const* m_mesh;
+  unsigned int m_index;
+  TriangleMesh const *m_mesh;
 };
-}
-
+} // namespace Discregrid
diff --git a/discregrid/src/cubic_lagrange_discrete_grid.cpp b/discregrid/src/cubic_lagrange_discrete_grid.cpp
index a6396a7..ff19cbd 100755
--- a/discregrid/src/cubic_lagrange_discrete_grid.cpp
+++ b/discregrid/src/cubic_lagrange_discrete_grid.cpp
@@ -818,7 +818,7 @@ CubicLagrangeDiscreteGrid::addFunction(ContinuousFunction const &func, bool verb
 
 			if (verbose && (++counter == n_nodes || duration_cast<milliseconds>(high_resolution_clock::now() - t0).count() > 1000u))
 			{
-				std::async(std::launch::async, [&]() {
+				auto _discard = std::async(std::launch::async, [&]() {
 					mutex.lock();
 					t0 = high_resolution_clock::now();
 					std::cout << "\r"