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 e2bec5f1 Improve implementation of uniform bspline mesh_fem and add a unit test e2bec5f1 is described below commit e2bec5f19452db681c12ecf69c306369f3501aa9 Author: Konstantinos Poulios <logar...@gmail.com> AuthorDate: Sun Jan 15 23:50:10 2023 +0100 Improve implementation of uniform bspline mesh_fem and add a unit test --- interface/src/gf_mesh_fem.cc | 82 ++++++- interface/tests/python/Makefile.am | 2 + interface/tests/python/check_bspline_mesh_fem.py | 162 ++++++++++++ src/getfem/getfem_global_function.h | 12 +- src/getfem/getfem_mesh_fem_global_function.h | 44 +++- src/getfem_global_function.cc | 98 +++++++- src/getfem_mesh_fem_global_function.cc | 298 +++++++++++++++++++---- 7 files changed, 632 insertions(+), 66 deletions(-) diff --git a/interface/src/gf_mesh_fem.cc b/interface/src/gf_mesh_fem.cc index 6b245acf..93f33046 100644 --- a/interface/src/gf_mesh_fem.cc +++ b/interface/src/gf_mesh_fem.cc @@ -213,33 +213,95 @@ void gf_mesh_fem(getfemint::mexargs_in& m_in, if (in.remaining() && in.front().is_integer()) q_dim = in.pop().to_integer(1,256); - std::vector<getfem::pglobal_function> vfunc(size_type(in_gf.narg())); - for (size_type i = 0; i < vfunc.size(); ++i) { + std::vector<getfem::pglobal_function> vfuncs(size_type(in_gf.narg())); + for (auto &vfunc : vfuncs) { getfem::pxy_function s = to_global_function_object(in_gf.pop()); - vfunc[i] = getfem::global_function_on_level_set(*pls, s); + vfunc = getfem::global_function_on_level_set(*pls, s); } auto mfgf = std::make_shared<getfem::mesh_fem_global_function>(*mm); mfgf->set_qdim(dim_type(q_dim)); - mfgf->set_functions(vfunc); + mfgf->set_functions(vfuncs); mmf = mfgf; ); - /*@INIT MF = ('bspline', @tmesh m, @int NX, @int NY, @int order) - Create a @tmf on mesh `m`, whose basis functions are global functions + /*@INIT MF = ('bspline_uniform', @tmesh m, @int NX[, @int NY,] @int order[, @str bcX_low[, @str bcY_low[, @str bcX_high][, @str bcY_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 - that spans the entire bounding box of `m`. @*/ + (just NX in 1s) 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'. @*/ sub_command - ("bspline", 3, 4, 0, 1, + ("bspline_uniform", 3, 8, 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"); size_type NX = in.pop().to_integer(1,1000); - size_type NY = 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 order = in.pop().to_integer(3,5); + std::string bcx_low("free"); + std::string bcy_low("free"); + std::string bcx_high(""); + std::string bcy_high(""); + if (in.remaining()) bcx_low = in.pop().to_string(); + if (dim == 2 && in.remaining()) bcy_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"); + getfem::bspline_boundary bcX_low(getfem::bspline_boundary::FREE); + getfem::bspline_boundary bcY_low(getfem::bspline_boundary::FREE); + getfem::bspline_boundary bcX_high(getfem::bspline_boundary::FREE); + getfem::bspline_boundary bcY_high(getfem::bspline_boundary::FREE); + if (bcx_low == "periodic") + bcX_high = bcX_low = getfem::bspline_boundary::PERIODIC; + else if (bcx_low == "symmetry") + bcX_high = bcX_low = getfem::bspline_boundary::SYMMETRY; + else if (bcx_low != "free") + THROW_ERROR("Unknown boundary condition " << bcx_low); + + if (bcy_low == "periodic") + bcY_high = bcY_low = getfem::bspline_boundary::PERIODIC; + else if (bcy_low == "symmetry") + bcY_high = bcY_low = getfem::bspline_boundary::SYMMETRY; + else if (bcy_low != "free") + THROW_ERROR("Unknown boundary condition " << bcy_low); + + if (!bcx_high.empty()) { + if (bcx_high == "periodic") + bcX_high = getfem::bspline_boundary::PERIODIC; + else if (bcx_high == "symmetry") + bcX_high = getfem::bspline_boundary::SYMMETRY; + else if (bcx_high == "free") + bcX_high = getfem::bspline_boundary::FREE; + else + THROW_ERROR("Unknown boundary condition " << bcx_high); + } + + if (!bcy_high.empty()) { + if (bcy_high == "periodic") + bcY_high = getfem::bspline_boundary::PERIODIC; + else if (bcy_high == "symmetry") + bcY_high = getfem::bspline_boundary::SYMMETRY; + else if (bcy_high == "free") + bcY_high = getfem::bspline_boundary::FREE; + else + THROW_ERROR("Unknown boundary condition " << bcy_high); + } auto mfgf = std::make_shared<getfem::mesh_fem_global_function>(*mm); mfgf->set_qdim(1.); - define_bspline_basis_functions_for_mesh_fem(*mfgf, NX, NY, order); + if (dim == 1) + define_uniform_bspline_basis_functions_for_mesh_fem + (*mfgf, NX, order, bcX_low, bcX_high); + else + define_uniform_bspline_basis_functions_for_mesh_fem + (*mfgf, NX, NY, order, bcX_low, bcY_low, bcX_high, bcY_high); mmf = mfgf; ); diff --git a/interface/tests/python/Makefile.am b/interface/tests/python/Makefile.am index 61debcd2..753c59aa 100644 --- a/interface/tests/python/Makefile.am +++ b/interface/tests/python/Makefile.am @@ -29,6 +29,7 @@ EXTRA_DIST= \ check_global_functions.py \ check_levelset.py \ check_asm.py \ + check_bspline_mesh_fem.py \ check_secondary_domain.py \ check_mixed_mesh.py \ demo_crack.py \ @@ -71,6 +72,7 @@ TESTS = \ check_export_vtu.py \ check_global_functions.py \ check_asm.py \ + check_bspline_mesh_fem.py \ check_secondary_domain.py \ check_mixed_mesh.py \ demo_truss.py \ diff --git a/interface/tests/python/check_bspline_mesh_fem.py b/interface/tests/python/check_bspline_mesh_fem.py new file mode 100644 index 00000000..7857f13f --- /dev/null +++ b/interface/tests/python/check_bspline_mesh_fem.py @@ -0,0 +1,162 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# Python GetFEM interface +# +# Copyright (C) 2023-2023 Konstantinos Poulios +# +# This file is a part of GetFEM +# +# GetFEM is free software; you can redistribute it and/or modify it +# under the terms of the GNU Lesser General Public License as published +# by the Free Software Foundation; either version 2.1 of the License, or +# (at your option) any later version. +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +# License for more details. +# You should have received a copy of the GNU Lesser General Public License +# along with this program; if not, write to the Free Software Foundation, +# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. +# +############################################################################ +""" Test of the gf.MeshFem("bspline_uniform", ...) object. + + This program is used to check that Python-GetFEM interface, and more + generally GetFEM are working. It focuses on testing the creation of mesh_fem + with bspline basis functions in 1D and 2D. + + $Id$ +""" +import numpy as np +import getfem as gf + +NX = 8 # Number of bspline intervals in x direction +NY = 5 # Number of bspline intervals in y direction + +use_quad = True # Quadrilaterals or triangles +export_vtu = False + +if (use_quad): + m1 = gf.Mesh('import','structured','GT="GT_QK(1,1)";SIZES=[12];NOISED=0;NSUBDIV=[%d];' % (2*NX)) + m2 = gf.Mesh('import','structured','GT="GT_QK(2,1)";SIZES=[12,7];NOISED=0;NSUBDIV=[%d,%d];' % (2*NX,2*NY)) +else: + m1 = gf.Mesh('import','structured','GT="GT_PK(1,1)";SIZES=[12];NOISED=0;NSUBDIV=[%d];' % (2*NX)) + m2 = gf.Mesh('import','structured','GT="GT_PK(2,1)";SIZES=[12,7];NOISED=0;NSUBDIV=[%d,%d];' % (2*NX,2*NY)) + +mim1 = gf.MeshIm(m1, 5) +mim2 = gf.MeshIm(m2, 5) + +for order in [3, 4, 5]: + mf1_free_free = gf.MeshFem("bspline_uniform", m1, NX, order) + mf1_free_sym = gf.MeshFem("bspline_uniform", m1, NX, order, "free", "symmetry") + mf1_sym_free = gf.MeshFem("bspline_uniform", m1, NX, order, "symmetry", "free") + mf1_sym_sym = gf.MeshFem("bspline_uniform", m1, NX, order, "symmetry") + mf1_periodic = gf.MeshFem("bspline_uniform", m1, NX, order, "periodic") + + mf2_free_free_free_free = gf.MeshFem("bspline_uniform", m2, NX, NY, order) + mf2_free_sym_sym_free = gf.MeshFem("bspline_uniform", m2, NX, NY, order, "free", "symmetry", "symmetry", "free") + mf2_sym_sym_sym_sym = gf.MeshFem("bspline_uniform", m2, NX, NY, order, "symmetry", "symmetry") + mf2_periodic_free_periodic_free = gf.MeshFem("bspline_uniform", m2, NX, NY, order, "periodic", "free") + mf2_sym_periodic_free_periodic = gf.MeshFem("bspline_uniform", m2, NX, NY, order, "symmetry", "periodic", "free", "periodic") + + print("order %d" % order) + print("mf1_free_free.nbdof()=", mf1_free_free.nbdof()) + print("mf1_free_sym.nbdof()=", mf1_free_sym.nbdof()) + print("mf1_sym_free.nbdof()=", mf1_sym_free.nbdof()) + print("mf1_sym_sym.nbdof()=", mf1_sym_sym.nbdof()) + print("mf1_periodic.nbdof()=", mf1_periodic.nbdof()) + if export_vtu: + for dof in range(mf1_free_free.nbdof()): + mf1_free_free.export_to_vtu(f"basis_funcs_order{order}_free_free_{dof}.vtu", + (np.arange(mf1_free_free.nbdof())==dof).astype(float), "basis function") + for dof in range(mf1_free_sym.nbdof()): + mf1_free_sym.export_to_vtu(f"basis_funcs_order{order}_free_sym_{dof}.vtu", + (np.arange(mf1_free_sym.nbdof())==dof).astype(float), "basis function") + for dof in range(mf1_sym_free.nbdof()): + mf1_sym_free.export_to_vtu(f"basis_funcs_order{order}_sym_free_{dof}.vtu", + (np.arange(mf1_sym_free.nbdof())==dof).astype(float), "basis function") + for dof in range(mf1_sym_sym.nbdof()): + mf1_sym_sym.export_to_vtu(f"basis_funcs_order{order}_sym_sym_{dof}.vtu", + (np.arange(mf1_sym_sym.nbdof())==dof).astype(float), "basis function") + for dof in range(mf1_periodic.nbdof()): + mf1_periodic.export_to_vtu(f"basis_funcs_order{order}_periodic_{dof}.vtu", + (np.arange(mf1_periodic.nbdof())==dof).astype(float), "basis function") + + print("mf2_free_free_free_free.nbdof()=", mf2_free_free_free_free.nbdof()) + print("mf2_free_sym_sym_free.nbdof()=", mf2_free_sym_sym_free.nbdof()) + print("mf2_sym_sym_sym_sym.nbdof()=", mf2_sym_sym_sym_sym.nbdof()) + print("mf2_periodic_free_periodic_free.nbdof()=", mf2_periodic_free_periodic_free.nbdof()) + print("mf2_sym_periodic_free_periodic.nbdof()=", mf2_sym_periodic_free_periodic.nbdof()) + if export_vtu: + for dof in range(mf2_free_free_free_free.nbdof()): + mf2_free_free_free_free.export_to_vtu(f"basis_funcs_order{order}_free_free_free_free_{dof}.vtu", + (np.arange(mf2_free_free_free_free.nbdof())==dof).astype(float), "basis function") + for dof in range(mf2_free_sym_sym_free.nbdof()): + mf2_free_sym_sym_free.export_to_vtu(f"basis_funcs_order{order}_free_sym_sym_free_{dof}.vtu", + (np.arange(mf2_free_sym_sym_free.nbdof())==dof).astype(float), "basis function") + for dof in range(mf2_sym_sym_sym_sym.nbdof()): + mf2_sym_sym_sym_sym.export_to_vtu(f"basis_funcs_order{order}_sym_sym_sym_sym_{dof}.vtu", + (np.arange(mf2_sym_sym_sym_sym.nbdof())==dof).astype(float), "basis function") + for dof in range(mf2_periodic_free_periodic_free.nbdof()): + mf2_periodic_free_periodic_free.export_to_vtu(f"basis_funcs_order{order}_periodic_free_periodic_free_{dof}.vtu", + (np.arange(mf2_periodic_free_periodic_free.nbdof())==dof).astype(float), "basis function") + for dof in range(mf2_sym_periodic_free_periodic.nbdof()): + mf2_sym_periodic_free_periodic.export_to_vtu(f"basis_funcs_order{order}_sym_periodic_free_periodic_{dof}.vtu", + (np.arange(mf2_sym_periodic_free_periodic.nbdof())==dof).astype(float), "basis function") + + assert(mf1_free_free.nbdof() == (10,11,12)[order-3]) + assert(mf1_free_sym.nbdof() == (9,10,10)[order-3]) + assert(mf1_sym_free.nbdof() == (9,10,10)[order-3]) + assert(mf1_sym_sym.nbdof() == (8,9,8)[order-3]) + assert(mf1_periodic.nbdof() == (8,8,8)[order-3]) + + assert(mf2_free_free_free_free.nbdof() == (70,88,108)[order-3]); + assert(mf2_free_sym_sym_free.nbdof() == (54,70,70)[order-3]); + assert(mf2_sym_sym_sym_sym.nbdof() == (40,54,40)[order-3]); + assert(mf2_periodic_free_periodic_free.nbdof() == (56,64,72)[order-3]); + assert(mf2_sym_periodic_free_periodic.nbdof() == (45,50,50)[order-3]); + + # partition of unity check + L2error = gf.asm_generic(mim1, 0, "Norm_sqr(1-t)", -1, + "t", False, mf1_free_free, np.ones(mf1_free_free.nbdof())) + assert(L2error < 1e-16); + print("mf1_free_free partition of unity error:", L2error) + L2error = gf.asm_generic(mim1, 0, "Norm_sqr(1-t)", -1, + "t", False, mf1_free_sym, np.ones(mf1_free_sym.nbdof())) + assert(L2error < 1e-16); + print("mf1_free_sym partition of unity error:", L2error) + L2error = gf.asm_generic(mim1, 0, "Norm_sqr(1-t)", -1, + "t", False, mf1_sym_free, np.ones(mf1_sym_free.nbdof())) + assert(L2error < 1e-16); + print("mf1_sym_free partition of unity error:", L2error) + L2error = gf.asm_generic(mim1, 0, "Norm_sqr(1-t)", -1, + "t", False, mf1_sym_sym, np.ones(mf1_sym_sym.nbdof())) + assert(L2error < 1e-16); + print("mf1_sym_sym partition of unity error:", L2error) + L2error = gf.asm_generic(mim1, 0, "Norm_sqr(1-t)", -1, + "t", False, mf1_periodic, np.ones(mf1_periodic.nbdof())) + assert(L2error < 1e-16); + print("mf1_periodic partition of unity error:", L2error) + + L2error = gf.asm_generic(mim2, 0, "Norm_sqr(1-t)", -1, + "t", False, mf2_free_free_free_free, np.ones(mf2_free_free_free_free.nbdof())) + assert(L2error < 1e-16); + print("mf2_free_free_free_free partition of unity error:", L2error) + L2error = gf.asm_generic(mim2, 0, "Norm_sqr(1-t)", -1, + "t", False, mf2_free_sym_sym_free, np.ones(mf2_free_sym_sym_free.nbdof())) + assert(L2error < 1e-16); + print("mf2_free_sym_sym_free partition of unity error:", L2error) + L2error = gf.asm_generic(mim2, 0, "Norm_sqr(1-t)", -1, + "t", False, mf2_sym_sym_sym_sym, np.ones(mf2_sym_sym_sym_sym.nbdof())) + assert(L2error < 1e-16); + print("mf2_sym_sym_sym_sym partition of unity error:", L2error) + L2error = gf.asm_generic(mim2, 0, "Norm_sqr(1-t)", -1, + "t", False, mf2_periodic_free_periodic_free, np.ones(mf2_periodic_free_periodic_free.nbdof())) + assert(L2error < 1e-16); + print("mf2_periodic_free_periodic_free partition of unity error:", L2error) + L2error = gf.asm_generic(mim2, 0, "Norm_sqr(1-t)", -1, + "t", False, mf2_sym_periodic_free_periodic, np.ones(mf2_sym_periodic_free_periodic.nbdof())) + assert(L2error < 1e-16); + print("mf2_sym_periodic_free_periodic partition of unity error:", L2error) + +exit(0) diff --git a/src/getfem/getfem_global_function.h b/src/getfem/getfem_global_function.h index 75fe122a..b07434e9 100644 --- a/src/getfem/getfem_global_function.h +++ b/src/getfem/getfem_global_function.h @@ -321,10 +321,14 @@ namespace getfem { 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); + global_function_bspline(const scalar_type xmin, const scalar_type xmax, + const size_type order, const size_type xtype); + + pglobal_function + global_function_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); } /* 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 65982de4..a85fb9c6 100644 --- a/src/getfem/getfem_mesh_fem_global_function.h +++ b/src/getfem/getfem_mesh_fem_global_function.h @@ -62,18 +62,54 @@ namespace getfem { virtual ~mesh_fem_global_function() { clear(); } }; + enum class bspline_boundary { FREE=0, PERIODIC=1, SYMMETRY=2}; + + /** This function will generate bspline basis functions on NX uniform + elements along a line. The dimensions of the domain correspond to + the bounding interval of the 1d 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 order, + bspline_boundary bcX_low=bspline_boundary::FREE, + bspline_boundary bcX_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 order, + bspline_boundary bcX_low, const mesh_im &mim=dummy_mesh_im()) { + define_uniform_bspline_basis_functions_for_mesh_fem + (mf, NX, order, bcX_low, bcX_low, mim); + } + + /** 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. + box of the 2d 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. + determine the support of he basis functions more precisely. */ - void define_bspline_basis_functions_for_mesh_fem + void define_uniform_bspline_basis_functions_for_mesh_fem (mesh_fem_global_function &mf, size_type NX, size_type NY, size_type order, + bspline_boundary bcX_low=bspline_boundary::FREE, + bspline_boundary bcY_low=bspline_boundary::FREE, + bspline_boundary bcX_high=bspline_boundary::FREE, + bspline_boundary bcY_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 order, + bspline_boundary bcX_low, bspline_boundary bcY_low, + const mesh_im &mim=dummy_mesh_im()) { + define_uniform_bspline_basis_functions_for_mesh_fem + (mf, NX, NY, order, bcX_low, bcY_low, bcX_low, bcY_low, mim); + } + } /* end of namespace getfem. */ #endif diff --git a/src/getfem_global_function.cc b/src/getfem_global_function.cc index 400d4dd0..0ed103f4 100644 --- a/src/getfem_global_function.cc +++ b/src/getfem_global_function.cc @@ -1244,6 +1244,89 @@ namespace getfem { + class global_function_x_bspline_ + : public global_function_simple, public context_dependencies { + scalar_type xmin, xmax, xscale; + std::function<scalar_type(scalar_type)> fx, fx_der, fx_der2; + public: + void update_from_context() const {} + + virtual scalar_type val(const base_node &pt) const + { + return fx(xscale*(pt[0]-xmin)); + } + virtual void grad(const base_node &pt, base_small_vector &g) const + { + scalar_type dx = xscale*(pt[0]-xmin); + g.resize(dim_); + g[0] = xscale * fx_der(dx); + } + virtual void hess(const base_node &pt, base_matrix &h) const + { + scalar_type dx = xscale*(pt[0]-xmin); + h.resize(dim_, dim_); + gmm::clear(h); + h(0,0) = xscale*xscale * fx_der2(dx); + } + + virtual bool is_in_support(const base_node &pt) const { + scalar_type dx = pt[0]-(xmin+xmax)/2; + return (gmm::abs(dx)+1e-9 < gmm::abs(xmax-xmin)/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); + bmax[0] = std::max(xmin,xmax); + } + + global_function_x_bspline_(const scalar_type &xmin_, const scalar_type &xmax_, + const size_type &order, const size_type &xtype) + : global_function_simple(1), xmin(xmin_), xmax(xmax_), + xscale(scalar_type(xtype)/(xmax-xmin)) + { + GMM_ASSERT1(order >= 3 && order <= 5, "Only orders 3 to 5 are supported"); + GMM_ASSERT1(xtype >= 1 && xtype <= 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; + } + } 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; + } + } 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; + } + } + } + + virtual ~global_function_x_bspline_() + { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function x bspline"); } + }; + + + class global_function_xy_bspline_ : public global_function_simple, public context_dependencies { scalar_type xmin, ymin, xmax, ymax, xscale, yscale; @@ -1372,10 +1455,17 @@ namespace getfem { 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) { + global_function_bspline(const scalar_type xmin, const scalar_type xmax, + const size_type order, const size_type xtype) { + return std::make_shared<global_function_x_bspline_> + (xmin, xmax, order, xtype); + } + + pglobal_function + global_function_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) { return std::make_shared<global_function_xy_bspline_> (xmin, xmax, ymin, ymax, order, xtype, ytype); } diff --git a/src/getfem_mesh_fem_global_function.cc b/src/getfem_mesh_fem_global_function.cc index d595d0ae..6b05af2f 100644 --- a/src/getfem_mesh_fem_global_function.cc +++ b/src/getfem_mesh_fem_global_function.cc @@ -52,59 +52,269 @@ namespace getfem { } - void define_bspline_basis_functions_for_mesh_fem +// examples of bspline basis functions assigned to 8 elements in 1D + +// n=8,k=3, free-free --> n-k+1 + 2*(k-1) = n+k-1 = 8+3-1 = 10 +// 1 2 3 4 5 6 7 8 | +//1 * | +//2 * * | +//3 * * * | +//4 * * * | +//5 * * * | +//6 * * * | +//7 * * * | +//8 * * * | +//9 * * | +//10 * | + +// n=8,k=4, free-free --> n-k+1 + 2*(k-1) = n+k-1 = 8+4-1 = 11 +// 1 2 3 4 5 6 7 8 +//1 * | +//2 * * | +//3 * * * | +//4 * * * * | +//5 * * * * | +//6 * * * * | +//7 * * * * | +//8 * * * * | +//9 * * * | +//10 * * | +//11 * | + +// n=8,k=3, periodic --> n-k+1 + k-1 = n +// 1 2 3 4 5 6 7 8 +//1 * * * | +//2 * * * | +//3 * * * | +//4 * * * | +//5 * * * | +//6 * * * | +//7 * * * | +//8 * * * | + +// n=8,k=4, periodic +// 1 2 3 4 5 6 7 8 +//1 * * * * | +//2 * * * * | +//3 * * * * | +//4 * * * * | +//5 * * * * | +//6 * * * * | +//7 * * * * | +//8 * * * * | + +// n=8,k=3, symmetry-symmetry --> n-k+1 + 2*floor(k/2) = 6 + 2 = 8 +// 1 2 3 4 5 6 7 8 +//1 + * | +//2 * * * | +//3 * * * | +//4 * * * | +//5 * * * | +//6 * * * | +//7 * * * | +//8 * + | + +// n=8,k=4, symmetry-symmetry --> n-k+1 + 2*floor(k/2) = 5 + 4 = 9 +// 1 2 3 4 5 6 7 8 | +//1 * * | +//2 + * * | +//3 * * * * | +//4 * * * * | +//5 * * * * | +//6 * * * * | +//7 * * * * | +//8 * * + | +//9 * * | + +// n=8,k=5, symmetry-symmetry --> n-k+1 + 2*floor(k/2) = 4 + 4 = 8 +// 1 2 3 4 5 6 7 8 | +//1 + + * | +//2 + * * * | +//3 * * * * * | +//4 * * * * * | +//5 * * * * * | +//6 * * * * * | +//7 * * * + | +//8 * + + | + +// n=8,k=6, symmetry-symmetry --> n-k+1 + 2*floor(k/2) = 3 + 6 = 9 +// 1 2 3 4 5 6 7 8 | +//1 * * * | +//2 + + * * | +//3 + * * * * | +//4 * * * * * * | +//5 * * * * * * | +//6 * * * * * * | +//7 * * * * + | +//8 * * + + | +//9 * * * | + +// n=8,k=3, free-symmetry --> n-k+1 + k-1 + floor(k/2) = 6 + 2 + 1 = 9 +// 1 2 3 4 5 6 7 8 | +//1 * | +//2 + * | +//3 * * * | +//4 * * * | +//5 * * * | +//6 * * * | +//7 * * * | +//8 * * * | +//9 * + | + + void params_for_uniform_1d_bspline_basis_functions + (scalar_type x0, scalar_type x1, size_type N, size_type order, + bspline_boundary bc_low, bspline_boundary bc_high, + std::vector<scalar_type> &xmin, std::vector<scalar_type> &xmax, + std::vector<scalar_type> &xshift, std::vector<size_type> &xtype) { + + if (bc_low == bspline_boundary::PERIODIC || + bc_high == bspline_boundary::PERIODIC) + GMM_ASSERT1(bc_low == bc_high, + "Periodic BC has to be assigned to both matching sides"); + const scalar_type dx = (x1-x0)/scalar_type(N); + size_type n_low, n_mid, n_high; + n_low = (bc_low == bspline_boundary::PERIODIC) ? 0 : + (bc_low == bspline_boundary::SYMMETRY ? order/2 : + /* FREE */ order-1); + n_high = (bc_high == bspline_boundary::PERIODIC) ? order-1 : + (bc_high == bspline_boundary::SYMMETRY ? order/2 : + /* FREE */ order-1); + n_mid = N - order + 1; + size_type n = n_low + n_mid + n_high; // number of basis functions + + xmin.resize(n); + xmax.resize(n); + xshift.resize(n); + xtype.resize(n); + for (size_type i=0; i < n; ++i) { + xshift[i] = 0.; + if (bc_low == bspline_boundary::FREE && i < n_low) { + xtype[i] = i+1; + xmin[i] = x0; + xmax[i] = xmin[i] + scalar_type(xtype[i])*dx; + } else if (bc_high == bspline_boundary::FREE && i >= n_low+n_mid) { + xtype[i] = n-i; // safe unsigned + xmin[i] = x1; + xmax[i] = xmin[i] - scalar_type(xtype[i])*dx; // yes, xmax < xmin + } else if (bc_low == bspline_boundary::SYMMETRY && i < n_low) { + xtype[i] = order; + xmin[i] = x0 - scalar_type(n_low-i)*dx; + xmax[i] = xmin[i] + scalar_type(xtype[i])*dx; + xshift[i] = -(xmin[i]+xmax[i]-2*x0); // this is 0 for already symmetric basis functions + } else if (bc_high == bspline_boundary::SYMMETRY && i >= n_low+n_mid) { + xtype[i] = order; + xmin[i] = x0 + scalar_type(i-n_low)*dx; // safe unsigned + xmax[i] = xmin[i] + scalar_type(xtype[i])*dx; + xshift[i] = 2*x1-xmin[i]-xmax[i]; // this is 0 for already symmetric basis functions + } else { // mid functions for periodic, free-free or free-symmetry or symmetry-free + GMM_ASSERT1(i >= n_low, "Internal error"); + xtype[i] = order; + xmin[i] = x0 + scalar_type(i-n_low)*dx; // safe unsigned + xmax[i] = xmin[i] + scalar_type(xtype[i])*dx; + } +//if (order==5) // && bc_low == bspline_boundary::SYMMETRY && bc_high == bspline_boundary::FREE) +//std::cout<<i<<":"<<xmin[i]<<","<<xmax[i]<<std::endl; + + if (bc_low == bspline_boundary::PERIODIC && xmax[i] > x1) + xshift[i] = -(x1-x0); // this will apply to the last order-1 functions + } + } + + void define_uniform_bspline_basis_functions_for_mesh_fem + (mesh_fem_global_function &mf, size_type NX, size_type order, + bspline_boundary bcX_low, bspline_boundary bcX_high, + const mesh_im &mim) { + + GMM_ASSERT1(mf.linked_mesh().dim() == 1, + "This function expects a mesh_fem defined in 1d"); + + base_node Pmin, Pmax; + mf.linked_mesh().bounding_box(Pmin, Pmax); + const scalar_type x0=Pmin[0], x1=Pmax[0]; + + std::vector<scalar_type> xmin, xmax, xshift; + std::vector<size_type> xtype; + params_for_uniform_1d_bspline_basis_functions + (x0, x1, NX, order, bcX_low, bcX_high, // input + xmin, xmax, xshift, xtype); // output + + std::vector<pglobal_function> funcs(0); + for (size_type i=0; i < xtype.size(); ++i) { + if (gmm::abs(xshift[i]) < 1e-10) + funcs.push_back(global_function_bspline + (xmin[i], xmax[i], order, xtype[i])); + else { + std::vector<pglobal_function> sum; + sum.push_back(global_function_bspline + (xmin[i], xmax[i], order, xtype[i])); + sum.push_back(global_function_bspline + (xmin[i]+xshift[i], xmax[i]+xshift[i], + order, xtype[i])); + funcs.push_back(std::make_shared<getfem::global_function_sum>(sum)); + } + } + mf.set_functions(funcs, mim); + } + + void define_uniform_bspline_basis_functions_for_mesh_fem (mesh_fem_global_function &mf, size_type NX, size_type NY, size_type order, + bspline_boundary bcX_low, bspline_boundary bcY_low, + bspline_boundary bcX_high, bspline_boundary bcY_high, const mesh_im &mim) { + GMM_ASSERT1(mf.linked_mesh().dim() == 2, + "This function expects a mesh_fem defined in 2d"); + 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; + const scalar_type x0=Pmin[0], x1=Pmax[0], + y0=Pmin[1], y1=Pmax[1]; + + std::vector<scalar_type> xmin, xmax, xshift; + std::vector<size_type> xtype; + params_for_uniform_1d_bspline_basis_functions + (x0, x1, NX, order, bcX_low, bcX_high, // input + xmin, xmax, xshift, xtype); // output + std::vector<scalar_type> ymin, ymax, yshift; + std::vector<size_type> ytype; + params_for_uniform_1d_bspline_basis_functions + (y0, y1, NY, order, bcY_low, bcY_high, // input + ymin, ymax, yshift, ytype); // output + + std::vector<pglobal_function> funcs(0); + for (size_type i=0; i < xtype.size(); ++i) { + for (size_type j=0; j < ytype.size(); ++j) { + if (gmm::abs(xshift[i]) < 1e-10 && + gmm::abs(yshift[j]) < 1e-10) + funcs.push_back(global_function_bspline + (xmin[i], xmax[i], ymin[j], ymax[j], + order, xtype[i], ytype[j])); + else { + std::vector<pglobal_function> sum; + sum.push_back(global_function_bspline + (xmin[i], xmax[i], ymin[j], ymax[j], + order, xtype[i], ytype[j])); + if (gmm::abs(xshift[i]) >= 1e-10) + sum.push_back(global_function_bspline + (xmin[i]+xshift[i], xmax[i]+xshift[i], + ymin[j], ymax[j], + order, xtype[i], ytype[j])); + if (gmm::abs(yshift[j]) >= 1e-10) { + sum.push_back(global_function_bspline + (xmin[i], xmax[i], + ymin[j]+yshift[j], ymax[j]+yshift[j], + order, xtype[i], ytype[j])); + if (gmm::abs(xshift[i]) >= 1e-10) + sum.push_back(global_function_bspline + (xmin[i]+xshift[i], xmax[i]+xshift[i], + ymin[j]+yshift[j], ymax[j]+yshift[j], + order, xtype[i], ytype[j])); + } + funcs.push_back(std::make_shared<getfem::global_function_sum>(sum)); } - funcs[i*(NY+order-1)+j] = global_function_bspline - (xmin, xmax, ymin, ymax, order, xtype, ytype); } } - mf.set_functions(funcs, mim); }