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