GetFEM  5.4.2
getfem_mesh_fem_sum.cc
1 /*===========================================================================
2 
3  Copyright (C) 1999-2020 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 
24 
25 namespace getfem {
26 
27  void fem_sum::init() {
28  cvr = pfems[0]->ref_convex(cv);
29  dim_ = cvr->structure()->dim();
30  is_equiv = !smart_global_dof_linking_;
31  real_element_defined = true;
32  is_polycomp = is_pol = is_lag = is_standard_fem = false;
33  es_degree = 5; /* humm ... */
34  ntarget_dim = 1;
35 
36  std::stringstream nm;
37  nm << "FEM_SUM(" << pfems[0]->debug_name() << ",";
38  for (size_type i = 1; i < pfems.size(); ++i)
39  nm << pfems[i]->debug_name() << ",";
40  nm << " cv:" << cv << ")";
41  debug_name_ = nm.str();
42 
43  init_cvs_node();
44  for (size_type i = 0; i < pfems.size(); ++i) {
45  GMM_ASSERT1(pfems[i]->target_dim() == 1, "Vectorial fems not supported");
46 
47  for (size_type k = 0; k < pfems[i]->nb_dof(cv); ++k) {
48  add_node(pfems[i]->dof_types()[k], pfems[i]->node_of_dof(cv,k));
49  }
50  }
51  }
52 
53  /* used only when smart_global_dof_linking_ is on */
54  void fem_sum::mat_trans(base_matrix &M,
55  const base_matrix &G,
56  bgeot::pgeometric_trans pgt) const {
57 
58  // cerr << "fem_sum::mat_trans " << debug_name_
59  // << " / smart_global_dof_linking_ = " << smart_global_dof_linking_
60  // << "\n";
61  pdof_description gdof = 0, lagdof = lagrange_dof(dim());
62  std::vector<pdof_description> hermdof(dim());
63  for (dim_type id = 0; id < dim(); ++id)
64  hermdof[id] = derivative_dof(dim(), id);
65  if (pfems.size()) gdof = global_dof(dim());
66  gmm::copy(gmm::identity_matrix(), M);
67  base_vector val(1), val2(1);
68  base_matrix grad(1, dim());
69 
70  // very inefficient, to be optimized ...
71  for (size_type ifem1 = 0, i=0; ifem1 < pfems.size(); ++ifem1) {
72  pfem pf1 = pfems[ifem1];
73  /* find global dofs */
74  for (size_type idof1 = 0; idof1 < pf1->nb_dof(cv); ++idof1, ++i) {
75  if (pf1->dof_types()[idof1] == gdof) {
76  base_vector coeff(pfems[ifem1]->nb_dof(cv));
77  coeff[idof1] = 1.0;
78  fem_interpolation_context fic(pgt, pf1, base_node(dim()), G, cv);
79  for (size_type ifem2 = 0, j=0; ifem2 < pfems.size(); ++ifem2) {
80  pfem pf2 = pfems[ifem2];
81  fem_interpolation_context fic2(pgt, pf2, base_node(dim()), G, cv);
82  for (size_type idof2 = 0; idof2 < pf2->nb_dof(cv); ++idof2, ++j) {
83  pdof_description pdd = pf2->dof_types()[idof2];
84  bool found = false;
85 
86  base_vector coeff2(pfems[ifem2]->nb_dof(cv));
87  coeff2[idof2] = 1.0;
88 
89  if (pdd != gdof) {
90  fic.set_xref(pf2->node_of_dof(cv, idof2));
91  fic2.set_xref(pf2->node_of_dof(cv, idof2));
92  if (dof_weak_compatibility(pdd, lagdof) == 0) {
93  pfems[ifem1]->interpolation(fic, coeff, val, 1);
94 
95  pfems[ifem2]->interpolation(fic2, coeff2, val2, 1);
96 
97  M(i, j) = -val[0]*val2[0];
98  /*if (pf2->nb_dof(cv) > 4)
99  cout << "dof " << idof2 << " ("
100  << pf2->node_of_dof(cv,idof2)
101  << ") " << (void*)&(*pdd)
102  << " compatible with lagrange\n";*/
103  found = true;
104  }
105  else for (size_type id = 0; id < dim(); ++id) {
106  if (dof_weak_compatibility(pdd, hermdof[id]) == 0) {
107  pfems[ifem1]->interpolation_grad(fic, coeff, grad, 1);
108  M(i, j) = -grad(0, id);
109  cout << "dof " << idof2 << "compatible with hermite "
110  << id << "\n";
111  found = true;
112  }
113  }
114  GMM_ASSERT1(found,
115  "Sorry, smart_global_dof_linking not "
116  "compatible with this kind of dof");
117  }
118  }
119  //if (pf2->nb_dof(cv) > 4) cout << " M=" << M << "\n";
120  }
121  }
122  }
123  }
124 
125  // static int cnt=0;
126  // if(++cnt < 40) cout << "fem = " << debug_name_ << ", M = " << M << "\n";
127  }
128 
129  size_type fem_sum::index_of_global_dof(size_type , size_type j) const {
130  for (size_type i = 0; i < pfems.size(); ++i) {
131  size_type nb = pfems[i]->nb_dof(cv);
132  if (j >= nb) j -= pfems[i]->nb_dof(cv);
133  else return pfems[i]->index_of_global_dof(cv, j);
134  }
135  GMM_ASSERT1(false, "incoherent global dof.");
136  }
137 
138  void fem_sum::base_value(const base_node &,
139  base_tensor &) const
140  { GMM_ASSERT1(false, "No base values, real only element."); }
141  void fem_sum::grad_base_value(const base_node &,
142  base_tensor &) const
143  { GMM_ASSERT1(false, "No base values, real only element."); }
144  void fem_sum::hess_base_value(const base_node &,
145  base_tensor &) const
146  { GMM_ASSERT1(false, "No base values, real only element."); }
147 
148  void fem_sum::real_base_value(const fem_interpolation_context &c,
149  base_tensor &t,
150  bool withM) const {
151  bgeot::multi_index mi(2);
152  mi[1] = target_dim(); mi[0] = short_type(nb_dof(0));
153  t.adjust_sizes(mi);
154  base_tensor::iterator it = t.begin(), itf;
155 
156  fem_interpolation_context c0 = c;
157  std::vector<base_tensor> val_e(pfems.size());
158  for (size_type k = 0; k < pfems.size(); ++k) {
159  if (c0.have_pfp()) {
160  c0.set_pfp(fem_precomp(pfems[k], c0.pfp()->get_ppoint_tab(),
161  c0.pfp()));
162  } else { c0.set_pf(pfems[k]); }
163  c0.base_value(val_e[k]);
164  }
165 
166  for (dim_type q = 0; q < target_dim(); ++q) {
167  for (size_type k = 0; k < pfems.size(); ++k) {
168  itf = val_e[k].begin() + q * pfems[k]->nb_dof(cv);
169  for (size_type i = 0; i < pfems[k]->nb_dof(cv); ++i)
170  *it++ = *itf++;
171  }
172  }
173  assert(it == t.end());
174  if (!is_equivalent() && withM) {
175  base_tensor tt(t);
176  t.mat_transp_reduction(tt, c.M(), 0);
177  }
178  //cerr << "fem_sum::real_base_value(" << c.xreal() << ")\n";
179  }
180 
181  void fem_sum::real_grad_base_value(const fem_interpolation_context &c,
182  base_tensor &t, bool withM) const {
183  bgeot::multi_index mi(3);
184  mi[2] = short_type(c.N()); mi[1] = target_dim();
185  mi[0] = short_type(nb_dof(0));
186  t.adjust_sizes(mi);
187  base_tensor::iterator it = t.begin(), itf;
188 
189  fem_interpolation_context c0 = c;
190  std::vector<base_tensor> grad_e(pfems.size());
191  for (size_type k = 0; k < pfems.size(); ++k) {
192  if (c0.have_pfp()) {
193  c0.set_pfp(fem_precomp(pfems[k], c0.pfp()->get_ppoint_tab(),
194  c0.pfp()));
195  } else { c0.set_pf(pfems[k]); }
196  c0.grad_base_value(grad_e[k]);
197  }
198 
199  for (dim_type k = 0; k < c.N() ; ++k) {
200  for (dim_type q = 0; q < target_dim(); ++q) {
201  for (size_type f = 0; f < pfems.size(); ++f) {
202  itf = grad_e[f].begin()
203  + (k * target_dim() + q) * pfems[f]->nb_dof(cv);
204  for (size_type i = 0; i < pfems[f]->nb_dof(cv); ++i)
205  *it++ = *itf++;
206  }
207  }
208  }
209  assert(it == t.end());
210  if (!is_equivalent() && withM) {
211  base_tensor tt(t);
212  t.mat_transp_reduction(tt, c.M(), 0);
213  }
214  }
215 
216  void fem_sum::real_hess_base_value(const fem_interpolation_context &c,
217  base_tensor &t, bool withM) const {
218  t.adjust_sizes(nb_dof(0), target_dim(), gmm::sqr(c.N()));
219  base_tensor::iterator it = t.begin(), itf;
220 
221  fem_interpolation_context c0 = c;
222  std::vector<base_tensor> hess_e(pfems.size());
223  for (size_type k = 0; k < pfems.size(); ++k) {
224  if (c0.have_pfp()) {
225  c0.set_pfp(fem_precomp(pfems[k], c0.pfp()->get_ppoint_tab(),
226  c0.pfp()));
227  } else { c0.set_pf(pfems[k]); }
228  c0.hess_base_value(hess_e[k]);
229  }
230 
231  dim_type NNdim = dim_type(gmm::sqr(c.N())*target_dim());
232  for (dim_type jkq = 0; jkq < NNdim ; ++jkq) {
233  for (size_type f = 0; f < pfems.size(); ++f) {
234  itf = hess_e[f].begin() + (jkq * pfems[f]->nb_dof(cv));
235  for (size_type i = 0; i < pfems[f]->nb_dof(cv); ++i)
236  *it++ = *itf++;
237  }
238  }
239  assert(it == t.end());
240  if (!is_equivalent() && withM) {
241  base_tensor tt(t);
242  t.mat_transp_reduction(tt, c.M(), 0);
243  }
244  }
245 
246  void mesh_fem_sum::clear_build_methods() {
247  for (size_type i = 0; i < build_methods.size(); ++i)
248  del_stored_object(build_methods[i]);
249  build_methods.clear();
250  }
251  void mesh_fem_sum::clear(void) {
252  mesh_fem::clear();
253  clear_build_methods();
254  situations.clear();
255  is_adapted = false;
256  }
257 
258  DAL_SIMPLE_KEY(special_mflsum_key, pfem);
259 
260  void mesh_fem_sum::adapt(void) {
261  context_check();
262  clear();
263 
264  for (size_type i = 0; i < mfs.size(); ++i)
265  GMM_ASSERT1(!(mfs[i]->is_reduced()),
266  "Sorry fem_sum for reduced mesh_fem is not implemented");
267 
268  for (dal::bv_visitor i(linked_mesh().convex_index()); !i.finished(); ++i) {
269  std::vector<pfem> pfems;
270  bool is_cv_dep = false;
271  for (size_type j = 0; j < mfs.size(); ++j) {
272  if (mfs[j]->convex_index().is_in(i)) {
273  pfem pf = mfs[j]->fem_of_element(i);
274  if (pf->nb_dof(i)) {
275  pfems.push_back(pf);
276  if (pf->is_on_real_element()) is_cv_dep = true;
277  }
278  }
279  }
280  if (pfems.size() == 1) {
281  set_finite_element(i, pfems[0]);
282  }
283  else if (pfems.size() > 0) {
284  if (situations.find(pfems) == situations.end() || is_cv_dep) {
285  pfem pf = std::make_shared<fem_sum>(pfems, i,
286  smart_global_dof_linking_);
287  dal::pstatic_stored_object_key
288  pk = std::make_shared<special_mflsum_key>(pf);
289  dal::add_stored_object(pk, pf, pf->ref_convex(0), pf->node_tab(0));
290  build_methods.push_back(pf);
291  situations[pfems] = pf;
292  }
293  set_finite_element(i, situations[pfems]);
294  }
295  }
296  is_adapted = true; touch();
297  }
298 
299 
300 } /* end of namespace getfem. */
301 
getfem::mesh_fem::is_reduced
bool is_reduced() const
Return true if a reduction matrix is applied to the dofs.
Definition: getfem_mesh_fem.h:195
getfem::pdof_description
dof_description * pdof_description
Type representing a pointer on a dof_description.
Definition: getfem_fem.h:160
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
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::fem_precomp
pfem_precomp fem_precomp(pfem pf, bgeot::pstored_point_tab pspt, dal::pstatic_stored_object dep)
Handles precomputations for FEM.
Definition: getfem_fem.cc:4333
dal::del_stored_object
void del_stored_object(const pstatic_stored_object &o, bool ignore_unstored)
Delete an object and the object which depend on it.
Definition: dal_static_stored_objects.cc:369
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::mesh_fem::convex_index
const dal::bit_vector & convex_index() const
Get the set of convexes where a finite element has been assigned.
Definition: getfem_mesh_fem.h:191
getfem::derivative_dof
pdof_description derivative_dof(dim_type d, dim_type r)
Description of a unique dof of derivative type.
Definition: getfem_fem.cc:487
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
getfem::mesh_fem::linked_mesh
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Definition: getfem_mesh_fem.h:279
getfem::lagrange_dof
pdof_description lagrange_dof(dim_type d)
Description of a unique dof of lagrange type (value at the node).
Definition: getfem_fem.cc:390
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
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
getfem::context_dependencies::context_check
bool context_check() const
return true if update_from_context was called
Definition: getfem_context.h:126
getfem_mesh_fem_sum.h
Implement a special mesh_fem with merges the FEMs of two (or more) mesh_fems.
getfem::global_dof
pdof_description global_dof(dim_type d)
Description of a global dof, i.e.
Definition: getfem_fem.cc:542