GetFEM  5.4.2
getfem_fem_global_function.cc
1 /*===========================================================================
2 
3  Copyright (C) 2004-2020 Yves Renard
4  Copyright (C) 2016 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 
24 
25 namespace getfem {
26 
27 
28  void fem_global_function::init() {
29  is_pol = is_lag = is_standard_fem = false; es_degree = 5;
30  is_equiv = real_element_defined = true;
31  ntarget_dim = 1; // An extension for vectorial elements should be easy
32 
33  std::stringstream nm;
34  nm << "GLOBAL_FEM(" << (void*)this << ")";
35  debug_name_ = nm.str();
36 
37  GMM_ASSERT1(functions.size() > 0, "Empty list of base functions.");
38  dim_ = functions[0]->dim();
39  for (size_type i=1; i < functions.size(); ++i)
40  GMM_ASSERT1(dim_ == functions[i]->dim(),
41  "Incompatible function space dimensions.");
42 
44  }
45 
46 
47  fem_global_function::fem_global_function
48  (const std::vector<pglobal_function> &funcs, const mesh &m_)
49  : functions(funcs), m(m_), mim(dummy_mesh_im()), has_mesh_im(false) {
50 
51  DAL_STORED_OBJECT_DEBUG_CREATED(this, "Global function fem");
52  GMM_ASSERT1(&m != &dummy_mesh(), "A non-empty mesh object"
53  " is expected.");
54  this->add_dependency(m_);
55  init();
56  }
57 
58  fem_global_function::fem_global_function
59  (const std::vector<pglobal_function> &funcs, const mesh_im &mim_)
60  : functions(funcs), m(mim_.linked_mesh()), mim(mim_), has_mesh_im(true) {
61 
62  DAL_STORED_OBJECT_DEBUG_CREATED(this, "Global function fem");
63  GMM_ASSERT1(&mim != &dummy_mesh_im(), "A non-empty mesh_im object"
64  " is expected.");
65  this->add_dependency(mim_);
66  init();
67  }
68 
69  void fem_global_function::update_from_context() const {
70 
71  if (precomps) {
72  for (const auto &cv_precomps : *precomps)
73  for (const auto &keyval : cv_precomps)
74  dal::del_dependency(precomps, keyval.first);
75  precomps->clear();
76  } else {
77  precomps = std::make_shared<precomp_pool>();
78  dal::pstatic_stored_object_key pkey
79  = std::make_shared<precomp_pool_key>(precomps);
80  dal::add_stored_object(pkey, precomps);
81  }
82 
83  size_type nb_total_dof(functions.size());
84  base_node bmin(dim()), bmax(dim());
85  bgeot::rtree boxtree{1E-13};
86  std::map<size_type, std::vector<size_type>> box_to_convexes_map;
87  for (size_type i=0; i < nb_total_dof; ++i) {
88  functions[i]->bounding_box(bmin, bmax);
89  box_to_convexes_map[boxtree.add_box(bmin, bmax, i)].push_back(i);
90  }
91  boxtree.build_tree();
92 
93  size_type max_dof(0);
94  index_of_global_dof_.clear();
95  index_of_global_dof_.resize(m.nb_allocated_convex());
96  for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
97  GMM_ASSERT1(dim_ == m.structure_of_convex(cv)->dim(),
98  "Convexes of different dimension: to be done");
99  bgeot::pgeometric_trans pgt = m.trans_of_convex(cv);
100 
101  bounding_box(bmin, bmax, m.points_of_convex(cv), pgt);
102 
103  bgeot::rtree::pbox_set boxlst;
104  boxtree.find_intersecting_boxes(bmin, bmax, boxlst);
105  index_of_global_dof_[cv].clear();
106 
107  if (has_mesh_im) {
108  pintegration_method pim = mim.int_method_of_element(cv);
109  GMM_ASSERT1(pim->type() == IM_APPROX, "You have to use approximated "
110  "integration in connection to a fem with global functions");
111  papprox_integration pai = pim->approx_method();
112 
113  for (const auto &box : boxlst) {
114  for (auto candidate : box_to_convexes_map.at(box->id)) {
115  for (size_type k = 0; k < pai->nb_points(); ++k) {
116  base_node gpt = pgt->transform(pai->point(k),
117  m.points_of_convex(cv));
118  if (functions[candidate]->is_in_support(gpt)) {
119  index_of_global_dof_[cv].push_back(candidate);
120  break;
121  }
122  }
123  }
124  }
125  } else { // !has_mesh_im
126  for (const auto &box : boxlst) {
127  for (auto candidate : box_to_convexes_map.at(box->id))
128  index_of_global_dof_[cv].push_back(candidate);
129  }
130  }
131  max_dof = std::max(max_dof, index_of_global_dof_[cv].size());
132  }
133 
134  /** setup global dofs, with dummy coordinates */
135  base_node P(dim()); gmm::fill(P,1./20);
136  std::vector<base_node> node_tab_(max_dof, P);
137  pspt_override = bgeot::store_point_tab(node_tab_);
138  pspt_valid = false;
139  dof_types_.resize(nb_total_dof);
140  std::fill(dof_types_.begin(), dof_types_.end(),
141  global_dof(dim()));
142  }
143 
144  size_type fem_global_function::nb_dof(size_type cv) const {
145  //return functions.size();
146  context_check();
147  return (cv < index_of_global_dof_.size()) ? index_of_global_dof_[cv].size()
148  : size_type(0);
149  }
150 
151  size_type fem_global_function::index_of_global_dof
152  (size_type cv, size_type i) const {
153  //return i;
154  context_check();
155  return (cv < index_of_global_dof_.size() &&
156  i < index_of_global_dof_[cv].size()) ? index_of_global_dof_[cv][i]
157  : size_type(-1);
158  }
159 
160  bgeot::pconvex_ref fem_global_function::ref_convex(size_type cv) const {
161  if (has_mesh_im)
162  return mim.int_method_of_element(cv)->approx_method()->ref_convex();
163  else
164  return bgeot::basic_convex_ref(m.trans_of_convex(cv)->convex_ref());
165  }
166 
168  fem_global_function::node_convex(size_type cv) const {
169  if (m.convex_index().is_in(cv))
171  (dim(), nb_dof(cv), m.structure_of_convex(cv)->nb_faces()));
172  else GMM_ASSERT1(false, "Wrong convex number: " << cv);
173  }
174 
175  void fem_global_function::base_value(const base_node &, base_tensor &) const
176  { GMM_ASSERT1(false, "No base values, real only element."); }
177  void fem_global_function::grad_base_value(const base_node &,
178  base_tensor &) const
179  { GMM_ASSERT1(false, "No grad values, real only element."); }
180  void fem_global_function::hess_base_value(const base_node &,
181  base_tensor &) const
182  { GMM_ASSERT1(false, "No hess values, real only element."); }
183 
184  void fem_global_function::real_base_value(const fem_interpolation_context& c,
185  base_tensor &t, bool) const {
186  assert(target_dim() == 1);
187  size_type cv = c.convex_num();
188  size_type nbdof = nb_dof(cv);
189  t.adjust_sizes(nbdof, target_dim());
190  if (c.have_pfp() && c.ii() != size_type(-1)) {
191  GMM_ASSERT1(precomps, "Internal error");
192  if (precomps->size() == 0)
193  precomps->resize(m.nb_allocated_convex());
194  GMM_ASSERT1(precomps->size() == m.nb_allocated_convex(), "Internal error");
195  const bgeot::pstored_point_tab ptab = c.pfp()->get_ppoint_tab();
196  auto it = (*precomps)[cv].find(ptab);
197  if (it == (*precomps)[cv].end()) {
198  it = (*precomps)[cv].emplace(ptab, precomp_data()).first;
199  dal::add_dependency(precomps, ptab);
200  // we could have added the dependency to this->shared_from_this()
201  // instead, but there is a risk that this will shadow the same
202  // dependency through a different path, so that it becomes dangerous
203  // to delete the dependency later
204  }
205  if (it->second.val.size() == 0) {
206  it->second.val.resize(ptab->size());
207  base_matrix G;
208  bgeot::vectors_to_base_matrix(G, m.points_of_convex(cv));
209  for (size_type k = 0; k < ptab->size(); ++k) {
211  ctx(m.trans_of_convex(cv), shared_from_this(), (*ptab)[k], G, cv);
212  real_base_value(ctx, it->second.val[k]);
213  }
214  }
215  gmm::copy(it->second.val[c.ii()].as_vector(), t.as_vector());
216  } else
217  for (size_type i=0; i < nbdof; ++i) {
218  /*cerr << "fem_global_function: real_base_value(" << c.xreal() << ")\n";
219  if (c.have_G()) cerr << "G = " << c.G() << "\n";
220  else cerr << "no G\n";*/
221  t[i] = functions[index_of_global_dof_[cv][i]]->val(c);
222  }
223  }
224 
225  void fem_global_function::real_grad_base_value
226  (const fem_interpolation_context& c, base_tensor &t, bool) const {
227  assert(target_dim() == 1);
228  size_type cv = c.convex_num();
229  size_type nbdof = nb_dof(cv);
230  t.adjust_sizes(nbdof, target_dim(), dim());
231  if (c.have_pfp() && c.ii() != size_type(-1)) {
232  GMM_ASSERT1(precomps, "Internal error");
233  if (precomps->size() == 0)
234  precomps->resize(m.nb_allocated_convex());
235  GMM_ASSERT1(precomps->size() == m.nb_allocated_convex(), "Internal error");
236  const bgeot::pstored_point_tab ptab = c.pfp()->get_ppoint_tab();
237  auto it = (*precomps)[cv].find(ptab);
238  if (it == (*precomps)[cv].end()) {
239  it = (*precomps)[cv].emplace(ptab, precomp_data()).first;
240  dal::add_dependency(precomps, ptab);
241  }
242  if (it->second.grad.size() == 0) {
243  it->second.grad.resize(ptab->size());
244  base_matrix G;
245  bgeot::vectors_to_base_matrix(G, m.points_of_convex(cv));
246  for (size_type k = 0; k < ptab->size(); ++k) {
248  ctx(m.trans_of_convex(cv), shared_from_this(), (*ptab)[k], G, cv);
249  real_grad_base_value(ctx, it->second.grad[k]);
250  }
251  }
252  gmm::copy(it->second.grad[c.ii()].as_vector(), t.as_vector());
253  } else {
254  base_small_vector G(dim());
255  for (size_type i=0; i < nbdof; ++i) {
256  functions[index_of_global_dof_[cv][i]]->grad(c,G);
257  for (size_type j=0; j < dim(); ++j)
258  t[j*nbdof + i] = G[j];
259  }
260  }
261  }
262 
263  void fem_global_function::real_hess_base_value
264  (const fem_interpolation_context &c, base_tensor &t, bool) const {
265  assert(target_dim() == 1);
266  size_type cv = c.convex_num();
267  size_type nbdof = nb_dof(cv);
268  t.adjust_sizes(nbdof, target_dim(), gmm::sqr(dim()));
269  if (c.have_pfp() && c.ii() != size_type(-1)) {
270  GMM_ASSERT1(precomps, "Internal error");
271  if (precomps->size() == 0)
272  precomps->resize(m.nb_allocated_convex());
273  GMM_ASSERT1(precomps->size() == m.nb_allocated_convex(), "Internal error");
274  const bgeot::pstored_point_tab ptab = c.pfp()->get_ppoint_tab();
275  auto it = (*precomps)[cv].find(ptab);
276  if (it == (*precomps)[cv].end()) {
277  it = (*precomps)[cv].emplace(ptab, precomp_data()).first;
278  dal::add_dependency(precomps, ptab);
279  }
280  if (it->second.hess.size() == 0) {
281  it->second.hess.resize(ptab->size());
282  base_matrix G;
283  bgeot::vectors_to_base_matrix(G, m.points_of_convex(cv));
284  for (size_type k = 0; k < ptab->size(); ++k) {
286  ctx(m.trans_of_convex(cv), shared_from_this(), (*ptab)[k], G, cv);
287  real_hess_base_value(ctx, it->second.hess[k]);
288  }
289  }
290  gmm::copy(it->second.hess[c.ii()].as_vector(), t.as_vector());
291  } else {
292  base_matrix H(dim(),dim());
293  for (size_type i=0; i < nbdof; ++i) {
294  functions[index_of_global_dof_[cv][i]]->hess(c,H);
295  for (size_type jk=0; jk < size_type(dim()*dim()); ++jk)
296  t[jk*nbdof + i] = H[jk];
297  }
298  }
299  }
300 
301 
302  DAL_SIMPLE_KEY(special_fem_globf_key, pfem);
303 
304  pfem new_fem_global_function(const std::vector<pglobal_function> &funcs,
305  const mesh &m) {
306  pfem pf = std::make_shared<fem_global_function>(funcs, m);
307  dal::pstatic_stored_object_key
308  pk = std::make_shared<special_fem_globf_key>(pf);
309  dal::add_stored_object(pk, pf);
310  return pf;
311  }
312 
313  pfem new_fem_global_function(const std::vector<pglobal_function> &funcs,
314  const mesh_im &mim) {
315  pfem pf = std::make_shared<fem_global_function>(funcs, mim);
316  dal::pstatic_stored_object_key
317  pk = std::make_shared<special_fem_globf_key>(pf);
318  dal::add_stored_object(pk, pf);
319  return pf;
320  }
321 
322 }
323 
324 /* end of namespace getfem */
getfem::virtual_fem::dim
dim_type dim() const
dimension of the reference element.
Definition: getfem_fem.h:311
bgeot::basic_convex_ref
pconvex_ref basic_convex_ref(pconvex_ref cvr)
return the associated order 1 reference convex.
Definition: bgeot_convex_ref.h:138
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::mesh_im
Describe an integration method linked to a mesh.
Definition: getfem_mesh_im.h:47
bgeot::convex< base_node >
dal::del_dependency
bool del_dependency(pstatic_stored_object o1, pstatic_stored_object o2)
remove a dependency.
Definition: dal_static_stored_objects.cc:256
getfem::fem_interpolation_context::convex_num
size_type convex_num() const
get the current convex number
Definition: getfem_fem.cc:52
bgeot::rtree
Balanced tree of n-dimensional rectangles.
Definition: bgeot_rtree.h:98
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::fem_interpolation_context
structure passed as the argument of fem interpolation functions.
Definition: getfem_fem.h:749
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:304
getfem_fem_global_function.h
Define mesh_fem whose base functions are global function given by the user.
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
getfem::fem_interpolation_context::have_pfp
bool have_pfp() const
true if a fem_precomp_ has been supplied.
Definition: getfem_fem.h:760
getfem::fem_global_function::update_from_context
virtual void update_from_context() const
this function has to be defined and should update the object when the context is modified.
Definition: getfem_fem_global_function.cc:69
getfem::dummy_mesh_im
const mesh_im & dummy_mesh_im()
Dummy mesh_im for default parameter of functions.
Definition: getfem_mesh_im.cc:234
dal::add_stored_object
void add_stored_object(pstatic_stored_object_key k, pstatic_stored_object o, permanence perm)
Add an object with two optional dependencies.
Definition: dal_static_stored_objects.cc:284
getfem::fem_interpolation_context::pfp
pfem_precomp pfp() const
get the current fem_precomp_
Definition: getfem_fem.h:794
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
bgeot::generic_dummy_convex_ref
pconvex_ref generic_dummy_convex_ref(dim_type nc, size_type n, short_type nf)
generic convex with n global nodes
Definition: bgeot_convex_ref.cc:867
getfem::global_dof
pdof_description global_dof(dim_type d)
Description of a global dof, i.e.
Definition: getfem_fem.cc:542
dal::add_dependency
void add_dependency(pstatic_stored_object o1, pstatic_stored_object o2)
Add a dependency, object o1 will depend on object o2.
Definition: dal_static_stored_objects.cc:230