[Getfem-commits] [getfem-commits] branch master updated: Implement bspline basis functions for mesh_fem in 2D

[Konstantinos Poulios]

This is an automated email from the git hooks/post-receive script.

logari81 pushed a commit to branch master
in repository getfem.

The following commit(s) were added to refs/heads/master by this push:
     new e7f75a1f Implement bspline basis functions for mesh_fem in 2D
e7f75a1f is described below

commit e7f75a1f61aae4dbad4a1104bea142c3866535df
Author: Konstantinos Poulios <logari81@gmail.com>
AuthorDate: Mon Dec 19 01:02:37 2022 +0100

    Implement bspline basis functions for mesh_fem in 2D
---
 src/getfem/getfem_global_function.h          |   5 +
 src/getfem/getfem_mesh_fem_global_function.h |  11 +
 src/getfem_global_function.cc                | 556 ++++++++++++++++++++++++++-
 src/getfem_mesh_fem_global_function.cc       |  61 ++-
 4 files changed, 629 insertions(+), 4 deletions(-)

diff --git a/src/getfem/getfem_global_function.h 
b/src/getfem/getfem_global_function.h
index 45ddd562..75fe122a 100644
--- a/src/getfem/getfem_global_function.h
+++ b/src/getfem/getfem_global_function.h
@@ -320,6 +320,11 @@ namespace getfem {
   global_function_on_level_sets(const std::vector<level_set> &lsets,
                                 const pxy_function &fn);
 
+  pglobal_function
+  global_function_bspline(scalar_type &xmin, scalar_type &xmax,
+                          scalar_type &ymin, scalar_type &ymax,
+                          size_type &order,
+                          size_type &xtype, size_type &ytype);
 
 }  /* end of namespace getfem.                                            */
 
diff --git a/src/getfem/getfem_mesh_fem_global_function.h 
b/src/getfem/getfem_mesh_fem_global_function.h
index e4686c77..65982de4 100644
--- a/src/getfem/getfem_mesh_fem_global_function.h
+++ b/src/getfem/getfem_mesh_fem_global_function.h
@@ -62,6 +62,17 @@ namespace getfem {
     virtual ~mesh_fem_global_function() { clear(); }
   };
 
+  /** This function will generate bspline basis functions in an NX x NY
+      rectilinear grid. The generated basis spans the entire bounding
+      box of the mesh linked by mf. The function will finally set the
+      generated bspline basis functions as the basis of mf.
+      In case mim is provided, this integration method will be used to
+      determine the support of he basis functions.
+  */
+  void define_bspline_basis_functions_for_mesh_fem
+  (mesh_fem_global_function &mf,
+   size_type NX, size_type NY, size_type order,
+   const mesh_im &mim=dummy_mesh_im());
 
 }  /* end of namespace getfem.                                            */
 
diff --git a/src/getfem_global_function.cc b/src/getfem_global_function.cc
index d552e188..e96503c6 100644
--- a/src/getfem_global_function.cc
+++ b/src/getfem_global_function.cc
@@ -1,7 +1,7 @@
 /*===========================================================================
 
- Copyright (C) 2004-2020 Yves Renard
- Copyright (C) 2016-2020 Konstantinos Poulios
+ Copyright (C) 2004-2022 Yves Renard
+ Copyright (C) 2016-2022 Konstantinos Poulios
 
  This file is a part of GetFEM
 
@@ -831,6 +831,558 @@ namespace getfem {
     return std::make_shared<global_function_on_levelsets_2D_>(ls, fn);
   }
 
+
+
+
+  // Global function for bspline basis
+  const scalar_type eps(1e-12);
+
+  // order k = 3
+  scalar_type bsp3_01(scalar_type t) {
+    return (t >= -eps && t < 1) ? pow(1.-t,2)
+                                : 0;
+  }
+  scalar_type bsp3_01_der(scalar_type t) {
+    return (t >= -eps && t < 1) ? 2.*t-2.
+                                : 0;
+  }
+  scalar_type bsp3_01_der2(scalar_type t) {
+    return (t >= eps && t <= 1.-eps) ? 2.
+                                     : 0;
+  }
+  scalar_type bsp3_01_der2_with_hint(scalar_type t, scalar_type t_hint) {
+    return (t > -eps && t < 1.+eps && t_hint > 0 && t_hint < 1) ? 2.
+                                                                : 0;
+  }
+  scalar_type bsp3_02(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1)
+        return (-1.5*t+2.)*t;
+      else if (t < 2)
+        return 0.5*pow(2.-t,2);
+    }
+    return 0;
+  }
+  scalar_type bsp3_02_der(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1)
+        return -3.*t+2.;
+      else if (t < 2)
+        return t-2.;
+    }
+    return 0;
+  }
+  scalar_type bsp3_02_der2(scalar_type t) {
+    if (t >= eps) {
+      if (t < 1.-eps)
+        return -3.;
+      else if (t < 1.+eps)
+        return 0;
+      else if (t <= 2.-eps)
+        return 1.;
+    }
+    return 0;
+  }
+  scalar_type bsp3_03(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return 0.5*t*t;
+      } else if (t < 2) {
+        t -= 1.;
+        return t*(1-t)+0.5;
+      } else if (t < 3) {
+        t = 3.-t;
+        return 0.5*t*t;
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp3_03_der(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return t;
+      } else if (t < 2) {
+        t -= 1.;
+        return 1.-2.*t;
+      } else if (t < 3) {
+        return t-3.;
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp3_03_der2(scalar_type t) {
+    if (t >= -eps) {
+      if (t < eps)
+        return 0;
+      else if (t <= 1.-eps)
+        return 1.;
+      else if (t < 1.+eps)
+        return 0;
+      else if (t <= 2.-eps)
+        return -2.;
+      else if (t < 2.+eps)
+        return 0;
+      else if (t <= 3.-eps)
+        return 1.;
+      else if (t < 3.+eps)
+        return 0;
+    }
+    return 0;
+  }
+
+  // order k = 4
+  scalar_type bsp4_01(scalar_type t) {
+    return (t >= -eps && t < 1) ? pow(1.-t,3)
+                                : 0;
+  }
+  scalar_type bsp4_01_der(scalar_type t) {
+    return (t >= -eps && t < 1) ? -3*pow(1.-t,2)
+                                : 0;
+  }
+  scalar_type bsp4_01_der2(scalar_type t) {
+    return (t >= -eps && t < 1) ? 6*(1.-t)
+                                : 0;
+  }
+  scalar_type bsp4_02(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return ((7./4.*t-9./2.)*t+3.)*t;
+      } else if (t < 2) {
+        return 1./4.*pow(2.-t,3);
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp4_02_der(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return (21./4.*t-9.)*t+3.;
+      } else if (t < 2) {
+        return -3./4.*pow(2.-t,2);
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp4_02_der2(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return 21./2.*t-9.;
+      } else if (t < 2) {
+        return 3./2.*(2.-t);
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp4_03(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return (-11./12.*t+3./2.)*t*t;
+      } else if (t < 2) {
+        t -= 1;
+        return ((7./12.*t-5./4.)*t+1./4.)*t+7./12.;
+      } else if (t < 3) {
+        return 1./6.*pow(3.-t,3);
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp4_03_der(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return (-11./4.*t+3.)*t;
+      } else if (t < 2) {
+        t -= 1;
+        return (7./4.*t-5./2.)*t+1./4.;
+      } else if (t < 3) {
+        return -1./2.*pow(3.-t,2);
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp4_03_der2(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return -11./2.*t+3.;
+      } else if (t < 2) {
+        t -= 1;
+        return 7./2.*t-5./2.;
+      } else if (t < 3) {
+        return 3.-t;
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp4_04(scalar_type t) {
+    if (t > 2) {
+      t = 4.-t;
+    }
+    if (t >= -eps) {
+      if (t < 1) {
+        return 1./6.*pow(t,3);
+      } else if (t <= 2) {
+        t = t-1.;
+        return ((-1./2.*t+1./2.)*t+1./2.)*t+1./6.;
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp4_04_der(scalar_type t) {
+    scalar_type sgn(1);
+    if (t > 2) {
+      t = 4.-t;
+      sgn = -1.;
+    }
+    if (t >= -eps) {
+      if (t < 1) {
+        return 1./2.*pow(t,2)*sgn;
+      } else if (t < 2) {
+        t = t-1.;
+        return ((-3./2.*t+1.)*t+1./2.)*sgn;
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp4_04_der2(scalar_type t) {
+    if (t > 2) {
+      t = 4.-t;
+    }
+    if (t >= -eps) {
+      if (t < 1) {
+        return t;
+      } else if (t < 2) {
+        t = t-1.;
+        return -3.*t+1.;
+      }
+    }
+    return 0;
+  }
+
+
+  // order k = 5
+  scalar_type bsp5_01(scalar_type t) {
+    return (t >= -eps && t < 1) ? pow(1.-t,4)
+                                : 0;
+  }
+  scalar_type bsp5_01_der(scalar_type t) {
+    return (t >= -eps && t < 1) ? -4.*pow(1.-t,3)
+                                : 0;
+  }
+  scalar_type bsp5_01_der2(scalar_type t) {
+    return (t >= -eps && t < 1) ? 12.*pow(1.-t,2)
+                                : 0;
+  }
+  scalar_type bsp5_02(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return (((-15./8.*t+7.)*t-9.)*t+4.)*t;
+      } else if (t < 2) {
+        return 1./8.*pow(2.-t,4);
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp5_02_der(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return ((-15./2.*t+21.)*t-18.)*t+4.;
+      } else if (t < 2) {
+        return -1./2.*pow(2.-t,3);
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp5_02_der2(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return (-45./2.*t+42.)*t-18.;
+      } else if (t < 2) {
+        return 3./2.*pow(2.-t,2);
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp5_03(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return ((85./72.*t-11./3.)*t+3.)*t*t;
+      } else if (t < 2) {
+        t = 2.-t;
+        return (((-23./72*t+2./9.)*t+1./3.)*t+2./9.)*t+1./18.;
+      } else if (t < 3) {
+        return 1./18.*pow(3.-t,4);
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp5_03_der(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return ((85./18.*t-11.)*t+6.)*t;
+      } else if (t < 2) {
+        t = 2.-t;
+        return -(((-23./18*t+2./3.)*t+2./3.)*t+2./9.);
+      } else if (t < 3) {
+        return -2./9.*pow(3.-t,3);
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp5_03_der2(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return (85./6.*t-22.)*t+6.;
+      } else if (t < 2) {
+        t = 2.-t;
+        return (-23./6*t+4./3.)*t+2./3.;
+      } else if (t < 3) {
+        return 2./3.*pow(3.-t,2);
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp5_04(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return (-25./72.*t+2./3.)*t*t*t;
+      } else if (t < 2) {
+        t = 2.-t;
+        return (((23./72.*t-5./9.)*t-1./3.)*t+4./9.)*t+4./9.;
+      } else if (t < 3) {
+        t = t-2.;
+        return (((-13./72.*t+5./9.)*t-1./3.)*t-4./9.)*t+4./9.;
+      } else if (t < 4) {
+        return 1./24.*pow(4.-t,4);
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp5_04_der(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return (-25./18.*t+2.)*t*t;
+      } else if (t < 2) {
+        t = 2.-t;
+        return -(((23./18.*t-5./3.)*t-2./3.)*t+4./9.);
+      } else if (t < 3) {
+        t = t-2.;
+        return ((-13./18.*t+5./3.)*t-2./3.)*t-4./9.;
+      } else if (t < 4) {
+        return -1./6.*pow(4.-t,3);
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp5_04_der2(scalar_type t) {
+    if (t >= -eps) {
+      if (t < 1) {
+        return (-25./6.*t+4.)*t;
+      } else if (t < 2) {
+        t = 2.-t;
+        return (23./6.*t-10./3.)*t-2./3.;
+      } else if (t < 3) {
+        t = t-2.;
+        return -((-13./6.*t+10./3.)*t-2./3.);
+      } else if (t < 4) {
+        return 1./2.*pow(4.-t,2);
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp5_05(scalar_type t) {
+    if (t > scalar_type(2.5)) {
+      t = 5.-t;
+    }
+    if (t >= -eps) {
+      if (t < 1) {
+        return 1./24.*pow(t,4);
+      } else if (t < 2) {
+        t = t-1.;
+        return (((-1./6.*t+1./6.)*t+1./4.)*t+1./6.)*t+1./24.;
+      } else if (t < 3) {
+        t = pow(t-2.5,2);
+        return (1./4.*t-5./8.)*t+115./192.;
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp5_05_der(scalar_type t) {
+    scalar_type sgn(1);
+    if (t > scalar_type(2.5)) {
+      t = 5.-t;
+      sgn = -1.;
+    }
+    if (t >= -eps) {
+      if (t < 1) {
+        return 1./6.*pow(t,3)*sgn;
+      } else if (t < 2) {
+        t = t-1.;
+        return (((-2./3.*t+1./2.)*t+1./2)*t+1./6.)*sgn;
+      } else if (t < 3) {
+        t = t-2.5;
+        return (t*t-5./4.)*t*sgn;
+      }
+    }
+    return 0;
+  }
+  scalar_type bsp5_05_der2(scalar_type t) {
+    if (t > scalar_type(2.5)) {
+      t = 5.-t;
+    }
+    if (t >= -eps) {
+      if (t < 1) {
+        return 1./2.*pow(t,2);
+      } else if (t < 2) {
+        t = t-1.;
+        return (-2.*t+1.)*t+1./2;
+      } else if (t < 3) {
+        t = t-2.5;
+        return 3*t*t-5./4.;
+      }
+    }
+    return 0;
+  }
+
+
+
+  class global_function_xy_bspline_
+    : public global_function_simple, public context_dependencies {
+    scalar_type xmin, ymin, xmax, ymax, xscale, yscale;
+    std::function<scalar_type(scalar_type)>
+      fx, fy, fx_der, fy_der, fx_der2, fy_der2;
+  public:
+    void update_from_context() const {}
+
+    virtual scalar_type val(const base_node &pt) const
+    {
+      return fx(xscale*(pt[0]-xmin))
+             * fy(yscale*(pt[1]-ymin));
+    }
+    virtual void grad(const base_node &pt, base_small_vector &g) const
+    {
+      scalar_type dx = xscale*(pt[0]-xmin),
+                  dy = yscale*(pt[1]-ymin);
+      g.resize(dim_);
+      g[0] = xscale * fx_der(dx) * fy(dy);
+      g[1] = fx(dx) * yscale * fy_der(dy);
+    }
+    virtual void hess(const base_node &pt, base_matrix &h) const
+    {
+      scalar_type dx = xscale*(pt[0]-xmin),
+                  dy = yscale*(pt[1]-ymin);
+      h.resize(dim_, dim_);
+      gmm::clear(h);
+      h(0,0) = xscale*xscale * fx_der2(dx) * fy(dy);
+      h(0,1) = h(1,0) = xscale * fx_der(dx) * yscale * fy_der(dy);
+      h(1,1) = fx(dx) * yscale*yscale * fy_der2(dy);
+    }
+
+    virtual bool is_in_support(const base_node &pt) const {
+      scalar_type dx = pt[0]-(xmin+xmax)/2,
+                  dy = pt[1]-(ymin+ymax)/2;
+      return (gmm::abs(dx)+1e-9 < gmm::abs(xmax-xmin)/2) &&
+             (gmm::abs(dy)+1e-9 < gmm::abs(ymax-ymin)/2);
+    }
+
+    virtual void bounding_box(base_node &bmin, base_node &bmax) const {
+      GMM_ASSERT1(bmin.size() == dim_ && bmax.size() == dim_,
+                  "Wrong dimensions");
+      bmin[0] = std::min(xmin,xmax);
+      bmin[1] = std::min(ymin,ymax);
+      bmax[0] = std::max(xmin,xmax);
+      bmax[1] = std::max(ymin,ymax);
+    }
+
+    global_function_xy_bspline_(const scalar_type &xmin_, const scalar_type 
&xmax_,
+                                const scalar_type &ymin_, const scalar_type 
&ymax_,
+                                const size_type &order, 
+                                const size_type &xtype, const size_type &ytype)
+    : global_function_simple(2), xmin(xmin_), ymin(ymin_), xmax(xmax_), 
ymax(ymax_),
+      xscale(scalar_type(xtype)/(xmax-xmin)), 
yscale(scalar_type(ytype)/(ymax-ymin))
+    {
+      GMM_ASSERT1(order >= 3 && order <= 5, "Wrong input");
+      GMM_ASSERT1(xtype >= 1 && xtype <= order &&
+                  ytype >= 1 && ytype <= order, "Wrong input");
+      if (order == 3) {
+        if (xtype == 1) {
+          fx = bsp3_01;   fx_der = bsp3_01_der;   fx_der2 = bsp3_01_der2;
+        } else if (xtype == 2) {
+          fx = bsp3_02;   fx_der = bsp3_02_der;   fx_der2 = bsp3_02_der2;
+        } else if (xtype == 3) {
+          fx = bsp3_03;   fx_der = bsp3_03_der;   fx_der2 = bsp3_03_der2;
+        }
+
+        if (ytype == 1) {
+          fy = bsp3_01;   fy_der = bsp3_01_der;   fy_der2 = bsp3_01_der2;
+        } else if (ytype == 2) {
+          fy = bsp3_02;   fy_der = bsp3_02_der;   fy_der2 = bsp3_02_der2;
+        } else if (ytype == 3) {
+          fy = bsp3_03;   fy_der = bsp3_03_der;   fy_der2 = bsp3_03_der2;
+        }
+      } else if (order == 4) {
+        if (xtype == 1) {
+          fx = bsp4_01;   fx_der = bsp4_01_der;   fx_der2 = bsp4_01_der2;
+        } else if (xtype == 2) {
+          fx = bsp4_02;   fx_der = bsp4_02_der;   fx_der2 = bsp4_02_der2;
+        } else if (xtype == 3) {
+          fx = bsp4_03;   fx_der = bsp4_03_der;   fx_der2 = bsp4_03_der2;
+        } else if (xtype == 4) {
+          fx = bsp4_04;   fx_der = bsp4_04_der;   fx_der2 = bsp4_04_der2;
+        }
+
+        if (ytype == 1) {
+          fy = bsp4_01;   fy_der = bsp4_01_der;   fy_der2 = bsp4_01_der2;
+        } else if (ytype == 2) {
+          fy = bsp4_02;   fy_der = bsp4_02_der;   fy_der2 = bsp4_02_der2;
+        } else if (ytype == 3) {
+          fy = bsp4_03;   fy_der = bsp4_03_der;   fy_der2 = bsp4_03_der2;
+        } else if (ytype == 4) {
+          fy = bsp4_04;   fy_der = bsp4_04_der;   fy_der2 = bsp4_04_der2;
+        }
+
+      } else if (order == 5) {
+        if (xtype == 1) {
+          fx = bsp5_01;   fx_der = bsp5_01_der;   fx_der2 = bsp5_01_der2;
+        } else if (xtype == 2) {
+          fx = bsp5_02;   fx_der = bsp5_02_der;   fx_der2 = bsp5_02_der2;
+        } else if (xtype == 3) {
+          fx = bsp5_03;   fx_der = bsp5_03_der;   fx_der2 = bsp5_03_der2;
+        } else if (xtype == 4) {
+          fx = bsp5_04;   fx_der = bsp5_04_der;   fx_der2 = bsp5_04_der2;
+        } else if (xtype == 5) {
+          fx = bsp5_05;   fx_der = bsp5_05_der;   fx_der2 = bsp5_05_der2;
+        }
+
+        if (ytype == 1) {
+          fy = bsp5_01;   fy_der = bsp5_01_der;   fy_der2 = bsp5_01_der2;
+        } else if (ytype == 2) {
+          fy = bsp5_02;   fy_der = bsp5_02_der;   fy_der2 = bsp5_02_der2;
+        } else if (ytype == 3) {
+          fy = bsp5_03;   fy_der = bsp5_03_der;   fy_der2 = bsp5_03_der2;
+        } else if (ytype == 4) {
+          fy = bsp5_04;   fy_der = bsp5_04_der;   fy_der2 = bsp5_04_der2;
+        } else if (ytype == 5) {
+          fy = bsp5_05;   fy_der = bsp5_05_der;   fy_der2 = bsp5_05_der2;
+        }
+      }
+    }
+
+    virtual ~global_function_xy_bspline_()
+    { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function xy bspline"); }
+  };
+
+
+  pglobal_function
+  global_function_bspline(scalar_type &xmin, scalar_type &xmax,
+                          scalar_type &ymin, scalar_type &ymax,
+                          size_type &order,
+                          size_type &xtype, size_type &ytype) {
+    return std::make_shared<global_function_xy_bspline_>
+                           (xmin, xmax, ymin, ymax, order, xtype, ytype);
+  }
+
+
+
+
   // interpolator on mesh_fem
 
   interpolator_on_mesh_fem::interpolator_on_mesh_fem
diff --git a/src/getfem_mesh_fem_global_function.cc 
b/src/getfem_mesh_fem_global_function.cc
index cffb4785..d595d0ae 100644
--- a/src/getfem_mesh_fem_global_function.cc
+++ b/src/getfem_mesh_fem_global_function.cc
@@ -1,7 +1,7 @@
 /*===========================================================================
 
- Copyright (C) 2004-2020 Yves Renard
- Copyright (C) 2016-2020 Konstantinos Poulios
+ Copyright (C) 2004-2022 Yves Renard
+ Copyright (C) 2016-2022 Konstantinos Poulios
 
  This file is a part of GetFEM
 
@@ -51,6 +51,63 @@ namespace getfem {
     }
   }
 
+
+  void define_bspline_basis_functions_for_mesh_fem
+  (mesh_fem_global_function &mf,
+   size_type NX, size_type NY, size_type order,
+   const mesh_im &mim) {
+
+    base_node Pmin, Pmax;
+    mf.linked_mesh().bounding_box(Pmin, Pmax);
+    scalar_type x0=Pmin[0], y0=Pmin[1], x1=Pmax[0], y1=Pmax[1];
+
+    scalar_type xmin, xmax, ymin, ymax;
+    size_type xtype, ytype;
+    base_node pt(2);
+
+    std::vector<pglobal_function> funcs((NX+order-1)*(NY+order-1));
+    for (size_type i=0; i < NX+order-1; ++i) {
+      if (i < order-1) {
+        xmin = x0;
+        xmax = x0+scalar_type(i+1)*(x1-x0)/scalar_type(NX);
+        xtype = i+1;
+        pt[0] = (i == 0) ? xmin : (xmin+(xmax-xmin)/3);
+      } else if (i >= NX) {
+        xmin = x1;
+        xmax = 
x1+(scalar_type(i-NX)-scalar_type(order-1))*(x1-x0)/scalar_type(NX);
+        xtype = NX-i+order-1;
+        pt[0] = (i == NX+1) ? xmin : (xmin+(xmax-xmin)/3);
+      } else {
+        xmin = x0+scalar_type(i-order+1)*(x1-x0)/scalar_type(NX);
+        xmax = x0+scalar_type(i+1)*(x1-x0)/scalar_type(NX);
+        xtype = order;
+        pt[0] = (xmin+xmax)/2;
+      }
+      for (size_type j=0; j < NY+order-1; ++j) {
+        if (j < order-1) {
+          ymin = y0;
+          ymax = y0+scalar_type(j+1)*(y1-y0)/scalar_type(NY);
+          ytype = j+1;
+          pt[1] = (j == 0) ? ymin : (ymin+(ymax-ymin)/3);
+        } else if (j >= NY) {
+          ymin = y1;
+          ymax = 
y1+(scalar_type(j-NY)-scalar_type(order-1))*(y1-y0)/scalar_type(NY);
+          ytype = NY-j+order-1;
+          pt[1] = (j == NY+1) ? ymin : (ymin+(ymax-ymin)/3);
+        } else {
+          ymin = y0+scalar_type(j-order+1)*(y1-y0)/scalar_type(NY);
+          ymax = y0+scalar_type(j+1)*(y1-y0)/scalar_type(NY);
+          ytype = order;
+          pt[1] = (ymin+ymax)/2;
+        }
+        funcs[i*(NY+order-1)+j] = global_function_bspline
+                                  (xmin, xmax, ymin, ymax, order, xtype, 
ytype);
+      }
+    }
+
+    mf.set_functions(funcs, mim);
+  }
+
 }
 
 /* end of namespace getfem  */