[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] [getfem-commits] branch master updated: Improve impleme
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] [getfem-commits] branch master updated: Improve implementation of uniform bspline mesh_fem and add a unit test |
|
Date: |
Sun, 15 Jan 2023 18:14:03 -0500 |
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 <logari81@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);
}
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] [getfem-commits] branch master updated: Improve implementation of uniform bspline mesh_fem and add a unit test,
Konstantinos Poulios <=