GetFEM  5.4.3
getfem_mesh_fem_global_function.cc
1 /*===========================================================================
2 
3  Copyright (C) 2004-2022 Yves Renard
4  Copyright (C) 2016-2022 Konstantinos Poulios
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21 ===========================================================================*/
22 
23 #include <getfem/getfem_mesh_fem_global_function.h>
24 
25 namespace getfem {
26 
27 
28  void mesh_fem_global_function::set_functions
29  (const std::vector<pglobal_function>& funcs, const mesh_im &mim) {
30  GMM_ASSERT1(linked_mesh_ != 0, "Mesh fem need to be initialized with"
31  " a mesh first.");
32  clear();
33  if (&mim == &dummy_mesh_im())
34  fem_ = getfem::new_fem_global_function(funcs, linked_mesh());
35  else {
36  GMM_ASSERT1(&(mim.linked_mesh()) == linked_mesh_,
37  "The provided mesh_im has to be linked to the same mesh"
38  " as this mesh_fem.");
39  fem_ = getfem::new_fem_global_function(funcs, mim);
40  }
41  set_finite_element(fem_);
42  }
43 
44  // mesh_fem_global_function::size_type memsize() const;
45 
46  void mesh_fem_global_function::clear() {
47  mesh_fem::clear();
48  if (fem_.get()) {
49  getfem::del_fem_global_function(fem_);
50  fem_.reset();
51  }
52  }
53 
54 
55 // examples of bspline basis functions assigned to 8 elements in 1D
56 
57 // n=8,k=3, free-free --> n-k+1 + 2*(k-1) = n+k-1 = 8+3-1 = 10
58 // 1 2 3 4 5 6 7 8 |
59 //1 * |
60 //2 * * |
61 //3 * * * |
62 //4 * * * |
63 //5 * * * |
64 //6 * * * |
65 //7 * * * |
66 //8 * * * |
67 //9 * * |
68 //10 * |
69 
70 // n=8,k=4, free-free --> n-k+1 + 2*(k-1) = n+k-1 = 8+4-1 = 11
71 // 1 2 3 4 5 6 7 8
72 //1 * |
73 //2 * * |
74 //3 * * * |
75 //4 * * * * |
76 //5 * * * * |
77 //6 * * * * |
78 //7 * * * * |
79 //8 * * * * |
80 //9 * * * |
81 //10 * * |
82 //11 * |
83 
84 // n=8,k=3, periodic --> n-k+1 + k-1 = n
85 // 1 2 3 4 5 6 7 8
86 //1 * * * |
87 //2 * * * |
88 //3 * * * |
89 //4 * * * |
90 //5 * * * |
91 //6 * * * |
92 //7 * * * |
93 //8 * * * |
94 
95 // n=8,k=4, periodic
96 // 1 2 3 4 5 6 7 8
97 //1 * * * * |
98 //2 * * * * |
99 //3 * * * * |
100 //4 * * * * |
101 //5 * * * * |
102 //6 * * * * |
103 //7 * * * * |
104 //8 * * * * |
105 
106 // n=8,k=3, symmetry-symmetry --> n-k+1 + 2*floor(k/2) = 6 + 2 = 8
107 // 1 2 3 4 5 6 7 8
108 //1 + * |
109 //2 * * * |
110 //3 * * * |
111 //4 * * * |
112 //5 * * * |
113 //6 * * * |
114 //7 * * * |
115 //8 * + |
116 
117 // n=8,k=4, symmetry-symmetry --> n-k+1 + 2*floor(k/2) = 5 + 4 = 9
118 // 1 2 3 4 5 6 7 8 |
119 //1 * * |
120 //2 + * * |
121 //3 * * * * |
122 //4 * * * * |
123 //5 * * * * |
124 //6 * * * * |
125 //7 * * * * |
126 //8 * * + |
127 //9 * * |
128 
129 // n=8,k=5, symmetry-symmetry --> n-k+1 + 2*floor(k/2) = 4 + 4 = 8
130 // 1 2 3 4 5 6 7 8 |
131 //1 + + * |
132 //2 + * * * |
133 //3 * * * * * |
134 //4 * * * * * |
135 //5 * * * * * |
136 //6 * * * * * |
137 //7 * * * + |
138 //8 * + + |
139 
140 // n=8,k=6, symmetry-symmetry --> n-k+1 + 2*floor(k/2) = 3 + 6 = 9
141 // 1 2 3 4 5 6 7 8 |
142 //1 * * * |
143 //2 + + * * |
144 //3 + * * * * |
145 //4 * * * * * * |
146 //5 * * * * * * |
147 //6 * * * * * * |
148 //7 * * * * + |
149 //8 * * + + |
150 //9 * * * |
151 
152 // n=8,k=3, free-symmetry --> n-k+1 + k-1 + floor(k/2) = 6 + 2 + 1 = 9
153 // 1 2 3 4 5 6 7 8 |
154 //1 * |
155 //2 + * |
156 //3 * * * |
157 //4 * * * |
158 //5 * * * |
159 //6 * * * |
160 //7 * * * |
161 //8 * * * |
162 //9 * + |
163 
164  void params_for_uniform_1d_bspline_basis_functions
165  (scalar_type x0, scalar_type x1, size_type N, size_type order,
166  bspline_boundary bc_low, bspline_boundary bc_high,
167  std::vector<scalar_type> &xmin, std::vector<scalar_type> &xmax,
168  std::vector<scalar_type> &xshift, std::vector<size_type> &xtype) {
169 
170  if (bc_low == bspline_boundary::PERIODIC ||
171  bc_high == bspline_boundary::PERIODIC)
172  GMM_ASSERT1(bc_low == bc_high,
173  "Periodic BC has to be assigned to both matching sides");
174  const scalar_type dx = (x1-x0)/scalar_type(N);
175  size_type n_low, n_mid, n_high;
176  n_low = (bc_low == bspline_boundary::PERIODIC) ? 0 :
177  (bc_low == bspline_boundary::SYMMETRY ? order/2 :
178  /* FREE */ order-1);
179  n_high = (bc_high == bspline_boundary::PERIODIC) ? order-1 :
180  (bc_high == bspline_boundary::SYMMETRY ? order/2 :
181  /* FREE */ order-1);
182  n_mid = N - order + 1;
183  size_type n = n_low + n_mid + n_high; // number of basis functions
184 
185  xmin.resize(n);
186  xmax.resize(n);
187  xshift.resize(n);
188  xtype.resize(n);
189  for (size_type i=0; i < n; ++i) {
190  xshift[i] = 0.;
191  if (bc_low == bspline_boundary::FREE && i < n_low) {
192  xtype[i] = i+1;
193  xmin[i] = x0;
194  xmax[i] = xmin[i] + scalar_type(xtype[i])*dx;
195  } else if (bc_high == bspline_boundary::FREE && i >= n_low+n_mid) {
196  xtype[i] = n-i; // safe unsigned
197  xmin[i] = x1;
198  xmax[i] = xmin[i] - scalar_type(xtype[i])*dx; // yes, xmax < xmin
199  } else if (bc_low == bspline_boundary::SYMMETRY && i < n_low) {
200  xtype[i] = order;
201  xmin[i] = x0 - scalar_type(n_low-i)*dx;
202  xmax[i] = xmin[i] + scalar_type(xtype[i])*dx;
203  xshift[i] = -(xmin[i]+xmax[i]-2*x0); // this is 0 for already symmetric basis functions
204  } else if (bc_high == bspline_boundary::SYMMETRY && i >= n_low+n_mid) {
205  xtype[i] = order;
206  xmin[i] = x0 + scalar_type(i-n_low)*dx; // safe unsigned
207  xmax[i] = xmin[i] + scalar_type(xtype[i])*dx;
208  xshift[i] = 2*x1-xmin[i]-xmax[i]; // this is 0 for already symmetric basis functions
209  } else { // mid functions for periodic, free-free or free-symmetry or symmetry-free
210  GMM_ASSERT1(i >= n_low, "Internal error");
211  xtype[i] = order;
212  xmin[i] = x0 + scalar_type(i-n_low)*dx; // safe unsigned
213  xmax[i] = xmin[i] + scalar_type(xtype[i])*dx;
214  }
215 //if (order==5) // && bc_low == bspline_boundary::SYMMETRY && bc_high == bspline_boundary::FREE)
216 //std::cout<<i<<":"<<xmin[i]<<","<<xmax[i]<<std::endl;
217 
218  if (bc_low == bspline_boundary::PERIODIC && xmax[i] > x1)
219  xshift[i] = -(x1-x0); // this will apply to the last order-1 functions
220  }
221  }
222 
223  void define_uniform_bspline_basis_functions_for_mesh_fem
224  (mesh_fem_global_function &mf, size_type NX, size_type order,
225  bspline_boundary bcX_low, bspline_boundary bcX_high,
226  const mesh_im &mim) {
227 
228  GMM_ASSERT1(mf.linked_mesh().dim() == 1,
229  "This function expects a mesh_fem defined in 1d");
230 
231  base_node Pmin, Pmax;
232  mf.linked_mesh().bounding_box(Pmin, Pmax);
233  const scalar_type x0=Pmin[0], x1=Pmax[0];
234 
235  std::vector<scalar_type> xmin, xmax, xshift;
236  std::vector<size_type> xtype;
237  params_for_uniform_1d_bspline_basis_functions
238  (x0, x1, NX, order, bcX_low, bcX_high, // input
239  xmin, xmax, xshift, xtype); // output
240 
241  std::vector<pglobal_function> funcs(0);
242  for (size_type i=0; i < xtype.size(); ++i) {
243  if (gmm::abs(xshift[i]) < 1e-10)
244  funcs.push_back(global_function_bspline
245  (xmin[i], xmax[i], order, xtype[i]));
246  else {
247  std::vector<pglobal_function> sum;
248  sum.push_back(global_function_bspline
249  (xmin[i], xmax[i], order, xtype[i]));
250  sum.push_back(global_function_bspline
251  (xmin[i]+xshift[i], xmax[i]+xshift[i],
252  order, xtype[i]));
253  funcs.push_back(std::make_shared<getfem::global_function_sum>(sum));
254  }
255  }
256  mf.set_functions(funcs, mim);
257  }
258 
259  void define_uniform_bspline_basis_functions_for_mesh_fem
260  (mesh_fem_global_function &mf,
261  size_type NX, size_type NY, size_type order,
262  bspline_boundary bcX_low, bspline_boundary bcY_low,
263  bspline_boundary bcX_high, bspline_boundary bcY_high,
264  const mesh_im &mim) {
265 
266  GMM_ASSERT1(mf.linked_mesh().dim() == 2,
267  "This function expects a mesh_fem defined in 2d");
268 
269  base_node Pmin, Pmax;
270  mf.linked_mesh().bounding_box(Pmin, Pmax);
271  const scalar_type x0=Pmin[0], x1=Pmax[0],
272  y0=Pmin[1], y1=Pmax[1];
273 
274  std::vector<scalar_type> xmin, xmax, xshift;
275  std::vector<size_type> xtype;
276  params_for_uniform_1d_bspline_basis_functions
277  (x0, x1, NX, order, bcX_low, bcX_high, // input
278  xmin, xmax, xshift, xtype); // output
279  std::vector<scalar_type> ymin, ymax, yshift;
280  std::vector<size_type> ytype;
281  params_for_uniform_1d_bspline_basis_functions
282  (y0, y1, NY, order, bcY_low, bcY_high, // input
283  ymin, ymax, yshift, ytype); // output
284 
285  std::vector<pglobal_function> funcs(0);
286  for (size_type i=0; i < xtype.size(); ++i) {
287  for (size_type j=0; j < ytype.size(); ++j) {
288  if (gmm::abs(xshift[i]) < 1e-10 &&
289  gmm::abs(yshift[j]) < 1e-10)
290  funcs.push_back(global_function_bspline
291  (xmin[i], xmax[i], ymin[j], ymax[j],
292  order, xtype[i], ytype[j]));
293  else {
294  std::vector<pglobal_function> sum;
295  sum.push_back(global_function_bspline
296  (xmin[i], xmax[i], ymin[j], ymax[j],
297  order, xtype[i], ytype[j]));
298  if (gmm::abs(xshift[i]) >= 1e-10)
299  sum.push_back(global_function_bspline
300  (xmin[i]+xshift[i], xmax[i]+xshift[i],
301  ymin[j], ymax[j],
302  order, xtype[i], ytype[j]));
303  if (gmm::abs(yshift[j]) >= 1e-10) {
304  sum.push_back(global_function_bspline
305  (xmin[i], xmax[i],
306  ymin[j]+yshift[j], ymax[j]+yshift[j],
307  order, xtype[i], ytype[j]));
308  if (gmm::abs(xshift[i]) >= 1e-10)
309  sum.push_back(global_function_bspline
310  (xmin[i]+xshift[i], xmax[i]+xshift[i],
311  ymin[j]+yshift[j], ymax[j]+yshift[j],
312  order, xtype[i], ytype[j]));
313  }
314  funcs.push_back(std::make_shared<getfem::global_function_sum>(sum));
315  }
316  }
317  }
318  mf.set_functions(funcs, mim);
319  }
320 
321  void define_uniform_bspline_basis_functions_for_mesh_fem
322  (mesh_fem_global_function &mf,
323  size_type NX, size_type NY, size_type NZ, size_type order,
324  bspline_boundary bcX_low,
325  bspline_boundary bcY_low,
326  bspline_boundary bcZ_low,
327  bspline_boundary bcX_high,
328  bspline_boundary bcY_high,
329  bspline_boundary bcZ_high, const mesh_im &mim) {
330 
331  GMM_ASSERT1(mf.linked_mesh().dim() == 3,
332  "This function expects a mesh_fem defined in 3d");
333 
334  base_node Pmin, Pmax;
335  mf.linked_mesh().bounding_box(Pmin, Pmax);
336  const scalar_type x0=Pmin[0], x1=Pmax[0],
337  y0=Pmin[1], y1=Pmax[1],
338  z0=Pmin[2], z1=Pmax[2];
339 
340  std::vector<scalar_type> xmin, xmax, xshift;
341  std::vector<size_type> xtype;
342  params_for_uniform_1d_bspline_basis_functions
343  (x0, x1, NX, order, bcX_low, bcX_high, // input
344  xmin, xmax, xshift, xtype); // output
345  std::vector<scalar_type> ymin, ymax, yshift;
346  std::vector<size_type> ytype;
347  params_for_uniform_1d_bspline_basis_functions
348  (y0, y1, NY, order, bcY_low, bcY_high, // input
349  ymin, ymax, yshift, ytype); // output
350  std::vector<scalar_type> zmin, zmax, zshift;
351  std::vector<size_type> ztype;
352  params_for_uniform_1d_bspline_basis_functions
353  (z0, z1, NZ, order, bcZ_low, bcZ_high, // input
354  zmin, zmax, zshift, ztype); // output
355 
356  std::vector<pglobal_function> funcs(0);
357  for (size_type i=0; i < xtype.size(); ++i) {
358  for (size_type j=0; j < ytype.size(); ++j) {
359  for (size_type k=0; k < ztype.size(); ++k) {
360 
361  bool has_xshift = gmm::abs(xshift[i]) >= 1e-10;
362  bool has_yshift = gmm::abs(yshift[j]) >= 1e-10;
363  bool has_zshift = gmm::abs(zshift[k]) >= 1e-10;
364  if (not(has_xshift) && not(has_yshift) && not(has_yshift))
365  funcs.push_back(global_function_bspline
366  (xmin[i], xmax[i],
367  ymin[j], ymax[j],
368  zmin[k], zmax[k],
369  order, xtype[i], ytype[j], ztype[k])) ;
370  else {
371  std::vector<pglobal_function> sum;
372  sum.push_back(global_function_bspline
373  (xmin[i], xmax[i],
374  ymin[j], ymax[j],
375  zmin[k], zmax[k],
376  order, xtype[i], ytype[j], ztype[k]));
377  if (has_xshift) // xshift
378  sum.push_back(global_function_bspline
379  (xmin[i]+xshift[i], xmax[i]+xshift[i],
380  ymin[j], ymax[j],
381  zmin[k], zmax[k],
382  order, xtype[i], ytype[j], ztype[k]));
383  if (has_yshift) // yshift
384  sum.push_back(global_function_bspline
385  (xmin[i], xmax[i],
386  ymin[j]+yshift[j], ymax[j]+yshift[j],
387  zmin[k], zmax[k],
388  order, xtype[i], ytype[j], ztype[k]));
389  if (has_zshift) // zshift
390  sum.push_back(global_function_bspline
391  (xmin[i], xmax[i],
392  ymin[j], ymax[j],
393  zmin[k]+zshift[k], zmax[k]+zshift[k],
394  order, xtype[i], ytype[j], ztype[k]));
395  if (has_xshift && has_yshift) // xshift + yshift
396  sum.push_back(global_function_bspline
397  (xmin[i]+xshift[i], xmax[i]+xshift[i],
398  ymin[j]+yshift[j], ymax[j]+yshift[j],
399  zmin[k], zmax[k],
400  order, xtype[i], ytype[j], ztype[k]));
401  if (has_yshift && has_zshift) // yshift + zshift
402  sum.push_back(global_function_bspline
403  (xmin[i], xmax[i],
404  ymin[j]+yshift[j], ymax[j]+yshift[j],
405  zmin[k]+zshift[k], zmax[k]+zshift[k],
406  order, xtype[i], ytype[j], ztype[k]));
407  if (has_zshift && has_xshift) // zshift + xshift
408  sum.push_back(global_function_bspline
409  (xmin[i]+xshift[i], xmax[i]+xshift[i],
410  ymin[j], ymax[j],
411  zmin[k]+zshift[k], zmax[k]+zshift[k],
412  order, xtype[i], ytype[j], ztype[k]));
413  if (has_xshift && has_yshift && has_zshift) // xshift + yshift + zshift
414  sum.push_back(global_function_bspline
415  (xmin[i]+xshift[i], xmax[i]+xshift[i],
416  ymin[j]+yshift[j], ymax[j]+yshift[j],
417  zmin[k]+zshift[k], zmax[k]+zshift[k],
418  order, xtype[i], ytype[j], ztype[k]));
419  funcs.push_back(std::make_shared<getfem::global_function_sum>(sum));
420  }
421  } // k
422  } // j
423  } // i
424  mf.set_functions(funcs, mim);
425  }
426 
427 }
428 
429 /* end of namespace getfem */
getfem::mesh_fem_global_function
this is a convenience class for defining a mesh_fem with base functions which are global functions (f...
Definition: getfem_mesh_fem_global_function.h:50
getfem::mesh_fem::linked_mesh
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Definition: getfem_mesh_fem.h:283
getfem::mesh_fem::set_finite_element
void set_finite_element(size_type cv, pfem pf)
Set the finite element method of a convex.
Definition: getfem_mesh_fem.cc:127
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
getfem::mesh::bounding_box
void bounding_box(base_node &Pmin, base_node &Pmax) const
Return the bounding box [Pmin - Pmax] of the mesh.
Definition: getfem_mesh.cc:261
getfem_mesh_fem_global_function.h
Define a mesh_fem with base functions which are global functions given by the user.
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:47
getfem::del_fem_global_function
void del_fem_global_function(const pfem &pf)
release a global function FEM
Definition: getfem_fem_global_function.h:115
getfem::define_uniform_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 order, bspline_boundary bcX_low=bspline_boundary::FREE, bspline_boundary bcX_high=bspline_boundary::FREE, const mesh_im &mim=dummy_mesh_im())
This function will generate bspline basis functions on NX uniform elements along a line.
Definition: getfem_mesh_fem_global_function.cc:224
getfem::dummy_mesh_im
const mesh_im & dummy_mesh_im()
Dummy mesh_im for default parameter of functions.
Definition: getfem_mesh_im.cc:234
getfem::new_fem_global_function
pfem new_fem_global_function(const std::vector< pglobal_function > &funcs, const mesh &m)
create a new global function FEM.
Definition: getfem_fem_global_function.cc:310