From 657cc1af37e6532e8ca96ec89c4e4733c7f801ed Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@uliege.be>
Date: Sat, 16 Jan 2021 10:07:09 +0100
Subject: [PATCH] use nullptr

---
 Geo/CGNSWrite.cpp                    |   8 +-
 Geo/CGNSZone.cpp                     |   6 +-
 Geo/CGNSZoneStruct.cpp               |  16 ++--
 Geo/CGNSZoneUnstruct.cpp             |  10 +--
 Geo/CellComplex.cpp                  |   4 +-
 Geo/ChainComplex.cpp                 |  70 +++++++--------
 Geo/ChainComplex.h                   |  18 ++--
 Geo/ExtrudeParams.cpp                |   2 +-
 Geo/GEdge.cpp                        |  18 ++--
 Geo/GEdge.h                          |   4 +-
 Geo/GEdgeLoop.cpp                    |  10 +--
 Geo/GEntity.cpp                      |   6 +-
 Geo/GEntity.h                        |  13 +--
 Geo/GFace.cpp                        |  26 +++---
 Geo/GFace.h                          |  13 ++-
 Geo/GModel.cpp                       |  50 +++++------
 Geo/GModelCreateTopologyFromMesh.cpp |  10 +--
 Geo/GModelIO_ACTRAN.cpp              |   2 +-
 Geo/GModelIO_BDF.cpp                 |   2 +-
 Geo/GModelIO_CGNS.cpp                |   2 +-
 Geo/GModelIO_DIFF.cpp                |  18 ++--
 Geo/GModelIO_GEO.cpp                 |  48 +++++-----
 Geo/GModelIO_MED.cpp                 |   4 +-
 Geo/GModelIO_MSH2.cpp                |  20 ++---
 Geo/GModelIO_MSH3.cpp                |  12 +--
 Geo/GModelIO_MSH4.cpp                | 112 +++++++++++------------
 Geo/GModelIO_OCC.cpp                 |  42 ++++-----
 Geo/GModelIO_OCC.h                   |   6 +-
 Geo/GModelIO_SAMCEF.cpp              |   2 +-
 Geo/GModelIO_UNV.cpp                 |   2 +-
 Geo/GModelParametrize.cpp            |  12 +--
 Geo/GPoint.h                         |   3 +-
 Geo/GRegion.cpp                      |  20 ++---
 Geo/GRegion.h                        |   2 +-
 Geo/GVertex.cpp                      |   4 +-
 Geo/GVertex.h                        |   2 +-
 Geo/Geo.cpp                          | 130 +++++++++++++--------------
 Geo/Geo.h                            |   2 +-
 Geo/GeoInterpolation.cpp             |   4 +-
 Geo/GeomMeshMatcher.cpp              |  24 ++---
 Geo/Homology.cpp                     |  14 +--
 Geo/MEdge.cpp                        |  10 +--
 Geo/MEdge.h                          |   2 +-
 Geo/MElement.cpp                     |  16 ++--
 Geo/MElement.h                       |  39 ++++----
 Geo/MElementCut.cpp                  |  36 ++++----
 Geo/MElementCut.h                    |  78 ++++++++--------
 Geo/MElementOctree.cpp               |   4 +-
 Geo/MFace.h                          |   2 +-
 Geo/MSubElement.cpp                  |  32 +++----
 Geo/MSubElement.h                    |  87 +++++++++---------
 Geo/MTriangle.h                      |   2 +-
 Geo/MVertex.h                        |   7 +-
 Geo/MVertexBoundaryLayerData.cpp     |   2 +-
 Geo/MVertexRTree.h                   |   8 +-
 Geo/OCCEdge.cpp                      |   4 +-
 Geo/OCCEdge.h                        |   2 +-
 Geo/OCCFace.cpp                      |   4 +-
 Geo/OCCRegion.cpp                    |   6 +-
 Geo/SOrientedBoundingBox.cpp         |   4 +-
 Geo/boundaryLayersData.cpp           |   2 +-
 Geo/boundaryLayersData.h             |  10 +--
 Geo/closestVertex.cpp                |   4 +-
 Geo/discreteEdge.cpp                 |   4 +-
 Geo/discreteFace.cpp                 |   4 +-
 Geo/discreteFace.h                   |   4 +-
 Geo/ghostEdge.h                      |   2 +-
 Geo/gmshFace.cpp                     |   6 +-
 Geo/gmshLevelset.cpp                 |   4 +-
 Geo/gmshLevelset.h                   |   2 +-
 Geo/gmshSurface.cpp                  |   4 +-
 Geo/gmshSurface.h                    |   2 +-
 Geo/partitionEdge.h                  |   6 +-
 Geo/partitionFace.h                  |   6 +-
 Geo/partitionRegion.h                |   7 +-
 Geo/partitionVertex.h                |   7 +-
 Geo/xyEdge.h                         |   2 +-
 Geo/xyFace.h                         |   2 +-
 78 files changed, 601 insertions(+), 595 deletions(-)

diff --git a/Geo/CGNSWrite.cpp b/Geo/CGNSWrite.cpp
index a6eec7b7b8..54abc18666 100644
--- a/Geo/CGNSWrite.cpp
+++ b/Geo/CGNSWrite.cpp
@@ -89,7 +89,7 @@ namespace {
     // names
     std::ostringstream ossEnt;
     ossEnt << entTypeStr;
-    if(geParent != 0) ossEnt << "_" << geParent->tag();
+    if(geParent != nullptr) ossEnt << "_" << geParent->tag();
     ossEnt << "_" << ge->tag();
     entityName = model->getElementaryName(ge->dim(), ge->tag());
     if(!entityName.empty()) ossEnt << "_" << entityName;
@@ -290,7 +290,7 @@ int writeZone(GModel *model, bool saveAll, double scalingFactor, int meshDim,
     // get entities
     GEntity *ge = entities[i];
     GEntity *geGeom = ge->getParentEntity();
-    if(geGeom == 0) geGeom = ge;
+    if(geGeom == nullptr) geGeom = ge;
 
     // get or create the names for the entity
     std::string entityName, geomEntityName;
@@ -484,7 +484,7 @@ void getEntitiesInPartitions(const std::vector<GEntity *> &entities,
 {
   for(std::size_t j = 0; j < entities.size(); j++) {
     GEntity *ge = entities[j];
-    const std::vector<int> *parts = 0;
+    const std::vector<int> *parts = nullptr;
     switch(ge->geomType()) {
     case GEntity::PartitionVolume: {
       partitionRegion *pr = static_cast<partitionRegion *>(ge);
@@ -512,7 +512,7 @@ void getEntitiesInPartitions(const std::vector<GEntity *> &entities,
     } break;
     default: break;
     } // switch
-    if(parts != 0) {
+    if(parts != nullptr) {
       for(std::size_t iPart = 0; iPart < parts->size(); iPart++) {
         entitiesPart[(*parts)[iPart]].push_back(ge);
       }
diff --git a/Geo/CGNSZone.cpp b/Geo/CGNSZone.cpp
index 8aab0fa025..5e9784cc5a 100644
--- a/Geo/CGNSZone.cpp
+++ b/Geo/CGNSZone.cpp
@@ -22,7 +22,7 @@ CGNSZone::CGNSZone(int fileIndex, int baseIndex, int zoneIndex,
                    const Family2EltNodeTransfo &allEltNodeTransfo, int &err)
   : fileIndex_(fileIndex), baseIndex_(baseIndex), meshDim_(meshDim),
     zoneIndex_(zoneIndex), type_(type), startNode_(startNode),
-    eltNodeTransfo_(0), nbPerConnect_(0)
+    eltNodeTransfo_(nullptr), nbPerConnect_(0)
 {
   int cgnsErr;
 
@@ -351,7 +351,7 @@ int CGNSZone::readBoundaryConditionRange(int iZoneBC,
 
   std::vector<cgsize_t> bcData(indexDataSize(2));
   cgnsErr =
-    cg_boco_read(fileIndex(), baseIndex(), index(), iZoneBC, bcData.data(), 0);
+    cg_boco_read(fileIndex(), baseIndex(), index(), iZoneBC, bcData.data(), nullptr);
   if(cgnsErr != CG_OK) return cgnsError(__FILE__, __LINE__, fileIndex());
 
   // get list of elements from range data
@@ -368,7 +368,7 @@ int CGNSZone::readBoundaryConditionList(int iZoneBC, cgsize_t nbVal,
   // read data
   std::vector<cgsize_t> bcData(indexDataSize(nbVal));
   cgnsErr =
-    cg_boco_read(fileIndex(), baseIndex(), index(), iZoneBC, bcData.data(), 0);
+    cg_boco_read(fileIndex(), baseIndex(), index(), iZoneBC, bcData.data(), nullptr);
   if(cgnsErr != CG_OK) return cgnsError(__FILE__, __LINE__, fileIndex());
 
   // get list of elements from list data
diff --git a/Geo/CGNSZoneStruct.cpp b/Geo/CGNSZoneStruct.cpp
index 15d6aeafd4..a1bbeb544b 100644
--- a/Geo/CGNSZoneStruct.cpp
+++ b/Geo/CGNSZoneStruct.cpp
@@ -251,7 +251,7 @@ namespace {
     }
 
     // do no add element if it is part of an internal interface between blocks
-    if(isInternalInterface) return 0;
+    if(isInternalInterface) return nullptr;
 
     // create element
     MElementFactory factory;
@@ -336,7 +336,7 @@ namespace {
     }
 
     // do no add element if it is part of an internal interface between blocks
-    if(isInternalInterface) return 0;
+    if(isInternalInterface) return nullptr;
 
     // create element
     MElementFactory factory;
@@ -430,11 +430,11 @@ int CGNSZoneStruct<DIM>::readElements(
       cgsize_t ijk[3] = {0, j * order, k * order};
       MElement *me;
       me = makeBndElement(ijk, dir, order, firstDefaultEnt, allVert, allElt);
-      if(me != 0) zoneElt.push_back(me);
+      if(me != nullptr) zoneElt.push_back(me);
       ijk[0] = nbNodeIJK(0) - 1;
       me =
         makeBndElement(ijk, dir, order, firstDefaultEnt + 1, allVert, allElt);
-      if(me != 0) zoneElt.push_back(me);
+      if(me != nullptr) zoneElt.push_back(me);
     }
   }
 
@@ -446,11 +446,11 @@ int CGNSZoneStruct<DIM>::readElements(
       MElement *me;
       me =
         makeBndElement(ijk, dir, order, firstDefaultEnt + 2, allVert, allElt);
-      if(me != 0) zoneElt.push_back(me);
+      if(me != nullptr) zoneElt.push_back(me);
       ijk[1] = nbNodeIJK(1) - 1;
       me =
         makeBndElement(ijk, dir, order, firstDefaultEnt + 3, allVert, allElt);
-      if(me != 0) zoneElt.push_back(me);
+      if(me != nullptr) zoneElt.push_back(me);
     }
   }
 
@@ -463,11 +463,11 @@ int CGNSZoneStruct<DIM>::readElements(
         MElement *me;
         me =
           makeBndElement(ijk, dir, order, firstDefaultEnt + 4, allVert, allElt);
-        if(me != 0) zoneElt.push_back(me);
+        if(me != nullptr) zoneElt.push_back(me);
         ijk[2] = nbNodeIJK(2) - 1;
         me =
           makeBndElement(ijk, dir, order, firstDefaultEnt + 5, allVert, allElt);
-        if(me != 0) zoneElt.push_back(me);
+        if(me != nullptr) zoneElt.push_back(me);
       }
     }
   }
diff --git a/Geo/CGNSZoneUnstruct.cpp b/Geo/CGNSZoneUnstruct.cpp
index a6ed94530b..a677ac4a98 100644
--- a/Geo/CGNSZoneUnstruct.cpp
+++ b/Geo/CGNSZoneUnstruct.cpp
@@ -41,8 +41,8 @@ namespace {
 
     // element high-order node transformation if specified (CPEX0045), otherwise
     // use the classic CGNS order
-    const std::vector<int> *nodeTransfo = 0;
-    if((mshEltType != MSH_PNT) && (eltNodeTransfo != 0) &&
+    const std::vector<int> *nodeTransfo = nullptr;
+    if((mshEltType != MSH_PNT) && (eltNodeTransfo != nullptr) &&
        (eltNodeTransfo->size() > 0)) {
       nodeTransfo = &((*eltNodeTransfo)[mshEltType]);
     }
@@ -133,12 +133,12 @@ int CGNSZoneUnstruct::readSection(
   if(sectEltType == CGNS_ENUMV(MIXED)) {
 #if CGNS_VERSION >= 4000
     cgnsErr = cg_poly_elements_read(fileIndex(), baseIndex(), index(), iSect,
-                                    sectData.data(), offsetData.data(), 0);
+                                    sectData.data(), offsetData.data(), nullptr);
 #endif
   }
   else {
     cgnsErr = cg_elements_read(fileIndex(), baseIndex(), index(), iSect,
-                               sectData.data(), 0);
+                               sectData.data(), nullptr);
   }
   if(cgnsErr != CG_OK) return cgnsError(__FILE__, __LINE__, fileIndex());
 
@@ -168,7 +168,7 @@ int CGNSZoneUnstruct::readElements(
   // data structures for node coordinate transformation (CPEX0045)
   // std::vector<bool> nodeUpdated;
   std::vector<SPoint3> rawNode;
-  if(eltNodeTransfo() != 0) {
+  if(eltNodeTransfo() != nullptr) {
     // nodeUpdated = std::vector<bool>(nbNode(), false);
     rawNode.resize(nbNode());
     for(int iN = 0; iN < nbNode(); iN++) {
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 7a8e1d6cc2..d2ce3c3b98 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -911,7 +911,7 @@ Cell *CellComplex::getACell(int dim, int domain)
       Msg::Warning("%d cells in domain", num);
     else if(domain == 2)
       Msg::Warning("%d cells in subdomain", num);
-    return NULL;
+    return nullptr;
   }
 
   for(auto cit = firstCell(dim); cit != lastCell(dim); cit++) {
@@ -921,7 +921,7 @@ Cell *CellComplex::getACell(int dim, int domain)
       return cell;
   }
   Msg::Debug("Domain cell counts not in sync.");
-  return NULL;
+  return nullptr;
 }
 
 bool CellComplex::restoreComplex()
diff --git a/Geo/ChainComplex.cpp b/Geo/ChainComplex.cpp
index bde99ed532..d0a3d1900d 100644
--- a/Geo/ChainComplex.cpp
+++ b/Geo/ChainComplex.cpp
@@ -16,12 +16,12 @@ ChainComplex::ChainComplex(CellComplex *cellComplex, int domain)
   _cellComplex = cellComplex;
 
   for(int i = 0; i < 5; i++) {
-    _hMatrix[i] = NULL;
-    _kerH[i] = NULL;
-    _codH[i] = NULL;
-    _jMatrix[i] = NULL;
-    _qMatrix[i] = NULL;
-    _hbasis[i] = NULL;
+    _hMatrix[i] = nullptr;
+    _kerH[i] = nullptr;
+    _codH[i] = nullptr;
+    _jMatrix[i] = nullptr;
+    _qMatrix[i] = nullptr;
+    _hbasis[i] = nullptr;
   }
 
   int lastCols = 0;
@@ -49,7 +49,7 @@ ChainComplex::ChainComplex(CellComplex *cellComplex, int domain)
     lastCols = cols;
 
     if(cols == 0) { // no dim-cells, no map
-      _hMatrix[dim] = NULL;
+      _hMatrix[dim] = nullptr;
     }
     else if(rows == 0) { // no dim-1-cells, maps everything to zero
       _hMatrix[dim] = create_gmp_matrix_zero(1, cols);
@@ -96,11 +96,11 @@ ChainComplex::ChainComplex(CellComplex *cellComplex, int domain)
       }
       mpz_clear(elem);
     }
-    _kerH[dim] = NULL;
-    _codH[dim] = NULL;
-    _jMatrix[dim] = NULL;
-    _qMatrix[dim] = NULL;
-    _hbasis[dim] = NULL;
+    _kerH[dim] = nullptr;
+    _codH[dim] = nullptr;
+    _jMatrix[dim] = nullptr;
+    _qMatrix[dim] = nullptr;
+    _hbasis[dim] = nullptr;
   }
 }
 
@@ -119,17 +119,17 @@ ChainComplex::~ChainComplex()
 void ChainComplex::transposeHMatrices()
 {
   for(int i = 0; i < 5; i++)
-    if(_hMatrix[i] != NULL) gmp_matrix_transp(_hMatrix[i]);
+    if(_hMatrix[i] != nullptr) gmp_matrix_transp(_hMatrix[i]);
 }
 void ChainComplex::transposeHMatrix(int dim)
 {
-  if(dim > -1 && dim < 5 && _hMatrix[dim] != NULL)
+  if(dim > -1 && dim < 5 && _hMatrix[dim] != nullptr)
     gmp_matrix_transp(_hMatrix[dim]);
 }
 
 void ChainComplex::KerCod(int dim)
 {
-  if(dim < 0 || dim > 3 || _hMatrix[dim] == NULL) return;
+  if(dim < 0 || dim > 3 || _hMatrix[dim] == nullptr) return;
 
   gmp_matrix *HMatrix =
     copy_gmp_matrix(_hMatrix[dim], 1, 1, gmp_matrix_rows(_hMatrix[dim]),
@@ -175,7 +175,7 @@ void ChainComplex::KerCod(int dim)
 // j:B_k->Z_k
 void ChainComplex::Inclusion(int lowDim, int highDim)
 {
-  if(getKerHMatrix(lowDim) == NULL || getCodHMatrix(highDim) == NULL ||
+  if(getKerHMatrix(lowDim) == nullptr || getCodHMatrix(highDim) == nullptr ||
      abs(lowDim - highDim) != 1)
     return;
 
@@ -263,7 +263,7 @@ void ChainComplex::Inclusion(int lowDim, int highDim)
 
 void ChainComplex::Quotient(int dim, int setDim)
 {
-  if(dim < 0 || dim > 4 || _jMatrix[dim] == NULL) return;
+  if(dim < 0 || dim > 4 || _jMatrix[dim] == nullptr) return;
   if(setDim < 0 || setDim > 4) return;
 
   gmp_matrix *JMatrix =
@@ -329,22 +329,22 @@ void ChainComplex::computeHomology(bool dual)
 
     // 1) no edges, but zero cells
     if(lowDim == 0 && !dual && gmp_matrix_cols(getHMatrix(lowDim)) > 0 &&
-       getHMatrix(highDim) == NULL) {
+       getHMatrix(highDim) == nullptr) {
       setHbasis(setDim, create_gmp_matrix_identity(
                           gmp_matrix_cols(getHMatrix(lowDim))));
     }
     else if(highDim == 0 && dual && gmp_matrix_rows(getHMatrix(highDim)) > 0 &&
-            getHMatrix(lowDim) == NULL) {
+            getHMatrix(lowDim) == nullptr) {
       setHbasis(setDim, create_gmp_matrix_identity(
                           gmp_matrix_rows(getHMatrix(highDim))));
     }
 
     // 2) this dimension is empty
-    else if(getHMatrix(setDim) == NULL) {
-      setHbasis(setDim, NULL);
+    else if(getHMatrix(setDim) == nullptr) {
+      setHbasis(setDim, nullptr);
     }
     // 3) No higher dimension cells -> none of the cycles are boundaries
-    else if(getHMatrix(highDim) == NULL) {
+    else if(getHMatrix(highDim) == nullptr) {
       setHbasis(setDim,
                 copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
                                 gmp_matrix_rows(getKerHMatrix(lowDim)),
@@ -357,21 +357,21 @@ void ChainComplex::computeHomology(bool dual)
     //   3) find quotient Z/j(B)
     else {
       // 4) No lower dimension cells -> all chains are cycles
-      if(getHMatrix(lowDim) == NULL) {
+      if(getHMatrix(lowDim) == nullptr) {
         setKerHMatrix(lowDim, create_gmp_matrix_identity(
                                 gmp_matrix_rows(getHMatrix(highDim))));
       }
       Inclusion(lowDim, highDim);
       Quotient(lowDim, setDim);
 
-      if(getCodHMatrix(highDim) == NULL) {
+      if(getCodHMatrix(highDim) == nullptr) {
         setHbasis(setDim,
                   copy_gmp_matrix(getKerHMatrix(lowDim), 1, 1,
                                   gmp_matrix_rows(getKerHMatrix(lowDim)),
                                   gmp_matrix_cols(getKerHMatrix(lowDim))));
       }
-      else if(getJMatrix(lowDim) == NULL || getQMatrix(lowDim) == NULL) {
-        setHbasis(setDim, NULL);
+      else if(getJMatrix(lowDim) == nullptr || getQMatrix(lowDim) == nullptr) {
+        setHbasis(setDim, nullptr);
       }
       else {
         setHbasis(setDim,
@@ -386,8 +386,8 @@ void ChainComplex::computeHomology(bool dual)
     destroy_gmp_matrix(getJMatrix(lowDim));
     destroy_gmp_matrix(getQMatrix(lowDim));
 
-    setJMatrix(lowDim, NULL);
-    setQMatrix(lowDim, NULL);
+    setJMatrix(lowDim, nullptr);
+    setQMatrix(lowDim, nullptr);
   }
   return;
 }
@@ -413,7 +413,7 @@ std::vector<int> ChainComplex::getCoeffVector(int dim, int chainNumber)
   std::vector<int> coeffVector;
 
   if(dim < 0 || dim > 4) return coeffVector;
-  if(_hbasis[dim] == NULL || (int)gmp_matrix_cols(_hbasis[dim]) < chainNumber)
+  if(_hbasis[dim] == nullptr || (int)gmp_matrix_cols(_hbasis[dim]) < chainNumber)
     return coeffVector;
 
   int rows = gmp_matrix_rows(_hbasis[dim]);
@@ -439,13 +439,13 @@ gmp_matrix *ChainComplex::getBasis(int dim, int basis)
 {
   if(dim > -2 && dim < 4 && basis == 2) return _codH[dim + 1];
   if(dim < 0 || dim > 4)
-    return NULL;
+    return nullptr;
   else if(basis == 1)
     return _kerH[dim];
   else if(basis == 3)
     return _hbasis[dim];
   else
-    return NULL;
+    return nullptr;
 }
 
 void ChainComplex::getBasisChain(std::map<Cell *, int, CellPtrLessThan> &chain,
@@ -456,7 +456,7 @@ void ChainComplex::getBasisChain(std::map<Cell *, int, CellPtrLessThan> &chain,
 
   chain.clear();
   if(dim < 0 || dim > 3) return;
-  if(basisMatrix == NULL || (int)gmp_matrix_cols(basisMatrix) < num) { return; }
+  if(basisMatrix == nullptr || (int)gmp_matrix_cols(basisMatrix) < num) { return; }
 
   int elemi;
   long int elemli;
@@ -491,7 +491,7 @@ void ChainComplex::getBasisChain(std::map<Cell *, int, CellPtrLessThan> &chain,
 int ChainComplex::getBasisSize(int dim, int basis)
 {
   gmp_matrix *basisMatrix;
-  if(basis == 0 && _hMatrix[dim] != NULL) {
+  if(basis == 0 && _hMatrix[dim] != nullptr) {
     return gmp_matrix_cols(_hMatrix[dim]);
   }
   else if(basis == 1)
@@ -503,7 +503,7 @@ int ChainComplex::getBasisSize(int dim, int basis)
   else
     return 0;
 
-  if(basisMatrix != NULL && gmp_matrix_rows(basisMatrix) != 0)
+  if(basisMatrix != nullptr && gmp_matrix_rows(basisMatrix) != 0)
     return gmp_matrix_cols(basisMatrix);
   else
     return 0;
@@ -512,7 +512,7 @@ int ChainComplex::getBasisSize(int dim, int basis)
 int ChainComplex::getTorsion(int dim, int num)
 {
   if(dim < 0 || dim > 4) return 0;
-  if(_hbasis[dim] == NULL || (int)gmp_matrix_cols(_hbasis[dim]) < num) return 0;
+  if(_hbasis[dim] == nullptr || (int)gmp_matrix_cols(_hbasis[dim]) < num) return 0;
   if(_torsion[dim].empty() || (int)_torsion[dim].size() < num)
     return 1;
   else
diff --git a/Geo/ChainComplex.h b/Geo/ChainComplex.h
index 082949f077..5137bc5557 100644
--- a/Geo/ChainComplex.h
+++ b/Geo/ChainComplex.h
@@ -93,42 +93,42 @@ private:
     if(dim > -1 && dim < 5)
       return _hMatrix[dim];
     else
-      return NULL;
+      return nullptr;
   }
   gmp_matrix *getKerHMatrix(int dim) const
   {
     if(dim > -1 && dim < 5)
       return _kerH[dim];
     else
-      return NULL;
+      return nullptr;
   }
   gmp_matrix *getCodHMatrix(int dim) const
   {
     if(dim > -1 && dim < 5)
       return _codH[dim];
     else
-      return NULL;
+      return nullptr;
   }
   gmp_matrix *getJMatrix(int dim) const
   {
     if(dim > -1 && dim < 5)
       return _jMatrix[dim];
     else
-      return NULL;
+      return nullptr;
   }
   gmp_matrix *getQMatrix(int dim) const
   {
     if(dim > -1 && dim < 5)
       return _qMatrix[dim];
     else
-      return NULL;
+      return nullptr;
   }
   gmp_matrix *getHbasis(int dim) const
   {
     if(dim > -1 && dim < 5)
       return _hbasis[dim];
     else
-      return NULL;
+      return nullptr;
   }
 
   // local deformation tools for chains
@@ -160,7 +160,7 @@ public:
     if(dim > -1 && dim < 5)
       return _hMatrix[dim];
     else
-      return NULL;
+      return nullptr;
   }
 
   // compute basis for kernel and codomain of boundary operator matrix
@@ -189,7 +189,7 @@ public:
   // get the cell index
   int getCellIndex(Cell *cell)
   {
-    citer cit = _cellIndices[cell->getDim()].find(cell);
+    auto cit = _cellIndices[cell->getDim()].find(cell);
     if(cit != lastCell(cell->getDim()))
       return cit->second;
     else
@@ -212,7 +212,7 @@ public:
   // apply a transformation T to a basis (T should be unimodular)
   void transformBasis(gmp_matrix *T, int dim, int basis)
   {
-    if(basis == 3 && _hbasis[dim] != NULL) {
+    if(basis == 3 && _hbasis[dim] != nullptr) {
       gmp_matrix_right_mult(_hbasis[dim], T);
     }
   }
diff --git a/Geo/ExtrudeParams.cpp b/Geo/ExtrudeParams.cpp
index 89b8d3d530..84572dfa73 100644
--- a/Geo/ExtrudeParams.cpp
+++ b/Geo/ExtrudeParams.cpp
@@ -7,7 +7,7 @@
 #include "Geo.h"
 #include "ExtrudeParams.h"
 
-smooth_data *ExtrudeParams::normals[2] = {0, 0};
+smooth_data *ExtrudeParams::normals[2] = {nullptr, nullptr};
 std::vector<SPoint3> ExtrudeParams::normalsCoherence;
 
 // Scale last layer size locally If one section of the boundary layer index = 0
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 59ce6f93bf..911e125af5 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -24,8 +24,8 @@
 #endif
 
 GEdge::GEdge(GModel *model, int tag, GVertex *v0, GVertex *v1)
-  : GEntity(model, tag), _length(0.), _tooSmall(false), _cp(0), _v0(v0),
-    _v1(v1), masterOrientation(0), compoundCurve(NULL)
+  : GEntity(model, tag), _length(0.), _tooSmall(false), _cp(nullptr), _v0(v0),
+    _v1(v1), masterOrientation(0), compoundCurve(nullptr)
 {
   if(_v0) _v0->addEdge(this);
   if(_v1 && _v1 != _v0) _v1->addEdge(this);
@@ -34,8 +34,8 @@ GEdge::GEdge(GModel *model, int tag, GVertex *v0, GVertex *v1)
 }
 
 GEdge::GEdge(GModel *model, int tag)
-  : GEntity(model, tag), _length(0.), _tooSmall(false), _cp(0), _v0(0), _v1(0),
-    masterOrientation(0), compoundCurve(NULL)
+  : GEntity(model, tag), _length(0.), _tooSmall(false), _cp(nullptr), _v0(nullptr), _v1(nullptr),
+    masterOrientation(0), compoundCurve(nullptr)
 {
   meshStatistics.status = GEdge::PENDING;
   GEdge::resetMeshAttributes();
@@ -177,14 +177,14 @@ void GEdge::getNumMeshElements(unsigned *const c) const
 
 MElement *const *GEdge::getStartElementType(int type) const
 {
-  if(lines.empty()) return 0; // msvc would throw an exception
+  if(lines.empty()) return nullptr; // msvc would throw an exception
   return reinterpret_cast<MElement *const *>(&lines[0]);
 }
 
 MElement *GEdge::getMeshElement(std::size_t index) const
 {
   if(index < lines.size()) return lines[index];
-  return 0;
+  return nullptr;
 }
 
 MElement *GEdge::getMeshElementByType(const int familyType,
@@ -192,7 +192,7 @@ MElement *GEdge::getMeshElementByType(const int familyType,
 {
   if(familyType == TYPE_LIN) return lines[index];
 
-  return 0;
+  return nullptr;
 }
 
 void GEdge::resetMeshAttributes()
@@ -201,7 +201,7 @@ void GEdge::resetMeshAttributes()
   meshAttributes.coeffTransfinite = 0.;
   meshAttributes.nbPointsTransfinite = 0;
   meshAttributes.typeTransfinite = 0;
-  meshAttributes.extrude = 0;
+  meshAttributes.extrude = nullptr;
   meshAttributes.meshSize = MAX_LC;
   meshAttributes.meshSizeFactor = 1.;
   meshAttributes.minimumMeshSegments = 1;
@@ -448,7 +448,7 @@ double GEdge::curvature(double par) const
 
 double GEdge::length(const double &u0, const double &u1, const int nbQuadPoints)
 {
-  double *t = 0, *w = 0;
+  double *t = nullptr, *w = nullptr;
   gmshGaussLegendre1D(nbQuadPoints, &t, &w);
   if(!t) {
     Msg::Error("Gauss-Legendre integration returned no points");
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 1a4d76aec7..4b4ee8e141 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -86,7 +86,7 @@ public:
   virtual int dim() const { return 1; }
 
   // returns the parent entity for partitioned entities
-  virtual GEntity *getParentEntity() { return 0; }
+  virtual GEntity *getParentEntity() { return nullptr; }
 
   // get the list of vertices
   virtual std::vector<GVertex *> vertices() const;
@@ -211,7 +211,7 @@ public:
 
   // get bounds of parametric coordinate
   virtual Range<double> parBounds(int i) const = 0;
-  virtual Range<double> parBoundsOnFace(GFace *face = NULL) const
+  virtual Range<double> parBoundsOnFace(GFace *face = nullptr) const
   {
     return parBounds(0);
   }
diff --git a/Geo/GEdgeLoop.cpp b/Geo/GEdgeLoop.cpp
index 04dd91d0f5..c6a1b451fc 100644
--- a/Geo/GEdgeLoop.cpp
+++ b/Geo/GEdgeLoop.cpp
@@ -80,7 +80,7 @@ GEdgeSigned nextOne(GEdgeSigned *thisOne, std::list<GEdge *> &wire)
   }
 
   // should never end up here
-  return GEdgeSigned(0, 0);
+  return GEdgeSigned(0, nullptr);
 }
 
 int GEdgeLoop::count(GEdge *ge) const
@@ -112,14 +112,14 @@ void GEdgeLoop::getSigns(std::vector<int> &signs) const
 static void loopTheLoop(std::list<GEdge *> &wire, std::list<GEdgeSigned> &loop,
                         GEdge **degeneratedToInsert)
 {
-  GEdgeSigned *prevOne = 0;
-  GEdgeSigned ges(0, 0);
+  GEdgeSigned *prevOne = nullptr;
+  GEdgeSigned ges(0, nullptr);
 
   while(wire.size()) {
     if(prevOne && (*degeneratedToInsert) &&
        (*degeneratedToInsert)->getBeginVertex() == prevOne->getEndVertex()) {
       ges = GEdgeSigned(1, *degeneratedToInsert);
-      *degeneratedToInsert = 0;
+      *degeneratedToInsert = nullptr;
       // printf("second degenerated edge inserted\n");
     }
     else
@@ -141,7 +141,7 @@ GEdgeLoop::GEdgeLoop(const std::vector<GEdge *> &cwire)
   // gmsh
   std::list<GEdge *> wire;
   std::vector<GEdge *> degenerated;
-  GEdge *degeneratedToInsert = 0;
+  GEdge *degeneratedToInsert = nullptr;
   for(auto it = cwire.begin(); it != cwire.end(); ++it) {
     GEdge *ed = *it;
     if(ed->degenerate(0))
diff --git a/Geo/GEntity.cpp b/Geo/GEntity.cpp
index a6e22f1505..77a13498d3 100644
--- a/Geo/GEntity.cpp
+++ b/Geo/GEntity.cpp
@@ -18,7 +18,7 @@
 
 GEntity::GEntity(GModel *m, int t)
   : _model(m), _tag(t), _meshMaster(this), _visible(1), _selection(0),
-    _allElementsVisible(1), _obb(0), va_lines(0), va_triangles(0)
+    _allElementsVisible(1), _obb(nullptr), va_lines(nullptr), va_triangles(nullptr)
 {
   _color = CTX::instance()->packColor(0, 0, 255, 0);
 }
@@ -26,9 +26,9 @@ GEntity::GEntity(GModel *m, int t)
 void GEntity::deleteVertexArrays()
 {
   if(va_lines) delete va_lines;
-  va_lines = 0;
+  va_lines = nullptr;
   if(va_triangles) delete va_triangles;
-  va_triangles = 0;
+  va_triangles = nullptr;
 }
 
 char GEntity::getVisibility()
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index baef6682d5..053e01bf7b 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -196,7 +196,7 @@ public:
   virtual int dim() const { return -1; }
 
   // returns the parent entity for partitioned entities
-  virtual GEntity *getParentEntity() { return 0; }
+  virtual GEntity *getParentEntity() { return nullptr; }
 
   // regions that bound this entity or that this entity bounds.
   virtual std::list<GRegion *> regions() const
@@ -265,7 +265,7 @@ public:
   virtual ModelType getNativeType() const { return UnknownModel; }
 
   // get the native pointer of the particular representation
-  virtual void *getNativePtr() const { return 0; }
+  virtual void *getNativePtr() const { return nullptr; }
 
   // get the native id (int) of the particular representation
   virtual int getNativeInt() const { return 0; }
@@ -351,15 +351,18 @@ public:
   virtual void getNumMeshElements(unsigned *const c) const {}
 
   // get the start of the array of a type of element
-  virtual MElement *const *getStartElementType(int type) const { return 0; }
+  virtual MElement *const *getStartElementType(int type) const
+  {
+    return nullptr;
+  }
 
   // get the element at the given index
-  virtual MElement *getMeshElement(std::size_t index) const { return 0; }
+  virtual MElement *getMeshElement(std::size_t index) const { return nullptr; }
   // get the element at the given index for a given familyType
   virtual MElement *getMeshElementByType(const int familyType,
                                          const std::size_t index) const
   {
-    return 0;
+    return nullptr;
   }
 
   // get/set all mesh element visibility flag
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 8e34965d8c..afaf675514 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -34,7 +34,7 @@
 #endif
 
 GFace::GFace(GModel *model, int tag)
-  : GEntity(model, tag), r1(0), r2(0), va_geom_triangles(0), compoundSurface(0)
+  : GEntity(model, tag), r1(nullptr), r2(nullptr), va_geom_triangles(nullptr), compoundSurface(nullptr)
 {
   meshStatistics.status = GFace::PENDING;
   meshStatistics.refineAllEdges = false;
@@ -180,16 +180,16 @@ MElement *const *GFace::getStartElementType(int type) const
 {
   switch(type) {
   case 0:
-    if(triangles.empty()) return 0; // msvc would throw an exception
+    if(triangles.empty()) return nullptr; // msvc would throw an exception
     return reinterpret_cast<MElement *const *>(&triangles[0]);
   case 1:
-    if(quadrangles.empty()) return 0; // msvc would throw an exception
+    if(quadrangles.empty()) return nullptr; // msvc would throw an exception
     return reinterpret_cast<MElement *const *>(&quadrangles[0]);
   case 2:
-    if(polygons.empty()) return 0;
+    if(polygons.empty()) return nullptr;
     return reinterpret_cast<MElement *const *>(&polygons[0]);
   }
-  return 0;
+  return nullptr;
 }
 
 MElement *GFace::getMeshElement(std::size_t index) const
@@ -200,7 +200,7 @@ MElement *GFace::getMeshElement(std::size_t index) const
     return quadrangles[index - triangles.size()];
   else if(index < triangles.size() + quadrangles.size() + polygons.size())
     return polygons[index - triangles.size() - quadrangles.size()];
-  return 0;
+  return nullptr;
 }
 
 MElement *GFace::getMeshElementByType(const int familyType,
@@ -213,7 +213,7 @@ MElement *GFace::getMeshElementByType(const int familyType,
   else if(familyType == TYPE_POLYG)
     return polygons[index];
 
-  return 0;
+  return nullptr;
 }
 
 void GFace::resetMeshAttributes()
@@ -223,7 +223,7 @@ void GFace::resetMeshAttributes()
   meshAttributes.method = MESH_UNSTRUCTURED;
   meshAttributes.transfiniteArrangement = 0;
   meshAttributes.transfiniteSmoothing = -1;
-  meshAttributes.extrude = 0;
+  meshAttributes.extrude = nullptr;
   meshAttributes.reverseMesh = false;
   meshAttributes.meshSize = MAX_LC;
   meshAttributes.meshSizeFactor = 1.;
@@ -1103,7 +1103,7 @@ private:
 public:
   data_wrapper()
   {
-    gf = NULL;
+    gf = nullptr;
     point = SPoint3();
   }
   ~data_wrapper() {}
@@ -1193,7 +1193,7 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint,
     data_wrapper w;
     w.set_point(queryPoint);
     w.set_face(this);
-    minlbfgsoptimize(state, bfgs_callback, NULL, &w);
+    minlbfgsoptimize(state, bfgs_callback, nullptr, &w);
 
     // Get results
     alglib::minlbfgsreport rep;
@@ -1768,7 +1768,7 @@ void GFace::setMeshMaster(GFace *master, const std::vector<double> &tfo)
     SPoint3 xyzTfo((*mvIter)->x(), (*mvIter)->y(), (*mvIter)->z());
     xyzTfo.transform(tfo);
 
-    GVertex *l_vertex = NULL;
+    GVertex *l_vertex = nullptr;
 
     double dist_min = 1.e22;
     for(auto lvIter = l_vertices.begin(); lvIter != l_vertices.end();
@@ -1782,7 +1782,7 @@ void GFace::setMeshMaster(GFace *master, const std::vector<double> &tfo)
       }
     }
 
-    if(l_vertex == NULL) {
+    if(l_vertex == nullptr) {
       Msg::Error("No corresponding point %d for periodic connection of surface "
                  "%d to %d (min. distance = %g, tolerance = %g)",
                  m_vertex->tag(), master->tag(), tag(), dist_min,
@@ -1813,7 +1813,7 @@ void GFace::setMeshMaster(GFace *master, const std::vector<double> &tfo)
     std::pair<GVertex *, GVertex *> backward(forward.second, forward.first);
     int numb = m_vtxToEdge.count(backward);
     int sign = 0;
-    GEdge *masterEdge = 0;
+    GEdge *masterEdge = nullptr;
 
     // unique matches
     if(!masterEdge && numf == 1 &&
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 065bb30538..70390e6d3a 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -73,12 +73,12 @@ public:
   void delRegion(GRegion *r)
   {
     if(r1 == r) r1 = r2;
-    r2 = NULL;
+    r2 = nullptr;
   }
   GRegion *getRegion(int const num) const { return num == 0 ? r1 : r2; }
 
   // get number of regions
-  std::size_t numRegions() const { return (r1 != NULL) + (r2 != NULL); }
+  std::size_t numRegions() const { return (r1 != nullptr) + (r2 != nullptr); }
 
   std::list<GRegion *> regions() const
   {
@@ -109,8 +109,7 @@ public:
 
   bool containsEdge(int const iEdge) const
   {
-    for(std::vector<GEdge *>::const_iterator it = l_edges.begin();
-        it != l_edges.end(); ++it)
+    for(auto it = l_edges.begin(); it != l_edges.end(); ++it)
       if((*it)->tag() == iEdge) return true;
     return false;
   }
@@ -139,7 +138,7 @@ public:
   virtual int dim() const { return 2; }
 
   // returns the parent entity for partitioned entities
-  virtual GEntity *getParentEntity() { return 0; }
+  virtual GEntity *getParentEntity() { return nullptr; }
 
   // set visibility flag
   virtual void setVisibility(char val, bool recursive = false);
@@ -283,8 +282,8 @@ public:
   // add points (and optionally normals) in vectors so that two
   // points are at most maxDist apart
   bool fillPointCloud(double maxDist, std::vector<SPoint3> *points,
-                      std::vector<SPoint2> *uvpoints = 0,
-                      std::vector<SVector3> *normals = 0);
+                      std::vector<SPoint2> *uvpoints = nullptr,
+                      std::vector<SVector3> *normals = nullptr);
 
   // tells if it's a sphere, and if it is, returns parameters
   virtual bool isSphere(double &radius, SPoint3 &center) const { return false; }
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 94cd1d2f52..3ef69115e3 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -65,10 +65,10 @@ std::vector<GModel *> GModel::list;
 int GModel::_current = -1;
 
 GModel::GModel(const std::string &name)
-  : _destroying(false), _name(name), _visible(1), _elementOctree(0),
-    _geo_internals(0), _occ_internals(0), _acis_internals(0),
-    _parasolid_internals(0), _fields(0), _currentMeshEntity(0),
-    _numPartitions(0), normals(0)
+  : _destroying(false), _name(name), _visible(1), _elementOctree(nullptr),
+    _geo_internals(nullptr), _occ_internals(nullptr), _acis_internals(nullptr),
+    _parasolid_internals(nullptr), _fields(nullptr), _currentMeshEntity(nullptr),
+    _numPartitions(0), normals(nullptr)
 {
   _maxVertexNum = CTX::instance()->mesh.firstNodeTag - 1;
   _maxElementNum = CTX::instance()->mesh.firstElementTag - 1;
@@ -156,7 +156,7 @@ GModel *GModel::findByName(const std::string &name, const std::string &fileName)
     if(list[i]->getName() == name &&
        (fileName.empty() || !list[i]->hasFileName(fileName)))
       return list[i];
-  return 0;
+  return nullptr;
 }
 
 void GModel::destroy(bool keepName)
@@ -173,7 +173,7 @@ void GModel::destroy(bool keepName)
   _maxElementNum = CTX::instance()->mesh.firstElementTag - 1;
   _checkPointedMaxVertexNum = _maxVertexNum;
   _checkPointedMaxElementNum = _maxElementNum;
-  _currentMeshEntity = 0;
+  _currentMeshEntity = nullptr;
   _lastMeshEntityError.clear();
   _lastMeshVertexError.clear();
 
@@ -198,7 +198,7 @@ void GModel::destroy(bool keepName)
   resetOCCInternals();
 
   if(normals) delete normals;
-  normals = 0;
+  normals = nullptr;
 
 #if defined(HAVE_MESH)
   _fields->reset();
@@ -221,7 +221,7 @@ void GModel::destroyMeshCaches()
   _elementIndexCache.clear();
   std::map<int, int>().swap(_elementIndexCache);
   delete _elementOctree;
-  _elementOctree = 0;
+  _elementOctree = nullptr;
 }
 
 void GModel::deleteMesh()
@@ -231,7 +231,7 @@ void GModel::deleteMesh()
   for(auto it = firstEdge(); it != lastEdge(); ++it) (*it)->deleteMesh();
   for(auto it = firstVertex(); it != lastVertex(); ++it) (*it)->deleteMesh();
   destroyMeshCaches();
-  _currentMeshEntity = 0;
+  _currentMeshEntity = nullptr;
   _lastMeshEntityError.clear();
   _lastMeshVertexError.clear();
 }
@@ -282,7 +282,7 @@ void GModel::deleteMesh(const std::vector<GEntity *> &entities)
     }
   }
   destroyMeshCaches();
-  _currentMeshEntity = 0;
+  _currentMeshEntity = nullptr;
   _lastMeshEntityError.clear();
   _lastMeshVertexError.clear();
 }
@@ -311,7 +311,7 @@ GRegion *GModel::getRegionByTag(int n) const
   if(it != regions.end())
     return *it;
   else
-    return 0;
+    return nullptr;
 }
 
 GFace *GModel::getFaceByTag(int n) const
@@ -321,7 +321,7 @@ GFace *GModel::getFaceByTag(int n) const
   if(it != faces.end())
     return *it;
   else
-    return 0;
+    return nullptr;
 }
 
 GEdge *GModel::getEdgeByTag(int n) const
@@ -331,7 +331,7 @@ GEdge *GModel::getEdgeByTag(int n) const
   if(it != edges.end())
     return *it;
   else
-    return 0;
+    return nullptr;
 }
 
 GVertex *GModel::getVertexByTag(int n) const
@@ -341,7 +341,7 @@ GVertex *GModel::getVertexByTag(int n) const
   if(it != vertices.end())
     return *it;
   else
-    return 0;
+    return nullptr;
 }
 
 GEntity *GModel::getEntityByTag(int dim, int n) const
@@ -352,7 +352,7 @@ GEntity *GModel::getEntityByTag(int dim, int n) const
   case 2: return getFaceByTag(n);
   case 3: return getRegionByTag(n);
   }
-  return 0;
+  return nullptr;
 }
 
 std::vector<int> GModel::getTagsForPhysicalName(int dim,
@@ -1236,14 +1236,14 @@ int GModel::adaptMesh(std::vector<int> technique,
             laplaceSmoothing(*fit, CTX::instance()->mesh.nbSmoothing);
           }
           if(_elementOctree) delete _elementOctree;
-          _elementOctree = 0;
+          _elementOctree = nullptr;
         }
       }
       else if(getDim() == 3) {
         for(auto rit = firstRegion(); rit != lastRegion(); ++rit) {
           refineMeshMMG(*rit);
           if(_elementOctree) delete _elementOctree;
-          _elementOctree = 0;
+          _elementOctree = nullptr;
         }
       }
 
@@ -1650,7 +1650,7 @@ void GModel::rebuildMeshVertexCache(bool onlyIfNecessary)
     getEntities(entities);
     if(dense) {
       // numbering starts at 1
-      _vertexVectorCache.resize(_maxVertexNum + 1, (MVertex *)0);
+      _vertexVectorCache.resize(_maxVertexNum + 1, (MVertex *)nullptr);
       for(std::size_t i = 0; i < entities.size(); i++)
         for(std::size_t j = 0; j < entities[i]->mesh_vertices.size(); j++)
           _vertexVectorCache[entities[i]->mesh_vertices[j]->getNum()] =
@@ -1686,7 +1686,7 @@ void GModel::rebuildMeshElementCache(bool onlyIfNecessary)
     getEntities(entities);
     if(dense) {
       // numbering starts at 1
-      _elementVectorCache.resize(_maxElementNum + 1, (MElement *)0);
+      _elementVectorCache.resize(_maxElementNum + 1, (MElement *)nullptr);
       for(std::size_t i = 0; i < entities.size(); i++)
         for(std::size_t j = 0; j < entities[i]->getNumMeshElements(); j++) {
           MElement *e = entities[i]->getMeshElement(j);
@@ -2015,7 +2015,7 @@ void GModel::_storeElementsInEntities(
     case TYPE_LIN: {
       GEdge *e = getEdgeByTag(it->first);
       if(!e) {
-        e = new discreteEdge(this, it->first, 0, 0);
+        e = new discreteEdge(this, it->first, nullptr, nullptr);
         add(e);
       }
       _addElements(e->lines, it->second);
@@ -2237,7 +2237,7 @@ void GModel::_storeVerticesInEntities(std::vector<MVertex *> &vertices)
         ge->mesh_vertices.push_back(v);
       else {
         delete v; // we delete all unused vertices
-        vertices[i] = 0;
+        vertices[i] = nullptr;
       }
     }
   }
@@ -2253,7 +2253,7 @@ void GModel::pruneMeshVertexAssociations()
       MElement *e = entities[i]->getMeshElement(j);
       for(std::size_t k = 0; k < e->getNumVertices(); k++) {
         MVertex *v = e->getVertex(k);
-        v->setEntity(0);
+        v->setEntity(nullptr);
         vertSet.insert(v);
       }
     }
@@ -2289,7 +2289,7 @@ void GModel::_storePhysicalTagsInEntities(
 {
   std::map<int, std::map<int, std::string> >::const_iterator it = map.begin();
   for(; it != map.end(); ++it) {
-    GEntity *ge = 0;
+    GEntity *ge = nullptr;
     switch(dim) {
     case 0: ge = getVertexByTag(it->first); break;
     case 1: ge = getEdgeByTag(it->first); break;
@@ -2513,7 +2513,7 @@ void GModel::alignPeriodicBoundaries()
     GEdge *tgt = *it;
     GEdge *src = dynamic_cast<GEdge *>(tgt->getMeshMaster());
 
-    if(src != NULL && src != tgt) {
+    if(src != nullptr && src != tgt) {
       // compose a search list on master edge
 
       std::map<MEdge, MLine *, MEdgeLessThan> srcLines;
@@ -2589,7 +2589,7 @@ void GModel::alignPeriodicBoundaries()
   for(auto it = firstFace(); it != lastFace(); ++it) {
     GFace *tgt = *it;
     GFace *src = dynamic_cast<GFace *>(tgt->getMeshMaster());
-    if(src != NULL && src != tgt) {
+    if(src != nullptr && src != tgt) {
       std::map<MFace, MElement *, MFaceLessThan> srcElmts;
 
       for(std::size_t i = 0; i < src->getNumMeshElements(); ++i) {
diff --git a/Geo/GModelCreateTopologyFromMesh.cpp b/Geo/GModelCreateTopologyFromMesh.cpp
index 44ea328aa8..c7d6217abd 100644
--- a/Geo/GModelCreateTopologyFromMesh.cpp
+++ b/Geo/GModelCreateTopologyFromMesh.cpp
@@ -60,7 +60,7 @@ std::vector<GEdge *> ensureSimplyConnectedEdge(GEdge *ge)
       auto it = _conn.find(ge->lines[i]->getVertex(j));
       if(it == _conn.end())
         _conn[ge->lines[i]->getVertex(j)] =
-          std::make_pair(ge->lines[i], (MLine *)NULL);
+          std::make_pair(ge->lines[i], (MLine *)nullptr);
       else
         it->second.second = ge->lines[i];
     }
@@ -105,7 +105,7 @@ std::vector<GEdge *> ensureSimplyConnectedEdge(GEdge *ge)
       ge->lines = _parts[i];
     else {
       discreteEdge *newE = new discreteEdge(
-        ge->model(), ge->model()->getMaxElementaryNumber(1) + 1, NULL, NULL);
+        ge->model(), ge->model()->getMaxElementaryNumber(1) + 1, nullptr, nullptr);
       ge->model()->add(newE);
       newE->lines = _parts[i];
       _all.push_back(newE);
@@ -141,7 +141,7 @@ void ensureManifoldFace(GFace *gf)
       if(_nonManifold.find(ed) == _nonManifold.end()) {
         auto it = _pairs.find(ed);
         if(it == _pairs.end()) {
-          _pairs[ed] = std::make_pair(e, (MElement *)NULL);
+          _pairs[ed] = std::make_pair(e, (MElement *)nullptr);
         }
         else {
           if(it->second.second == NULL) { it->second.second = e; }
@@ -421,7 +421,7 @@ void createTopologyFromMesh2D(GModel *gm, int &num)
     auto gfIter = gFacesToGEdge.find(gfaces);
     if(gfIter == gFacesToGEdge.end()) {
       discreteEdge *de =
-        new discreteEdge(gm, gm->getMaxElementaryNumber(1) + 1, NULL, NULL);
+        new discreteEdge(gm, gm->getMaxElementaryNumber(1) + 1, nullptr, nullptr);
       num++;
       gm->add(de);
       auto gfIter = gfaces.begin();
@@ -579,7 +579,7 @@ void createTopologyFromMesh3D(GModel *gm, int &num)
         else {
           auto frIter = tFaceToGRegionPair.find(f);
           if(frIter == tFaceToGRegionPair.end()) {
-            tFaceToGRegionPair[f] = std::make_pair((GRegion *)NULL, *it);
+            tFaceToGRegionPair[f] = std::make_pair((GRegion *)nullptr, *it);
           }
           else
             frIter->second.first = gr;
diff --git a/Geo/GModelIO_ACTRAN.cpp b/Geo/GModelIO_ACTRAN.cpp
index 9dfa1046c1..e18b4e74ca 100644
--- a/Geo/GModelIO_ACTRAN.cpp
+++ b/Geo/GModelIO_ACTRAN.cpp
@@ -77,7 +77,7 @@ int GModel::readACTRAN(const std::string &name)
           sscanf(buffer, "%d %lf %lf %lf", &num, &x, &y, &z);
         else
           sscanf(buffer, "%d %lf %lf", &num, &x, &y);
-        _vertexMapCache[num] = new MVertex(x, y, z, 0, num);
+        _vertexMapCache[num] = new MVertex(x, y, z, nullptr, num);
       }
     }
     else if(!strcmp(str, "BEGIN") && !strcmp(str2, "ELEMENT")) {
diff --git a/Geo/GModelIO_BDF.cpp b/Geo/GModelIO_BDF.cpp
index c522368cc2..5b8bc6c0a5 100644
--- a/Geo/GModelIO_BDF.cpp
+++ b/Geo/GModelIO_BDF.cpp
@@ -228,7 +228,7 @@ int GModel::readBDF(const std::string &name)
         int num;
         double x, y, z;
         if(!readVertexBDF(fp, buffer, 4, &num, &x, &y, &z)) break;
-        vertexMap[num] = new MVertex(x, y, z, 0, num);
+        vertexMap[num] = new MVertex(x, y, z, nullptr, num);
       }
     }
   }
diff --git a/Geo/GModelIO_CGNS.cpp b/Geo/GModelIO_CGNS.cpp
index c75fe9aa9c..342fb44a20 100644
--- a/Geo/GModelIO_CGNS.cpp
+++ b/Geo/GModelIO_CGNS.cpp
@@ -114,7 +114,7 @@ int GModel::readCGNS(const std::string &name,
                              geomName2Phys);
   // destroy all zones
   for(std::size_t iZone = 0; iZone < allZones.size(); iZone++) {
-    if(allZones[iZone] != 0) delete allZones[iZone];
+    if(allZones[iZone] != nullptr) delete allZones[iZone];
   }
 
   // reconstruct geometrical topology if required
diff --git a/Geo/GModelIO_DIFF.cpp b/Geo/GModelIO_DIFF.cpp
index 91ddba72d6..eaa525892c 100644
--- a/Geo/GModelIO_DIFF.cpp
+++ b/Geo/GModelIO_DIFF.cpp
@@ -54,7 +54,7 @@ int GModel::readDIFF(const std::string &name)
   std::vector<MVertex *> vertexVector;
 
   {
-    while(strstr(str, "Number of space dim. =") == NULL) {
+    while(strstr(str, "Number of space dim. =") == nullptr) {
       if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
     }
 
@@ -70,7 +70,7 @@ int GModel::readDIFF(const std::string &name)
       fclose(fp);
       return 0;
     }
-    while(strstr(str, "Number of elements   =") == NULL) {
+    while(strstr(str, "Number of elements   =") == nullptr) {
       if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
     }
     if(sscanf(str, "%*s %*s %*s %*s %d", &numElements) != 1) {
@@ -84,7 +84,7 @@ int GModel::readDIFF(const std::string &name)
       fclose(fp);
       return 0;
     }
-    while(strstr(str, "Number of nodes      =") == NULL) {
+    while(strstr(str, "Number of nodes      =") == nullptr) {
       if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
     }
     if(sscanf(str, "%*s %*s %*s %*s %d", &numVertices) != 1) {
@@ -98,7 +98,7 @@ int GModel::readDIFF(const std::string &name)
       fclose(fp);
       return 0;
     }
-    while(strstr(str, "Max number of nodes in an element:") == NULL) {
+    while(strstr(str, "Max number of nodes in an element:") == nullptr) {
       if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
     }
     if(sscanf(str, "%*s %*s %*s %*s %*s %*s %*s %d", &numVerticesPerElement) !=
@@ -131,8 +131,8 @@ int GModel::readDIFF(const std::string &name)
       fclose(fp);
       return 0;
     }
-    while(strstr(str, "Boundary indicators:") == NULL &&
-          strstr(str, "boundary indicators:") == NULL) {
+    while(strstr(str, "Boundary indicators:") == nullptr &&
+          strstr(str, "boundary indicators:") == nullptr) {
       if(!fgets(str, sizeof(str), fp) || feof(fp)) break;
     }
     if(sscanf(str, "%d %*s %*s", &nbi) != 1) {
@@ -185,7 +185,7 @@ int GModel::readDIFF(const std::string &name)
       if(vertexMap.count(num))
         Msg::Warning("Skipping duplicate node %d", num);
       else
-        vertexMap[num] = new MVertex(xyz[0], xyz[1], xyz[2], 0, num);
+        vertexMap[num] = new MVertex(xyz[0], xyz[1], xyz[2], nullptr, num);
       if(numVertices > 100000) Msg::ProgressMeter(i + 1, true, "Reading nodes");
       // If the vertex numbering is dense, tranfer the map into a
       // vector to speed up element creation
@@ -195,9 +195,9 @@ int GModel::readDIFF(const std::string &name)
         Msg::Info("Vertex numbering is dense");
         vertexVector.resize(vertexMap.size() + 1);
         if(minVertex == 1)
-          vertexVector[0] = 0;
+          vertexVector[0] = nullptr;
         else
-          vertexVector[numVertices] = 0;
+          vertexVector[numVertices] = nullptr;
         std::map<int, MVertex *>::const_iterator it = vertexMap.begin();
         for(; it != vertexMap.end(); ++it) vertexVector[it->first] = it->second;
         vertexMap.clear();
diff --git a/Geo/GModelIO_GEO.cpp b/Geo/GModelIO_GEO.cpp
index df3f789b73..298f63d8a8 100644
--- a/Geo/GModelIO_GEO.cpp
+++ b/Geo/GModelIO_GEO.cpp
@@ -163,7 +163,7 @@ bool GEO_Internals::addLine(int &tag, const std::vector<int> &pointTags)
     List_Add(tmp, &t);
   }
   bool ok = true;
-  Curve *c = CreateCurve(tag, MSH_SEGM_LINE, 1, tmp, NULL, -1, -1, 0., 1., ok);
+  Curve *c = CreateCurve(tag, MSH_SEGM_LINE, 1, tmp, nullptr, -1, -1, 0., 1., ok);
   Tree_Add(Curves, &c);
   CreateReversedCurve(c);
   List_Delete(tmp);
@@ -184,7 +184,7 @@ bool GEO_Internals::addCircleArc(int &tag, int startTag, int centerTag,
   List_Add(tmp, &centerTag);
   List_Add(tmp, &endTag);
   bool ok = true;
-  Curve *c = CreateCurve(tag, MSH_SEGM_CIRC, 2, tmp, NULL, -1, -1, 0., 1., ok);
+  Curve *c = CreateCurve(tag, MSH_SEGM_CIRC, 2, tmp, nullptr, -1, -1, 0., 1., ok);
   if(nx || ny || nz) {
     c->Circle.n[0] = nx;
     c->Circle.n[1] = ny;
@@ -219,7 +219,7 @@ bool GEO_Internals::addEllipseArc(int &tag, int startTag, int centerTag,
   List_Add(tmp, &majorTag);
   List_Add(tmp, &endTag);
   bool ok = true;
-  Curve *c = CreateCurve(tag, MSH_SEGM_ELLI, 2, tmp, NULL, -1, -1, 0., 1., ok);
+  Curve *c = CreateCurve(tag, MSH_SEGM_ELLI, 2, tmp, nullptr, -1, -1, 0., 1., ok);
   if(nx || ny || nz) {
     c->Circle.n[0] = nx;
     c->Circle.n[1] = ny;
@@ -256,7 +256,7 @@ bool GEO_Internals::addSpline(int &tag, const std::vector<int> &pointTags)
     List_Add(tmp, &t);
   }
   bool ok = true;
-  Curve *c = CreateCurve(tag, MSH_SEGM_SPLN, 3, tmp, NULL, -1, -1, 0., 1., ok);
+  Curve *c = CreateCurve(tag, MSH_SEGM_SPLN, 3, tmp, nullptr, -1, -1, 0., 1., ok);
   Tree_Add(Curves, &c);
   CreateReversedCurve(c);
   List_Delete(tmp);
@@ -282,7 +282,7 @@ bool GEO_Internals::addBezier(int &tag, const std::vector<int> &pointTags)
   }
   bool ok = true;
   Curve *c =
-    CreateCurve(tag, MSH_SEGM_BEZIER, 2, tmp, NULL, -1, -1, 0., 1., ok);
+    CreateCurve(tag, MSH_SEGM_BEZIER, 2, tmp, nullptr, -1, -1, 0., 1., ok);
   Tree_Add(Curves, &c);
   CreateReversedCurve(c);
   List_Delete(tmp);
@@ -307,10 +307,10 @@ bool GEO_Internals::addBSpline(int &tag, const std::vector<int> &pointTags,
     int t = pointTags[i];
     List_Add(tmp, &t);
   }
-  Curve *c = NULL;
+  Curve *c = nullptr;
   bool ok = true;
   if(seqknots.empty()) {
-    c = CreateCurve(tag, MSH_SEGM_BSPLN, 2, tmp, NULL, -1, -1, 0., 1., ok);
+    c = CreateCurve(tag, MSH_SEGM_BSPLN, 2, tmp, nullptr, -1, -1, 0., 1., ok);
   }
   else {
     int order = seqknots.size() - pointTags.size() - 1;
@@ -370,9 +370,9 @@ bool GEO_Internals::_addCompoundSpline(int &tag,
   bool ok = true;
   Curve *c;
   if(bspline)
-    c = CreateCurve(tag, MSH_SEGM_BSPLN, 2, tmp, NULL, -1, -1, 0., 1., ok);
+    c = CreateCurve(tag, MSH_SEGM_BSPLN, 2, tmp, nullptr, -1, -1, 0., 1., ok);
   else // often too oscillatory for non-uniform distribution of control points
-    c = CreateCurve(tag, MSH_SEGM_SPLN, 3, tmp, NULL, -1, -1, 0., 1., ok);
+    c = CreateCurve(tag, MSH_SEGM_SPLN, 3, tmp, nullptr, -1, -1, 0., 1., ok);
   Tree_Add(Curves, &c);
   CreateReversedCurve(c);
   List_Delete(tmp);
@@ -448,7 +448,7 @@ static bool SortCurvesConsecutive(const std::vector<Curve *> &e,
 
     auto it0 = c.find(v0), it1 = c.find(v1);
     if(it0 == c.end())
-      c[v0] = std::make_pair(v1, (Vertex *)NULL);
+      c[v0] = std::make_pair(v1, (Vertex *)nullptr);
     else {
       if(it0->second.second == NULL) { it0->second.second = v1; }
       else {
@@ -457,7 +457,7 @@ static bool SortCurvesConsecutive(const std::vector<Curve *> &e,
       }
     }
     if(it1 == c.end())
-      c[v1] = std::make_pair(v0, (Vertex *)NULL);
+      c[v1] = std::make_pair(v0, (Vertex *)nullptr);
     else {
       if(it1->second.second == NULL) { it1->second.second = v0; }
       else {
@@ -471,7 +471,7 @@ static bool SortCurvesConsecutive(const std::vector<Curve *> &e,
 
   while(!c.empty()) {
     std::vector<Vertex *> v;
-    Vertex *start = NULL;
+    Vertex *start = nullptr;
     {
       auto it = c.begin();
       start = it->first;
@@ -496,7 +496,7 @@ static bool SortCurvesConsecutive(const std::vector<Curve *> &e,
       }
       v.push_back(current);
       auto it = c.find(current);
-      if(it == c.end() || it->first == NULL) {
+      if(it == c.end() || it->first == nullptr) {
         Msg::Error("Impossible to find point %d", current->Num);
         return false;
       }
@@ -513,7 +513,7 @@ static bool SortCurvesConsecutive(const std::vector<Curve *> &e,
       }
       prev = temp;
       if(current == start) { v.push_back(current); }
-    } while(current != start && current != NULL);
+    } while(current != start && current != nullptr);
     if(v.size() > 2 && v[v.size() - 2] == v[v.size() - 1]) {
       v.erase(v.begin() + v.size() - 1);
     }
@@ -1432,7 +1432,7 @@ void GEO_Internals::synchronize(GModel *model, bool resetMeshAttributes)
           model->add(e);
         }
         else if(!e) {
-          e = new gmshEdge(model, c, 0, 0);
+          e = new gmshEdge(model, c, nullptr, nullptr);
           model->add(e);
         }
         else {
@@ -1442,7 +1442,7 @@ void GEO_Internals::synchronize(GModel *model, bool resetMeshAttributes)
                 ->resetNativePtr(c, model->getVertexByTag(c->beg->Num),
                                  model->getVertexByTag(c->end->Num));
             else
-              ((gmshEdge *)e)->resetNativePtr(c, 0, 0);
+              ((gmshEdge *)e)->resetNativePtr(c, nullptr, nullptr);
           }
           if(resetMeshAttributes) e->resetMeshAttributes();
         }
@@ -1502,7 +1502,7 @@ void GEO_Internals::synchronize(GModel *model, bool resetMeshAttributes)
     for(int j = 0; j < List_Nbr(p->Entities); j++) {
       int num;
       List_Read(p->Entities, j, &num);
-      GEntity *ge = 0;
+      GEntity *ge = nullptr;
       int tag = CTX::instance()->geom.orientedPhysicals ? abs(num) : num;
       switch(p->Typ) {
       case MSH_PHYSICAL_POINT: ge = model->getVertexByTag(tag); break;
@@ -1527,7 +1527,7 @@ void GEO_Internals::synchronize(GModel *model, bool resetMeshAttributes)
     std::vector<GEntity *> ents;
     for(std::size_t i = 0; i < compound.size(); i++) {
       int tag = compound[i];
-      GEntity *ent = NULL;
+      GEntity *ent = nullptr;
       switch(dim) {
       case 1: ent = model->getEdgeByTag(tag); break;
       case 2: ent = model->getFaceByTag(tag); break;
@@ -1569,12 +1569,12 @@ gmshSurface *GEO_Internals::newGeometrySphere(int tag, int centerTag,
   Vertex *v1 = FindPoint(centerTag);
   if(!v1) {
     Msg::Error("Unknown sphere center point %d", centerTag);
-    return 0;
+    return nullptr;
   }
   Vertex *v2 = FindPoint(pointTag);
   if(!v2) {
     Msg::Error("Unknown sphere point %d", pointTag);
-    return 0;
+    return nullptr;
   }
   return gmshSphere::NewSphere(
     tag, v1->Pos.X, v1->Pos.Y, v1->Pos.Z,
@@ -1589,12 +1589,12 @@ gmshSurface *GEO_Internals::newGeometryPolarSphere(int tag, int centerTag,
   Vertex *v1 = FindPoint(centerTag);
   if(!v1) {
     Msg::Error("Unknown polar sphere center point %d", centerTag);
-    return 0;
+    return nullptr;
   }
   Vertex *v2 = FindPoint(pointTag);
   if(!v2) {
     Msg::Error("Unknown polar sphere point %d", pointTag);
-    return 0;
+    return nullptr;
   }
   return gmshPolarSphere::NewPolarSphere(
     tag, v1->Pos.X, v1->Pos.Y, v1->Pos.Z,
@@ -1610,7 +1610,7 @@ void GModel::createGEOInternals() { _geo_internals = new GEO_Internals; }
 void GModel::deleteGEOInternals()
 {
   if(_geo_internals) delete _geo_internals;
-  _geo_internals = 0;
+  _geo_internals = nullptr;
 }
 
 #if defined(HAVE_MESH)
@@ -1812,7 +1812,7 @@ int GModel::exportDiscreteGEOInternals()
   for(auto it = firstEdge(); it != lastEdge(); it++) {
     if((*it)->geomType() == GEntity::DiscreteCurve) {
       bool ok = true;
-      Curve *c = CreateCurve((*it)->tag(), MSH_SEGM_DISCRETE, 1, NULL, NULL, -1,
+      Curve *c = CreateCurve((*it)->tag(), MSH_SEGM_DISCRETE, 1, nullptr, nullptr, -1,
                              -1, 0., 1., ok);
       c->Control_Points = List_Create(2, 1, sizeof(Vertex *));
       GVertex *gvb = (*it)->getBeginVertex();
diff --git a/Geo/GModelIO_MED.cpp b/Geo/GModelIO_MED.cpp
index 8e582e9a47..00e538276f 100644
--- a/Geo/GModelIO_MED.cpp
+++ b/Geo/GModelIO_MED.cpp
@@ -392,7 +392,7 @@ int GModel::readMED(const std::string &name, int meshIndex)
   for(int i = 0; i < numNodes; i++)
     verts[i] = new MVertex(coord[spaceDim * i],
                            (spaceDim > 1) ? coord[spaceDim * i + 1] : 0.,
-                           (spaceDim > 2) ? coord[spaceDim * i + 2] : 0., 0,
+                           (spaceDim > 2) ? coord[spaceDim * i + 2] : 0., nullptr,
                            nodeTags.empty() ? 0 : nodeTags[i]);
 
   std::vector<med_int> nodeFamily(numNodes, 0);
@@ -637,7 +637,7 @@ static void writeElementsMED(med_idt &fid, char *meshName,
 #if(MED_MAJOR_NUM >= 3)
   if(MEDmeshElementWr(fid, meshName, MED_NO_DT, MED_NO_IT, 0., MED_CELL, type,
                       MED_NODAL, MED_FULL_INTERLACE, (med_int)fam.size(),
-                      &conn[0], MED_FALSE, 0, MED_TRUE, &tags[0], MED_TRUE,
+                      &conn[0], MED_FALSE, nullptr, MED_TRUE, &tags[0], MED_TRUE,
                       &fam[0]) < 0)
 #else
   if(MEDelementsEcr(fid, meshName, (med_int)3, &conn[0], MED_FULL_INTERLACE, 0,
diff --git a/Geo/GModelIO_MSH2.cpp b/Geo/GModelIO_MSH2.cpp
index 4c88402aa3..39ffd0f126 100644
--- a/Geo/GModelIO_MSH2.cpp
+++ b/Geo/GModelIO_MSH2.cpp
@@ -71,8 +71,8 @@ createElementMSH2(GModel *m, int num, int typeMSH, int physical, int reg,
                   unsigned int part, std::vector<MVertex *> &v,
                   std::map<int, std::vector<MElement *> > elements[10],
                   std::map<int, std::map<int, std::string> > physicals[4],
-                  bool owner = false, MElement *parent = 0, MElement *d1 = 0,
-                  MElement *d2 = 0)
+                  bool owner = false, MElement *parent = nullptr, MElement *d1 = nullptr,
+                  MElement *d2 = nullptr)
 {
   if(CTX::instance()->mesh.switchElementTags) {
     int tmp = reg;
@@ -85,7 +85,7 @@ createElementMSH2(GModel *m, int num, int typeMSH, int physical, int reg,
 
   if(!e) {
     Msg::Error("Unknown type of element %d", typeMSH);
-    return NULL;
+    return nullptr;
   }
 
   int dim = e->getDim();
@@ -197,7 +197,7 @@ int GModel::_readMSH2(const std::string &name)
       for(int i = 0; i < numVertices; i++) {
         int num;
         double xyz[3], uv[2];
-        MVertex *newVertex = 0;
+        MVertex *newVertex = nullptr;
         if(!parametric) {
           if(!binary) {
             if(fscanf(fp, "%d %lf %lf %lf", &num, &xyz[0], &xyz[1], &xyz[2]) !=
@@ -218,7 +218,7 @@ int GModel::_readMSH2(const std::string &name)
             }
             if(swap) SwapBytes((char *)xyz, sizeof(double), 3);
           }
-          newVertex = new MVertex(xyz[0], xyz[1], xyz[2], 0, num);
+          newVertex = new MVertex(xyz[0], xyz[1], xyz[2], nullptr, num);
         }
         else {
           int iClasDim, iClasTag;
@@ -313,9 +313,9 @@ int GModel::_readMSH2(const std::string &name)
         Msg::Debug("Vertex numbering is dense");
         vertexVector.resize(vertexMap.size() + 1);
         if(minVertex == 1)
-          vertexVector[0] = 0;
+          vertexVector[0] = nullptr;
         else
-          vertexVector[numVertices] = 0;
+          vertexVector[numVertices] = nullptr;
         std::map<int, MVertex *>::const_iterator it = vertexMap.begin();
         for(; it != vertexMap.end(); ++it) vertexVector[it->first] = it->second;
         vertexMap.clear();
@@ -428,7 +428,7 @@ int GModel::_readMSH2(const std::string &name)
               return 0;
             }
           }
-          MElement *p = NULL;
+          MElement *p = nullptr;
           bool own = false;
 
           // search parent element
@@ -451,7 +451,7 @@ int GModel::_readMSH2(const std::string &name)
           }
 
           // search domains
-          MElement *doms[2] = {NULL, NULL};
+          MElement *doms[2] = {nullptr, nullptr};
           if(dom1) {
             auto ite = elems.find(dom1);
             if(ite == elems.end())
@@ -537,7 +537,7 @@ int GModel::_readMSH2(const std::string &name)
                 return 0;
               }
             }
-            MElement *p = NULL;
+            MElement *p = nullptr;
             bool own = false;
             if(parent) {
               auto ite = elems.find(parent);
diff --git a/Geo/GModelIO_MSH3.cpp b/Geo/GModelIO_MSH3.cpp
index 2c2c2fea41..75ff216140 100644
--- a/Geo/GModelIO_MSH3.cpp
+++ b/Geo/GModelIO_MSH3.cpp
@@ -57,7 +57,7 @@ static void readMSHEntities(FILE *fp, GModel *gm)
     if(fscanf(fp, "%d %d", &tag, &n) != 2) return;
     GEdge *ge = gm->getEdgeByTag(tag);
     if(!ge) {
-      GVertex *v1 = 0, *v2 = 0;
+      GVertex *v1 = nullptr, *v2 = nullptr;
       for(int j = 0; j < n; j++) {
         int tagv;
         if(fscanf(fp, "%d", &tagv) != 1) {
@@ -131,7 +131,7 @@ void readMSHPeriodicNodes(FILE *fp, GModel *gm) // also used in MSH2
 
     if(fscanf(fp, "%d %d %d", &dim, &slave, &master) != 3) continue;
 
-    GEntity *s = 0, *m = 0;
+    GEntity *s = nullptr, *m = nullptr;
     switch(dim) {
     case 0:
       s = gm->getVertexByTag(slave);
@@ -295,7 +295,7 @@ int GModel::_readMSH3(const std::string &name)
       for(int i = 0; i < numVertices; i++) {
         int num, entity, dim;
         double xyz[3];
-        MVertex *vertex = 0;
+        MVertex *vertex = nullptr;
         if(!binary) {
           if(fscanf(fp, "%d %lf %lf %lf %d", &num, &xyz[0], &xyz[1], &xyz[2],
                     &entity) != 5) {
@@ -320,7 +320,7 @@ int GModel::_readMSH3(const std::string &name)
           }
           if(swap) SwapBytes((char *)&entity, sizeof(int), 1);
         }
-        if(!entity) { vertex = new MVertex(xyz[0], xyz[1], xyz[2], 0, num); }
+        if(!entity) { vertex = new MVertex(xyz[0], xyz[1], xyz[2], nullptr, num); }
         else {
           if(!binary) {
             if(fscanf(fp, "%d", &dim) != 1) {
@@ -420,9 +420,9 @@ int GModel::_readMSH3(const std::string &name)
         Msg::Debug("Vertex numbering is dense");
         _vertexVectorCache.resize(_vertexMapCache.size() + 1);
         if(minVertex == 1)
-          _vertexVectorCache[0] = 0;
+          _vertexVectorCache[0] = nullptr;
         else
-          _vertexVectorCache[numVertices] = 0;
+          _vertexVectorCache[numVertices] = nullptr;
         for(std::map<int, MVertex *>::const_iterator it =
               _vertexMapCache.begin();
             it != _vertexMapCache.end(); ++it)
diff --git a/Geo/GModelIO_MSH4.cpp b/Geo/GModelIO_MSH4.cpp
index 0d9eee5808..441c0d123a 100644
--- a/Geo/GModelIO_MSH4.cpp
+++ b/Geo/GModelIO_MSH4.cpp
@@ -386,13 +386,13 @@ static bool readMSH4Entities(GModel *const model, FILE *fp, bool partition,
         GEdge *ge = model->getEdgeByTag(tag);
         if(!ge) {
           if(partition) {
-            ge = new partitionEdge(model, tag, 0, 0, partitions);
+            ge = new partitionEdge(model, tag, nullptr, nullptr, partitions);
             if(parentTag)
               static_cast<partitionEdge *>(ge)->setParentEntity(
                 model->getEntityByTag(parentDim, parentTag));
           }
           else {
-            ge = new discreteEdge(model, tag, 0, 0);
+            ge = new discreteEdge(model, tag, nullptr, nullptr);
           }
           model->add(ge);
         }
@@ -469,7 +469,7 @@ readMSH4Nodes(GModel *const model, FILE *fp, bool binary, bool &dense,
 
   if(binary) {
     std::size_t data[4];
-    if(fread(data, sizeof(std::size_t), 4, fp) != 4) { return 0; }
+    if(fread(data, sizeof(std::size_t), 4, fp) != 4) { return nullptr; }
     if(swap) SwapBytes((char *)data, sizeof(std::size_t), 4);
     numBlock = data[0];
     totalNumNodes = data[1];
@@ -480,11 +480,11 @@ readMSH4Nodes(GModel *const model, FILE *fp, bool binary, bool &dense,
     if(version >= 4.1) {
       if(fscanf(fp, "%lu %lu %lu %lu", &numBlock, &totalNumNodes, &minTag,
                 &maxTag) != 4) {
-        return 0;
+        return nullptr;
       }
     }
     else {
-      if(fscanf(fp, "%lu %lu", &numBlock, &totalNumNodes) != 2) { return 0; }
+      if(fscanf(fp, "%lu %lu", &numBlock, &totalNumNodes) != 2) { return nullptr; }
     }
   }
 
@@ -506,7 +506,7 @@ readMSH4Nodes(GModel *const model, FILE *fp, bool binary, bool &dense,
       int data[3];
       if(fread(data, sizeof(int), 3, fp) != 3) {
         delete[] vertexCache;
-        return 0;
+        return nullptr;
       }
       if(swap) SwapBytes((char *)data, sizeof(int), 3);
       entityDim = data[0];
@@ -515,7 +515,7 @@ readMSH4Nodes(GModel *const model, FILE *fp, bool binary, bool &dense,
 
       if(fread(&numNodes, sizeof(std::size_t), 1, fp) != 1) {
         delete[] vertexCache;
-        return 0;
+        return nullptr;
       }
       if(swap) SwapBytes((char *)&numNodes, sizeof(std::size_t), 1);
     }
@@ -524,14 +524,14 @@ readMSH4Nodes(GModel *const model, FILE *fp, bool binary, bool &dense,
         if(fscanf(fp, "%d %d %d %lu", &entityDim, &entityTag, &parametric,
                   &numNodes) != 4) {
           delete[] vertexCache;
-          return 0;
+          return nullptr;
         }
       }
       else {
         if(fscanf(fp, "%d %d %d %lu", &entityTag, &entityDim, &parametric,
                   &numNodes) != 4) {
           delete[] vertexCache;
-          return 0;
+          return nullptr;
         }
       }
     }
@@ -548,7 +548,7 @@ readMSH4Nodes(GModel *const model, FILE *fp, bool binary, bool &dense,
       }
       case 1: {
         Msg::Info("Creating discrete curve %d", entityTag);
-        GEdge *ge = new discreteEdge(model, entityTag, 0, 0);
+        GEdge *ge = new discreteEdge(model, entityTag, nullptr, nullptr);
         GModel::current()->add(ge);
         entity = ge;
         break;
@@ -570,7 +570,7 @@ readMSH4Nodes(GModel *const model, FILE *fp, bool binary, bool &dense,
       default: {
         Msg::Error("Invalid dimension %d to create discrete entity", entityDim);
         delete[] vertexCache;
-        return 0;
+        return nullptr;
       }
       }
     }
@@ -582,18 +582,18 @@ readMSH4Nodes(GModel *const model, FILE *fp, bool binary, bool &dense,
     if(binary) {
       if(fread(&tags[0], sizeof(std::size_t), numNodes, fp) != numNodes) {
         delete[] vertexCache;
-        return 0;
+        return nullptr;
       }
       if(swap) SwapBytes((char *)&tags[0], sizeof(std::size_t), numNodes);
       std::vector<double> coord(n * numNodes);
       if(fread(&coord[0], sizeof(double), n * numNodes, fp) != n * numNodes) {
         delete[] vertexCache;
-        return 0;
+        return nullptr;
       }
       if(swap) SwapBytes((char *)&coord[0], sizeof(double), n * numNodes);
       std::size_t k = 0;
       for(std::size_t j = 0; j < numNodes; j++) {
-        MVertex *mv = 0;
+        MVertex *mv = nullptr;
         std::size_t tagNode = tags[j];
         if(n == 5) {
           mv = new MFaceVertex(coord[k], coord[k + 1], coord[k + 2], entity,
@@ -623,7 +623,7 @@ readMSH4Nodes(GModel *const model, FILE *fp, bool binary, bool &dense,
         for(std::size_t j = 0; j < numNodes; j++) {
           if(fscanf(fp, "%lu", &tags[j]) != 1) {
             delete[] vertexCache;
-            return 0;
+            return nullptr;
           }
         }
       }
@@ -633,15 +633,15 @@ readMSH4Nodes(GModel *const model, FILE *fp, bool binary, bool &dense,
         else {
           if(fscanf(fp, "%lu", &tagNode) != 1) {
             delete[] vertexCache;
-            return 0;
+            return nullptr;
           }
         }
-        MVertex *mv = 0;
+        MVertex *mv = nullptr;
         if(n == 5) {
           double x, y, z, u, v;
           if(fscanf(fp, "%lf %lf %lf %lf %lf", &x, &y, &z, &u, &v) != 5) {
             delete[] vertexCache;
-            return 0;
+            return nullptr;
           }
           mv = new MFaceVertex(x, y, z, entity, u, v, tagNode);
         }
@@ -649,7 +649,7 @@ readMSH4Nodes(GModel *const model, FILE *fp, bool binary, bool &dense,
           double x, y, z, u;
           if(fscanf(fp, "%lf %lf %lf %lf", &x, &y, &z, &u) != 4) {
             delete[] vertexCache;
-            return 0;
+            return nullptr;
           }
           mv = new MEdgeVertex(x, y, z, entity, u, tagNode);
         }
@@ -657,14 +657,14 @@ readMSH4Nodes(GModel *const model, FILE *fp, bool binary, bool &dense,
           double x, y, z;
           if(fscanf(fp, "%lf %lf %lf", &x, &y, &z) != 3) {
             delete[] vertexCache;
-            return 0;
+            return nullptr;
           }
           // discard extra parametric coordinates, as Gmsh does not use them
           for(std::size_t k = 3; k < n; k++) {
             double dummy;
             if(fscanf(fp, "%lf", &dummy) != 1) {
               delete[] vertexCache;
-              return 0;
+              return nullptr;
             }
           }
           mv = new MVertex(x, y, z, entity, tagNode);
@@ -719,7 +719,7 @@ readMSH4Elements(GModel *const model, FILE *fp, bool binary, bool &dense,
 
   if(binary) {
     std::size_t data[4];
-    if(fread(data, sizeof(std::size_t), 4, fp) != 4) { return 0; }
+    if(fread(data, sizeof(std::size_t), 4, fp) != 4) { return nullptr; }
     if(swap) SwapBytes((char *)data, sizeof(std::size_t), 4);
     numBlock = data[0];
     totalNumElements = data[1];
@@ -730,11 +730,11 @@ readMSH4Elements(GModel *const model, FILE *fp, bool binary, bool &dense,
     if(version >= 4.1) {
       if(fscanf(fp, "%lu %lu %lu %lu", &numBlock, &totalNumElements, &minTag,
                 &maxTag) != 4) {
-        return 0;
+        return nullptr;
       }
     }
     else {
-      if(fscanf(fp, "%lu %lu", &numBlock, &totalNumElements) != 2) { return 0; }
+      if(fscanf(fp, "%lu %lu", &numBlock, &totalNumElements) != 2) { return nullptr; }
     }
   }
 
@@ -754,7 +754,7 @@ readMSH4Elements(GModel *const model, FILE *fp, bool binary, bool &dense,
       int data[3];
       if(fread(data, sizeof(int), 3, fp) != 3) {
         delete[] elementCache;
-        return 0;
+        return nullptr;
       }
       if(swap) SwapBytes((char *)data, sizeof(int), 3);
       entityDim = data[0];
@@ -763,7 +763,7 @@ readMSH4Elements(GModel *const model, FILE *fp, bool binary, bool &dense,
 
       if(fread(&numElements, sizeof(std::size_t), 1, fp) != 1) {
         delete[] elementCache;
-        return 0;
+        return nullptr;
       }
       if(swap) SwapBytes((char *)&numElements, sizeof(std::size_t), 1);
     }
@@ -772,14 +772,14 @@ readMSH4Elements(GModel *const model, FILE *fp, bool binary, bool &dense,
         if(fscanf(fp, "%d %d %d %lu", &entityDim, &entityTag, &elmType,
                   &numElements) != 4) {
           delete[] elementCache;
-          return 0;
+          return nullptr;
         }
       }
       else {
         if(fscanf(fp, "%d %d %d %lu", &entityTag, &entityDim, &elmType,
                   &numElements) != 4) {
           delete[] elementCache;
-          return 0;
+          return nullptr;
         }
       }
     }
@@ -788,7 +788,7 @@ readMSH4Elements(GModel *const model, FILE *fp, bool binary, bool &dense,
     if(!entity) {
       Msg::Error("Unknown entity %d of dimension %d", entityTag, entityDim);
       delete[] elementCache;
-      return 0;
+      return nullptr;
     }
     if(entity->geomType() == GEntity::GhostCurve) {
       static_cast<ghostEdge *>(entity)->haveMesh(true);
@@ -807,12 +807,12 @@ readMSH4Elements(GModel *const model, FILE *fp, bool binary, bool &dense,
       if(fread(&data[0], sizeof(std::size_t), numElements * n, fp) !=
          numElements * n) {
         delete[] elementCache;
-        return 0;
+        return nullptr;
       }
       if(swap)
         SwapBytes((char *)&data[0], sizeof(std::size_t), numElements * n);
 
-      std::vector<MVertex *> vertices(numVertPerElm, (MVertex *)0);
+      std::vector<MVertex *> vertices(numVertPerElm, (MVertex *)nullptr);
       for(std::size_t j = 0; j < numElements * n; j += n) {
         for(int k = 0; k < numVertPerElm; k++) {
           vertices[k] = model->getMeshVertexByTag(data[j + k + 1]);
@@ -820,18 +820,18 @@ readMSH4Elements(GModel *const model, FILE *fp, bool binary, bool &dense,
             Msg::Error("Unknown node %lu in element %lu", data[j + k + 1],
                        data[j]);
             delete[] elementCache;
-            return 0;
+            return nullptr;
           }
         }
 
         MElementFactory elementFactory;
         MElement *element = elementFactory.create(elmType, vertices, data[j], 0,
-                                                  false, 0, 0, 0, 0);
+                                                  false, 0, nullptr, nullptr, nullptr);
         if(!element) {
           Msg::Error("Could not create element %lu of type %d", data[j],
                      elmType);
           delete[] elementCache;
-          return 0;
+          return nullptr;
         }
         if(entity->geomType() != GEntity::GhostCurve &&
            entity->geomType() != GEntity::GhostSurface &&
@@ -855,27 +855,27 @@ readMSH4Elements(GModel *const model, FILE *fp, bool binary, bool &dense,
         std::size_t elmTag = 0;
         if(fscanf(fp, "%lu", &elmTag) != 1) {
           delete[] elementCache;
-          return 0;
+          return nullptr;
         }
         if(!fgets(str, sizeof(str), fp)) {
           delete[] elementCache;
-          return 0;
+          return nullptr;
         }
 
-        std::vector<MVertex *> vertices(numVertPerElm, (MVertex *)0);
+        std::vector<MVertex *> vertices(numVertPerElm, (MVertex *)nullptr);
 
         for(int k = 0; k < numVertPerElm; k++) {
           std::size_t vertexTag = 0;
           if(k != numVertPerElm - 1) {
             if(sscanf(str, "%lu %[0-9- ]", &vertexTag, str) != 2) {
               delete[] elementCache;
-              return 0;
+              return nullptr;
             }
           }
           else {
             if(sscanf(str, "%lu", &vertexTag) != 1) {
               delete[] elementCache;
-              return 0;
+              return nullptr;
             }
           }
 
@@ -883,18 +883,18 @@ readMSH4Elements(GModel *const model, FILE *fp, bool binary, bool &dense,
           if(!vertices[k]) {
             Msg::Error("Unknown node %lu in element %lu", vertexTag, elmTag);
             delete[] elementCache;
-            return 0;
+            return nullptr;
           }
         }
 
         MElementFactory elementFactory;
         MElement *element = elementFactory.create(elmType, vertices, elmTag, 0,
-                                                  false, 0, 0, 0, 0);
+                                                  false, 0, nullptr, nullptr, nullptr);
         if(!element) {
           Msg::Error("Could not create element %lu of type %d", elmTag,
                      elmType);
           delete[] elementCache;
-          return 0;
+          return nullptr;
         }
         if(entity->geomType() != GEntity::GhostCurve &&
            entity->geomType() != GEntity::GhostSurface &&
@@ -964,7 +964,7 @@ static bool readMSH4PeriodicNodes(GModel *const model, FILE *fp, bool binary,
       }
     }
 
-    GEntity *slave = 0, *master = 0;
+    GEntity *slave = nullptr, *master = nullptr;
     switch(slaveDim) {
     case 0:
       slave = model->getVertexByTag(slaveTag);
@@ -1148,7 +1148,7 @@ static bool readMSH4GhostElements(GModel *const model, FILE *fp, bool binary,
     }
   }
 
-  std::vector<GEntity *> ghostEntities(model->getNumPartitions() + 1, 0);
+  std::vector<GEntity *> ghostEntities(model->getNumPartitions() + 1, nullptr);
   std::vector<GEntity *> entities;
   model->getEntities(entities);
   for(std::size_t i = 0; i < entities.size(); i++) {
@@ -1365,7 +1365,7 @@ int GModel::_readMSH4(const std::string &name)
         return false;
       }
       if(dense) {
-        _vertexVectorCache.resize(maxNodeNum + 1, 0);
+        _vertexVectorCache.resize(maxNodeNum + 1, nullptr);
         for(std::size_t i = 0; i < totalNumNodes; i++) {
           if(!_vertexVectorCache[vertexCache[i].first]) {
             _vertexVectorCache[vertexCache[i].first] = vertexCache[i].second;
@@ -1400,7 +1400,7 @@ int GModel::_readMSH4(const std::string &name)
         return 0;
       }
       if(dense) {
-        _elementVectorCache.resize(maxElementNum + 1, (MElement *)0);
+        _elementVectorCache.resize(maxElementNum + 1, (MElement *)nullptr);
         for(std::size_t i = 0; i < totalNumElements; i++) {
           if(!_elementVectorCache[elementCache[i].first]) {
             _elementVectorCache[elementCache[i].first] = elementCache[i].second;
@@ -2745,7 +2745,7 @@ int GModel::_writeMSH4(const std::string &name, double version, bool binary,
                        bool saveAll, bool saveParametric, double scalingFactor,
                        bool append)
 {
-  FILE *fp = 0;
+  FILE *fp = nullptr;
   if(append)
     fp = Fopen(name.c_str(), binary ? "ab" : "a");
   else
@@ -2846,7 +2846,7 @@ static void associateVertices(GModel *model)
     for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) {
       for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices();
           k++) {
-        (*it)->getMeshElement(j)->getVertex(k)->setEntity(0);
+        (*it)->getMeshElement(j)->getVertex(k)->setEntity(nullptr);
       }
     }
     (*it)->mesh_vertices.clear();
@@ -2855,7 +2855,7 @@ static void associateVertices(GModel *model)
     for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) {
       for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices();
           k++) {
-        (*it)->getMeshElement(j)->getVertex(k)->setEntity(0);
+        (*it)->getMeshElement(j)->getVertex(k)->setEntity(nullptr);
       }
     }
     (*it)->mesh_vertices.clear();
@@ -2864,7 +2864,7 @@ static void associateVertices(GModel *model)
     for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) {
       for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices();
           k++) {
-        (*it)->getMeshElement(j)->getVertex(k)->setEntity(0);
+        (*it)->getMeshElement(j)->getVertex(k)->setEntity(nullptr);
       }
     }
     (*it)->mesh_vertices.clear();
@@ -2873,7 +2873,7 @@ static void associateVertices(GModel *model)
     for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) {
       for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices();
           k++) {
-        (*it)->getMeshElement(j)->getVertex(k)->setEntity(0);
+        (*it)->getMeshElement(j)->getVertex(k)->setEntity(nullptr);
       }
     }
     (*it)->mesh_vertices.clear();
@@ -2883,7 +2883,7 @@ static void associateVertices(GModel *model)
     for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) {
       for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices();
           k++) {
-        if((*it)->getMeshElement(j)->getVertex(k)->onWhat() == 0) {
+        if((*it)->getMeshElement(j)->getVertex(k)->onWhat() == nullptr) {
           (*it)->getMeshElement(j)->getVertex(k)->setEntity(*it);
           (*it)->mesh_vertices.push_back(
             (*it)->getMeshElement(j)->getVertex(k));
@@ -2895,7 +2895,7 @@ static void associateVertices(GModel *model)
     for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) {
       for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices();
           k++) {
-        if((*it)->getMeshElement(j)->getVertex(k)->onWhat() == 0) {
+        if((*it)->getMeshElement(j)->getVertex(k)->onWhat() == nullptr) {
           (*it)->getMeshElement(j)->getVertex(k)->setEntity(*it);
           (*it)->mesh_vertices.push_back(
             (*it)->getMeshElement(j)->getVertex(k));
@@ -2907,7 +2907,7 @@ static void associateVertices(GModel *model)
     for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) {
       for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices();
           k++) {
-        if((*it)->getMeshElement(j)->getVertex(k)->onWhat() == 0) {
+        if((*it)->getMeshElement(j)->getVertex(k)->onWhat() == nullptr) {
           (*it)->getMeshElement(j)->getVertex(k)->setEntity(*it);
           (*it)->mesh_vertices.push_back(
             (*it)->getMeshElement(j)->getVertex(k));
@@ -2919,7 +2919,7 @@ static void associateVertices(GModel *model)
     for(std::size_t j = 0; j < (*it)->getNumMeshElements(); j++) {
       for(std::size_t k = 0; k < (*it)->getMeshElement(j)->getNumVertices();
           k++) {
-        if((*it)->getMeshElement(j)->getVertex(k)->onWhat() == 0) {
+        if((*it)->getMeshElement(j)->getVertex(k)->onWhat() == nullptr) {
           (*it)->getMeshElement(j)->getVertex(k)->setEntity(*it);
           (*it)->mesh_vertices.push_back(
             (*it)->getMeshElement(j)->getVertex(k));
@@ -2951,7 +2951,7 @@ int GModel::_writePartitionedMSH4(const std::string &baseName, double version,
 
   for(std::size_t i = 1; i <= getNumPartitions(); i++) {
     std::set<GEntity *> entitiesSet;
-    GEntity *ghostEntity = 0;
+    GEntity *ghostEntity = nullptr;
 
     for(std::size_t j = 0; j < entities.size(); j++) {
       GEntity *ge = entities[j];
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 804fac3611..783e61efbd 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -2776,7 +2776,7 @@ void OCC_Internals::_copyExtrudedAttributes(TopoDS_Edge edge, GEdge *ge)
     while(++recur < CTX::instance()->mesh.maxRetries) {
       if(ge->tag() == p->geo.Source) {
         Msg::Info("Extrusion layer cycle detected for curve %d", ge->tag());
-        e = 0;
+        e = nullptr;
         break;
       }
       GEdge *src = ge->model()->getEdgeByTag(p->geo.Source);
@@ -2810,7 +2810,7 @@ void OCC_Internals::_copyExtrudedAttributes(TopoDS_Face face, GFace *gf)
     while(++recur < CTX::instance()->mesh.maxRetries) {
       if(gf->tag() == p->geo.Source) {
         Msg::Info("Extrusion layer cycle detected for surface %d", gf->tag());
-        e = 0;
+        e = nullptr;
         break;
       }
       GFace *src = gf->model()->getFaceByTag(p->geo.Source);
@@ -2911,7 +2911,7 @@ bool OCC_Internals::_extrudePerDim(
       }
       result = p.Shape();
       const BRepSweep_Prism &prism(p.Prism());
-      _setExtrudedAttributes(c, (BRepSweep_Prism *)&prism, 0, e, 0., 0., 0., dx,
+      _setExtrudedAttributes(c, (BRepSweep_Prism *)&prism, nullptr, e, 0., 0., 0., dx,
                              dy, dz, 0., 0., 0., 0.);
       dim = getReturnedShapes(c, (BRepSweep_Prism *)&prism, top, body, lateral);
     }
@@ -2925,7 +2925,7 @@ bool OCC_Internals::_extrudePerDim(
       }
       result = r.Shape();
       const BRepSweep_Revol &revol(r.Revol());
-      _setExtrudedAttributes(c, 0, (BRepSweep_Revol *)&revol, e, x, y, z, 0.,
+      _setExtrudedAttributes(c, nullptr, (BRepSweep_Revol *)&revol, e, x, y, z, 0.,
                              0., 0., ax, ay, az, angle);
       dim = getReturnedShapes(c, (BRepSweep_Revol *)&revol, top, body, lateral);
     }
@@ -3684,7 +3684,7 @@ bool OCC_Internals::translate(
     gp_Trsf t;
     t.SetTranslation(gp_Pnt(0, 0, 0), gp_Pnt(dx, dy, dz));
     BRepBuilderAPI_Transform tfo(t);
-    return _transform(inDimTags, &tfo, 0);
+    return _transform(inDimTags, &tfo, nullptr);
   } catch(Standard_Failure &err) {
     Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
     return false;
@@ -3700,7 +3700,7 @@ bool OCC_Internals::rotate(const std::vector<std::pair<int, int> > &inDimTags,
     gp_Ax1 axisOfRevolution(gp_Pnt(x, y, z), gp_Dir(ax, ay, az));
     t.SetRotation(axisOfRevolution, angle);
     BRepBuilderAPI_Transform tfo(t);
-    return _transform(inDimTags, &tfo, 0);
+    return _transform(inDimTags, &tfo, nullptr);
   } catch(Standard_Failure &err) {
     Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
     return false;
@@ -3716,7 +3716,7 @@ bool OCC_Internals::dilate(const std::vector<std::pair<int, int> > &inDimTags,
     gt.SetVectorialPart(gp_Mat(a, 0, 0, 0, b, 0, 0, 0, c));
     gt.SetTranslationPart(gp_XYZ(x * (1 - a), y * (1 - b), z * (1 - c)));
     BRepBuilderAPI_GTransform gtfo(gt);
-    return _transform(inDimTags, 0, &gtfo);
+    return _transform(inDimTags, nullptr, &gtfo);
   } catch(Standard_Failure &err) {
     Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
     return false;
@@ -3736,7 +3736,7 @@ bool OCC_Internals::symmetry(const std::vector<std::pair<int, int> > &inDimTags,
                                1. + c * c * f));
     gt.SetTranslationPart(gp_XYZ(a * d * f, b * d * f, c * d * f));
     BRepBuilderAPI_GTransform gtfo(gt);
-    return _transform(inDimTags, 0, &gtfo);
+    return _transform(inDimTags, nullptr, &gtfo);
   } catch(Standard_Failure &err) {
     Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
     return false;
@@ -3757,7 +3757,7 @@ bool OCC_Internals::affine(const std::vector<std::pair<int, int> > &inDimTags,
       gp_Mat(a[0], a[1], a[2], a[4], a[5], a[6], a[8], a[9], a[10]));
     gt.SetTranslationPart(gp_XYZ(a[3], a[7], a[11]));
     BRepBuilderAPI_GTransform gtfo(gt);
-    return _transform(inDimTags, 0, &gtfo);
+    return _transform(inDimTags, nullptr, &gtfo);
   } catch(Standard_Failure &err) {
     Msg::Error("OpenCASCADE exception %s", err.GetMessageString());
     return false;
@@ -4453,7 +4453,7 @@ GVertex *OCC_Internals::getVertexForOCCShape(GModel *model,
 {
   if(_vertexTag.IsBound(toFind))
     return model->getVertexByTag(_vertexTag.Find(toFind));
-  return 0;
+  return nullptr;
 }
 
 GEdge *OCC_Internals::getEdgeForOCCShape(GModel *model,
@@ -4461,7 +4461,7 @@ GEdge *OCC_Internals::getEdgeForOCCShape(GModel *model,
 {
   if(_edgeTag.IsBound(toFind))
     return model->getEdgeByTag(_edgeTag.Find(toFind));
-  return 0;
+  return nullptr;
 }
 
 GFace *OCC_Internals::getFaceForOCCShape(GModel *model,
@@ -4469,7 +4469,7 @@ GFace *OCC_Internals::getFaceForOCCShape(GModel *model,
 {
   if(_faceTag.IsBound(toFind))
     return model->getFaceByTag(_faceTag.Find(toFind));
-  return 0;
+  return nullptr;
 }
 
 GRegion *OCC_Internals::getRegionForOCCShape(GModel *model,
@@ -4477,7 +4477,7 @@ GRegion *OCC_Internals::getRegionForOCCShape(GModel *model,
 {
   if(_solidTag.IsBound(toFind))
     return model->getRegionByTag(_solidTag.Find(toFind));
-  return 0;
+  return nullptr;
 }
 
 void OCC_Internals::_addShapeToMaps(const TopoDS_Shape &shape)
@@ -5137,7 +5137,7 @@ bool OCC_Internals::makeFaceSTL(const TopoDS_Face &s,
                                 std::vector<SPoint2> &vertices_uv,
                                 std::vector<int> &triangles)
 {
-  return makeSTL(s, &vertices_uv, 0, 0, triangles);
+  return makeSTL(s, &vertices_uv, nullptr, nullptr, triangles);
 }
 
 bool OCC_Internals::makeFaceSTL(const TopoDS_Face &s,
@@ -5154,7 +5154,7 @@ bool OCC_Internals::makeFaceSTL(const TopoDS_Face &s,
                                 std::vector<SVector3> &normals,
                                 std::vector<int> &triangles)
 {
-  return makeSTL(s, 0, &vertices, &normals, triangles);
+  return makeSTL(s, nullptr, &vertices, &normals, triangles);
 }
 
 bool OCC_Internals::_makeSTL(const TopoDS_Shape &s,
@@ -5166,7 +5166,7 @@ bool OCC_Internals::_makeSTL(const TopoDS_Shape &s,
   TopExp_Explorer exp0;
   for(exp0.Init(s, TopAbs_FACE); exp0.More(); exp0.Next()) {
     TopoDS_Face face = TopoDS::Face(exp0.Current());
-    bool tmp = makeSTL(TopoDS::Face(face.Oriented(TopAbs_FORWARD)), 0,
+    bool tmp = makeSTL(TopoDS::Face(face.Oriented(TopAbs_FORWARD)), nullptr,
                        &vertices, &normals, triangles);
     if(!tmp) ret = false;
   }
@@ -5303,7 +5303,7 @@ void GModel::createOCCInternals()
 void GModel::deleteOCCInternals()
 {
   if(_occ_internals) delete _occ_internals;
-  _occ_internals = 0;
+  _occ_internals = nullptr;
 }
 
 void GModel::resetOCCInternals()
@@ -5376,7 +5376,7 @@ int GModel::importOCCShape(const void *shape)
 
 GVertex *GModel::getVertexForOCCShape(const void *shape)
 {
-  if(!_occ_internals) return 0;
+  if(!_occ_internals) return nullptr;
 #if defined(HAVE_OCC)
   return _occ_internals->getVertexForOCCShape(this, *(TopoDS_Vertex *)shape);
 #else
@@ -5386,7 +5386,7 @@ GVertex *GModel::getVertexForOCCShape(const void *shape)
 
 GEdge *GModel::getEdgeForOCCShape(const void *shape)
 {
-  if(!_occ_internals) return 0;
+  if(!_occ_internals) return nullptr;
 #if defined(HAVE_OCC)
   return _occ_internals->getEdgeForOCCShape(this, *(TopoDS_Edge *)shape);
 #else
@@ -5396,7 +5396,7 @@ GEdge *GModel::getEdgeForOCCShape(const void *shape)
 
 GFace *GModel::getFaceForOCCShape(const void *shape)
 {
-  if(!_occ_internals) return 0;
+  if(!_occ_internals) return nullptr;
 #if defined(HAVE_OCC)
   return _occ_internals->getFaceForOCCShape(this, *(TopoDS_Face *)shape);
 #else
@@ -5406,7 +5406,7 @@ GFace *GModel::getFaceForOCCShape(const void *shape)
 
 GRegion *GModel::getRegionForOCCShape(const void *shape)
 {
-  if(!_occ_internals) return 0;
+  if(!_occ_internals) return nullptr;
 #if defined(HAVE_OCC)
   return _occ_internals->getRegionForOCCShape(this, *(TopoDS_Solid *)shape);
 #else
diff --git a/Geo/GModelIO_OCC.h b/Geo/GModelIO_OCC.h
index 60540e92ff..33e2bf25ea 100644
--- a/Geo/GModelIO_OCC.h
+++ b/Geo/GModelIO_OCC.h
@@ -135,7 +135,7 @@ private:
                 double x, double y, double z, double dx, double dy, double dz,
                 double ax, double ay, double az, double angle, int wireTag,
                 std::vector<std::pair<int, int> > &outDimTags,
-                ExtrudeParams *e = 0, const std::string &trihedron = "");
+                ExtrudeParams *e = nullptr, const std::string &trihedron = "");
 
   // apply fillet-like operations
   bool _fillet(int mode, const std::vector<int> &volumeTags,
@@ -272,11 +272,11 @@ public:
   bool extrude(const std::vector<std::pair<int, int> > &inDimTags, double dx,
                double dy, double dz,
                std::vector<std::pair<int, int> > &outDimTags,
-               ExtrudeParams *e = 0);
+               ExtrudeParams *e = nullptr);
   bool revolve(const std::vector<std::pair<int, int> > &inDimTags, double x,
                double y, double z, double ax, double ay, double az,
                double angle, std::vector<std::pair<int, int> > &outDimTags,
-               ExtrudeParams *e = 0);
+               ExtrudeParams *e = nullptr);
   bool addPipe(const std::vector<std::pair<int, int> > &inDimTags, int wireTag,
                std::vector<std::pair<int, int> > &outDimTags,
                const std::string &trihedron = "");
diff --git a/Geo/GModelIO_SAMCEF.cpp b/Geo/GModelIO_SAMCEF.cpp
index 2880cb8acf..138d92f5a7 100644
--- a/Geo/GModelIO_SAMCEF.cpp
+++ b/Geo/GModelIO_SAMCEF.cpp
@@ -54,7 +54,7 @@ int GModel::readSAMCEF(const std::string &name)
         if(sscanf(buffer, "%s %d %s %lf %s %lf %s %lf", dummy, &num, dummy, &x,
                   dummy, &y, dummy, &z) != 8)
           return 0;
-        _vertexMapCache[num] = new MVertex(x, y, z, 0, num);
+        _vertexMapCache[num] = new MVertex(x, y, z, nullptr, num);
       }
       Msg::Info("Read %d mesh nodes", (int)_vertexMapCache.size());
     }
diff --git a/Geo/GModelIO_UNV.cpp b/Geo/GModelIO_UNV.cpp
index 23d3c87a2e..2cb6c9ccda 100644
--- a/Geo/GModelIO_UNV.cpp
+++ b/Geo/GModelIO_UNV.cpp
@@ -66,7 +66,7 @@ int GModel::readUNV(const std::string &name)
           for(std::size_t i = 0; i < strlen(buffer); i++)
             if(buffer[i] == 'D') buffer[i] = 'E';
           if(sscanf(buffer, "%lf %lf %lf", &x, &y, &z) != 3) break;
-          _vertexMapCache[num] = new MVertex(x, y, z, 0, num);
+          _vertexMapCache[num] = new MVertex(x, y, z, nullptr, num);
         }
       }
       else if(record == 2412) { // elements
diff --git a/Geo/GModelParametrize.cpp b/Geo/GModelParametrize.cpp
index 50f2258d9b..d52e28b387 100644
--- a/Geo/GModelParametrize.cpp
+++ b/Geo/GModelParametrize.cpp
@@ -45,7 +45,7 @@ getModelEdge(GModel *gm, std::vector<GFace *> &gfs,
              std::vector<std::pair<GEdge *, std::vector<GFace *> > > &newEdges,
              size_t &MAX1)
 {
-  if(gfs.size() == 2 && gfs[0] == gfs[1]) return NULL;
+  if(gfs.size() == 2 && gfs[0] == gfs[1]) return nullptr;
   for(size_t i = 0; i < newEdges.size(); i++) {
     if(gfs.size() == newEdges[i].second.size()) {
       bool found = true;
@@ -59,7 +59,7 @@ getModelEdge(GModel *gm, std::vector<GFace *> &gfs,
     }
   }
 
-  discreteEdge *ge = new discreteEdge(gm, (MAX1++) + 1, 0, 0);
+  discreteEdge *ge = new discreteEdge(gm, (MAX1++) + 1, nullptr, nullptr);
   newEdges.push_back(std::make_pair(ge, gfs));
   return ge;
 }
@@ -171,7 +171,7 @@ void classifyFaces(GModel *gm, double curveAngleThreshold)
   // reset classification of all mesh nodes
   for(auto it = touched.begin(); it != touched.end(); it++) {
     for(std::size_t j = 0; j < (*it)->getNumVertices(); j++)
-      (*it)->getVertex(j)->setEntity(0);
+      (*it)->getVertex(j)->setEntity(nullptr);
   }
 
   std::map<MTriangle *, GFace *> reverse;
@@ -414,7 +414,7 @@ void classifyFaces(GModel *gm, double angleThreshold, bool includeBoundary,
                     (*it)->quadrangles.end());
   }
 
-  discreteEdge *edge = new discreteEdge(gm, (MAX1++) + 1, 0, 0);
+  discreteEdge *edge = new discreteEdge(gm, (MAX1++) + 1, nullptr, nullptr);
   gm->add(edge);
 
   e2t_cont adj;
@@ -789,12 +789,12 @@ void makeMLinesUnique(std::vector<MLine *> &v)
 class twoT {
 public:
   MTriangle *t1, *t2;
-  twoT(MTriangle *t) : t1(t), t2(NULL) {}
+  twoT(MTriangle *t) : t1(t), t2(nullptr) {}
   MTriangle *other(MTriangle *t) const
   {
     if(t == t1) return t2;
     if(t == t2) return t1;
-    return NULL;
+    return nullptr;
   }
 };
 
diff --git a/Geo/GPoint.h b/Geo/GPoint.h
index 06f1845e52..54c48dc9ef 100644
--- a/Geo/GPoint.h
+++ b/Geo/GPoint.h
@@ -27,7 +27,8 @@ public:
   inline double u() const { return par[0]; }
   inline double v() const { return par[1]; }
   inline const GEntity *g() const { return e; }
-  GPoint(double _x = 0, double _y = 0, double _z = 0, const GEntity *onwhat = 0)
+  GPoint(double _x = 0, double _y = 0, double _z = 0,
+         const GEntity *onwhat = nullptr)
     : X(_x), Y(_y), Z(_z), e(onwhat), success(true)
   {
     par[0] = -1.;
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index b17aa9829a..1be089c78b 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -100,25 +100,25 @@ MElement *const *GRegion::getStartElementType(int type) const
 {
   switch(type) {
   case 0:
-    if(tetrahedra.empty()) return 0; // msvc would throw an exception
+    if(tetrahedra.empty()) return nullptr; // msvc would throw an exception
     return reinterpret_cast<MElement *const *>(&tetrahedra[0]);
   case 1:
-    if(hexahedra.empty()) return 0; // msvc would throw an exception
+    if(hexahedra.empty()) return nullptr; // msvc would throw an exception
     return reinterpret_cast<MElement *const *>(&hexahedra[0]);
   case 2:
-    if(prisms.empty()) return 0; // msvc would throw an exception
+    if(prisms.empty()) return nullptr; // msvc would throw an exception
     return reinterpret_cast<MElement *const *>(&prisms[0]);
   case 3:
-    if(pyramids.empty()) return 0; // msvc would throw an exception
+    if(pyramids.empty()) return nullptr; // msvc would throw an exception
     return reinterpret_cast<MElement *const *>(&pyramids[0]);
   case 4:
-    if(trihedra.empty()) return 0;
+    if(trihedra.empty()) return nullptr;
     return reinterpret_cast<MElement *const *>(&trihedra[0]);
   case 5:
-    if(polyhedra.empty()) return 0;
+    if(polyhedra.empty()) return nullptr;
     return reinterpret_cast<MElement *const *>(&polyhedra[0]);
   }
-  return 0;
+  return nullptr;
 }
 
 MElement *GRegion::getMeshElement(std::size_t index) const
@@ -142,7 +142,7 @@ MElement *GRegion::getMeshElement(std::size_t index) const
     return polyhedra[index - tetrahedra.size() - hexahedra.size() -
                      prisms.size() - pyramids.size() - trihedra.size()];
 
-  return 0;
+  return nullptr;
 }
 
 MElement *GRegion::getMeshElementByType(const int familyType,
@@ -161,14 +161,14 @@ MElement *GRegion::getMeshElementByType(const int familyType,
   else if(familyType == TYPE_POLYH)
     return polyhedra[index];
 
-  return 0;
+  return nullptr;
 }
 
 void GRegion::resetMeshAttributes()
 {
   meshAttributes.recombine3D = 0;
   meshAttributes.method = MESH_UNSTRUCTURED;
-  meshAttributes.extrude = 0;
+  meshAttributes.extrude = nullptr;
   meshAttributes.QuadTri = NO_QUADTRI;
   meshAttributes.meshSize = MAX_LC;
 }
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index 7dbc85662a..027857db11 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -45,7 +45,7 @@ public:
   virtual int dim() const { return 3; }
 
   // returns the parent entity for partitioned entities
-  virtual GEntity *getParentEntity() { return 0; }
+  virtual GEntity *getParentEntity() { return nullptr; }
 
   // set the visibility flag
   virtual void setVisibility(char val, bool recursive = false);
diff --git a/Geo/GVertex.cpp b/Geo/GVertex.cpp
index dc9f801bd0..c9e3a7e27f 100644
--- a/Geo/GVertex.cpp
+++ b/Geo/GVertex.cpp
@@ -110,7 +110,7 @@ void GVertex::getNumMeshElements(unsigned *const c) const
 MElement *GVertex::getMeshElement(std::size_t index) const
 {
   if(index < points.size()) return points[index];
-  return 0;
+  return nullptr;
 }
 
 MElement *GVertex::getMeshElementByType(const int familyType,
@@ -118,7 +118,7 @@ MElement *GVertex::getMeshElementByType(const int familyType,
 {
   if(familyType == TYPE_PNT) return points[index];
 
-  return 0;
+  return nullptr;
 }
 
 bool GVertex::isOnSeam(const GFace *gf) const
diff --git a/Geo/GVertex.h b/Geo/GVertex.h
index b2f950483c..9c875ce82c 100644
--- a/Geo/GVertex.h
+++ b/Geo/GVertex.h
@@ -63,7 +63,7 @@ public:
   virtual int dim() const { return 0; }
 
   // returns the parent entity for partitioned entities
-  virtual GEntity *getParentEntity() { return 0; }
+  virtual GEntity *getParentEntity() { return nullptr; }
 
   // get the geometric type of the vertex
   virtual GeomType geomType() const { return Point; }
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index b3bbb33f24..45ab6eb000 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -25,7 +25,7 @@
 #include "Parser.h"
 #endif
 
-static List_T *ListOfTransformedPoints = NULL;
+static List_T *ListOfTransformedPoints = nullptr;
 
 int CompareVertex(const void *a, const void *b)
 {
@@ -111,7 +111,7 @@ Vertex *CreateVertex(int Num, double X, double Y, double Z, double lc, double u)
   GModel::current()->getGEOInternals()->setMaxTag(
     0, std::max(GModel::current()->getGEOInternals()->getMaxTag(0), Num));
   pV->u = u;
-  pV->geometry = 0;
+  pV->geometry = nullptr;
   return pV;
 }
 
@@ -135,7 +135,7 @@ void FreeVertex(void *a, void *b)
   Vertex *v = *(Vertex **)a;
   if(v) {
     delete v;
-    v = NULL;
+    v = nullptr;
   }
 }
 
@@ -161,7 +161,7 @@ void FreePhysicalGroup(void *a, void *b)
   if(p) {
     List_Delete(p->Entities);
     delete p;
-    p = NULL;
+    p = nullptr;
   }
 }
 
@@ -186,7 +186,7 @@ void FreeEdgeLoop(void *a, void *b)
   if(l) {
     List_Delete(l->Curves);
     delete l;
-    l = NULL;
+    l = nullptr;
   }
 }
 
@@ -211,7 +211,7 @@ void FreeSurfaceLoop(void *a, void *b)
   if(l) {
     List_Delete(l->Surfaces);
     delete l;
-    l = NULL;
+    l = nullptr;
   }
 }
 
@@ -236,7 +236,7 @@ bool EndCurve(Curve *c)
     for(int i = 1; i < NN; i++) {
       List_Read(c->Control_Points, i, &pV);
       if(c->geometry != pV->geometry) {
-        c->geometry = 0;
+        c->geometry = nullptr;
         break;
       }
     }
@@ -254,7 +254,7 @@ bool EndCurve(Curve *c)
     if(List_Nbr(c->Control_Points) == 4)
       List_Read(c->Control_Points, 2, &v[3]);
     else
-      v[3] = NULL;
+      v[3] = nullptr;
     if(c->Typ == MSH_SEGM_CIRC_INV || c->Typ == MSH_SEGM_ELLI_INV) {
       List_Read(c->Control_Points, 0, &v[2]);
       List_Read(c->Control_Points, 1, &v[1]);
@@ -436,7 +436,7 @@ void EndSurface(Surface *s)
     for(int i = 1; i < NN; i++) {
       List_Read(s->Generatrices, i, &c);
       if(c->geometry != s->geometry) {
-        s->geometry = 0;
+        s->geometry = nullptr;
         break;
       }
     }
@@ -457,7 +457,7 @@ Curve *CreateCurve(int Num, int Typ, int Order, List_T *Liste, List_T *Knots,
     {-1, 3, -3, 1}, {3, -6, 3, 0}, {-3, 3, 0, 0}, {1, 0, 0, 0}};
 
   Curve *pC = new Curve;
-  pC->Extrude = NULL;
+  pC->Extrude = nullptr;
   pC->Typ = Typ;
   pC->Num = Num;
   GModel::current()->getGEOInternals()->setMaxTag(
@@ -467,14 +467,14 @@ Curve *CreateCurve(int Num, int Typ, int Order, List_T *Liste, List_T *Knots,
   pC->Circle.n[0] = 0.0;
   pC->Circle.n[1] = 0.0;
   pC->Circle.n[2] = 1.0;
-  pC->geometry = 0;
+  pC->geometry = nullptr;
   pC->nbPointsTransfinite = 0;
   pC->typeTransfinite = 0;
   pC->coeffTransfinite = 0.;
   pC->ReverseMesh = 0;
-  pC->beg = NULL;
-  pC->end = NULL;
-  pC->Control_Points = NULL;
+  pC->beg = nullptr;
+  pC->end = nullptr;
+  pC->Control_Points = nullptr;
   pC->degenerated = false;
 
   if(Typ == MSH_SEGM_SPLN) {
@@ -507,7 +507,7 @@ Curve *CreateCurve(int Num, int Typ, int Order, List_T *Liste, List_T *Knots,
     }
   }
   else
-    pC->k = NULL;
+    pC->k = nullptr;
 
   if(List_Nbr(Liste)) {
     pC->Control_Points = List_Create(List_Nbr(Liste), 1, sizeof(Vertex *));
@@ -562,7 +562,7 @@ void FreeCurve(void *a, void *b)
     List_Delete(pC->Control_Points);
     if(pC->Extrude) delete pC->Extrude;
     delete pC;
-    pC = NULL;
+    pC = nullptr;
   }
 }
 
@@ -572,8 +572,8 @@ Surface *CreateSurface(int Num, int Typ)
   pS->Num = Num;
   GModel::current()->getGEOInternals()->setMaxTag(
     2, std::max(GModel::current()->getGEOInternals()->getMaxTag(2), Num));
-  pS->geometry = 0;
-  pS->InSphereCenter = 0;
+  pS->geometry = nullptr;
+  pS->InSphereCenter = nullptr;
   pS->Typ = Typ;
   pS->Method = MESH_UNSTRUCTURED;
   pS->Recombine = 0;
@@ -581,10 +581,10 @@ Surface *CreateSurface(int Num, int Typ)
   pS->Recombine_Dir = -1;
   pS->TransfiniteSmoothing = -1;
   pS->TrsfPoints = List_Create(4, 4, sizeof(Vertex *));
-  pS->Generatrices = NULL;
-  pS->GeneratricesByTag = NULL;
-  pS->Extrude = NULL;
-  pS->geometry = NULL;
+  pS->Generatrices = nullptr;
+  pS->GeneratricesByTag = nullptr;
+  pS->Extrude = nullptr;
+  pS->geometry = nullptr;
   pS->ReverseMesh = 0;
   pS->MeshAlgorithm = 0;
   pS->MeshSizeFromBoundary = -1;
@@ -600,7 +600,7 @@ void FreeSurface(void *a, void *b)
     List_Delete(pS->GeneratricesByTag);
     if(pS->Extrude) delete pS->Extrude;
     delete pS;
-    pS = NULL;
+    pS = nullptr;
   }
 }
 
@@ -618,7 +618,7 @@ Volume *CreateVolume(int Num, int Typ)
   pV->Surfaces = List_Create(1, 2, sizeof(Surface *));
   pV->SurfacesOrientations = List_Create(1, 2, sizeof(int));
   pV->SurfacesByTag = List_Create(1, 2, sizeof(int));
-  pV->Extrude = NULL;
+  pV->Extrude = nullptr;
   return pV;
 }
 
@@ -632,7 +632,7 @@ void FreeVolume(void *a, void *b)
     List_Delete(pV->SurfacesByTag);
     if(pV->Extrude) delete pV->Extrude;
     delete pV;
-    pV = NULL;
+    pV = nullptr;
   }
 }
 
@@ -673,7 +673,7 @@ static Vertex *FindPoint(int inum, Tree_T *t)
   pc = &C;
   pc->Num = inum;
   if(Tree_Query(t, &pc)) { return pc; }
-  return NULL;
+  return nullptr;
 }
 
 Vertex *FindPoint(int inum)
@@ -687,7 +687,7 @@ static Curve *FindCurve(int inum, Tree_T *t)
   pc = &C;
   pc->Num = inum;
   if(Tree_Query(t, &pc)) { return pc; }
-  return NULL;
+  return nullptr;
 }
 
 Curve *FindCurve(int inum)
@@ -701,7 +701,7 @@ static Surface *FindSurface(int inum, Tree_T *t)
   ps = &S;
   ps->Num = inum;
   if(Tree_Query(t, &ps)) { return ps; }
-  return NULL;
+  return nullptr;
 }
 
 Surface *FindSurface(int inum)
@@ -717,7 +717,7 @@ Volume *FindVolume(int inum)
   if(Tree_Query(GModel::current()->getGEOInternals()->Volumes, &pv)) {
     return pv;
   }
-  return NULL;
+  return nullptr;
 }
 
 EdgeLoop *FindEdgeLoop(int inum)
@@ -728,7 +728,7 @@ EdgeLoop *FindEdgeLoop(int inum)
   if(Tree_Query(GModel::current()->getGEOInternals()->EdgeLoops, &ps)) {
     return ps;
   }
-  return NULL;
+  return nullptr;
 }
 
 SurfaceLoop *FindSurfaceLoop(int inum)
@@ -739,7 +739,7 @@ SurfaceLoop *FindSurfaceLoop(int inum)
   if(Tree_Query(GModel::current()->getGEOInternals()->SurfaceLoops, &ps)) {
     return ps;
   }
-  return NULL;
+  return nullptr;
 }
 
 PhysicalGroup *FindPhysicalGroup(int num, int type)
@@ -753,7 +753,7 @@ PhysicalGroup *FindPhysicalGroup(int num, int type)
         ComparePhysicalGroup))) {
     return *ppp;
   }
-  return NULL;
+  return nullptr;
 }
 
 static void CopyVertex(Vertex *v, Vertex *vv)
@@ -767,7 +767,7 @@ static void CopyVertex(Vertex *v, Vertex *vv)
 
 Vertex *DuplicateVertex(Vertex *v)
 {
-  if(!v) return NULL;
+  if(!v) return nullptr;
   Vertex *pv = CreateVertex(NEWPOINT(), 0, 0, 0, 0, 0);
   CopyVertex(v, pv);
   Tree_Insert(GModel::current()->getGEOInternals()->Points, &pv);
@@ -807,7 +807,7 @@ static void CopyCurve(Curve *c, Curve *cc)
 Curve *DuplicateCurve(Curve *c)
 {
   bool ok = true;
-  Curve *pc = CreateCurve(NEWLINE(), 0, 1, NULL, NULL, -1, -1, 0., 1., ok);
+  Curve *pc = CreateCurve(NEWLINE(), 0, 1, nullptr, nullptr, -1, -1, 0., 1., ok);
   CopyCurve(c, pc);
   Tree_Insert(GModel::current()->getGEOInternals()->Curves, &pc);
   for(int i = 0; i < List_Nbr(c->Control_Points); i++) {
@@ -842,7 +842,7 @@ static void CopySurface(Surface *s, Surface *ss)
     List_Create(List_Nbr(s->GeneratricesByTag) + 1, 1, sizeof(int));
   // after copy (and subsequent transformation), the sphere center does not make
   // sense anymore
-  ss->InSphereCenter = 0;
+  ss->InSphereCenter = nullptr;
   List_Copy(s->Generatrices, ss->Generatrices);
   List_Copy(s->GeneratricesByTag, ss->GeneratricesByTag);
   EndSurface(ss);
@@ -1086,7 +1086,7 @@ void DeletePhysicalVolume(int num)
 Curve *CreateReversedCurve(Curve *c)
 {
   bool ok = true;
-  Curve *newc = CreateCurve(-c->Num, c->Typ, 1, NULL, NULL, -1, -1, 0., 1., ok);
+  Curve *newc = CreateCurve(-c->Num, c->Typ, 1, nullptr, nullptr, -1, -1, 0., 1., ok);
 
   if(List_Nbr(c->Control_Points)) {
     newc->Control_Points =
@@ -1137,7 +1137,7 @@ Curve *CreateReversedCurve(Curve *c)
   Curve **pc;
   if((pc = (Curve **)Tree_PQuery(GModel::current()->getGEOInternals()->Curves,
                                  &newc))) {
-    FreeCurve(&newc, NULL);
+    FreeCurve(&newc, nullptr);
     return *pc;
   }
   else {
@@ -1784,7 +1784,7 @@ static void ReplaceDuplicatePointsNew(double tol = -1.)
   Msg::Info("Done new Coherence (removed %d additional points)", start - end);
 }
 
-static void ReplaceDuplicatePoints(std::map<int, int> *v_report = 0)
+static void ReplaceDuplicatePoints(std::map<int, int> *v_report = nullptr)
 {
   // This routine is logically wrong: the CompareTwoPoints() function used in
   // the avl tree is not an appropriate comparison function. Fixing the routine
@@ -1961,7 +1961,7 @@ static void ReplaceDuplicatePoints(std::map<int, int> *v_report = 0)
   if(!success) ReplaceDuplicatePointsNew();
 }
 
-static void ReplaceDuplicateCurves(std::map<int, int> *c_report = 0)
+static void ReplaceDuplicateCurves(std::map<int, int> *c_report = nullptr)
 {
   Curve *c, *c2, **pc, **pc2;
   Surface *s;
@@ -2283,7 +2283,7 @@ bool Surface::degenerate() const
   return Nd == 0;
 }
 
-static void ReplaceDuplicateSurfaces(std::map<int, int> *s_report = 0)
+static void ReplaceDuplicateSurfaces(std::map<int, int> *s_report = nullptr)
 {
   Surface *s, *s2, **ps, **ps2;
   Volume *vol;
@@ -2411,9 +2411,9 @@ static void ReplaceDuplicateSurfaces(std::map<int, int> *s_report = 0)
 
 static void ReplaceAllDuplicates(std::vector<std::map<int, int> > &report)
 {
-  std::map<int, int> *vertex_report = 0;
-  std::map<int, int> *curve_report = 0;
-  std::map<int, int> *surface_report = 0;
+  std::map<int, int> *vertex_report = nullptr;
+  std::map<int, int> *curve_report = nullptr;
+  std::map<int, int> *surface_report = nullptr;
   if(report.size() >= 1 && report[0].size()) vertex_report = &(report[0]);
   if(report.size() >= 2 && report[1].size()) curve_report = &(report[1]);
   if(report.size() >= 3 && report[2].size()) surface_report = &(report[2]);
@@ -2440,8 +2440,8 @@ void ReplaceAllDuplicatesNew(double tol)
 {
   if(tol < 0) tol = CTX::instance()->geom.tolerance * CTX::instance()->lc;
   ReplaceDuplicatePointsNew(tol);
-  ReplaceDuplicateCurves(NULL);
-  ReplaceDuplicateSurfaces(NULL);
+  ReplaceDuplicateCurves(nullptr);
+  ReplaceDuplicateSurfaces(nullptr);
 }
 
 void ProtudeXYZ(double &x, double &y, double &z, ExtrudeParams *e)
@@ -2478,14 +2478,14 @@ int ExtrudePoint(int type, int ip, double T0, double T1, double T2, double A0,
                  ExtrudeParams *e)
 {
   double matrix[4][4], T[3], Ax[3], d;
-  Vertex V, *newp = NULL;
-  Curve *c = NULL;
+  Vertex V, *newp = nullptr;
+  Curve *c = nullptr;
   int i;
   bool ok = true;
 
   Vertex *pv = &V;
   pv->Num = ip;
-  *pc = *prc = NULL;
+  *pc = *prc = nullptr;
   if(!Tree_Query(GModel::current()->getGEOInternals()->Points, &pv)) return 0;
 
   Msg::Debug("Extrude Point %d", ip);
@@ -2502,7 +2502,7 @@ int ExtrudePoint(int type, int ip, double T0, double T1, double T2, double A0,
     ApplyTransformationToPoint(matrix, chapeau);
     if(!ComparePosition(&pv, &chapeau)) return pv->Num;
     c =
-      CreateCurve(NEWLINE(), MSH_SEGM_LINE, 1, NULL, NULL, -1, -1, 0., 1., ok);
+      CreateCurve(NEWLINE(), MSH_SEGM_LINE, 1, nullptr, nullptr, -1, -1, 0., 1., ok);
     c->Control_Points = List_Create(2, 1, sizeof(Vertex *));
     c->Extrude = new ExtrudeParams;
     c->Extrude->fill(type, T0, T1, T2, A0, A1, A2, X0, X1, X2, alpha);
@@ -2515,7 +2515,7 @@ int ExtrudePoint(int type, int ip, double T0, double T1, double T2, double A0,
   case BOUNDARY_LAYER:
     chapeau->Typ = MSH_POINT_BND_LAYER;
     if(e) chapeau->boundaryLayerIndex = e->mesh.BoundaryLayerIndex;
-    c = CreateCurve(NEWLINE(), MSH_SEGM_BND_LAYER, 1, NULL, NULL, -1, -1, 0.,
+    c = CreateCurve(NEWLINE(), MSH_SEGM_BND_LAYER, 1, nullptr, nullptr, -1, -1, 0.,
                     1., ok);
     c->Control_Points = List_Create(2, 1, sizeof(Vertex *));
     c->Extrude = new ExtrudeParams;
@@ -2547,7 +2547,7 @@ int ExtrudePoint(int type, int ip, double T0, double T1, double T2, double A0,
     ApplyTransformationToPoint(matrix, chapeau);
     if(!ComparePosition(&pv, &chapeau)) return pv->Num;
     c =
-      CreateCurve(NEWLINE(), MSH_SEGM_CIRC, 1, NULL, NULL, -1, -1, 0., 1., ok);
+      CreateCurve(NEWLINE(), MSH_SEGM_CIRC, 1, nullptr, nullptr, -1, -1, 0., 1., ok);
     c->Control_Points = List_Create(3, 1, sizeof(Vertex *));
     c->Extrude = new ExtrudeParams;
     c->Extrude->fill(type, T0, T1, T2, A0, A1, A2, X0, X1, X2, alpha);
@@ -2575,7 +2575,7 @@ int ExtrudePoint(int type, int ip, double T0, double T1, double T2, double A0,
     d = CTX::instance()->geom.extrudeSplinePoints;
     d = d ? d : 1;
     c =
-      CreateCurve(NEWLINE(), MSH_SEGM_SPLN, 1, NULL, NULL, -1, -1, 0., 1., ok);
+      CreateCurve(NEWLINE(), MSH_SEGM_SPLN, 1, nullptr, nullptr, -1, -1, 0., 1., ok);
     c->Control_Points = List_Create(
       CTX::instance()->geom.extrudeSplinePoints + 1, 1, sizeof(Vertex *));
     c->Extrude = new ExtrudeParams;
@@ -2637,7 +2637,7 @@ int ExtrudePoint(int type, int ip, double T0, double T1, double T2, double A0,
       chap_num = report[0][chap_num];
     else
       chap_num = 0;
-    if(report[1][body_num] != body_num) *pc = *prc = NULL;
+    if(report[1][body_num] != body_num) *pc = *prc = nullptr;
   }
   return chap_num;
 }
@@ -2654,7 +2654,7 @@ int ExtrudeCurve(int type, int ic, double T0, double T1, double T2, double A0,
 
   pc = FindCurve(ic);
   revpc = FindCurve(-ic);
-  *ps = NULL;
+  *ps = nullptr;
 
   if(!pc || !revpc) { return 0; }
 
@@ -2823,7 +2823,7 @@ int ExtrudeCurve(int type, int ic, double T0, double T1, double T2, double A0,
       chap_num = report[1][chap_num];
     else
       chap_num = 0;
-    if(report[2][body_num] != body_num) *ps = NULL;
+    if(report[2][body_num] != body_num) *ps = nullptr;
   }
 
   return chap_num;
@@ -2838,7 +2838,7 @@ int ExtrudeSurface(int type, int is, double T0, double T1, double T2, double A0,
   int i;
   Surface *s, *ps, *chapeau;
 
-  *pv = NULL;
+  *pv = nullptr;
 
   // 'is' can be negative, to signify that the surface orientation
   // should be reversed. This orientation information is only used at
@@ -3032,7 +3032,7 @@ void ExtrudeShapes(int type, List_T *list_in, double T0, double T1, double T2,
     List_Read(list_in, i, &shape);
     switch(shape.Type) {
     case MSH_POINT: {
-      Curve *pc = 0, *prc = 0;
+      Curve *pc = nullptr, *prc = nullptr;
       Shape top;
       top.Num = ExtrudePoint(type, shape.Num, T0, T1, T2, A0, A1, A2, X0, X1,
                              X2, alpha, &pc, &prc, 1, e);
@@ -3054,7 +3054,7 @@ void ExtrudeShapes(int type, List_T *list_in, double T0, double T1, double T2,
     case MSH_SEGM_ELLI:
     case MSH_SEGM_ELLI_INV:
     case MSH_SEGM_NURBS: {
-      Surface *ps = 0;
+      Surface *ps = nullptr;
       Shape top;
       top.Num = ExtrudeCurve(type, shape.Num, T0, T1, T2, A0, A1, A2, X0, X1,
                              X2, alpha, &ps, 1, e);
@@ -3084,7 +3084,7 @@ void ExtrudeShapes(int type, List_T *list_in, double T0, double T1, double T2,
     case MSH_SURF_TRIC:
     case MSH_SURF_PLAN:
     case MSH_SURF_DISCRETE: {
-      Volume *pv = 0;
+      Volume *pv = nullptr;
       Shape top;
       top.Num = ExtrudeSurface(type, shape.Num, T0, T1, T2, A0, A1, A2, X0, X1,
                                X2, alpha, &pv, e);
@@ -3214,21 +3214,21 @@ static Curve *_create_splitted_curve(Curve *c, List_T *nodes)
   List_Read(nodes, 0, &beg);
   List_Read(nodes, List_Nbr(nodes) - 1, &end);
   int id = NEWLINE();
-  Curve *cnew = NULL;
+  Curve *cnew = nullptr;
   bool ok = true;
   switch(c->Typ) {
   case MSH_SEGM_LINE:
-    cnew = CreateCurve(id, c->Typ, 1, nodes, NULL, -1, -1, 0., 1., ok);
+    cnew = CreateCurve(id, c->Typ, 1, nodes, nullptr, -1, -1, 0., 1., ok);
     break;
   case MSH_SEGM_SPLN:
-    cnew = CreateCurve(id, c->Typ, 3, nodes, NULL, -1, -1, 0., 1., ok);
+    cnew = CreateCurve(id, c->Typ, 3, nodes, nullptr, -1, -1, 0., 1., ok);
     break;
   case MSH_SEGM_BSPLN:
-    cnew = CreateCurve(id, c->Typ, 2, nodes, NULL, -1, -1, 0., 1., ok);
+    cnew = CreateCurve(id, c->Typ, 2, nodes, nullptr, -1, -1, 0., 1., ok);
     break;
   default: // should never reach this point...
     Msg::Error("Cannot split a curve with type %i", c->Typ);
-    return NULL;
+    return nullptr;
   }
   Tree_Add(GModel::current()->getGEOInternals()->Curves, &cnew);
   CreateReversedCurve(cnew);
@@ -3492,7 +3492,7 @@ bool SetSurfaceGeneratrices(Surface *s, List_T *loops)
     if(!(el = FindEdgeLoop(abs(iLoop)))) {
       Msg::Error("Unknown curve loop %d in GEO surface %d", iLoop, s->Num);
       List_Delete(s->Generatrices);
-      s->Generatrices = NULL;
+      s->Generatrices = nullptr;
       return false;
     }
     else {
diff --git a/Geo/Geo.h b/Geo/Geo.h
index 7eed104d81..c08ed704fb 100644
--- a/Geo/Geo.h
+++ b/Geo/Geo.h
@@ -41,7 +41,7 @@ public:
   int boundaryLayerIndex;
   Vertex(double X = 0., double Y = 0., double Z = 0., double l = 1.,
          double W = 1.)
-    : Num(0), lc(l), u(0.), w(W), geometry(0), boundaryLayerIndex(0)
+    : Num(0), lc(l), u(0.), w(W), geometry(nullptr), boundaryLayerIndex(0)
   {
     Typ = MSH_POINT;
     Pos.X = X;
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index e8ae4c84e9..9acd51fa83 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -805,7 +805,7 @@ static void TransfiniteSph(Vertex S, Vertex center, Vertex *T)
 
 static Vertex InterpolateRuledSurface(Surface *s, double u, double v)
 {
-  Curve *C[4] = {0, 0, 0, 0};
+  Curve *C[4] = {nullptr, nullptr, nullptr, nullptr};
 
   if(!List_Nbr(s->Generatrices)) {
     Msg::Error("No curves on boundary of ruled surface");
@@ -826,7 +826,7 @@ static Vertex InterpolateRuledSurface(Surface *s, double u, double v)
     }
   }
 
-  Vertex *O = 0;
+  Vertex *O = nullptr;
   bool isSphere = true;
 
   // Ugly hack: "fix" transfinite interpolation if we have a sphere
diff --git a/Geo/GeomMeshMatcher.cpp b/Geo/GeomMeshMatcher.cpp
index 9edf86d0a4..0060c688b9 100644
--- a/Geo/GeomMeshMatcher.cpp
+++ b/Geo/GeomMeshMatcher.cpp
@@ -32,7 +32,7 @@
 
 #include "closestVertex.h"
 
-GeomMeshMatcher *GeomMeshMatcher::_gmm_instance = 0;
+GeomMeshMatcher *GeomMeshMatcher::_gmm_instance = nullptr;
 
 template <class T, class container>
 void getIntersection(std::vector<T> &res, std::vector<container> &lists)
@@ -88,7 +88,7 @@ GeomMeshMatcher::matchVertices(GModel *m1, GModel *m2, bool &ok)
     // FIXME: need a *much* better way to fix the tolerance...
     double tol = CTX::instance()->geom.matchMeshTolerance;
 
-    discreteVertex *choice = 0;
+    discreteVertex *choice = nullptr;
     double best_score = DBL_MAX;
 
     for(auto vit2 = m2->firstVertex(); vit2 != m2->lastVertex(); vit2++) {
@@ -172,7 +172,7 @@ GeomMeshMatcher::matchEdges(GModel *m1, GModel *m2,
       if(ok1 && ok2) getIntersection<GEdge *>(common_edges, lists);
     }
 
-    GEdge *choice = 0;
+    GEdge *choice = nullptr;
     if(common_edges.size() == 0) continue;
     if(common_edges.size() == 1) { choice = common_edges[0]; }
     else {
@@ -247,7 +247,7 @@ GeomMeshMatcher::matchFaces(GModel *m1, GModel *m2,
     }
     std::vector<GFace *> common_faces;
     getIntersection<GFace *>(common_faces, lists);
-    GFace *choice = 0;
+    GFace *choice = nullptr;
 
     if(common_faces.size() == 0) {
       Msg::Debug("Could not match face %i (geom).", f1->tag());
@@ -363,7 +363,7 @@ GeomMeshMatcher::matchRegions(GModel *m1, GModel *m2,
       // Then, compute the minimal bounding box
       SOrientedBoundingBox geo_obb = ((GRegion *)*entity1)->getOBB();
 
-      GRegion *choice = 0;
+      GRegion *choice = nullptr;
       double best_score = DBL_MAX;
       // Next, let's iterate over the mesh entities.
       for(auto candidate = common_regions.begin();
@@ -411,7 +411,7 @@ void GeomMeshMatcher::destroy()
 
 static GVertex *getGVertex(MVertex *v1, GModel *gm, const double TOL)
 {
-  GVertex *best = 0;
+  GVertex *best = nullptr;
   double bestScore = TOL;
   for(auto it = gm->firstEdge(); it != gm->lastEdge(); ++it) {
     {
@@ -526,7 +526,7 @@ int GeomMeshMatcher::forceTomatch(GModel *geom, GModel *mesh, const double TOL)
       MVertex *v2 =
         geom->getMeshVertexByTag((*it)->lines[i]->getVertex(1)->getNum());
       if(v1 && v2) {
-        GEdge *ge = 0;
+        GEdge *ge = nullptr;
         if(v1->onWhat()->dim() == 1) ge = (GEdge *)v1->onWhat();
         if(v2->onWhat()->dim() == 1) ge = (GEdge *)v2->onWhat();
         if(ge) {
@@ -610,7 +610,7 @@ static void copy_periodicity(std::vector<Pair<GEType *, GEType *> > &eCor,
     GEType *oldTgt = srcIter->first;
     GEType *oldSrc = dynamic_cast<GEType *>(oldTgt->getMeshMaster());
 
-    if(oldSrc != NULL && oldSrc != oldTgt) {
+    if(oldSrc != nullptr && oldSrc != oldTgt) {
       GEType *newTgt = srcIter->second;
       auto tgtIter = eMap.find(oldSrc);
       if(tgtIter == eMap.end()) {
@@ -861,10 +861,10 @@ int GeomMeshMatcher::match(GModel *geom, GModel *mesh)
 
   bool ok = true;
 
-  std::vector<Pair<GVertex *, GVertex *> > *coresp_v(NULL);
-  std::vector<Pair<GEdge *, GEdge *> > *coresp_e(NULL);
-  std::vector<Pair<GFace *, GFace *> > *coresp_f(NULL);
-  std::vector<Pair<GRegion *, GRegion *> > *coresp_r(NULL);
+  std::vector<Pair<GVertex *, GVertex *> > *coresp_v(nullptr);
+  std::vector<Pair<GEdge *, GEdge *> > *coresp_e(nullptr);
+  std::vector<Pair<GFace *, GFace *> > *coresp_f(nullptr);
+  std::vector<Pair<GRegion *, GRegion *> > *coresp_r(nullptr);
 
   coresp_v = matchVertices(geom, mesh, ok);
   if(ok) {
diff --git a/Geo/Homology.cpp b/Geo/Homology.cpp
index ea6267dc46..f191cfe65c 100644
--- a/Geo/Homology.cpp
+++ b/Geo/Homology.cpp
@@ -16,7 +16,7 @@ Homology::Homology(GModel *model, const std::vector<int> &physicalDomain,
                    int combine, bool omit, bool smoothen, int heuristic)
   : _model(model), _domain(physicalDomain), _subdomain(physicalSubdomain),
     _imdomain(physicalImdomain), _saveOrig(saveOrig), _combine(combine),
-    _omit(omit), _smoothen(smoothen), _heuristic(heuristic), _cellComplex(NULL)
+    _omit(omit), _smoothen(smoothen), _heuristic(heuristic), _cellComplex(nullptr)
 {
   _fileName = "";
 
@@ -107,7 +107,7 @@ void Homology::_createCellComplex()
   _getElements(_nonsubdomainEntities, nonsubdomainElements);
   _getElements(_immuneEntities, immuneElements);
 
-  if(_cellComplex != NULL) delete _cellComplex;
+  if(_cellComplex != nullptr) delete _cellComplex;
   _cellComplex = new CellComplex(_model, domainElements, subdomainElements,
                                  nondomainElements, nonsubdomainElements,
                                  immuneElements, _saveOrig);
@@ -151,7 +151,7 @@ void Homology::_deleteCochains(std::vector<int> dim)
 
 Homology::~Homology()
 {
-  if(_cellComplex != NULL) delete _cellComplex;
+  if(_cellComplex != nullptr) delete _cellComplex;
   _deleteChains();
   _deleteCochains();
 }
@@ -168,7 +168,7 @@ void Homology::findHomologyBasis(std::vector<int> dim)
     return;
   }
 
-  if(_cellComplex == NULL) _createCellComplex();
+  if(_cellComplex == nullptr) _createCellComplex();
   if(_cellComplex->isReduced()) _cellComplex->restoreComplex();
 
   Msg::StatusBar(true, "Reducing cell complex...");
@@ -257,7 +257,7 @@ void Homology::findCohomologyBasis(std::vector<int> dim)
     return;
   }
 
-  if(_cellComplex == NULL) _createCellComplex();
+  if(_cellComplex == nullptr) _createCellComplex();
   if(_cellComplex->isReduced()) _cellComplex->restoreComplex();
 
   Msg::StatusBar(true, "Reducing cell complex...");
@@ -531,7 +531,7 @@ void Homology::getCohomologyBasis(int dim, std::vector<Chain<int> > &coh)
 void Homology::findBettiNumbers()
 {
   if(!isBettiComputed()) {
-    if(_cellComplex == NULL) _createCellComplex();
+    if(_cellComplex == nullptr) _createCellComplex();
     if(_cellComplex->isReduced()) _cellComplex->restoreComplex();
 
     Msg::StatusBar(true, "Reducing cell complex...");
@@ -586,7 +586,7 @@ int Homology::betti(int dim)
 
 int Homology::eulerCharacteristic()
 {
-  if(_cellComplex == NULL) _createCellComplex();
+  if(_cellComplex == nullptr) _createCellComplex();
   return _cellComplex->eulerCharacteristic();
 }
 
diff --git a/Geo/MEdge.cpp b/Geo/MEdge.cpp
index 5b01876814..7a1d743388 100644
--- a/Geo/MEdge.cpp
+++ b/Geo/MEdge.cpp
@@ -78,7 +78,7 @@ bool SortEdgeConsecutive(const std::vector<MEdge> &e,
 
     auto it0 = c.find(v0), it1 = c.find(v1);
     if(it0 == c.end())
-      c[v0] = std::make_pair(v1, (MVertex *)NULL);
+      c[v0] = std::make_pair(v1, (MVertex *)nullptr);
     else {
       if(it0->second.second == NULL) { it0->second.second = v1; }
       else {
@@ -87,7 +87,7 @@ bool SortEdgeConsecutive(const std::vector<MEdge> &e,
       }
     }
     if(it1 == c.end())
-      c[v1] = std::make_pair(v0, (MVertex *)NULL);
+      c[v1] = std::make_pair(v0, (MVertex *)nullptr);
     else {
       if(it1->second.second == NULL) { it1->second.second = v0; }
       else {
@@ -102,7 +102,7 @@ bool SortEdgeConsecutive(const std::vector<MEdge> &e,
 
   while(!c.empty()) {
     std::vector<MVertex *> v;
-    MVertex *start = NULL;
+    MVertex *start = nullptr;
     {
       auto it = c.begin();
       start = it->first;
@@ -127,7 +127,7 @@ bool SortEdgeConsecutive(const std::vector<MEdge> &e,
       }
       v.push_back(current);
       auto it = c.find(current);
-      if(it == c.end() || it->first == NULL) {
+      if(it == c.end() || it->first == nullptr) {
         Msg::Error("Impossible to find %d", current->getNum());
         return false;
       }
@@ -144,7 +144,7 @@ bool SortEdgeConsecutive(const std::vector<MEdge> &e,
       }
       prev = temp;
       if(current == start) { v.push_back(current); }
-    } while(current != start && current != NULL);
+    } while(current != start && current != nullptr);
     if(v.size() > 2 && v[v.size() - 2] == v[v.size() - 1]) {
       v.erase(v.begin() + v.size() - 1);
     }
diff --git a/Geo/MEdge.h b/Geo/MEdge.h
index d8a3f39b7e..05d2cd5ef0 100644
--- a/Geo/MEdge.h
+++ b/Geo/MEdge.h
@@ -19,7 +19,7 @@ private:
 public:
   MEdge()
   {
-    _v[0] = _v[1] = 0;
+    _v[0] = _v[1] = nullptr;
     _si[0] = _si[1] = 0;
   }
   MEdge(MVertex *v0, MVertex *v1)
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 42f16a11bc..6cfda67c8d 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -666,7 +666,7 @@ const nodalBasis *MElement::getFunctionSpace(int order, bool serendip) const
 {
   if(order == -1) return BasisFactory::getNodalBasis(getTypeForMSH());
   int type = ElementType::getType(getType(), order, serendip);
-  return type ? BasisFactory::getNodalBasis(type) : NULL;
+  return type ? BasisFactory::getNodalBasis(type) : nullptr;
 }
 
 const FuncSpaceData MElement::getFuncSpaceData(int order, bool serendip) const
@@ -679,7 +679,7 @@ const JacobianBasis *MElement::getJacobianFuncSpace(int orderElement) const
 {
   if(orderElement == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
   int tag = ElementType::getType(getType(), orderElement);
-  return tag ? BasisFactory::getJacobianBasis(tag) : NULL;
+  return tag ? BasisFactory::getJacobianBasis(tag) : nullptr;
 }
 
 const FuncSpaceData MElement::getJacobianFuncSpaceData(int orderElement) const
@@ -2372,7 +2372,7 @@ MElement *MElement::copy(std::map<int, MVertex *> &vertexMap,
       if(vertexMap.count(numV))
         vmv.push_back(vertexMap[numV]);
       else {
-        MVertex *mv = new MVertex(v->x(), v->y(), v->z(), 0, numV);
+        MVertex *mv = new MVertex(v->x(), v->y(), v->z(), nullptr, numV);
         vmv.push_back(mv);
         vertexMap[numV] = mv;
       }
@@ -2386,7 +2386,7 @@ MElement *MElement::copy(std::map<int, MVertex *> &vertexMap,
         if(vertexMap.count(numV))
           vmv.push_back(vertexMap[numV]);
         else {
-          MVertex *mv = new MVertex(v->x(), v->y(), v->z(), 0, numV);
+          MVertex *mv = new MVertex(v->x(), v->y(), v->z(), nullptr, numV);
           vmv.push_back(mv);
           vertexMap[numV] = mv;
         }
@@ -2394,7 +2394,7 @@ MElement *MElement::copy(std::map<int, MVertex *> &vertexMap,
     }
   }
 
-  MElement *parent = 0;
+  MElement *parent = nullptr;
   if(eParent && !getDomain(0) && !getDomain(1)) {
     auto it = newParents.find(eParent);
     MElement *newParent;
@@ -2557,7 +2557,7 @@ MElement *MElementFactory::create(int type, std::vector<MVertex *> &v,
   case MSH_PYR_285: return new MPyramidN(v, 8, num, part);
   case MSH_PYR_385: return new MPyramidN(v, 9, num, part);
   case MSH_TRIH_4: return new MTrihedron(v, num, part);
-  default: return 0;
+  default: return nullptr;
   }
 }
 
@@ -2582,13 +2582,13 @@ MElement *MElementFactory::create(int num, int type,
       if(v) { vertices[i] = v; }
       else {
         Msg::Error("Unknown node %d in element %d", numVertex, num);
-        return 0;
+        return nullptr;
       }
     }
   }
   else {
     Msg::Error("Missing data in element %d", num);
-    return 0;
+    return nullptr;
   }
 
   unsigned int part = 0;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 76087621a7..1130700b34 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -218,11 +218,11 @@ public:
   }
 
   // get and set parent and children for hierarchial grids
-  virtual MElement *getParent() const { return NULL; }
+  virtual MElement *getParent() const { return nullptr; }
   virtual void setParent(MElement *p, bool owner = false) {}
   virtual void updateParent(GModel *gm) {} // used only for reading msh files
   virtual int getNumChildren() const { return 0; }
-  virtual MElement *getChild(int i) const { return NULL; }
+  virtual MElement *getChild(int i) const { return nullptr; }
   virtual bool ownsParent() const { return false; }
 
   // get base element in case of MSubElement
@@ -230,7 +230,7 @@ public:
   virtual MElement *getBaseElement() { return this; }
 
   // get and set domain for borders
-  virtual MElement *getDomain(int i) const { return NULL; }
+  virtual MElement *getDomain(int i) const { return nullptr; }
   virtual void setDomain(MElement *e, int i) {}
 
   // get the type of the element
@@ -270,10 +270,10 @@ public:
   double minScaledJacobian(bool knownValid = false, bool reversedOk = false);
   virtual double angleShapeMeasure() { return 1.0; }
   virtual void scaledJacRange(double &jmin, double &jmax,
-                              GEntity *ge = 0) const;
-  virtual void idealJacRange(double &jmin, double &jmax, GEntity *ge = 0);
+                              GEntity *ge = nullptr) const;
+  virtual void idealJacRange(double &jmin, double &jmax, GEntity *ge = nullptr);
   virtual void signedInvCondNumRange(double &iCNMin, double &iCNMax,
-                                     GEntity *ge = 0);
+                                     GEntity *ge = nullptr);
   virtual void signedInvGradErrorRange(double &minSIGE, double &maxSIGE);
 
   // get the radius of the inscribed circle/sphere if it exists, otherwise get
@@ -411,7 +411,8 @@ public:
   double interpolate(double val[], double u, double v, double w, int stride = 1,
                      int order = -1);
   void interpolateGrad(double val[], double u, double v, double w, double f[],
-                       int stride = 1, double invjac[3][3] = 0, int order = -1);
+                       int stride = 1, double invjac[3][3] = nullptr,
+                       int order = -1);
   void interpolateCurl(double val[], double u, double v, double w, double f[],
                        int stride = 3, int order = -1);
   double interpolateDiv(double val[], double u, double v, double w,
@@ -423,7 +424,7 @@ public:
     Msg::Error("No integration points defined for this type of element: %d",
                this->getType());
     *npts = 0;
-    *pts = 0;
+    *pts = nullptr;
   }
   double integrate(double val[], int pOrder, int stride = 1, int order = -1);
   // val[] must contain interpolation data for face/edge vertices of given
@@ -435,9 +436,9 @@ public:
   virtual void writeMSH2(FILE *fp, double version = 1.0, bool binary = false,
                          int num = 0, int elementary = 1, int physical = 1,
                          int parentNum = 0, int dom1Num = 0, int dom2Num = 0,
-                         std::vector<short> *ghosts = 0);
+                         std::vector<short> *ghosts = nullptr);
   virtual void writeMSH3(FILE *fp, bool binary = false, int elementary = 1,
-                         std::vector<short> *ghosts = 0);
+                         std::vector<short> *ghosts = nullptr);
   virtual void writePOS(FILE *fp, bool printElementary, bool printElementNumber,
                         bool printSICN, bool printSIGE, bool printGamma,
                         bool printDisto, double scalingFactor = 1.0,
@@ -472,16 +473,16 @@ public:
   virtual int getTypeForMSH() const { return 0; }
   virtual int getTypeForUNV() const { return 0; }
   virtual int getTypeForVTK() const { return 0; }
-  virtual const char *getStringForTOCHNOG() const { return 0; }
-  virtual const char *getStringForPOS() const { return 0; }
-  virtual const char *getStringForBDF() const { return 0; }
-  virtual const char *getStringForDIFF() const { return 0; }
-  virtual const char *getStringForINP() const { return 0; }
-  virtual const char *getStringForKEY() const { return 0; }
+  virtual const char *getStringForTOCHNOG() const { return nullptr; }
+  virtual const char *getStringForPOS() const { return nullptr; }
+  virtual const char *getStringForBDF() const { return nullptr; }
+  virtual const char *getStringForDIFF() const { return nullptr; }
+  virtual const char *getStringForINP() const { return nullptr; }
+  virtual const char *getStringForKEY() const { return nullptr; }
 
   // return the number of vertices, as well as the element name if 'name' != 0
   static unsigned int getInfoMSH(const int typeMSH,
-                                 const char **const name = 0);
+                                 const char **const name = nullptr);
   std::string getName();
   virtual std::size_t getNumVerticesForMSH() { return getNumVertices(); }
   virtual void getVerticesIdForMSH(std::vector<int> &verts);
@@ -501,8 +502,8 @@ class MElementFactory {
 public:
   MElement *create(int type, std::vector<MVertex *> &v, std::size_t num = 0,
                    int part = 0, bool owner = false, int parent = 0,
-                   MElement *parent_ptr = NULL, MElement *d1 = 0,
-                   MElement *d2 = 0);
+                   MElement *parent_ptr = nullptr, MElement *d1 = nullptr,
+                   MElement *d2 = nullptr);
   MElement *create(int num, int type, const std::vector<int> &data,
                    GModel *model);
 };
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index cc62fc6eff..c28fd9ff6d 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -979,12 +979,12 @@ static void elementCutMesh(
         isCut = (isCut || iC);
       }
     }
-    MPolyhedron *p1 = NULL, *p2 = NULL;
+    MPolyhedron *p1 = nullptr, *p2 = nullptr;
     if(isCut) {
       std::vector<MTetrahedron *> poly[2];
 
       for(std::size_t i = nbTe; i < tetras.size(); i++) {
-        MVertex *mv[4] = {NULL, NULL, NULL, NULL};
+        MVertex *mv[4] = {nullptr, nullptr, nullptr, nullptr};
         for(int j = 0; j < 4; j++) {
           int numV = getElementVertexNum(tetras[i]->pt(j), e);
           if(numV == -1) {
@@ -999,7 +999,7 @@ static void elementCutMesh(
             auto it = vertexMap.find(numV);
             if(it == vertexMap.end()) {
               mv[j] = new MVertex(tetras[i]->x(j), tetras[i]->y(j),
-                                  tetras[i]->z(j), 0, numV);
+                                  tetras[i]->z(j), nullptr, numV);
               vertexMap[numV] = mv[j];
             }
             else
@@ -1053,7 +1053,7 @@ static void elementCutMesh(
       borders[1].erase(itr.first, itr.second);
       for(std::size_t i = 0; i < bords.size(); i++) borders[1].insert(bords[i]);
       if(eParent) {
-        copy->setParent(NULL, false);
+        copy->setParent(nullptr, false);
         delete copy;
       }
     }
@@ -1080,7 +1080,7 @@ static void elementCutMesh(
     }
 
     for(std::size_t i = nbTr; i < triangles.size(); i++) {
-      MVertex *mv[3] = {NULL, NULL, NULL};
+      MVertex *mv[3] = {nullptr, nullptr, nullptr};
       for(int j = 0; j < 3; j++) {
         int numV = getElementVertexNum(triangles[i]->pt(j), e);
         if(numV == -1) {
@@ -1095,7 +1095,7 @@ static void elementCutMesh(
           auto it = vertexMap.find(numV);
           if(it == vertexMap.end()) {
             mv[j] = new MVertex(triangles[i]->x(j), triangles[i]->y(j),
-                                triangles[i]->z(j), 0, numV);
+                                triangles[i]->z(j), nullptr, numV);
             vertexMap[numV] = mv[j];
           }
           else
@@ -1179,12 +1179,12 @@ static void elementCutMesh(
       }
     }
 
-    MPolygon *p1 = NULL, *p2 = NULL;
+    MPolygon *p1 = nullptr, *p2 = nullptr;
     if(isCut) {
       std::vector<MTriangle *> poly[2];
 
       for(std::size_t i = nbTr; i < triangles.size(); i++) {
-        MVertex *mv[3] = {NULL, NULL, NULL};
+        MVertex *mv[3] = {nullptr, nullptr, nullptr};
         for(int j = 0; j < 3; j++) {
           int numV = getElementVertexNum(triangles[i]->pt(j), e);
           if(numV == -1) {
@@ -1199,7 +1199,7 @@ static void elementCutMesh(
             auto it = vertexMap.find(numV);
             if(it == vertexMap.end()) {
               mv[j] = new MVertex(triangles[i]->x(j), triangles[i]->y(j),
-                                  triangles[i]->z(j), 0, numV);
+                                  triangles[i]->z(j), nullptr, numV);
               vertexMap[numV] = mv[j];
             }
             else
@@ -1214,7 +1214,7 @@ static void elementCutMesh(
       }
       // if quads
       for(std::size_t i = nbQ; i < quads.size(); i++) {
-        MVertex *mv[4] = {NULL, NULL, NULL, NULL};
+        MVertex *mv[4] = {nullptr, nullptr, nullptr, nullptr};
         for(int j = 0; j < 4; j++) {
           int numV = getElementVertexNum(quads[i]->pt(j), e);
           if(numV == -1) {
@@ -1229,7 +1229,7 @@ static void elementCutMesh(
             auto it = vertexMap.find(numV);
             if(it == vertexMap.end()) {
               mv[j] = new MVertex(quads[i]->x(j), quads[i]->y(j),
-                                  quads[i]->z(j), 0, numV);
+                                  quads[i]->z(j), nullptr, numV);
               vertexMap[numV] = mv[j];
             }
             else
@@ -1304,7 +1304,7 @@ static void elementCutMesh(
       borders[0].erase(itr.first, itr.second);
       for(std::size_t i = 0; i < bords.size(); i++) borders[0].insert(bords[i]);
       if(eParent) {
-        copy->setParent(NULL, false);
+        copy->setParent(nullptr, false);
         delete copy;
       }
     }
@@ -1331,7 +1331,7 @@ static void elementCutMesh(
     }
 
     for(std::size_t i = nbL; i < lines.size(); i++) {
-      MVertex *mv[2] = {NULL, NULL};
+      MVertex *mv[2] = {nullptr, nullptr};
       for(int j = 0; j < 2; j++) {
         int numV = getElementVertexNum(lines[i]->pt(j), e);
         if(numV == -1) {
@@ -1346,7 +1346,7 @@ static void elementCutMesh(
           auto it = vertexMap.find(numV);
           if(it == vertexMap.end()) {
             mv[j] = new MVertex(lines[i]->x(j), lines[i]->y(j), lines[i]->z(j),
-                                0, numV);
+                                nullptr, numV);
             vertexMap[numV] = mv[j];
           }
           else
@@ -1396,7 +1396,7 @@ static void elementCutMesh(
     if(isCut) {
       bool own = (eParent && !e->ownsParent()) ? false : true;
       for(std::size_t i = nbL; i < lines.size(); i++) {
-        MVertex *mv[2] = {NULL, NULL};
+        MVertex *mv[2] = {nullptr, nullptr};
         for(int j = 0; j < 2; j++) {
           int numV = getElementVertexNum(lines[i]->pt(j), e);
           if(numV == -1) {
@@ -1411,7 +1411,7 @@ static void elementCutMesh(
             auto it = vertexMap.find(numV);
             if(it == vertexMap.end()) {
               mv[j] = new MVertex(lines[i]->x(j), lines[i]->y(j),
-                                  lines[i]->z(j), 0, numV);
+                                  lines[i]->z(j), nullptr, numV);
               vertexMap[numV] = mv[j];
             }
             else
@@ -1439,7 +1439,7 @@ static void elementCutMesh(
               std::pair<MElement *, MElement *>(ml->getDomain(i), ml));
       }
       if(eParent) {
-        copy->setParent(NULL, false);
+        copy->setParent(nullptr, false);
         delete copy;
       }
     }
@@ -1532,7 +1532,7 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
       if(primS > 1)
         verticesLs(k, vi->getIndex()) = (*ls)(vi->x(), vi->y(), vi->z());
 
-      MVertex *vn = new MVertex(vi->x(), vi->y(), vi->z(), 0, vi->getNum());
+      MVertex *vn = new MVertex(vi->x(), vi->y(), vi->z(), nullptr, vi->getNum());
       vertexMap[vi->getNum()] = vn;
     }
   }
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index 757b2a5fa5..0b3befc518 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -32,8 +32,8 @@ protected:
 
 public:
   MPolyhedron(std::vector<MVertex *> v, int num = 0, int part = 0,
-              bool owner = false, MElement *orig = NULL)
-    : MElement(num, part), _owner(owner), _orig(orig), _intpt(0)
+              bool owner = false, MElement *orig = nullptr)
+    : MElement(num, part), _owner(owner), _orig(orig), _intpt(nullptr)
   {
     if(v.size() % 4) {
       Msg::Error("Got %d nodes for polyhedron", (int)v.size());
@@ -44,8 +44,8 @@ public:
     _init();
   }
   MPolyhedron(std::vector<MTetrahedron *> vT, int num = 0, int part = 0,
-              bool owner = false, MElement *orig = NULL)
-    : MElement(num, part), _owner(owner), _orig(orig), _intpt(0)
+              bool owner = false, MElement *orig = nullptr)
+    : MElement(num, part), _owner(owner), _orig(orig), _intpt(nullptr)
   {
     for(std::size_t i = 0; i < vT.size(); i++) _parts.push_back(vT[i]);
     _init();
@@ -129,11 +129,11 @@ public:
   virtual const nodalBasis *getFunctionSpace(int order = -1,
                                              bool serendip = false) const
   {
-    return (_orig ? _orig->getFunctionSpace(order, serendip) : 0);
+    return (_orig ? _orig->getFunctionSpace(order, serendip) : nullptr);
   }
   virtual const JacobianBasis *getJacobianFuncSpace(int order = -1) const
   {
-    return (_orig ? _orig->getJacobianFuncSpace(order) : 0);
+    return (_orig ? _orig->getJacobianFuncSpace(order) : nullptr);
   }
   virtual void getShapeFunctions(double u, double v, double w, double s[],
                                  int o) const
@@ -160,11 +160,11 @@ public:
   }
   virtual const MVertex *getShapeFunctionNode(int i) const
   {
-    return (_orig ? _orig->getShapeFunctionNode(i) : 0);
+    return (_orig ? _orig->getShapeFunctionNode(i) : nullptr);
   }
   virtual MVertex *getShapeFunctionNode(int i)
   {
-    return (_orig ? _orig->getShapeFunctionNode(i) : 0);
+    return (_orig ? _orig->getShapeFunctionNode(i) : nullptr);
   }
 
   // the parametric coordinates of the polyhedron are
@@ -210,16 +210,16 @@ protected:
 
 public:
   MPolygon(std::vector<MVertex *> v, int num = 0, int part = 0,
-           bool owner = false, MElement *orig = NULL)
-    : MElement(num, part), _owner(owner), _orig(orig), _intpt(0)
+           bool owner = false, MElement *orig = nullptr)
+    : MElement(num, part), _owner(owner), _orig(orig), _intpt(nullptr)
   {
     for(std::size_t i = 0; i < v.size() / 3; i++)
       _parts.push_back(new MTriangle(v[i * 3], v[i * 3 + 1], v[i * 3 + 2]));
     _initVertices();
   }
   MPolygon(std::vector<MTriangle *> vT, int num = 0, int part = 0,
-           bool owner = false, MElement *orig = NULL)
-    : MElement(num, part), _owner(owner), _orig(orig), _intpt(0)
+           bool owner = false, MElement *orig = nullptr)
+    : MElement(num, part), _owner(owner), _orig(orig), _intpt(nullptr)
   {
     for(std::size_t i = 0; i < vT.size(); i++) {
       MTriangle *t = (MTriangle *)vT[i];
@@ -304,11 +304,11 @@ public:
   virtual const nodalBasis *getFunctionSpace(int order = -1,
                                              bool serendip = false) const
   {
-    return (_orig ? _orig->getFunctionSpace(order, serendip) : 0);
+    return (_orig ? _orig->getFunctionSpace(order, serendip) : nullptr);
   }
   virtual const JacobianBasis *getJacobianFuncSpace(int order = -1) const
   {
-    return (_orig ? _orig->getJacobianFuncSpace(order) : 0);
+    return (_orig ? _orig->getJacobianFuncSpace(order) : nullptr);
   }
   virtual void getShapeFunctions(double u, double v, double w, double s[],
                                  int o) const
@@ -335,11 +335,11 @@ public:
   }
   virtual const MVertex *getShapeFunctionNode(int i) const
   {
-    return (_orig ? _orig->getShapeFunctionNode(i) : 0);
+    return (_orig ? _orig->getShapeFunctionNode(i) : nullptr);
   }
   virtual MVertex *getShapeFunctionNode(int i)
   {
-    return (_orig ? _orig->getShapeFunctionNode(i) : 0);
+    return (_orig ? _orig->getShapeFunctionNode(i) : nullptr);
   }
 
   // the parametric coordinates of the polygon are
@@ -371,13 +371,13 @@ protected:
 
 public:
   MLineChild(MVertex *v0, MVertex *v1, int num = 0, int part = 0,
-             bool owner = false, MElement *orig = NULL)
-    : MLine(v0, v1, num, part), _owner(owner), _orig(orig), _intpt(0)
+             bool owner = false, MElement *orig = nullptr)
+    : MLine(v0, v1, num, part), _owner(owner), _orig(orig), _intpt(nullptr)
   {
   }
   MLineChild(const std::vector<MVertex *> &v, int num = 0, int part = 0,
-             bool owner = false, MElement *orig = NULL)
-    : MLine(v, num, part), _owner(owner), _orig(orig), _intpt(0)
+             bool owner = false, MElement *orig = nullptr)
+    : MLine(v, num, part), _owner(owner), _orig(orig), _intpt(nullptr)
   {
   }
   ~MLineChild()
@@ -389,12 +389,12 @@ public:
                                              bool serendip = false) const
   {
     if(_orig) return _orig->getFunctionSpace(order, serendip);
-    return 0;
+    return nullptr;
   }
   virtual const JacobianBasis *getJacobianFuncSpace(int order = -1) const
   {
     if(_orig) return _orig->getJacobianFuncSpace(order);
-    return 0;
+    return nullptr;
   }
   virtual void getShapeFunctions(double u, double v, double w, double s[],
                                  int o) const
@@ -433,15 +433,15 @@ protected:
 
 public:
   MTriangleBorder(MVertex *v0, MVertex *v1, MVertex *v2, int num = 0,
-                  int part = 0, MElement *d1 = NULL, MElement *d2 = NULL)
-    : MTriangle(v0, v1, v2, num, part), _intpt(0)
+                  int part = 0, MElement *d1 = nullptr, MElement *d2 = nullptr)
+    : MTriangle(v0, v1, v2, num, part), _intpt(nullptr)
   {
     _domains[0] = d1;
     _domains[1] = d2;
   }
   MTriangleBorder(const std::vector<MVertex *> &v, int num = 0, int part = 0,
-                  MElement *d1 = NULL, MElement *d2 = NULL)
-    : MTriangle(v, num, part), _intpt(0)
+                  MElement *d1 = nullptr, MElement *d2 = nullptr)
+    : MTriangle(v, num, part), _intpt(nullptr)
   {
     _domains[0] = d1;
     _domains[1] = d2;
@@ -453,7 +453,7 @@ public:
   {
     if(_domains[0]) return _domains[0]->getParent();
     if(_domains[1]) return _domains[1]->getParent();
-    return NULL;
+    return nullptr;
   }
   virtual int getTypeForMSH() const { return MSH_TRI_B; }
   virtual bool isInside(double u, double v, double w) const;
@@ -469,17 +469,17 @@ protected:
 
 public:
   MPolygonBorder(const std::vector<MTriangle *> &v, int num = 0, int part = 0,
-                 bool own = false, MElement *p = NULL, MElement *d1 = NULL,
-                 MElement *d2 = NULL)
-    : MPolygon(v, num, part, own, p), _intpt(0)
+                 bool own = false, MElement *p = nullptr,
+                 MElement *d1 = nullptr, MElement *d2 = nullptr)
+    : MPolygon(v, num, part, own, p), _intpt(nullptr)
   {
     _domains[0] = d1;
     _domains[1] = d2;
   }
   MPolygonBorder(const std::vector<MVertex *> &v, int num = 0, int part = 0,
-                 bool own = false, MElement *p = NULL, MElement *d1 = NULL,
-                 MElement *d2 = NULL)
-    : MPolygon(v, num, part, own, p), _intpt(0)
+                 bool own = false, MElement *p = nullptr,
+                 MElement *d1 = nullptr, MElement *d2 = nullptr)
+    : MPolygon(v, num, part, own, p), _intpt(nullptr)
   {
     _domains[0] = d1;
     _domains[1] = d2;
@@ -491,7 +491,7 @@ public:
   {
     if(_domains[0]) return _domains[0]->getParent();
     if(_domains[1]) return _domains[1]->getParent();
-    return NULL;
+    return nullptr;
   }
   virtual int getTypeForMSH() const { return MSH_POLYG_B; }
 };
@@ -503,15 +503,15 @@ protected:
 
 public:
   MLineBorder(MVertex *v0, MVertex *v1, int num = 0, int part = 0,
-              MElement *d1 = NULL, MElement *d2 = NULL)
-    : MLine(v0, v1, num, part), _intpt(0)
+              MElement *d1 = nullptr, MElement *d2 = nullptr)
+    : MLine(v0, v1, num, part), _intpt(nullptr)
   {
     _domains[0] = d1;
     _domains[1] = d2;
   }
   MLineBorder(const std::vector<MVertex *> &v, int num = 0, int part = 0,
-              MElement *d1 = NULL, MElement *d2 = NULL)
-    : MLine(v, num, part), _intpt(0)
+              MElement *d1 = nullptr, MElement *d2 = nullptr)
+    : MLine(v, num, part), _intpt(nullptr)
   {
     _domains[0] = d1;
     _domains[1] = d2;
@@ -523,7 +523,7 @@ public:
   {
     if(_domains[0]) return _domains[0]->getParent();
     if(_domains[1]) return _domains[1]->getParent();
-    return NULL;
+    return nullptr;
   }
   virtual int getTypeForMSH() const { return MSH_LIN_B; }
   virtual bool isInside(double u, double v, double w) const;
diff --git a/Geo/MElementOctree.cpp b/Geo/MElementOctree.cpp
index 339ae31f62..36dd85f23d 100644
--- a/Geo/MElementOctree.cpp
+++ b/Geo/MElementOctree.cpp
@@ -116,7 +116,7 @@ MElementOctree::MElementOctree(GModel *m) : _gm(m)
 }
 
 MElementOctree::MElementOctree(const std::vector<MElement *> &v)
-  : _gm(0), _elems(v)
+  : _gm(nullptr), _elems(v)
 {
   SBoundingBox3d bb;
   for(std::size_t i = 0; i < v.size(); i++) {
@@ -258,5 +258,5 @@ MElement *MElementOctree::find(double x, double y, double z, int dim,
     MElement::setTolerance(initialTol);
     // Msg::Warning("Point %g %g %g not found",x,y,z);
   }
-  return NULL;
+  return nullptr;
 }
diff --git a/Geo/MFace.h b/Geo/MFace.h
index 625385d335..26fc3e1877 100644
--- a/Geo/MFace.h
+++ b/Geo/MFace.h
@@ -24,7 +24,7 @@ private:
 
 public:
   MFace() {}
-  MFace(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3 = 0);
+  MFace(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3 = nullptr);
   MFace(const std::vector<MVertex *> &v);
   std::size_t getNumVertices() const { return _v.size(); }
   MVertex *getVertex(std::size_t i) const { return _v[i]; }
diff --git a/Geo/MSubElement.cpp b/Geo/MSubElement.cpp
index 1249398752..8b11536053 100644
--- a/Geo/MSubElement.cpp
+++ b/Geo/MSubElement.cpp
@@ -27,13 +27,13 @@ const nodalBasis *MSubTetrahedron::getFunctionSpace(int order,
                                                     bool serendip) const
 {
   if(_orig) return _orig->getFunctionSpace(order, serendip);
-  return 0;
+  return nullptr;
 }
 
 const JacobianBasis *MSubTetrahedron::getJacobianFuncSpace(int order) const
 {
   if(_orig) return _orig->getJacobianFuncSpace(order);
-  return 0;
+  return nullptr;
 }
 
 void MSubTetrahedron::getShapeFunctions(double u, double v, double w,
@@ -102,13 +102,13 @@ std::size_t MSubTetrahedron::getNumPrimaryShapeFunctions() const
 const MVertex *MSubTetrahedron::getShapeFunctionNode(int i) const
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
-  return 0;
+  return nullptr;
 }
 
 MVertex *MSubTetrahedron::getShapeFunctionNode(int i)
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
-  return 0;
+  return nullptr;
 }
 
 void MSubTetrahedron::xyz2uvw(double xyz[3], double uvw[3]) const
@@ -236,13 +236,13 @@ void MSubTriangle::updateParent(GModel *gm)
 const nodalBasis *MSubTriangle::getFunctionSpace(int order, bool serendip) const
 {
   if(_orig) return _orig->getFunctionSpace(order, serendip);
-  return 0;
+  return nullptr;
 }
 
 const JacobianBasis *MSubTriangle::getJacobianFuncSpace(int order) const
 {
   if(_orig) return _orig->getJacobianFuncSpace(order);
-  return 0;
+  return nullptr;
 }
 
 void MSubTriangle::getShapeFunctions(double u, double v, double w, double s[],
@@ -363,13 +363,13 @@ std::size_t MSubTriangle::getNumPrimaryShapeFunctions() const
 const MVertex *MSubTriangle::getShapeFunctionNode(int i) const
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
-  return 0;
+  return nullptr;
 }
 
 MVertex *MSubTriangle::getShapeFunctionNode(int i)
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
-  return 0;
+  return nullptr;
 }
 
 void MSubTriangle::xyz2uvw(double xyz[3], double uvw[3]) const
@@ -494,13 +494,13 @@ void MSubLine::updateParent(GModel *gm)
 const nodalBasis *MSubLine::getFunctionSpace(int order, bool serendip) const
 {
   if(_orig) return _orig->getFunctionSpace(order, serendip);
-  return 0;
+  return nullptr;
 }
 
 const JacobianBasis *MSubLine::getJacobianFuncSpace(int order) const
 {
   if(_orig) return _orig->getJacobianFuncSpace(order);
-  return 0;
+  return nullptr;
 }
 
 void MSubLine::getShapeFunctions(double u, double v, double w, double s[],
@@ -610,13 +610,13 @@ std::size_t MSubLine::getNumPrimaryShapeFunctions() const
 const MVertex *MSubLine::getShapeFunctionNode(int i) const
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
-  return 0;
+  return nullptr;
 }
 
 MVertex *MSubLine::getShapeFunctionNode(int i)
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
-  return 0;
+  return nullptr;
 }
 
 void MSubLine::xyz2uvw(double xyz[3], double uvw[3]) const
@@ -742,13 +742,13 @@ void MSubPoint::updateParent(GModel *gm)
 const nodalBasis *MSubPoint::getFunctionSpace(int order, bool serendip) const
 {
   if(_orig) return _orig->getFunctionSpace(order, serendip);
-  return 0;
+  return nullptr;
 }
 
 const JacobianBasis *MSubPoint::getJacobianFuncSpace(int order) const
 {
   if(_orig) return _orig->getJacobianFuncSpace(order);
-  return 0;
+  return nullptr;
 }
 
 void MSubPoint::getShapeFunctions(double u, double v, double w, double s[],
@@ -816,13 +816,13 @@ std::size_t MSubPoint::getNumPrimaryShapeFunctions() const
 const MVertex *MSubPoint::getShapeFunctionNode(int i) const
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
-  return 0;
+  return nullptr;
 }
 
 MVertex *MSubPoint::getShapeFunctionNode(int i)
 {
   if(_orig) return _orig->getShapeFunctionNode(i);
-  return 0;
+  return nullptr;
 }
 
 void MSubPoint::xyz2uvw(double xyz[3], double uvw[3]) const
diff --git a/Geo/MSubElement.h b/Geo/MSubElement.h
index 9481f4dc23..79b2c0e934 100644
--- a/Geo/MSubElement.h
+++ b/Geo/MSubElement.h
@@ -36,8 +36,8 @@ protected:
 
   MSubTetrahedron(const std::vector<MVertex *> &v, int num, int part,
                   bool owner, int orig)
-    : MTetrahedron(v, num, part), _owner(owner), _orig_N(orig), _base(0),
-      _pOrder(-1), _npts(0), _pts(0)
+    : MTetrahedron(v, num, part), _owner(owner), _orig_N(orig), _base(nullptr),
+      _pOrder(-1), _npts(0), _pts(nullptr)
   {
   }
   virtual void updateParent(
@@ -46,21 +46,21 @@ protected:
 public:
   MSubTetrahedron(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3,
                   int num = 0, int part = 0, bool owner = false,
-                  MElement *orig = NULL)
+                  MElement *orig = nullptr)
     : MTetrahedron(v0, v1, v2, v3, num, part), _owner(owner), _orig(orig),
-      _base(0), _pOrder(-1), _npts(0), _pts(0)
+      _base(nullptr), _pOrder(-1), _npts(0), _pts(nullptr)
   {
   }
   MSubTetrahedron(const std::vector<MVertex *> &v, int num = 0, int part = 0,
-                  bool owner = false, MElement *orig = NULL)
-    : MTetrahedron(v, num, part), _owner(owner), _orig(orig), _base(0),
-      _pOrder(-1), _npts(0), _pts(0)
+                  bool owner = false, MElement *orig = nullptr)
+    : MTetrahedron(v, num, part), _owner(owner), _orig(orig), _base(nullptr),
+      _pOrder(-1), _npts(0), _pts(nullptr)
   {
   }
   MSubTetrahedron(const MTetrahedron &tet, bool owner = false,
-                  MElement *orig = NULL)
-    : MTetrahedron(tet), _owner(owner), _orig(orig), _base(0), _pOrder(-1),
-      _npts(0), _pts(0)
+                  MElement *orig = nullptr)
+    : MTetrahedron(tet), _owner(owner), _orig(orig), _base(nullptr),
+      _pOrder(-1), _npts(0), _pts(nullptr)
   {
   }
   ~MSubTetrahedron();
@@ -146,28 +146,29 @@ protected:
 
   MSubTriangle(const std::vector<MVertex *> &v, int num, int part, bool owner,
                int orig)
-    : MTriangle(v, num, part), _owner(owner), _orig_N(orig), _base(0),
-      _pOrder(-1), _npts(0), _pts(0)
+    : MTriangle(v, num, part), _owner(owner), _orig_N(orig), _base(nullptr),
+      _pOrder(-1), _npts(0), _pts(nullptr)
   {
   }
   virtual void updateParent(
     GModel *gm); // NEVER ever use this ! (except for reading msh files !)
 public:
   MSubTriangle(MVertex *v0, MVertex *v1, MVertex *v2, int num = 0, int part = 0,
-               bool owner = false, MElement *orig = NULL)
-    : MTriangle(v0, v1, v2, num, part), _owner(owner), _orig(orig), _base(0),
-      _pOrder(-1), _npts(0), _pts(0)
+               bool owner = false, MElement *orig = nullptr)
+    : MTriangle(v0, v1, v2, num, part), _owner(owner), _orig(orig),
+      _base(nullptr), _pOrder(-1), _npts(0), _pts(nullptr)
   {
   }
   MSubTriangle(const std::vector<MVertex *> &v, int num = 0, int part = 0,
-               bool owner = false, MElement *orig = NULL)
-    : MTriangle(v, num, part), _owner(owner), _orig(orig), _base(0),
-      _pOrder(-1), _npts(0), _pts(0)
+               bool owner = false, MElement *orig = nullptr)
+    : MTriangle(v, num, part), _owner(owner), _orig(orig), _base(nullptr),
+      _pOrder(-1), _npts(0), _pts(nullptr)
   {
   }
-  MSubTriangle(const MTriangle &tri, bool owner = false, MElement *orig = NULL)
-    : MTriangle(tri), _owner(owner), _orig(orig), _base(0), _pOrder(-1),
-      _npts(0), _pts(0)
+  MSubTriangle(const MTriangle &tri, bool owner = false,
+               MElement *orig = nullptr)
+    : MTriangle(tri), _owner(owner), _orig(orig), _base(nullptr), _pOrder(-1),
+      _npts(0), _pts(nullptr)
   {
   }
   ~MSubTriangle();
@@ -253,28 +254,28 @@ protected:
 
   MSubLine(const std::vector<MVertex *> &v, int num, int part, bool owner,
            int orig)
-    : MLine(v, num, part), _owner(owner), _orig_N(orig), _base(0), _pOrder(-1),
-      _npts(0), _pts(0)
+    : MLine(v, num, part), _owner(owner), _orig_N(orig), _base(nullptr),
+      _pOrder(-1), _npts(0), _pts(nullptr)
   {
   }
   virtual void updateParent(
     GModel *gm); // NEVER ever use this ! (except for reading msh files !)
 public:
   MSubLine(MVertex *v0, MVertex *v1, int num = 0, int part = 0,
-           bool owner = false, MElement *orig = NULL)
-    : MLine(v0, v1, num, part), _owner(owner), _orig(orig), _base(0),
-      _pOrder(-1), _npts(0), _pts(0)
+           bool owner = false, MElement *orig = nullptr)
+    : MLine(v0, v1, num, part), _owner(owner), _orig(orig), _base(nullptr),
+      _pOrder(-1), _npts(0), _pts(nullptr)
   {
   }
   MSubLine(const std::vector<MVertex *> &v, int num = 0, int part = 0,
-           bool owner = false, MElement *orig = NULL)
-    : MLine(v, num, part), _owner(owner), _orig(orig), _base(0), _pOrder(-1),
-      _npts(0), _pts(0)
+           bool owner = false, MElement *orig = nullptr)
+    : MLine(v, num, part), _owner(owner), _orig(orig), _base(nullptr),
+      _pOrder(-1), _npts(0), _pts(nullptr)
   {
   }
-  MSubLine(const MLine &lin, bool owner = false, MElement *orig = NULL)
-    : MLine(lin), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0),
-      _pts(0)
+  MSubLine(const MLine &lin, bool owner = false, MElement *orig = nullptr)
+    : MLine(lin), _owner(owner), _orig(orig), _base(nullptr), _pOrder(-1),
+      _npts(0), _pts(nullptr)
   {
   }
   ~MSubLine();
@@ -360,28 +361,28 @@ protected:
 
   MSubPoint(const std::vector<MVertex *> &v, int num, int part, bool owner,
             int orig)
-    : MPoint(v, num, part), _owner(owner), _orig_N(orig), _base(0), _pOrder(-1),
-      _npts(0), _pts(0)
+    : MPoint(v, num, part), _owner(owner), _orig_N(orig), _base(nullptr),
+      _pOrder(-1), _npts(0), _pts(nullptr)
   {
   }
   virtual void updateParent(
     GModel *gm); // NEVER ever use this ! (except for reading msh files !)
 public:
   MSubPoint(MVertex *v0, int num = 0, int part = 0, bool owner = false,
-            MElement *orig = NULL)
-    : MPoint(v0, num, part), _owner(owner), _orig(orig), _base(0), _pOrder(-1),
-      _npts(0), _pts(0)
+            MElement *orig = nullptr)
+    : MPoint(v0, num, part), _owner(owner), _orig(orig), _base(nullptr),
+      _pOrder(-1), _npts(0), _pts(nullptr)
   {
   }
   MSubPoint(const std::vector<MVertex *> &v, int num = 0, int part = 0,
-            bool owner = false, MElement *orig = NULL)
-    : MPoint(v, num, part), _owner(owner), _orig(orig), _base(0), _pOrder(-1),
-      _npts(0), _pts(0)
+            bool owner = false, MElement *orig = nullptr)
+    : MPoint(v, num, part), _owner(owner), _orig(orig), _base(nullptr),
+      _pOrder(-1), _npts(0), _pts(nullptr)
   {
   }
-  MSubPoint(const MPoint &pt, bool owner = false, MElement *orig = NULL)
-    : MPoint(pt), _owner(owner), _orig(orig), _base(0), _pOrder(-1), _npts(0),
-      _pts(0)
+  MSubPoint(const MPoint &pt, bool owner = false, MElement *orig = nullptr)
+    : MPoint(pt), _owner(owner), _orig(orig), _base(nullptr), _pOrder(-1),
+      _npts(0), _pts(nullptr)
   {
   }
   ~MSubPoint();
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 0fa0009be0..65fdd22271 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -68,7 +68,7 @@ public:
     if(_v[0] != v1 && _v[0] != v2) return _v[0];
     if(_v[1] != v1 && _v[1] != v2) return _v[1];
     if(_v[2] != v1 && _v[2] != v2) return _v[2];
-    return 0;
+    return nullptr;
   }
   virtual int getNumEdges() const { return 3; }
   virtual MEdge getEdge(int num) const
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 2486f46848..661e48fb3f 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -43,7 +43,8 @@ protected:
   GEntity *_ge;
 
 public:
-  MVertex(double x, double y, double z, GEntity *ge = 0, std::size_t num = 0);
+  MVertex(double x, double y, double z, GEntity *ge = nullptr,
+          std::size_t num = 0);
   virtual ~MVertex() {}
   void deleteLast();
 
@@ -134,7 +135,7 @@ public:
 
   MEdgeVertex(double x, double y, double z, GEntity *ge, double u,
               std::size_t num = 0, double lc = -1.0)
-    : MVertex(x, y, z, ge, num), _u(u), _lc(lc), bl_data(0)
+    : MVertex(x, y, z, ge, num), _u(u), _lc(lc), bl_data(nullptr)
   {
   }
   virtual ~MEdgeVertex()
@@ -165,7 +166,7 @@ public:
 
   MFaceVertex(double x, double y, double z, GEntity *ge, double u, double v,
               std::size_t num = 0)
-    : MVertex(x, y, z, ge, num), _u(u), _v(v), bl_data(0)
+    : MVertex(x, y, z, ge, num), _u(u), _v(v), bl_data(nullptr)
   {
   }
   virtual ~MFaceVertex()
diff --git a/Geo/MVertexBoundaryLayerData.cpp b/Geo/MVertexBoundaryLayerData.cpp
index 0592187c6c..262cab084f 100644
--- a/Geo/MVertexBoundaryLayerData.cpp
+++ b/Geo/MVertexBoundaryLayerData.cpp
@@ -9,7 +9,7 @@ std::vector<MVertex *> *MVertexBoundaryLayerData::getChildren(int i)
 {
   if(i < (int)this->children.size() && i >= 0) { return &(children[i]); }
   else {
-    return 0;
+    return nullptr;
   }
 }
 
diff --git a/Geo/MVertexRTree.h b/Geo/MVertexRTree.h
index 2ac23facda..730d7259ef 100644
--- a/Geo/MVertexRTree.h
+++ b/Geo/MVertexRTree.h
@@ -36,14 +36,14 @@ public:
     delete _rtree;
   }
   MVertex *insert(MVertex *v, bool warnIfExists = false,
-                  std::set<MVertex *, MVertexPtrLessThan> *duplicates = 0)
+                  std::set<MVertex *, MVertexPtrLessThan> *duplicates = nullptr)
   {
     MVertex *out;
     double _min[3] = {v->x() - _tol, v->y() - _tol, v->z() - _tol};
     double _max[3] = {v->x() + _tol, v->y() + _tol, v->z() + _tol};
     if(!_rtree->Search(_min, _max, rtree_callback, &out)) {
       _rtree->Insert(_min, _max, v);
-      return 0;
+      return nullptr;
     }
     else {
       if(duplicates) {
@@ -60,7 +60,7 @@ public:
     }
   }
   int insert(std::vector<MVertex *> &v, bool warnIfExists = false,
-             std::set<MVertex *, MVertexPtrLessThan> *duplicates = 0)
+             std::set<MVertex *, MVertexPtrLessThan> *duplicates = nullptr)
   {
     int num = 0;
     for(std::size_t i = 0; i < v.size(); i++)
@@ -73,7 +73,7 @@ public:
     double _min[3] = {x - _tol, y - _tol, z - _tol};
     double _max[3] = {x + _tol, y + _tol, z + _tol};
     if(_rtree->Search(_min, _max, rtree_callback, &out)) return out;
-    return 0;
+    return nullptr;
   }
   std::size_t size() { return _rtree->Count(); }
 };
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index acd4e0d111..b16da5cc4c 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -37,7 +37,7 @@
 #include <BOPTools_AlgoTools.hxx>
 
 OCCEdge::OCCEdge(GModel *m, TopoDS_Edge c, int num, GVertex *v1, GVertex *v2)
-  : GEdge(m, num, v1, v2), _c(c), _trimmed(0)
+  : GEdge(m, num, v1, v2), _c(c), _trimmed(nullptr)
 {
   // force orientation of internal/external edges: otherwise reverse will not
   // produce the expected result
@@ -100,7 +100,7 @@ void OCCEdge::setTrimmed(OCCFace *f)
     _trimmed = f;
     const TopoDS_Face *s = (TopoDS_Face *)_trimmed->getNativePtr();
     _curve2d = BRep_Tool::CurveOnSurface(_c, *s, _s0, _s1);
-    if(_curve2d.IsNull()) _trimmed = 0;
+    if(_curve2d.IsNull()) _trimmed = nullptr;
   }
 }
 
diff --git a/Geo/OCCEdge.h b/Geo/OCCEdge.h
index 149a345865..136400e952 100644
--- a/Geo/OCCEdge.h
+++ b/Geo/OCCEdge.h
@@ -35,7 +35,7 @@ public:
   virtual ~OCCEdge();
   virtual SBoundingBox3d bounds(bool fast = false);
   virtual Range<double> parBounds(int i) const;
-  virtual Range<double> parBoundsOnFace(GFace *face = NULL) const;
+  virtual Range<double> parBoundsOnFace(GFace *face = nullptr) const;
   virtual GeomType geomType() const;
   virtual bool degenerate(int) const { return BRep_Tool::Degenerated(_c); }
   virtual GPoint point(double p) const;
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 500e8ca54a..1511603aed 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -71,7 +71,7 @@ void OCCFace::_setup()
     std::vector<GEdge *> l_wire;
     for(exp3.Init(wire, TopAbs_EDGE); exp3.More(); exp3.Next()) {
       TopoDS_Edge edge = TopoDS::Edge(exp3.Current());
-      GEdge *e = 0;
+      GEdge *e = nullptr;
       if(model()->getOCCInternals())
         e = model()->getOCCInternals()->getEdgeForOCCShape(model(), edge);
       if(!e) { Msg::Error("Unknown curve in surface %d", tag()); }
@@ -113,7 +113,7 @@ void OCCFace::_setup()
   for(exp2.Init(_s.Oriented(TopAbs_FORWARD), TopAbs_VERTEX, TopAbs_EDGE);
       exp2.More(); exp2.Next()) {
     TopoDS_Vertex vertex = TopoDS::Vertex(exp2.Current());
-    GVertex *v = 0;
+    GVertex *v = nullptr;
     if(model()->getOCCInternals())
       v = model()->getOCCInternals()->getVertexForOCCShape(model(), vertex);
     if(!v) { Msg::Error("Unknown point in surface %d", tag()); }
diff --git a/Geo/OCCRegion.cpp b/Geo/OCCRegion.cpp
index 233cdd7e6d..5a31b6708a 100644
--- a/Geo/OCCRegion.cpp
+++ b/Geo/OCCRegion.cpp
@@ -45,7 +45,7 @@ void OCCRegion::_setup()
     Msg::Debug("OCC volume %d - new shell", tag());
     for(exp3.Init(shell, TopAbs_FACE); exp3.More(); exp3.Next()) {
       TopoDS_Face face = TopoDS::Face(exp3.Current());
-      GFace *f = 0;
+      GFace *f = nullptr;
       if(model()->getOCCInternals())
         f = model()->getOCCInternals()->getFaceForOCCShape(model(), face);
       if(!f) { Msg::Error("Unknown surface in volume %d", tag()); }
@@ -62,7 +62,7 @@ void OCCRegion::_setup()
 
   for(exp3.Init(_s, TopAbs_EDGE, TopAbs_FACE); exp3.More(); exp3.Next()) {
     TopoDS_Edge edge = TopoDS::Edge(exp3.Current());
-    GEdge *e = 0;
+    GEdge *e = nullptr;
     if(model()->getOCCInternals())
       e = model()->getOCCInternals()->getEdgeForOCCShape(model(), edge);
     if(!e) { Msg::Error("Unknown curve in volume %d", tag()); }
@@ -76,7 +76,7 @@ void OCCRegion::_setup()
 
   for(exp3.Init(_s, TopAbs_VERTEX, TopAbs_FACE); exp3.More(); exp3.Next()) {
     TopoDS_Vertex vertex = TopoDS::Vertex(exp3.Current());
-    GVertex *v = 0;
+    GVertex *v = nullptr;
     if(model()->getOCCInternals())
       v = model()->getOCCInternals()->getVertexForOCCShape(model(), vertex);
     if(!v) { Msg::Error("Unknown point in volume %d", tag()); }
diff --git a/Geo/SOrientedBoundingBox.cpp b/Geo/SOrientedBoundingBox.cpp
index 5a0003d7b1..f7db28a88e 100644
--- a/Geo/SOrientedBoundingBox.cpp
+++ b/Geo/SOrientedBoundingBox.cpp
@@ -292,7 +292,7 @@ SOrientedBoundingBox::buildOBB(std::vector<SPoint3> &vertices)
   // Find the convex hull from a delaunay triangulation of the points
   DocRecord record(points.size());
   record.numPoints = points.size();
-  srand((unsigned)time(0));
+  srand((unsigned)time(nullptr));
   for(std::size_t i = 0; i < points.size(); i++) {
     record.points[i].where.h =
       points[i]->x() + (10e-6) * sizes[smallest_comp == 0 ? 1 : 0] *
@@ -300,7 +300,7 @@ SOrientedBoundingBox::buildOBB(std::vector<SPoint3> &vertices)
     record.points[i].where.v =
       points[i]->y() + (10e-6) * sizes[smallest_comp == 2 ? 1 : 0] *
                          (-0.5 + ((double)rand()) / RAND_MAX);
-    record.points[i].adjacent = NULL;
+    record.points[i].adjacent = nullptr;
   }
 
   try {
diff --git a/Geo/boundaryLayersData.cpp b/Geo/boundaryLayersData.cpp
index f19b807325..252dee36d0 100644
--- a/Geo/boundaryLayersData.cpp
+++ b/Geo/boundaryLayersData.cpp
@@ -425,7 +425,7 @@ bool buildAdditionalPoints2D(GFace *gf)
   for(int i = 0; i < nBL; ++i) {
     // GET THE FIELD THAT DEFINES THE DISTANCE FUNCTION
     Field *bl_field = fields->get(fields->getBoundaryLayerField(i));
-    if(bl_field == NULL) continue;
+    if(bl_field == nullptr) continue;
     BoundaryLayerField *blf = dynamic_cast<BoundaryLayerField *>(bl_field);
 
     if(!blf->setupFor2d(gf->tag())) continue;
diff --git a/Geo/boundaryLayersData.h b/Geo/boundaryLayersData.h
index 25458dad46..ab7380eb39 100644
--- a/Geo/boundaryLayersData.h
+++ b/Geo/boundaryLayersData.h
@@ -94,13 +94,13 @@ public:
   }
   inline const BoundaryLayerFan *getFan(MVertex *v) const
   {
-    std::map<MVertex *, BoundaryLayerFan>::const_iterator it = _fans.find(v);
+    auto it = _fans.find(v);
     if(it != _fans.end()) { return &it->second; }
-    return 0;
+    return nullptr;
   }
   inline const BoundaryLayerData &getColumn(MVertex *v, MEdge e) const
   {
-    std::map<MVertex *, BoundaryLayerFan>::const_iterator it = _fans.find(v);
+    auto it = _fans.find(v);
     int N = getNbColumns(v);
     if(N == 1) return getColumn(v, 0);
     MEdgeEqual aaa;
@@ -119,9 +119,7 @@ public:
   inline const BoundaryLayerData &getColumn(MVertex *v, int iColumn) const
   {
     int count = 0;
-    for(std::multimap<MVertex *, BoundaryLayerData>::const_iterator itm =
-          _data.lower_bound(v);
-        itm != _data.upper_bound(v); ++itm) {
+    for(auto itm = _data.lower_bound(v); itm != _data.upper_bound(v); ++itm) {
       if(count++ == iColumn) return itm->second;
     }
     static BoundaryLayerData error;
diff --git a/Geo/closestVertex.cpp b/Geo/closestVertex.cpp
index 438257ec68..69d1f8cf8d 100644
--- a/Geo/closestVertex.cpp
+++ b/Geo/closestVertex.cpp
@@ -55,7 +55,7 @@ MVertex *closestVertexFinder ::operator()(const SPoint3 &p)
 {
 #if defined(HAVE_ANN)
 
-  if(nbVtcs == 0) return NULL;
+  if(nbVtcs == 0) return nullptr;
 
   double xyz[3] = {p.x(), p.y(), p.z()};
   kdtree->annkSearch(xyz, 1, index, dist);
@@ -70,7 +70,7 @@ MVertex *closestVertexFinder ::operator()(const SPoint3 &p,
 {
 #if defined(HAVE_ANN)
 
-  if(nbVtcs == 0) return NULL;
+  if(nbVtcs == 0) return nullptr;
 
   double ori[4] = {p.x(), p.y(), p.z(), 1};
   double xyz[4] = {0, 0, 0, 0};
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 32a0548bad..7859f868a8 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -21,7 +21,7 @@ discreteEdge::discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1)
   : GEdge(model, num, _v0, _v1)
 {
   bool ok;
-  Curve *c = CreateCurve(num, MSH_SEGM_DISCRETE, 0, 0, 0, -1, -1, 0., 1., ok);
+  Curve *c = CreateCurve(num, MSH_SEGM_DISCRETE, 0, nullptr, nullptr, -1, -1, 0., 1., ok);
   Tree_Add(model->getGEOInternals()->Curves, &c);
   CreateReversedCurve(c);
 }
@@ -29,7 +29,7 @@ discreteEdge::discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1)
 discreteEdge::discreteEdge(GModel *model, int num) : GEdge(model, num)
 {
   bool ok;
-  Curve *c = CreateCurve(num, MSH_SEGM_DISCRETE, 0, 0, 0, -1, -1, 0., 1., ok);
+  Curve *c = CreateCurve(num, MSH_SEGM_DISCRETE, 0, nullptr, nullptr, -1, -1, 0., 1., ok);
   Tree_Add(model->getGEOInternals()->Curves, &c);
   CreateReversedCurve(c);
 }
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 3e60dacf37..d3b2ecfb6b 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -166,7 +166,7 @@ public:
   double _distance;
   SPoint3 _closestPoint;
   MTriangle *_t3d, *_t2d;
-  dfWrapper(const SPoint3 &p) : _p(p), _distance(1.e22), _t3d(NULL), _t2d(NULL)
+  dfWrapper(const SPoint3 &p) : _p(p), _distance(1.e22), _t3d(nullptr), _t2d(nullptr)
   {
   }
 };
@@ -566,7 +566,7 @@ GPoint discreteFace::intersectionWithCircle(const SVector3 &n1,
   if(_param.empty()) return 0.;
 
   MTriangle *t2d = (MTriangle *)_param.oct->find(uv[0], uv[1], 0.0, -1, true);
-  MTriangle *t3d = NULL;
+  MTriangle *t3d = nullptr;
   if(t2d) {
     int position = (int)(t2d - &_param.t2d[0]);
     t3d = &_param.t3d[position];
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index 00fd8d0b4d..240ab11aa7 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -26,7 +26,7 @@ private:
     std::vector<MTriangle> t3d;
     std::vector<SVector3> CURV;
     double umin, umax, vmin, vmax;
-    param() : oct(NULL), umin(-1), umax(1), vmin(-1), vmax(1) {}
+    param() : oct(nullptr), umin(-1), umax(1), vmin(-1), vmax(1) {}
     ~param();
     bool empty() const { return t2d.empty(); }
     void clear();
@@ -47,7 +47,7 @@ public:
   Range<double> parBounds(int i) const;
   bool containsParam(const SPoint2 &pt);
   GPoint closestPoint(const SPoint3 &queryPoint, double maxDistance,
-                      SVector3 *normal = NULL) const;
+                      SVector3 *normal = nullptr) const;
   GPoint closestPoint(const SPoint3 &queryPoint,
                       const double initialGuess[2]) const;
   SVector3 normal(const SPoint2 &param) const;
diff --git a/Geo/ghostEdge.h b/Geo/ghostEdge.h
index 031e7771c0..1e142bb5e8 100644
--- a/Geo/ghostEdge.h
+++ b/Geo/ghostEdge.h
@@ -21,7 +21,7 @@ private:
 
 public:
   ghostEdge(GModel *model, const int num, const int partition)
-    : discreteEdge(model, num, NULL, NULL), _partition(partition),
+    : discreteEdge(model, num, nullptr, nullptr), _partition(partition),
       _ghostCells(), _saveMesh(false), _haveMesh(false)
   {
   }
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 3243965e3d..ced2ab6531 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -78,7 +78,7 @@ void gmshFace::resetNativePtr(Surface *s)
   std::vector<GEdge *> l_wire;
   l_wire.reserve(eds.size());
 
-  GVertex *first = 0;
+  GVertex *first = nullptr;
   for(std::size_t i = 0; i < eds.size(); i++) {
     GEdge *e = eds[i];
     int num = nums[i];
@@ -89,7 +89,7 @@ void gmshFace::resetNativePtr(Surface *s)
     if(next == first) {
       edgeLoops.push_back(GEdgeLoop(l_wire));
       l_wire.clear();
-      first = 0;
+      first = nullptr;
     }
     l_edges.push_back(e);
     e->addFace(this);
@@ -283,7 +283,7 @@ GPoint gmshFace::closestPoint(const SPoint3 &qp,
   v.Pos.Z = qp.z();
   double u[2] = {initialGuess[0], initialGuess[1]};
   bool result = ProjectPointOnSurface(_s, v, u);
-  if(!result) return GPoint(-1.e22, -1.e22, -1.e22, 0, u);
+  if(!result) return GPoint(-1.e22, -1.e22, -1.e22, nullptr, u);
   return GPoint(v.Pos.X, v.Pos.Y, v.Pos.Z, this, u);
 }
 
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index ae34337bd3..2a3f2ff708 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -37,7 +37,7 @@ gLevelset *gLevelset::find(int tag)
 {
   gLevelset l(tag);
   auto it = all_.find(&l);
-  if(it == all_.end()) return 0;
+  if(it == all_.end()) return nullptr;
   return *it;
 }
 
@@ -1114,7 +1114,7 @@ gLevelsetPostView::gLevelsetPostView(int index, int tag)
   }
   else {
     Msg::Error("Unknown View[%d] in PostView levelset", _viewIndex);
-    _octree = 0;
+    _octree = nullptr;
   }
 }
 
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index 295852c3d4..e52776569b 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -84,7 +84,7 @@ public:
   virtual ~gLevelset() {}
   static gLevelset *find(int tag);
   static void add(gLevelset *l);
-  virtual gLevelset *clone() const { return 0; }
+  virtual gLevelset *clone() const { return nullptr; }
   virtual double operator()(double x, double y, double z) const { return 0.; }
   bool isInsideDomain(const double &x, const double &y, const double &z) const
   {
diff --git a/Geo/gmshSurface.cpp b/Geo/gmshSurface.cpp
index 68939c5c3b..c731aa05ef 100644
--- a/Geo/gmshSurface.cpp
+++ b/Geo/gmshSurface.cpp
@@ -53,7 +53,7 @@ gmshSurface *gmshSurface::getSurface(int iSurface)
   auto it = allGmshSurfaces.find(iSurface);
   if(it == allGmshSurfaces.end()) {
     Msg::Error("gmshSurface %d does not exist", iSurface);
-    return 0;
+    return nullptr;
   }
   return it->second;
 }
@@ -120,7 +120,7 @@ gmshParametricSurface::gmshParametricSurface(char *valX, char *valY, char *valZ)
   _f = new mathEvaluator(expressions, variables);
   if(expressions.empty()) {
     delete _f;
-    _f = 0;
+    _f = nullptr;
   }
 }
 
diff --git a/Geo/gmshSurface.h b/Geo/gmshSurface.h
index d11369838f..e8eb59fefe 100644
--- a/Geo/gmshSurface.h
+++ b/Geo/gmshSurface.h
@@ -26,7 +26,7 @@ public:
   virtual ~gmshSurface() {}
   static void reset()
   {
-    std::map<int, gmshSurface *>::iterator it = allGmshSurfaces.begin();
+    auto it = allGmshSurfaces.begin();
     for(; it != allGmshSurfaces.end(); ++it) {
       if(!it->second->vertex_defined_on_surface) delete it->second;
     }
diff --git a/Geo/partitionEdge.h b/Geo/partitionEdge.h
index 7af06b2312..45e67e5c49 100644
--- a/Geo/partitionEdge.h
+++ b/Geo/partitionEdge.h
@@ -18,15 +18,15 @@ public:
   partitionEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1,
                 const std::vector<int> &partitions)
     : discreteEdge(model, num, _v0, _v1), _partitions(partitions),
-      _parentEntity(NULL)
+      _parentEntity(nullptr)
   {
   }
   partitionEdge(GModel *model, int num, const std::vector<int> &partitions)
-    : discreteEdge(model, num), _partitions(partitions), _parentEntity(NULL)
+    : discreteEdge(model, num), _partitions(partitions), _parentEntity(nullptr)
   {
   }
   partitionEdge(GModel *model, const std::vector<int> &partitions)
-    : discreteEdge(model), _partitions(partitions), _parentEntity(NULL)
+    : discreteEdge(model), _partitions(partitions), _parentEntity(nullptr)
   {
   }
   virtual ~partitionEdge() {}
diff --git a/Geo/partitionFace.h b/Geo/partitionFace.h
index 3418722144..aadf3b7c4b 100644
--- a/Geo/partitionFace.h
+++ b/Geo/partitionFace.h
@@ -16,15 +16,15 @@ private:
 
 public:
   partitionFace(GModel *model, int num, const std::vector<int> &partitions)
-    : discreteFace(model, num), _partitions(partitions), _parentEntity(NULL)
+    : discreteFace(model, num), _partitions(partitions), _parentEntity(nullptr)
   {
   }
   partitionFace(GModel *model, int num)
-    : discreteFace(model, num), _partitions(), _parentEntity(NULL)
+    : discreteFace(model, num), _partitions(), _parentEntity(nullptr)
   {
   }
   partitionFace(GModel *model, const std::vector<int> &partitions)
-    : discreteFace(model), _partitions(partitions), _parentEntity(NULL)
+    : discreteFace(model), _partitions(partitions), _parentEntity(nullptr)
   {
   }
   virtual ~partitionFace() {}
diff --git a/Geo/partitionRegion.h b/Geo/partitionRegion.h
index 45b39fe41b..4effacb3f6 100644
--- a/Geo/partitionRegion.h
+++ b/Geo/partitionRegion.h
@@ -16,15 +16,16 @@ private:
 
 public:
   partitionRegion(GModel *model, int num, const std::vector<int> &partitions)
-    : discreteRegion(model, num), _partitions(partitions), _parentEntity(NULL)
+    : discreteRegion(model, num), _partitions(partitions),
+      _parentEntity(nullptr)
   {
   }
   partitionRegion(GModel *model, int num)
-    : discreteRegion(model, num), _partitions(), _parentEntity(NULL)
+    : discreteRegion(model, num), _partitions(), _parentEntity(nullptr)
   {
   }
   partitionRegion(GModel *model, const std::vector<int> &partitions)
-    : discreteRegion(model), _partitions(partitions), _parentEntity(NULL)
+    : discreteRegion(model), _partitions(partitions), _parentEntity(nullptr)
   {
   }
   virtual ~partitionRegion() {}
diff --git a/Geo/partitionVertex.h b/Geo/partitionVertex.h
index 6f697c0141..11b5629553 100644
--- a/Geo/partitionVertex.h
+++ b/Geo/partitionVertex.h
@@ -16,15 +16,16 @@ private:
 
 public:
   partitionVertex(GModel *model, int num, const std::vector<int> &partitions)
-    : discreteVertex(model, num), _partitions(partitions), _parentEntity(NULL)
+    : discreteVertex(model, num), _partitions(partitions),
+      _parentEntity(nullptr)
   {
   }
   partitionVertex(GModel *model, int num)
-    : discreteVertex(model, num), _partitions(), _parentEntity(NULL)
+    : discreteVertex(model, num), _partitions(), _parentEntity(nullptr)
   {
   }
   partitionVertex(GModel *model, const std::vector<int> &partitions)
-    : discreteVertex(model), _partitions(partitions), _parentEntity(NULL)
+    : discreteVertex(model), _partitions(partitions), _parentEntity(nullptr)
   {
   }
   virtual ~partitionVertex() {}
diff --git a/Geo/xyEdge.h b/Geo/xyEdge.h
index 1a0aac0ffe..fb5bb53435 100644
--- a/Geo/xyEdge.h
+++ b/Geo/xyEdge.h
@@ -18,7 +18,7 @@ public:
   virtual SVector3 firstDer(double par) const { return SVector3(); }
   virtual SVector3 secondDer(double par) const { return SVector3(); }
   ModelType getNativeType() const { return GmshModel; }
-  void *getNativePtr() const { return NULL; }
+  void *getNativePtr() const { return nullptr; }
   virtual SPoint2 reparamOnFace(const GFace *face, double epar, int dir) const
   {
     return SPoint2();
diff --git a/Geo/xyFace.h b/Geo/xyFace.h
index 4f5006a177..7a3bf9e75e 100644
--- a/Geo/xyFace.h
+++ b/Geo/xyFace.h
@@ -49,7 +49,7 @@ public:
   }
   virtual GEntity::GeomType geomType() const { return GEntity::Plane; }
   ModelType getNativeType() const { return GmshModel; }
-  void *getNativePtr() const { return NULL; }
+  void *getNativePtr() const { return nullptr; }
   virtual SPoint2 parFromPoint(const SPoint3 &p, bool onSurface = true) const
   {
     return SPoint2(p.x(), p.y());
-- 
GitLab