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); + }