diff --git a/CHANGELOG.txt b/CHANGELOG.txt
index 50b21bd19371728f5c19e863af74d1eedcfd901a..136b6eac9544ea173c489c9457836940059551e8 100644
--- a/CHANGELOG.txt
+++ b/CHANGELOG.txt
@@ -1,7 +1,9 @@
 4.9.3 (Work-in-progress): improved handling of degenerate curves in periodic
-surfaces and boundary layer extrusion; moved all kernel sources to src/
-subdirectory; renamed demos/ as examples/ and tutorial/ as tutorials/; extended
-Mesh.SaveGroupsOfElements capabilities for INP export; small bug fixes.
+surfaces and boundary layer extrusion; extended Mesh.SaveGroupsOfElements
+capabilities for INP export; extended Mesh.MeshSizeExtendFromBoundary + new
+"Extend" mesh size field to enable alternative mesh size extensions from
+boundary; moved all kernel sources to src/ subdirectory; renamed demos/ as
+examples/ and tutorial/ as tutorials/; small bug fixes.
 
 4.9.2 (December 23, 2021): faster projection on OCC entities; extended
 Mesh.SaveGroupsOfNodes capabilities for INP export; improved transfinite meshing
diff --git a/contrib/MathEx/mathex.cpp b/contrib/MathEx/mathex.cpp
index 81a609d825b6c6bbf54baad6a77457f5aa85bc93..76522c355343f141911cafb27adf8c68e6f07a0d 100644
--- a/contrib/MathEx/mathex.cpp
+++ b/contrib/MathEx/mathex.cpp
@@ -182,7 +182,7 @@
          // in future, put more precisery checking, for overflow ?
          {
             if (y == 0)
-               throw mathex::error("Error [binary_divide()]: divisin by zero");
+               throw mathex::error("Error [binary_divide()]: division by zero");
             else
                return (x/y);
          } // binary_divide()
diff --git a/doc/texinfo/gmsh.texi b/doc/texinfo/gmsh.texi
index 981778d01ce88b2673e32434cc4f7f4cb85388e9..52b17375bf6b36cc88c7c23524bb336479b34b2d 100644
--- a/doc/texinfo/gmsh.texi
+++ b/doc/texinfo/gmsh.texi
@@ -2882,10 +2882,9 @@ the smallest element size is selected at any given point.  Using the
 Gmsh API, the resulting value can be further modified using a C++, C,
 Python or Julia mesh size callback function provided via
 @code{gmsh/model/mesh/setSizeCallback} (@pxref{Namespace
-gmsh/model/mesh}). In addition, boundary mesh sizes (on curves or
-surfaces) are interpolated inside the enclosed entity (surface or
-volume, respectively) if the option
-@code{Mesh.MeshSizeExtendFromBoundary} is set (it is by default).
+gmsh/model/mesh}). In addition, boundary mesh sizes are interpolated
+inside surfaces and/or volumes depending on the value of
+@code{Mesh.MeshSizeExtendFromBoundary}.
 
 All element sizes are further constrained in the interval [
 @code{Mesh.MeshSizeMin}, @code{Mesh.MeshSizeMax} ] (which can also be
diff --git a/doc/texinfo/opt_fields.texi b/doc/texinfo/opt_fields.texi
index d6580e73cb1c009679bd46bdca5e188b6856a6ba..4100c7e6925923ed6517d4d19f174da6914bea0f 100644
--- a/doc/texinfo/opt_fields.texi
+++ b/doc/texinfo/opt_fields.texi
@@ -356,6 +356,37 @@ type: list@*
 default value: @code{@{@}}
 @end table
 
+@item Extend
+Compute an extension of the mesh sizes from the given boundary curves (resp. surfaces) inside the surfaces (resp. volumes) being meshed. If the mesh size on the boundary, computed as the local average length of the edges of the boundary elements, is denoted by SizeBnd, the extension is computed as:@*
+@*
+  F = f * SizeBnd + (1 - f) * SizeMax, if d < DistMax@*
+@*
+  F = SizeMax, if d >= DistMax@*
+where d denotes the distance to the boundary and f = ((DistMax - d) / DistMax)^Power.@*
+Options:@*
+@table @code
+@item CurvesList
+Tags of curves in the geometric model@*
+type: list@*
+default value: @code{@{@}}
+@item DistMax
+Maximum distance of the size extension@*
+type: float@*
+default value: @code{1}
+@item Power
+Power exponent used to interpolate the mesh size@*
+type: float@*
+default value: @code{1}
+@item SizeMax
+Mesh size outside DistMax@*
+type: float@*
+default value: @code{1}
+@item SurfacesList
+Tags of surfaces in the geometric model@*
+type: list@*
+default value: @code{@{@}}
+@end table
+
 @item ExternalProcess
 **This Field is experimental**@*
 Call an external process that received coordinates triple (x,y,z) as binary double precision numbers on stdin and is supposed to write the field value on stdout as a binary double precision number.@*
@@ -800,15 +831,15 @@ Return F = SizeMin if Field[InField] <= DistMin, F = SizeMax if Field[InField] >
 Options:@*
 @table @code
 @item DistMax
-Distance from entity after which the mesh size will be SizeMax@*
+Value after which the mesh size will be SizeMax@*
 type: float@*
 default value: @code{10}
 @item DistMin
-Distance from entity up to which the mesh size will be SizeMin@*
+Value up to which the mesh size will be SizeMin@*
 type: float@*
 default value: @code{1}
 @item InField
-Tag of the field to evaluate@*
+Tag of the field computing the input value, usually a distance@*
 type: integer@*
 default value: @code{0}
 @item Sigmoid
@@ -816,11 +847,11 @@ True to interpolate between SizeMin and LcMax using a sigmoid, false to interpol
 type: boolean@*
 default value: @code{0}
 @item SizeMax
-Mesh size outside DistMax@*
+Mesh size when value > DistMax@*
 type: float@*
 default value: @code{1}
 @item SizeMin
-Mesh size inside DistMin@*
+Mesh size when value < DistMin@*
 type: float@*
 default value: @code{0.1}
 @item StopAtDistMax
diff --git a/doc/texinfo/opt_mesh.texi b/doc/texinfo/opt_mesh.texi
index 5ed24b958800f20622884faa6c31f39d31e9cabe..25ca5c6525f5e25305d0c80ef40469cc6cdffe78 100644
--- a/doc/texinfo/opt_mesh.texi
+++ b/doc/texinfo/opt_mesh.texi
@@ -327,7 +327,7 @@ Default value: @code{0}@*
 Saved in: @code{General.OptionsFileName}
 
 @item Mesh.MeshSizeExtendFromBoundary
-Extend computation of mesh element sizes from the boundaries into the interior (for 3D Delaunay, use 1: longest or 2: shortest surface edge length)@*
+Extend computation of mesh element sizes from the boundaries into the interior (0: never; 1: for surfaces and volumes; 2: for surfaces and volumes, but use smallest surface element edge length instead of longest length in 3D Delaunay; -2: only for surfaces; -3: only for volumes)@*
 Default value: @code{1}@*
 Saved in: @code{General.OptionsFileName}
 
diff --git a/examples/boolean/extend_field.geo b/examples/boolean/extend_field.geo
new file mode 100644
index 0000000000000000000000000000000000000000..80a3fe627dd61efcc7df887ebd6e9886b35d248e
--- /dev/null
+++ b/examples/boolean/extend_field.geo
@@ -0,0 +1,38 @@
+SetFactory("OpenCASCADE");
+Box(1) = {0,0,0, 1,1,0.5};
+Sphere(2) = {1,1,0.5,0.4};
+BooleanDifference(3) = { Volume{1}; Delete; }{ Volume{2}; Delete; };
+
+Box(10) = {0.3,0.3,0.5, 0.1,0.1,0.05};
+Box(11) = {0.5,0.5,0.5, 0.1,0.1,0.05};
+BooleanFragments{ Volume{10, 11}; Delete; }{ Volume{3}; Delete; }
+
+MeshSize{ : } = 0.04;
+//MeshSize{ : } = 0.002;
+MeshSize{ PointsOf{ Volume{10, 11}; } } = 0.002;
+Mesh.Algorithm3D = 10;
+
+DefineConstant[
+  use_field = {1, Choices{0, 1},
+    Name "Parameters/0Use Extend field?"}
+  size_max = {0.04, Min 0.002, Max 0.2, Step 0.001,
+    Name "Parameters/Mesh size after Extend", Visible use_field}
+  dist_max = {0.1, Min 0.01, Max 1, Step 0.01,
+    Name "Parameters/Extend distance", Visible use_field}
+  power = {2, Min 0.1, Max 10, Step 0.1,
+    Name "Parameters/Extend power", Visible use_field}
+];
+
+If(use_field == 0)
+  Mesh.MeshSizeExtendFromBoundary = 1;
+Else
+  // set to -2 or -3 to only extend in 2D or 3D
+  Mesh.MeshSizeExtendFromBoundary = 0;
+  Field[1] = Extend;
+  Field[1].SurfacesList = {Surface{:}};
+  Field[1].CurvesList = {Curve{:}};
+  Field[1].DistMax = dist_max;
+  Field[1].SizeMax = size_max;
+  Field[1].Power = power;
+  Background Field = 1;
+EndIf
diff --git a/src/common/DefaultOptions.h b/src/common/DefaultOptions.h
index 46c9a96413955c0fdb9825b424fcfdc05c39de79..30b443e7e9ec2eefaaaeaa363ad93d57f21b6248 100644
--- a/src/common/DefaultOptions.h
+++ b/src/common/DefaultOptions.h
@@ -1246,7 +1246,9 @@ StringXNumber MeshOptions_Number[] = {
     "Mesh only entities that have no existing mesh" },
   { F|O, "MeshSizeExtendFromBoundary" , opt_mesh_lc_extend_from_boundary, 1. ,
     "Extend computation of mesh element sizes from the boundaries into the interior "
-    "(for 3D Delaunay, use 1: longest or 2: shortest surface edge length)"},
+    "(0: never; 1: for surfaces and volumes; 2: for surfaces and volumes, but use "
+    "smallest surface element edge length instead of longest length in 3D Delaunay; "
+    "-2: only for surfaces; -3: only for volumes)"},
   { F|O, "MeshSizeFactor" , opt_mesh_lc_factor , 1.0 ,
     "Factor applied to all mesh element sizes" },
   { F|O, "MeshSizeMin" , opt_mesh_lc_min, 0.0 ,
diff --git a/src/mesh/BackgroundMeshTools.cpp b/src/mesh/BackgroundMeshTools.cpp
index 1d19d77f1971a7f6c72df105c3b3c42ab5f0461b..3f51f05117aa2a95d2a64c949117423f8d3a550b 100644
--- a/src/mesh/BackgroundMeshTools.cpp
+++ b/src/mesh/BackgroundMeshTools.cpp
@@ -338,12 +338,14 @@ SMetric3 BGM_MeshMetric(GEntity *ge, double U, double V, double X, double Y,
 
 bool Extend1dMeshIn2dSurfaces(GFace *gf)
 {
-  return gf->getMeshSizeFromBoundary();
+  int val = gf->getMeshSizeFromBoundary();
+  return (val > 0 || val == -2);
 }
 
 bool Extend2dMeshIn3dVolumes()
 {
-  return CTX::instance()->mesh.lcExtendFromBoundary;
+  int val = CTX::instance()->mesh.lcExtendFromBoundary;
+  return (val > 0 || val == -3);
 }
 
 SMetric3 max_edge_curvature_metric(const GVertex *gv)
diff --git a/src/mesh/Field.cpp b/src/mesh/Field.cpp
index 04cd871a576ef85f9aeebe372621ab310d8bf50d..a5ca502ce38d7050a99d52115ad68a773300ce0a 100644
--- a/src/mesh/Field.cpp
+++ b/src/mesh/Field.cpp
@@ -633,15 +633,16 @@ public:
     _stopAtDistMax = false;
 
     options["InField"] =
-      new FieldOptionInt(_inField, "Tag of the field to evaluate");
+      new FieldOptionInt(_inField, "Tag of the field computing the input "
+                         "value, usually a distance");
     options["DistMin"] = new FieldOptionDouble(
-      _dMin, "Distance from entity up to which the mesh size will be SizeMin");
+      _dMin, "Value up to which the mesh size will be SizeMin");
     options["DistMax"] = new FieldOptionDouble(
-      _dMax, "Distance from entity after which the mesh size will be SizeMax");
+      _dMax, "Value after which the mesh size will be SizeMax");
     options["SizeMin"] =
-      new FieldOptionDouble(_lcMin, "Mesh size inside DistMin");
+      new FieldOptionDouble(_lcMin, "Mesh size when value < DistMin");
     options["SizeMax"] =
-      new FieldOptionDouble(_lcMax, "Mesh size outside DistMax");
+      new FieldOptionDouble(_lcMax, "Mesh size when value > DistMax");
     options["Sigmoid"] = new FieldOptionBool(
       _sigmoid,
       "True to interpolate between SizeMin and LcMax using a sigmoid, "
@@ -2718,83 +2719,231 @@ public:
   }
   void update()
   {
-    _infos.clear();
-    _pc.pts.clear();
-    if(_kdtree) delete _kdtree;
+    if(updateNeeded) {
+      _infos.clear();
+      _pc.pts.clear();
+      if(_kdtree) delete _kdtree;
 
-    for(auto it = _pointTags.begin(); it != _pointTags.end(); ++it) {
-      GVertex *gv = GModel::current()->getVertexByTag(*it);
-      if(gv) {
-        _pc.pts.push_back(SPoint3(gv->x(), gv->y(), gv->z()));
-        _infos.push_back(AttractorInfo(*it, 0, 0, 0));
-      }
-      else {
-        Msg::Warning("Unknown point %d", *it);
+      for(auto it = _pointTags.begin(); it != _pointTags.end(); ++it) {
+        GVertex *gv = GModel::current()->getVertexByTag(*it);
+        if(gv) {
+          _pc.pts.push_back(SPoint3(gv->x(), gv->y(), gv->z()));
+          _infos.push_back(AttractorInfo(*it, 0, 0, 0));
+        }
+        else {
+          Msg::Warning("Unknown point %d", *it);
+        }
       }
-    }
 
-    for(auto it = _curveTags.begin(); it != _curveTags.end(); ++it) {
-      GEdge *e = GModel::current()->getEdgeByTag(*it);
-      if(e) {
-        if(e->mesh_vertices.size()) {
-          for(std::size_t i = 0; i < e->mesh_vertices.size(); i++) {
-            _pc.pts.push_back(SPoint3(e->mesh_vertices[i]->x(),
-                                      e->mesh_vertices[i]->y(),
-                                      e->mesh_vertices[i]->z()));
-            double t = 0.;
-            e->mesh_vertices[i]->getParameter(0, t);
+      for(auto it = _curveTags.begin(); it != _curveTags.end(); ++it) {
+        GEdge *e = GModel::current()->getEdgeByTag(*it);
+        if(e) {
+          if(e->mesh_vertices.size()) {
+            for(std::size_t i = 0; i < e->mesh_vertices.size(); i++) {
+              _pc.pts.push_back(SPoint3(e->mesh_vertices[i]->x(),
+                                        e->mesh_vertices[i]->y(),
+                                        e->mesh_vertices[i]->z()));
+              double t = 0.;
+              e->mesh_vertices[i]->getParameter(0, t);
+              _infos.push_back(AttractorInfo(*it, 1, t, 0));
+            }
+          }
+          int NNN = _sampling - e->mesh_vertices.size();
+          for(int i = 1; i < NNN - 1; i++) {
+            double u = (double)i / (NNN - 1);
+            Range<double> b = e->parBounds(0);
+            double t = b.low() + u * (b.high() - b.low());
+            GPoint gp = e->point(t);
+            _pc.pts.push_back(SPoint3(gp.x(), gp.y(), gp.z()));
             _infos.push_back(AttractorInfo(*it, 1, t, 0));
           }
         }
-        int NNN = _sampling - e->mesh_vertices.size();
-        for(int i = 1; i < NNN - 1; i++) {
-          double u = (double)i / (NNN - 1);
-          Range<double> b = e->parBounds(0);
-          double t = b.low() + u * (b.high() - b.low());
-          GPoint gp = e->point(t);
-          _pc.pts.push_back(SPoint3(gp.x(), gp.y(), gp.z()));
-          _infos.push_back(AttractorInfo(*it, 1, t, 0));
+        else {
+          Msg::Warning("Unknown curve %d", *it);
         }
       }
-      else {
-        Msg::Warning("Unknown curve %d", *it);
+
+      for(auto it = _surfaceTags.begin(); it != _surfaceTags.end(); ++it) {
+        GFace *f = GModel::current()->getFaceByTag(*it);
+        if(f) {
+          double maxDist = f->bounds().diag() / _sampling;
+          std::vector<SPoint2> uvpoints;
+          f->fillPointCloud(maxDist, &_pc.pts, &uvpoints);
+          for(std::size_t i = 0; i < uvpoints.size(); i++)
+            _infos.push_back
+              (AttractorInfo(*it, 2, uvpoints[i].x(), uvpoints[i].y()));
+        }
+        else {
+          Msg::Warning("Unknown surface %d", *it);
+        }
       }
+
+      // construct a kd-tree index:
+      _kdtree = new SPoint3KDTree(3, _pc2kdtree,
+                                  nanoflann::KDTreeSingleIndexAdaptorParams(10));
+      _kdtree->buildIndex();
+      updateNeeded = false;
     }
+  }
+  using Field::operator();
+  virtual double operator()(double X, double Y, double Z, GEntity *ge = nullptr)
+  {
+    if(!_kdtree) return MAX_LC;
+    double pt[3] = {X, Y, Z};
+    nanoflann::KNNResultSet<double> res(1);
+    double outDistSqr;
+    res.init(&_outIndex, &outDistSqr);
+    _kdtree->findNeighbors(res, &pt[0], nanoflann::SearchParams(10));
+    return sqrt(outDistSqr);
+  }
+};
 
-    for(auto it = _surfaceTags.begin(); it != _surfaceTags.end(); ++it) {
-      GFace *f = GModel::current()->getFaceByTag(*it);
-      if(f) {
-        double maxDist = f->bounds().diag() / _sampling;
-        std::vector<SPoint2> uvpoints;
-        f->fillPointCloud(maxDist, &_pc.pts, &uvpoints);
-        for(std::size_t i = 0; i < uvpoints.size(); i++)
-          _infos.push_back
-            (AttractorInfo(*it, 2, uvpoints[i].x(), uvpoints[i].y()));
+class ExtendField : public Field {
+  std::list<int> _tagCurves, _tagSurfaces;
+  std::vector<double> _sizeCurves, _sizeSurfaces;
+  SPoint3Cloud _pcCurves, _pcSurfaces;
+  SPoint3CloudAdaptor<SPoint3Cloud> _pc2kdtreeCurves, _pc2kdtreeSurfaces;
+  SPoint3KDTree *_kdtreeCurves, *_kdtreeSurfaces;
+  double _distMax, _sizeMax, _power;
+public:
+  ExtendField()
+    : _pc2kdtreeCurves(_pcCurves), _pc2kdtreeSurfaces(_pcSurfaces),
+      _kdtreeCurves(nullptr), _kdtreeSurfaces(nullptr)
+  {
+    _distMax = 1.;
+    _sizeMax = 1.;
+    _power = 1.;
+    options["CurvesList"] = new FieldOptionList(
+      _tagCurves, "Tags of curves in the geometric model", &updateNeeded);
+    options["SurfacesList"] = new FieldOptionList(
+      _tagSurfaces, "Tags of surfaces in the geometric model", &updateNeeded);
+    options["DistMax"] = new FieldOptionDouble(
+      _distMax, "Maximum distance of the size extension");
+    options["SizeMax"] = new FieldOptionDouble(
+      _sizeMax, "Mesh size outside DistMax");
+    options["Power"] = new FieldOptionDouble(
+      _power, "Power exponent used to interpolate the mesh size");
+  }
+  ~ExtendField()
+  {
+    if(_kdtreeCurves) delete _kdtreeCurves;
+    if(_kdtreeSurfaces) delete _kdtreeSurfaces;
+  }
+  const char *getName() { return "Extend"; }
+  std::string getDescription()
+  {
+    return "Compute an extension of the mesh sizes from the given boundary "
+           "curves (resp. surfaces) inside the surfaces (resp. volumes) "
+           "being meshed. If the mesh size on the boundary, computed "
+           "as the local average length of the edges of the boundary elements, "
+           "is denoted by SizeBnd, the extension is computed as:\n\n"
+           "  F = f * SizeBnd + (1 - f) * SizeMax, if d < DistMax\n\n"
+           "  F = SizeMax, if d >= DistMax\n"
+           "where d denotes the distance to the boundary and "
+           "f = ((DistMax - d) / DistMax)^Power.";
+  }
+  void recomputeCurves()
+  {
+    _sizeCurves.clear();
+    _pcCurves.pts.clear();
+    if(_kdtreeCurves) delete _kdtreeCurves;
+    for(auto t : _tagCurves) {
+      GEdge *ge = GModel::current()->getEdgeByTag(t);
+      if(ge) {
+        for(auto e : ge->lines) {
+          _pcCurves.pts.push_back(e->barycenter());
+          _sizeCurves.push_back(e->getEdge(0).length());
+        }
       }
       else {
-        Msg::Warning("Unknown surface %d", *it);
+        Msg::Warning("Unknown curve %d", t);
       }
     }
-
-    // construct a kd-tree index:
-    _kdtree = new SPoint3KDTree(3, _pc2kdtree,
-                                nanoflann::KDTreeSingleIndexAdaptorParams(10));
-    _kdtree->buildIndex();
-    updateNeeded = false;
+    if(_sizeCurves.empty()) {
+      _kdtreeCurves = nullptr;
+    }
+    else {
+      _kdtreeCurves = new SPoint3KDTree(3, _pc2kdtreeCurves,
+                                        nanoflann::KDTreeSingleIndexAdaptorParams(10));
+      _kdtreeCurves->buildIndex();
+    }
+  }
+  void recomputeSurfaces()
+  {
+    _sizeSurfaces.clear();
+    _pcSurfaces.pts.clear();
+    if(_kdtreeSurfaces) delete _kdtreeSurfaces;
+    for(auto t : _tagSurfaces) {
+      GFace *gf = GModel::current()->getFaceByTag(t);
+      if(gf) {
+        for(auto e : gf->triangles) {
+          _pcSurfaces.pts.push_back(e->barycenter());
+          double s = e->getEdge(0).length() + e->getEdge(1).length() +
+            e->getEdge(2).length();
+          _sizeSurfaces.push_back(s / 3.);
+        }
+      }
+      else {
+        Msg::Warning("Unknown surface %d", t);
+      }
+    }
+    if(_sizeSurfaces.empty()) {
+      _kdtreeSurfaces = nullptr;
+    }
+    else {
+      _kdtreeSurfaces = new SPoint3KDTree(3, _pc2kdtreeSurfaces,
+                                          nanoflann::KDTreeSingleIndexAdaptorParams(10));
+      _kdtreeSurfaces->buildIndex();
+    }
   }
   using Field::operator();
   virtual double operator()(double X, double Y, double Z, GEntity *ge = nullptr)
   {
+    if(!ge) return MAX_LC;
+    if(ge->dim() != 2 && ge->dim() != 3) return MAX_LC;
+    if(ge->dim() == 2 && _tagCurves.empty()) return MAX_LC;
+    if(ge->dim() == 3 && _tagSurfaces.empty()) return MAX_LC;
 #pragma omp critical
-    {
-      if(!_kdtree || updateNeeded) update();
+    if(updateNeeded ||
+       (ge->dim() == 2 && _tagCurves.size() && _sizeCurves.empty())) {
+      // we are meshing our first surface; recompute distance to the elements on
+      // curves, and invalidate the distance to surfaces
+      recomputeCurves();
+      _sizeSurfaces.clear();
+      updateNeeded = false;
     }
-    double query_pt[3] = {X, Y, Z};
-    nanoflann::KNNResultSet<double> resultSet(1);
-    double outDistSqr;
-    resultSet.init(&_outIndex, &outDistSqr);
-    _kdtree->findNeighbors(resultSet, &query_pt[0], nanoflann::SearchParams(10));
-    return sqrt(outDistSqr);
+#pragma omp critical
+    if(updateNeeded ||
+       (ge->dim() == 3 && _tagSurfaces.size() && _sizeSurfaces.empty())) {
+      // we are meshing our first volume; recompute distance to the elements on
+      // surfaces, and invalidate the distance to curves (to be ready for
+      // subsequent surface meshing pass)
+      recomputeSurfaces();
+      _sizeCurves.clear();
+      updateNeeded = false;
+    }
+    double pt[3] = {X, Y, Z};
+    nanoflann::KNNResultSet<double> res(1);
+    std::size_t index = 0;
+    double dist2 = 0., sbnd = 0.;
+    res.init(&index, &dist2);
+    if(ge->dim() == 2 && _kdtreeCurves) {
+      _kdtreeCurves->findNeighbors(res, &pt[0], nanoflann::SearchParams(10));
+      sbnd = _sizeCurves[index];
+    }
+    else if(ge->dim() == 3 && _kdtreeSurfaces) {
+      _kdtreeSurfaces->findNeighbors(res, &pt[0], nanoflann::SearchParams(10));
+      sbnd = _sizeSurfaces[index];
+    }
+    if(sbnd <= 0 || _distMax <= 0 || _sizeMax <= 0) return MAX_LC;
+    double d = sqrt(dist2);
+    if(d >= _distMax) return _sizeMax;
+    double f = (_distMax - d) / _distMax;
+    if(_power != 1.) {
+      f = std::pow(f, _power);
+    }
+    double s = f * sbnd + (1. - f) * _sizeMax;
+    return s;
   }
 };
 
@@ -3225,6 +3374,7 @@ FieldManager::FieldManager()
   mapTypeName["Gradient"] = new FieldFactoryT<GradientField>();
   mapTypeName["Octree"] = new FieldFactoryT<OctreeField>();
   mapTypeName["Distance"] = new FieldFactoryT<DistanceField>();
+  mapTypeName["Extend"] = new FieldFactoryT<ExtendField>();
   mapTypeName["Restrict"] = new FieldFactoryT<RestrictField>();
   mapTypeName["Constant"] = new FieldFactoryT<ConstantField>();
   mapTypeName["Min"] = new FieldFactoryT<MinField>();
diff --git a/src/mesh/meshGRegionHxt.cpp b/src/mesh/meshGRegionHxt.cpp
index b4d6bda9734ce2a4386628b8757b100a1e55f147..ab87c8a36c342e5c060684c51c92293d249d1380 100644
--- a/src/mesh/meshGRegionHxt.cpp
+++ b/src/mesh/meshGRegionHxt.cpp
@@ -41,7 +41,7 @@ static HXTStatus nodalSizesCallBack(double *pts, uint32_t *volume,
   std::vector<GRegion *> *allGR = (std::vector<GRegion *> *)userData;
 
   double lcGlob = CTX::instance()->lc;
-  int useInterpolatedSize = CTX::instance()->mesh.lcExtendFromBoundary;
+  int useInterpolatedSize = Extend2dMeshIn3dVolumes();
 
   HXT_INFO("Computing %smesh sizes...", useInterpolatedSize ? "interpolated " : "");
 
diff --git a/src/mesh/meshGRegionMMG.cpp b/src/mesh/meshGRegionMMG.cpp
index d5995068e227fadfae5d1b50182a80e4b8eba20f..852079f33cd40fd220867c8080755926480ef404 100644
--- a/src/mesh/meshGRegionMMG.cpp
+++ b/src/mesh/meshGRegionMMG.cpp
@@ -172,7 +172,7 @@ static void gmsh2MMG(GRegion *gr, MMG5_pMesh mmg, MMG5_pSol sol,
     auto itv = LCS.find(v);
     if(itv != LCS.end()) {
       mmg2gmsh[(*it)->getNum()] = *it;
-      // if (CTX::instance()->mesh.lcExtendFromBoundary){
+      // if (Extend2dMeshIn3dVolumes()){
       double LL = itv->second.first / itv->second.second;
       SMetric3 l4(1. / (LL * LL));
       SMetric3 MM = intersection_conserve_mostaniso(l4, m);
@@ -218,7 +218,7 @@ static void updateSizes(GRegion *gr, MMG5_pMesh mmg, MMG5_pSol sol,
   std::vector<GFace *> f = gr->faces();
 
   std::map<MVertex *, std::pair<double, int> > LCS;
-  // if (CTX::instance()->mesh.lcExtendFromBoundary){
+  // if (Extend2dMeshIn3dVolumes()){
   for(auto it = f.begin(); it != f.end(); ++it) {
     for(unsigned int i = 0; i < (*it)->triangles.size(); i++) {
       MTriangle *t = (*it)->triangles[i];
@@ -250,7 +250,7 @@ static void updateSizes(GRegion *gr, MMG5_pMesh mmg, MMG5_pSol sol,
 
     auto it = mmg2gmsh.find(k);
 
-    if(it != mmg2gmsh.end() && CTX::instance()->mesh.lcExtendFromBoundary) {
+    if(it != mmg2gmsh.end() && Extend2dMeshIn3dVolumes()) {
       auto itv = LCS.find(it->second);
       if(itv != LCS.end()) {
         double LL = itv->second.first / itv->second.second;
diff --git a/src/parser/Gmsh.tab.cpp b/src/parser/Gmsh.tab.cpp
index 7fd34f5822a2d64d642f7667735ac7ffb7e77c68..877816585260338e0be59f1185c9433ac7913a6b 100644
--- a/src/parser/Gmsh.tab.cpp
+++ b/src/parser/Gmsh.tab.cpp
@@ -12243,8 +12243,8 @@ yyreduce:
   case 376:
 #line 4771 "Gmsh.y"
     {
-      // lcExtendFromBoundary onstraints are stored in GEO internals in addition
-      // to GModel, as they can be copied around during GEO operations
+      // mesh size from boundary onstraints are stored in GEO internals in
+      // addition to GModel, as they can be copied around during GEO operations
       if(GModel::current()->getOCCInternals() &&
          GModel::current()->getOCCInternals()->getChanged())
         GModel::current()->getOCCInternals()->synchronize(GModel::current());
diff --git a/src/parser/Gmsh.y b/src/parser/Gmsh.y
index 900cd2644b05b1f4e9836f5076259a240aa7490c..04e3a429721d687bb2213794fd0d3762e73e4211 100644
--- a/src/parser/Gmsh.y
+++ b/src/parser/Gmsh.y
@@ -4769,8 +4769,8 @@ Constraints :
     }
   | tMeshSizeFromBoundary tSurface '{' RecursiveListOfDouble '}' tAFFECT FExpr tEND
     {
-      // lcExtendFromBoundary onstraints are stored in GEO internals in addition
-      // to GModel, as they can be copied around during GEO operations
+      // mesh size from boundary onstraints are stored in GEO internals in
+      // addition to GModel, as they can be copied around during GEO operations
       if(GModel::current()->getOCCInternals() &&
          GModel::current()->getOCCInternals()->getChanged())
         GModel::current()->getOCCInternals()->synchronize(GModel::current());
diff --git a/tutorials/c++/t10.cpp b/tutorials/c++/t10.cpp
index d475662abe1459de2adeb39b4345ca92e4b30e64..4c0d1061c2ef7e260dadd25fef92a4cf11858b52 100644
--- a/tutorials/c++/t10.cpp
+++ b/tutorials/c++/t10.cpp
@@ -129,9 +129,9 @@ int main(int argc, char **argv)
   // The value can then be further modified by the mesh size callback, if any,
   // before being constrained in the interval [`Mesh.MeshSizeMin',
   // `Mesh.MeshSizeMax'] and multiplied by `Mesh.MeshSizeFactor'. In addition,
-  // boundary mesh sizes (on curves or surfaces) are interpolated inside the
-  // enclosed entity (surface or volume, respectively) if the option
-  // `Mesh.MeshSizeExtendFromBoundary' is set (which is the case by default).
+  // boundary mesh sizes are interpolated inside surfaces and/or volumes
+  // depending on the value of `Mesh.MeshSizeExtendFromBoundary' (which is set
+  // by default).
   //
   // When the element size is fully specified by a background mesh (as it is in
   // this example), it is thus often desirable to set
diff --git a/tutorials/c/t2.c b/tutorials/c/t2.c
index 28cda7ef5a125d1b5856b52cd1a9f4684e8b8a15..1df37d3122ec42ce9d5b037432da24900109a9c4 100644
--- a/tutorials/c/t2.c
+++ b/tutorials/c/t2.c
@@ -172,8 +172,8 @@ int main(int argc, char **argv)
 
   int *ov2;
   size_t ov2_n;
-  int numElements[] = {};
-  double heights[] = {};
+  int *numElements;
+  double *heights;
   size_t numElements_n = 0;
   size_t heights_n = 0;
   int recombine = 0;
diff --git a/tutorials/python/t10.py b/tutorials/python/t10.py
index b11d8f6a876a343e0b039adcafc826871555b710..c92b27756511fe6413078e7ce11c46004d83e42b 100644
--- a/tutorials/python/t10.py
+++ b/tutorials/python/t10.py
@@ -121,9 +121,8 @@ gmsh.model.mesh.setSizeCallback(meshSizeCallback)
 # The value can then be further modified by the mesh size callback, if any,
 # before being constrained in the interval [`Mesh.MeshSizeMin',
 # `Mesh.MeshSizeMax'] and multiplied by `Mesh.MeshSizeFactor'.  In addition,
-# boundary mesh sizes (on curves or surfaces) are interpolated inside the
-# enclosed entity (surface or volume, respectively) if the option
-# `Mesh.MeshSizeExtendFromBoundary' is set (which is the case by default).
+# boundary mesh sizes are interpolated inside surfaces and/or volumes depending
+# on the value of `Mesh.MeshSizeExtendFromBoundary' (which is set by default).
 #
 # When the element size is fully specified by a background mesh (as it is in
 # this example), it is thus often desirable to set
diff --git a/tutorials/t10.geo b/tutorials/t10.geo
index cf92f2dad31f105c6234ff98cd4fc12ce285d0a6..8e39ea3ba7e4828ce58201b91113e15a1e512978 100644
--- a/tutorials/t10.geo
+++ b/tutorials/t10.geo
@@ -96,10 +96,9 @@ Background Field = 7;
 // 5) any per-entity mesh size constraint.
 //
 // This value is then constrained in the interval [`Mesh.MeshSizeMin',
-// `Mesh.MeshSizeMax'] and multiplied by `Mesh.MeshSizeFactor'.  In addition,
-// boundary mesh sizes (on curves or surfaces) are interpolated inside the
-// enclosed entity (surface or volume, respectively) if the option
-// `Mesh.MeshSizeExtendFromBoundary' is set (which is the case by default).
+// `Mesh.MeshSizeMax'] and multiplied by `Mesh.MeshSizeFactor'. In addition,
+// boundary mesh sizes are interpolated inside surfaces and/or volumes depending
+// on the value of `Mesh.MeshSizeExtendFromBoundary' (which is set by default).
 //
 // When the element size is fully specified by a background mesh size field (as
 // it is in this example), it is thus often desirable to set