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 6da9ff14 Add 3D version of bspline basis functions on uniform
rectilinear grid
6da9ff14 is described below
commit 6da9ff14cc0ccdea3c26781f7004d1789f67a7bd
Author: Konstantinos Poulios <logar...@gmail.com>
AuthorDate: Mon May 6 15:50:55 2024 +0200
Add 3D version of bspline basis functions on uniform rectilinear grid
---
interface/src/gf_mesh_fem.cc | 60 ++++++---
src/getfem/getfem_global_function.h | 7 +
src/getfem/getfem_mesh_fem_global_function.h | 30 +++++
src/getfem_global_function.cc | 192 +++++++++++++++++++++++++++
4 files changed, 274 insertions(+), 15 deletions(-)
diff --git a/interface/src/gf_mesh_fem.cc b/interface/src/gf_mesh_fem.cc
index cf675d81..ee26fbfe 100644
--- a/interface/src/gf_mesh_fem.cc
+++ b/interface/src/gf_mesh_fem.cc
@@ -253,38 +253,46 @@ void gf_mesh_fem(getfemint::mexargs_in& m_in,
);
- /*@INIT MF = ('bspline_uniform', @tmesh m, @int NX[, @int NY,] @int
order[, @str bcX_low[, @str bcY_low[, @str bcX_high][, @str bcY_high]]])
+ /*@INIT MF = ('bspline_uniform', @tmesh m, @int NX[, @int NY[, @int NZ]],
@int order[, @str bcX_low[, @str bcY_low[, @str bcZ_low]][, @str bcX_high[,
@str bcY_high[, @str bcZ_high]]]])
Create a @tmf on mesh `m`, whose base functions are global functions
- corresponding to bspline basis of order `order`, in an NX x NY grid
- (just NX in 1s) that spans the entire bounding box of `m`.
+ corresponding to bspline basis of order `order`, in an NX x NY x NZ
+ grid (just NX in 1D or NX x NY in 2D) that spans the entire bounding
+ box of `m`.
Optionally boundary conditions at the edges of the domain can be
- defined with `bcX_low`, `bcY_low`, `bcX_high`, abd `bcY_high` set to
- 'free' (default) or 'periodic' or 'symmetry'. @*/
+ defined with `bcX_low`, `bcY_low`, `bcZ_low`, `bcX_high`, `bcY_high`,
+ and `bcZ_high` set to 'free' (default) or 'periodic' or 'symmetry'. @*/
sub_command
- ("bspline_uniform", 3, 8, 0, 1,
+ ("bspline_uniform", 3, 11, 0, 1,
mm = extract_mesh_object(in.pop());
dim_type dim = mm->dim();
- if (dim > 2)
- THROW_ERROR("Uniform bspline only supported for dim = 1 or 2");
+ if (dim > 3)
+ THROW_ERROR("Uniform bspline only supported for dim = 1,2,3");
size_type NX = in.pop().to_integer(1,1000);
size_type NY = (dim >= 2) ? in.pop().to_integer(1,1000) : 0;
- if (dim == 2 && (!in.remaining() || !in.front().is_integer()))
- THROW_ERROR("In 2d, 3 integers are expected for NX,NY,order");
+ size_type NZ = (dim == 3) ? in.pop().to_integer(1,1000) : 0;
+ if (!in.remaining() || !in.front().is_integer())
+ THROW_ERROR("One integer was expected for bspline order");
size_type order = in.pop().to_integer(3,5);
std::string bcx_low("free");
std::string bcy_low("free");
+ std::string bcz_low("free");
std::string bcx_high("");
std::string bcy_high("");
+ std::string bcz_high("");
if (in.remaining()) bcx_low = in.pop().to_string();
- if (dim == 2 && in.remaining()) bcy_low = in.pop().to_string();
+ if (dim >= 2 && in.remaining()) bcy_low = in.pop().to_string();
+ if (dim == 3 && in.remaining()) bcz_low = in.pop().to_string();
if (in.remaining()) bcx_high = in.pop().to_string();
- if (dim == 2 && in.remaining()) bcy_high = in.pop().to_string();
- if (dim == 1 && in.remaining())
- THROW_ERROR("Too many arguments for 1d bspline");
+ if (dim >= 2 && in.remaining()) bcy_high = in.pop().to_string();
+ if (dim == 3 && in.remaining()) bcz_high = in.pop().to_string();
+ if (in.remaining())
+ THROW_ERROR("Too many arguments for bspline mesh_fem");
getfem::bspline_boundary bcX_low(getfem::bspline_boundary::FREE);
getfem::bspline_boundary bcY_low(getfem::bspline_boundary::FREE);
+ getfem::bspline_boundary bcZ_low(getfem::bspline_boundary::FREE);
getfem::bspline_boundary bcX_high(getfem::bspline_boundary::FREE);
getfem::bspline_boundary bcY_high(getfem::bspline_boundary::FREE);
+ getfem::bspline_boundary bcZ_high(getfem::bspline_boundary::FREE);
if (bcx_low == "periodic")
bcX_high = bcX_low = getfem::bspline_boundary::PERIODIC;
else if (bcx_low == "symmetry")
@@ -299,6 +307,13 @@ void gf_mesh_fem(getfemint::mexargs_in& m_in,
else if (bcy_low != "free")
THROW_ERROR("Unknown boundary condition " << bcy_low);
+ if (bcz_low == "periodic")
+ bcZ_high = bcZ_low = getfem::bspline_boundary::PERIODIC;
+ else if (bcz_low == "symmetry")
+ bcZ_high = bcZ_low = getfem::bspline_boundary::SYMMETRY;
+ else if (bcz_low != "free")
+ THROW_ERROR("Unknown boundary condition " << bcz_low);
+
if (!bcx_high.empty()) {
if (bcx_high == "periodic")
bcX_high = getfem::bspline_boundary::PERIODIC;
@@ -321,14 +336,29 @@ void gf_mesh_fem(getfemint::mexargs_in& m_in,
THROW_ERROR("Unknown boundary condition " << bcy_high);
}
+ if (!bcz_high.empty()) {
+ if (bcz_high == "periodic")
+ bcZ_high = getfem::bspline_boundary::PERIODIC;
+ else if (bcz_high == "symmetry")
+ bcZ_high = getfem::bspline_boundary::SYMMETRY;
+ else if (bcz_high == "free")
+ bcZ_high = getfem::bspline_boundary::FREE;
+ else
+ THROW_ERROR("Unknown boundary condition " << bcz_high);
+ }
+
auto mfgf = std::make_shared<getfem::mesh_fem_global_function>(*mm);
mfgf->set_qdim(1.);
if (dim == 1)
define_uniform_bspline_basis_functions_for_mesh_fem
(*mfgf, NX, order, bcX_low, bcX_high);
- else
+ else if (dim == 2)
define_uniform_bspline_basis_functions_for_mesh_fem
(*mfgf, NX, NY, order, bcX_low, bcY_low, bcX_high, bcY_high);
+ else
+ define_uniform_bspline_basis_functions_for_mesh_fem
+ (*mfgf, NX, NY, NZ, order, bcX_low, bcY_low, bcZ_low,
+ bcX_high, bcY_high, bcZ_high);
mmf = mfgf;
);
diff --git a/src/getfem/getfem_global_function.h
b/src/getfem/getfem_global_function.h
index e40ea564..7f805fa7 100644
--- a/src/getfem/getfem_global_function.h
+++ b/src/getfem/getfem_global_function.h
@@ -331,6 +331,13 @@ namespace getfem {
const size_type order,
const size_type xtype, const size_type ytype);
+ pglobal_function
+ global_function_bspline(const scalar_type xmin, const scalar_type xmax,
+ const scalar_type ymin, const scalar_type ymax,
+ const scalar_type zmin, const scalar_type zmax,
+ const size_type order,
+ const size_type xtype, const size_type ytype,
+ const size_type ztype);
} /* end of namespace getfem. */
#endif
diff --git a/src/getfem/getfem_mesh_fem_global_function.h
b/src/getfem/getfem_mesh_fem_global_function.h
index a85fb9c6..91b25702 100644
--- a/src/getfem/getfem_mesh_fem_global_function.h
+++ b/src/getfem/getfem_mesh_fem_global_function.h
@@ -110,6 +110,36 @@ namespace getfem {
(mf, NX, NY, order, bcX_low, bcY_low, bcX_low, bcY_low, mim);
}
+
+ /** This function will generate bspline basis functions in an
+ NX x NY x NZ rectilinear grid. The generated basis spans the
+ entire bounding box of the 3d mesh linked by mf. The generated
+ bspline basis functions are then set as the basis of mf.
+ In case mim is provided, this integration method will be used to
+ determine the support of he basis functions more precisely.
+ */
+ void define_uniform_bspline_basis_functions_for_mesh_fem
+ (mesh_fem_global_function &mf,
+ size_type NX, size_type NY, size_type NZ, size_type order,
+ bspline_boundary bcX_low=bspline_boundary::FREE,
+ bspline_boundary bcY_low=bspline_boundary::FREE,
+ bspline_boundary bcZ_low=bspline_boundary::FREE,
+ bspline_boundary bcX_high=bspline_boundary::FREE,
+ bspline_boundary bcY_high=bspline_boundary::FREE,
+ bspline_boundary bcZ_high=bspline_boundary::FREE,
+ const mesh_im &mim=dummy_mesh_im());
+
+ inline void define_uniform_bspline_basis_functions_for_mesh_fem
+ (mesh_fem_global_function &mf,
+ size_type NX, size_type NY, size_type NZ, size_type order,
+ bspline_boundary bcX_low, bspline_boundary bcY_low,
+ bspline_boundary bcZ_low, const mesh_im &mim=dummy_mesh_im()) {
+ define_uniform_bspline_basis_functions_for_mesh_fem
+ (mf, NX, NY, NZ, order, bcX_low, bcY_low, bcZ_low,
+ bcX_low, bcY_low, bcZ_low, mim);
+ }
+
+
} /* end of namespace getfem. */
#endif
diff --git a/src/getfem_global_function.cc b/src/getfem_global_function.cc
index 80bdc66d..31077a57 100644
--- a/src/getfem_global_function.cc
+++ b/src/getfem_global_function.cc
@@ -1455,6 +1455,186 @@ namespace getfem {
};
+ class global_function_xyz_bspline_
+ : public global_function_simple, public context_dependencies {
+ scalar_type xmin, ymin, zmin, xmax, ymax, zmax, xscale, yscale, zscale;
+ std::function<scalar_type(scalar_type)>
+ fx, fy, fz, fx_der, fy_der, fz_der, fx_der2, fy_der2, fz_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))
+ * fz(zscale*(pt[2]-zmin));
+ }
+ virtual void grad(const base_node &pt, base_small_vector &g) const
+ {
+ scalar_type dx = xscale*(pt[0]-xmin),
+ dy = yscale*(pt[1]-ymin),
+ dz = zscale*(pt[2]-zmin);
+ g.resize(dim_);
+ g[0] = xscale * fx_der(dx) * fy(dy) * fz(dz);
+ g[1] = fx(dx) * yscale * fy_der(dy) * fz(dz);
+ g[2] = fx(dx) * fy(dy) * zscale * fz_der(dz);
+ }
+ virtual void hess(const base_node &pt, base_matrix &h) const
+ {
+ scalar_type dx = xscale*(pt[0]-xmin),
+ dy = yscale*(pt[1]-ymin),
+ dz = zscale*(pt[2]-zmin);
+ h.resize(dim_, dim_);
+ gmm::clear(h);
+ h(0,0) = xscale*xscale * fx_der2(dx) * fy(dy) * fz(dz);
+ h(1,1) = fx(dx) * yscale*yscale * fy_der2(dy) * fz(dz);
+ h(2,2) = fx(dx) * fy(dy) * zscale*zscale * fz_der2(dz);
+ h(0,1) = h(1,0) = xscale * fx_der(dx) * yscale * fy_der(dy) * fz(dz);
+ h(0,2) = h(2,0) = xscale * fx_der(dx) * fy(dy) * zscale * fz_der(dz);
+ h(1,2) = h(2,1) = fx(dx) * yscale * fy_der(dy) * zscale * fz_der(dz);
+ }
+
+ virtual bool is_in_support(const base_node &pt) const {
+ scalar_type dx = pt[0]-(xmin+xmax)/2,
+ dy = pt[1]-(ymin+ymax)/2,
+ dz = pt[2]-(zmin+zmax)/2;
+ return (gmm::abs(dx)+1e-9 < gmm::abs(xmax-xmin)/2) &&
+ (gmm::abs(dy)+1e-9 < gmm::abs(ymax-ymin)/2) &&
+ (gmm::abs(dz)+1e-9 < gmm::abs(zmax-zmin)/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);
+ bmin[2] = std::min(zmin,zmax);
+ bmax[0] = std::max(xmin,xmax);
+ bmax[1] = std::max(ymin,ymax);
+ bmax[2] = std::max(zmin,zmax);
+ }
+
+ global_function_xyz_bspline_(const scalar_type &xmin_, const scalar_type
&xmax_,
+ const scalar_type &ymin_, const scalar_type
&ymax_,
+ const scalar_type &zmin_, const scalar_type
&zmax_,
+ const size_type &order,
+ const size_type &xtype,
+ const size_type &ytype,
+ const size_type &ztype)
+ : global_function_simple(3), xmin(xmin_), ymin(ymin_), zmin(zmin_),
+ xmax(xmax_), ymax(ymax_), zmax(zmax_),
+ xscale(scalar_type(xtype)/(xmax-xmin)),
+ yscale(scalar_type(ytype)/(ymax-ymin)),
+ zscale(scalar_type(ztype)/(zmax-zmin))
+ {
+ GMM_ASSERT1(order >= 3 && order <= 5, "Wrong input");
+ GMM_ASSERT1(xtype >= 1 && xtype <= order &&
+ ytype >= 1 && ytype <= order &&
+ ztype >= 1 && ztype <= 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;
+ }
+
+ if (ztype == 1) {
+ fz = bsp3_01; fz_der = bsp3_01_der; fz_der2 = bsp3_01_der2;
+ } else if (ztype == 2) {
+ fz = bsp3_02; fz_der = bsp3_02_der; fz_der2 = bsp3_02_der2;
+ } else if (ztype == 3) {
+ fz = bsp3_03; fz_der = bsp3_03_der; fz_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;
+ }
+
+ if (ztype == 1) {
+ fz = bsp4_01; fz_der = bsp4_01_der; fz_der2 = bsp4_01_der2;
+ } else if (ztype == 2) {
+ fz = bsp4_02; fz_der = bsp4_02_der; fz_der2 = bsp4_02_der2;
+ } else if (ztype == 3) {
+ fz = bsp4_03; fz_der = bsp4_03_der; fz_der2 = bsp4_03_der2;
+ } else if (ztype == 4) {
+ fz = bsp4_04; fz_der = bsp4_04_der; fz_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;
+ }
+
+ if (ztype == 1) {
+ fz = bsp5_01; fz_der = bsp5_01_der; fz_der2 = bsp5_01_der2;
+ } else if (ztype == 2) {
+ fz = bsp5_02; fz_der = bsp5_02_der; fz_der2 = bsp5_02_der2;
+ } else if (ztype == 3) {
+ fz = bsp5_03; fz_der = bsp5_03_der; fz_der2 = bsp5_03_der2;
+ } else if (ztype == 4) {
+ fz = bsp5_04; fz_der = bsp5_04_der; fz_der2 = bsp5_04_der2;
+ } else if (ztype == 5) {
+ fz = bsp5_05; fz_der = bsp5_05_der; fz_der2 = bsp5_05_der2;
+ }
+
+ }
+ }
+
+ virtual ~global_function_xyz_bspline_()
+ { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function xyz bspline"); }
+ };
+
+
pglobal_function
global_function_bspline(const scalar_type xmin, const scalar_type xmax,
const size_type order, const size_type xtype) {
@@ -1471,6 +1651,18 @@ namespace getfem {
(xmin, xmax, ymin, ymax, order, xtype, ytype);
}
+ pglobal_function
+ global_function_bspline(const scalar_type xmin, const scalar_type xmax,
+ const scalar_type ymin, const scalar_type ymax,
+ const scalar_type zmin, const scalar_type zmax,
+ const size_type order,
+ const size_type xtype,
+ const size_type ytype,
+ const size_type ztype) {
+ return std::make_shared<global_function_xyz_bspline_>
+ (xmin, xmax, ymin, ymax, zmin, zmax, order,
+ xtype, ytype, ztype);
+ }
