GetFEM  5.4.2
getfem_fem_level_set.cc
Go to the documentation of this file.
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 /** \file getfem_fem_level_set.cc
22  \brief an FEM which should be used with getfem::mesh_fem_level_set.
23 */
25 
26 namespace getfem {
27 
28  void fem_level_set::init() {
29  cvr = bfem->ref_convex(0);
30  dim_ = cvr->structure()->dim();
31  is_equiv = true; real_element_defined = true;
32  is_polycomp = is_pol = is_lag = is_standard_fem = false;
33  es_degree = 5; /* humm ... */
34  ntarget_dim = bfem->target_dim();
35 
36  std::stringstream nm;
37  nm << "FEM_LEVEL_SET(" << bfem->debug_name() << ")";
38  debug_name_ = nm.str();
39 
40  ls_index.sup(0, mls.nb_level_sets());
41 
42  common_ls_zones.resize(mls.nb_level_sets());
43  /* fill ls_index .. */
44  for (size_type i=0; i < mls.nb_level_sets(); ++i) {
45  char c = '*';
46  for (size_type k=0; k < bfem->nb_dof(0); ++k) {
47  const mesh_level_set::zoneset *ze = dofzones[k];
48  if (ze) {
49  for (mesh_level_set::zoneset::const_iterator itz = ze->begin();
50  itz != ze->end(); ++itz) {
51  const mesh_level_set::zone *z = *itz;
52  for (mesh_level_set::zone::const_iterator it = z->begin();
53  it != z->end(); ++it) {
54  assert((**it).size() == mls.nb_level_sets());
55  char d = (*(*it))[i];
56  if (c == '*') c = d;
57  else if (c != d) { ls_index.add(i); break; }
58  }
59  }
60  }
61  }
62  common_ls_zones[i] = c;
63  }
64 
65  init_cvs_node();
66  for (size_type k = 0; k < bfem->nb_dof(0); ++k) {
67  if (!dofzones[k]) {
68  add_node(bfem->dof_types()[k], bfem->node_of_dof(0,k));
69  } else {
70  for (size_type j = 0; j < dofzones[k]->size(); ++j) {
71  // cout << " -> +dof: '" << *(*dofzones[k])[j] << "'\n";
72  add_node(xfem_dof(bfem->dof_types()[k], j+xfem_index),
73  bfem->node_of_dof(0,k));
74  }
75  }
76  }
77  }
78 
79  void fem_level_set::base_value(const base_node &, base_tensor &) const
80  { GMM_ASSERT1(false, "No base values, real only element."); }
81  void fem_level_set::grad_base_value(const base_node &,
82  base_tensor &) const
83  { GMM_ASSERT1(false, "No base values, real only element."); }
84  void fem_level_set::hess_base_value(const base_node &,
85  base_tensor &) const
86  { GMM_ASSERT1(false, "No base values, real only element."); }
87 
88  static bool are_zones_compatible_(const std::string a, const std::string b) {
89  if (a.size() != b.size()) return false;
90  for (size_type i = 0; i < a.size(); ++i)
91  if (a[i] != '0' && a[i] != b[i]) return false;
92  return true;
93  }
94 
95  void fem_level_set::find_zone_id(const fem_interpolation_context &c,
96  std::vector<bool> &ids, int side) const {
97  size_type s = 0, cv = c.convex_num();
98  for (size_type i = 0; i < dofzones.size(); ++i)
99  if (dofzones[i]) s += dofzones[i]->size();
100  ids.resize(0); ids.resize(s, false);
101  std::string z(common_ls_zones);
102  base_vector coeff(32);
103 
104  mesher_level_set eval;
105  size_type iclosest = size_type(-1); scalar_type vclosest = 1E100;
106  for (dal::bv_visitor i(ls_index); !i.finished(); ++i) {
107  const level_set *ls = mls.get_level_set(i);
108  const mesh_fem &mf = ls->get_mesh_fem();
109  slice_vector_on_basic_dof_of_element(mf, ls->values(), cv, coeff);
110  eval.init_base(mf.fem_of_element(cv), coeff);
111  eval.set_shift(ls->get_shift()); // Deprecated
112 
113  // mesher_level_set eval = mls.get_level_set(i)->mls_of_convex(cv);
114 
115  scalar_type v = eval(c.xref());
116  if (side != 0) {
117  if (gmm::abs(v) < vclosest) { vclosest = gmm::abs(v); iclosest = i; }
118  }
119  z[size_type(i)] = (v > 0.) ? '+' : '-';
120  }
121 
122  if (side != 0 && iclosest != size_type(-1)) // Forces the side of the
123  // closest level-set (in order to compute jumps).
124  z[iclosest] = (side > 0) ? '+' : '-';
125 
126 
127  unsigned cnt = 0;
128  for (unsigned d = 0; d < dofzones.size(); ++d) {
129  if (!dofzones[d]) continue;
130  for (mesh_level_set::zoneset::const_iterator it = dofzones[d]->begin();
131  it != dofzones[d]->end(); ++it, ++cnt) {
132  ids[cnt] = false;
133  for (mesh_level_set::zone::const_iterator it2 = (*it)->begin();
134  it2 != (*it)->end(); ++it2) {
135  if (are_zones_compatible_(z,*(*it2))) { ids[cnt] = true; break; }
136  }
137  }
138  }
139  }
140 
142  base_tensor &t, bool) const {
143  // bgeot::multi_index mi(2);
144  // mi[1] = target_dim(); mi[0] = short_type(nb_base(0));
145  t.adjust_sizes(nb_base(0), target_dim());
146  base_tensor::iterator it = t.begin();
148  if (c0.have_pfp())
149  c0.set_pfp(fem_precomp(bfem, c0.pfp()->get_ppoint_tab(), c0.pfp()));
150  else c0.set_pf(bfem);
151  base_tensor tt; c0.base_value(tt);
152  base_tensor::const_iterator itf = tt.begin();
153 
154  std::vector<bool> zid;
155  find_zone_id(c, zid, c.xfem_side());
156  for (dim_type q = 0; q < target_dim(); ++q) {
157  unsigned cnt = 0;
158  for (size_type d = 0; d < bfem->nb_dof(0); ++d, ++itf) {
159  if (dofzones[d]) { /* enriched dof ? */
160  for (size_type k = 0; k < dofzones[d]->size(); ++k, ++cnt)
161  *it++ = zid[cnt] ? *itf : 0;
162  } else *it++ = *itf;
163  }
164  }
165  assert(it == t.end());
166  }
167 
169  base_tensor &t, bool) const {
170  // bgeot::multi_index mi(3);
171  // mi[2] = short_type(c.N()); mi[1] = target_dim();
172  // mi[0] = short_type(nb_base(0));
173  t.adjust_sizes(nb_base(0), target_dim(), c.N());
175  if (c0.have_pfp())
176  c0.set_pfp(fem_precomp(bfem, c0.pfp()->get_ppoint_tab(), c0.pfp()));
177  else c0.set_pf(bfem);
178  base_tensor tt; c0.grad_base_value(tt);
179 
180  base_tensor::iterator it = t.begin();
181  base_tensor::const_iterator itf = tt.begin();
182 
183  std::vector<bool> zid;
184  find_zone_id(c, zid, c.xfem_side());
185 
186  for (dim_type i = 0; i < c.N() ; ++i) {
187  for (dim_type q = 0; q < target_dim(); ++q) {
188  unsigned cnt = 0;
189  for (size_type d = 0; d < bfem->nb_dof(0); ++d, ++itf) {
190  if (dofzones[d]) { /* enriched dof ? */
191  for (size_type k = 0; k < dofzones[d]->size(); ++k, ++cnt)
192  *it++ = zid[cnt] ? *itf : 0;
193  } else *it++ = *itf;
194  }
195  }
196  }
197  assert(it == t.end());
198  }
199 
201  base_tensor &t, bool) const {
202  t.adjust_sizes(nb_base(0), target_dim(), gmm::sqr(c.N()));
204  if (c0.have_pfp())
205  c0.set_pfp(fem_precomp(bfem, c0.pfp()->get_ppoint_tab(), c0.pfp()));
206  else c0.set_pf(bfem);
207  base_tensor tt; c0.hess_base_value(tt);
208 
209  base_tensor::iterator it = t.begin();
210  base_tensor::const_iterator itf = tt.begin();
211 
212  std::vector<bool> zid;
213  find_zone_id(c, zid, c.xfem_side());
214 
215  dim_type NNdim = dim_type(gmm::sqr(c.N())*target_dim());
216  for (dim_type ijq = 0; ijq < NNdim ; ++ijq) {
217  unsigned cnt = 0;
218  for (size_type d = 0; d < bfem->nb_dof(0); ++d, ++itf) {
219  if (dofzones[d]) /* enriched dof ? */
220  for (size_type k = 0; k < dofzones[d]->size(); ++k, ++cnt)
221  *it++ = zid[cnt] ? *itf : 0;
222  else
223  *it++ = *itf;
224  }
225  }
226  assert(it == t.end());
227  }
228 
229 
230 } /* end of namespace getfem. */
231 
getfem::fem_level_set::hess_base_value
void hess_base_value(const base_node &x, base_tensor &t) const
Give the value of all hessians (on ref.
Definition: getfem_fem_level_set.cc:84
getfem_fem_level_set.h
FEM associated with getfem::mesh_fem_level_set objects.
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::fem_level_set::real_grad_base_value
void real_grad_base_value(const fem_interpolation_context &c, base_tensor &t, bool=true) const
Give the gradient of all components of the base functions at the current point of the fem_interpolati...
Definition: getfem_fem_level_set.cc:168
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
getfem::fem_interpolation_context::hess_base_value
void hess_base_value(base_tensor &t, bool withM=true) const
fill the tensor with the hessian of the base functions (taken at point this->xref())
Definition: getfem_fem.cc:255
getfem::mesh_level_set::nb_level_sets
size_type nb_level_sets(void) const
Get number of level-sets referenced in this object.
Definition: getfem_mesh_level_set.h:86
getfem::fem_level_set::base_value
void base_value(const base_node &x, base_tensor &t) const
Give the value of all components of the base functions at the point x of the reference element.
Definition: getfem_fem_level_set.cc:79
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::virtual_fem::add_node
void add_node(const pdof_description &d, const base_node &pt, const dal::bit_vector &faces)
internal function adding a node to an element for the creation of a finite element method.
Definition: getfem_fem.cc:648
getfem::fem_level_set::grad_base_value
void grad_base_value(const base_node &x, base_tensor &t) const
Give the value of all gradients (on ref.
Definition: getfem_fem_level_set.cc:81
getfem::virtual_fem::target_dim
dim_type target_dim() const
dimension of the target space.
Definition: getfem_fem.h:314
getfem::virtual_fem::nb_base
virtual size_type nb_base(size_type cv) const
Number of basis functions.
Definition: getfem_fem.h:299
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_interpolation_context::pfp
pfem_precomp pfp() const
get the current fem_precomp_
Definition: getfem_fem.h:794
getfem::slice_vector_on_basic_dof_of_element
void slice_vector_on_basic_dof_of_element(const mesh_fem &mf, const VEC1 &vec, size_type cv, VEC2 &coeff, size_type qmult1=size_type(-1), size_type qmult2=size_type(-1))
Given a mesh_fem.
Definition: getfem_mesh_fem.h:659
getfem::fem_level_set::real_hess_base_value
void real_hess_base_value(const fem_interpolation_context &c, base_tensor &t, bool=true) const
Give the hessian of all components of the base functions at the current point of the fem_interpolatio...
Definition: getfem_fem_level_set.cc:200
getfem::fem_interpolation_context::grad_base_value
void grad_base_value(base_tensor &t, bool withM=true) const
fill the tensor with the gradient of the base functions (taken at point this->xref())
Definition: getfem_fem.cc:202
getfem::xfem_dof
pdof_description xfem_dof(pdof_description p, size_type ind)
Description of a special dof for Xfem.
Definition: getfem_fem.cc:433
getfem::fem_interpolation_context::base_value
void base_value(base_tensor &t, bool withM=true) const
fill the tensor with the values of the base functions (taken at point this->xref())
Definition: getfem_fem.cc:120
getfem::fem_level_set::real_base_value
void real_base_value(const fem_interpolation_context &c, base_tensor &t, bool=true) const
Give the value of all components of the base functions at the current point of the fem_interpolation_...
Definition: getfem_fem_level_set.cc:141