GetFEM  5.4.2
getfem_partial_mesh_fem.cc
1 /*===========================================================================
2 
3  Copyright (C) 2006-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 
23 
24 namespace getfem {
25 
26 
27  void partial_mesh_fem::clear(void)
28  { mesh_fem::clear(); is_adapted = false; }
29 
30  partial_mesh_fem::partial_mesh_fem(const mesh_fem &mef)
31  : mesh_fem(mef.linked_mesh()), mf(mef)
32  { is_adapted = false; }
33 
34  static getfem::mesh void_mesh__;
35 
36  partial_mesh_fem::partial_mesh_fem(const mesh_fem *mef)
37  : mesh_fem(*(mef ? &(mef->linked_mesh()) : &(void_mesh__))), mf(*mef)
38  { is_adapted = false; }
39 
40  DAL_SIMPLE_KEY(special_partialmf_key, pfem);
41  void partial_mesh_fem::adapt(const dal::bit_vector &kept_dofs,
42  const dal::bit_vector &rejected_elt) {
43  mf.context_check();
44 
45  if (!(mi.is_equal(mf.get_qdims()))) {
46  mi = mf.get_qdims();
47  Qdim = mf.get_qdim();
48  dof_enumeration_made = false; touch(); v_num = act_counter();
49  }
50 
51  fe_convex = mf.convex_index();
52  fe_convex.setminus(rejected_elt);
53 
54  gmm::row_matrix<gmm::rsvector<scalar_type> >
55  RR(kept_dofs.card(), mf.nb_dof());
56  size_type j = 0;
57  for (dal::bv_visitor i(kept_dofs); !i.finished(); ++i, ++j)
58  RR(j, i) = scalar_type(1);
59 
60  R_ = REDUCTION_MATRIX(kept_dofs.card(), mf.nb_basic_dof());
61  E_ = EXTENSION_MATRIX(mf.nb_basic_dof(), kept_dofs.card());
62 
63  if (mf.is_reduced()) {
64  gmm::row_matrix<gmm::rsvector<scalar_type> >
65  A(kept_dofs.card(), mf.nb_basic_dof());
66  gmm::mult(RR, mf.reduction_matrix(), A);
67  gmm::copy(A, R_);
68  gmm::row_matrix<gmm::rsvector<scalar_type> >
69  B(mf.nb_basic_dof(), kept_dofs.card());
70  gmm::mult(mf.extension_matrix(), gmm::transposed(RR), B);
71  gmm::copy(B, E_);
72  }
73  else {
74  gmm::copy(RR, R_); gmm::copy(gmm::transposed(RR), E_);
75  }
76  use_reduction = true;
77 
78  is_adapted = true; touch(); v_num = act_counter();
79  }
80 
81  // invalid function for a mesh change.
82  // dal::bit_vector partial_mesh_fem::retrieve_kept_dofs() const
83  // {
84  // base_vector full(nb_basic_dof());
85  // for (size_type i = 0; i < full.size(); ++i) full[i] = i;
86  // base_vector reduced(nb_dof());
87  //
88  // if (R_.ncols() > 0) gmm::mult(R_, full, reduced);
89  // else reduced = full;
90  //
91  // dal::bit_vector kept_dofs;
92  // for (size_type i=0; i < reduced.size(); ++i) kept_dofs.add(reduced[i]);
93  //
94  // return kept_dofs;
95  // }
96 
97  void partial_mesh_fem::write_to_file(std::ostream &ost) const
98  { context_check(); mf.context_check();
99  gmm::stream_standard_locale sl(ost);
100  ost << '\n' << "BEGIN MESH_FEM" << '\n' << '\n';
101  mf.write_basic_to_file(ost);
102  write_reduction_matrices_to_file(ost);
103  ost << "END MESH_FEM" << '\n';
104  }
105 
106  void partial_mesh_fem::write_to_file(const std::string &name,
107  bool with_mesh) const {
108  std::ofstream o(name.c_str());
109  GMM_ASSERT1(o, "impossible to open file '" << name << "'");
110  o << "% GETFEM MESH_FEM FILE " << '\n';
111  o << "% GETFEM VERSION " << GETFEM_VERSION << '\n' << '\n' << '\n';
112  if (with_mesh) mf.linked_mesh().write_to_file(o);
113  write_to_file(o);
114  }
115 
116  dal::bit_vector select_dofs_from_im(const mesh_fem &mf, const mesh_im &mim,
117  unsigned P) {
118  const mesh &m = mf.linked_mesh();
119  unsigned N = m.dim();
120  if (P == unsigned(-1)) P = N;
121  base_matrix G;
122  bgeot::pgeometric_trans pgt_old = 0;
123  bgeot::pgeotrans_precomp pgp2 = 0;
124  getfem::pfem pf_old = 0;
125  getfem::pfem_precomp pfp = 0;
126  pintegration_method pim1 = 0;
127 
128  std::vector<scalar_type> areas(mf.nb_basic_dof());
129  std::vector<scalar_type> area_supports(mf.nb_basic_dof());
130  dal::bit_vector kept_dofs;
131 
132  for (dal::bv_visitor cv(mim.convex_index()); !cv.finished(); ++cv) {
133  m.points_of_convex(cv, G);
134  bgeot::pgeometric_trans pgt = m.trans_of_convex(cv);
135  pintegration_method pim = mim.int_method_of_element(cv);
136  if (pim == im_none()) continue;
137  getfem::pfem pf = mf.fem_of_element(cv);
138  GMM_ASSERT1(pim->type() == IM_APPROX,
139  "Works only with approximate integration");
140  papprox_integration pai2= pim->approx_method();
141  static papprox_integration pai2_old = 0;
142  if (pgt_old != pgt || pai2 != pai2_old) {
143  pim1 = getfem::classical_approx_im(pgt, 2);
144  pgp2 = bgeot::geotrans_precomp(pgt, pai2->pintegration_points(),pim);
145  }
146  if (pai2 != pai2_old || pf != pf_old) {
147  pf_old = pf;
148  pfp = getfem::fem_precomp(pf, pai2->pintegration_points(), pim);
149  }
150  pai2_old = pai2;
151  pgt_old = pgt;
152 
154  scalar_type area1 = convex_area_estimate(pgt, G, pim1);
155 
156  size_type tdim = mf.get_qdim() / pf->target_dim();
157 
158  for (size_type i = 0; i < pai2->nb_points_on_convex(); ++i) {
159  for (unsigned d = 0; d < pf->nb_dof(cv); ++d) {
160  for (size_type j = 0; j < tdim; ++j) {
161  if (i == 0)
162  areas[mf.ind_basic_dof_of_element(cv)[d*tdim+j]] += area1;
163  c2.set_ii(i);
164  area_supports[mf.ind_basic_dof_of_element(cv)[d*tdim+j]]
165  += pai2->coeff(i) * c2.J() * gmm::sqr(pfp->val(i)[d]);
166  }
167  // * ((gmm::abs(pfp->val(i)[d]) < 1e-10) ? 0.0 : 1.0);
168  }
169  }
170  }
171 
172 
173  std::vector<scalar_type> areas2(mf.nb_dof());
174  std::vector<scalar_type> area_supports2(mf.nb_dof());
175 
176  if (mf.is_reduced()) {
177  gmm::mult(gmm::transposed(mf.extension_matrix()), areas, areas2);
178  gmm::mult(gmm::transposed(mf.extension_matrix()), area_supports,
179  area_supports2);
180  }
181  else {
182  gmm::copy(areas, areas2);
183  gmm::copy(area_supports, area_supports2);
184  }
185 
186  for (size_type i = 0; i < mf.nb_dof(); ++i) {
187  if (area_supports2[i] > pow(1e-14 * areas2[i], scalar_type(P) / N))
188  kept_dofs.add(i);
189  }
190 
191  return kept_dofs;
192  }
193 
194 
195 
196 } /* end of namespace getfem. */
197 
getfem::partial_mesh_fem::adapt
void adapt(const dal::bit_vector &kept_dof, const dal::bit_vector &rejected_elt=dal::bit_vector())
build the mesh_fem keeping only the dof of the original mesh_fem which are listed in kept_dof.
Definition: getfem_partial_mesh_fem.cc:41
getfem::mesh_fem::get_qdim
virtual dim_type get_qdim() const
Return the Q dimension.
Definition: getfem_mesh_fem.h:312
getfem::mesh::write_to_file
void write_to_file(const std::string &name) const
Write the mesh to a file.
Definition: getfem_mesh.cc:709
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::partial_mesh_fem::write_to_file
void write_to_file(std::ostream &ost) const
Write the mesh_fem to a stream.
Definition: getfem_partial_mesh_fem.cc:97
getfem::mesh_im::convex_index
const dal::bit_vector & convex_index(void) const
Get the set of convexes where an integration method has been assigned.
Definition: getfem_mesh_im.h:73
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
getfem::convex_area_estimate
scalar_type APIDECL convex_area_estimate(bgeot::pgeometric_trans pgt, const base_matrix &pts, pintegration_method pim)
rough estimate of the convex area.
Definition: getfem_mesh.cc:745
bgeot::geotrans_interpolation_context::J
scalar_type J() const
get the Jacobian of the geometric trans (taken at point xref() )
Definition: bgeot_geometric_trans.h:453
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::mesh_fem::fem_of_element
virtual pfem fem_of_element(size_type cv) const
Return the basic fem associated with an element (if no fem is associated, the function will crash!...
Definition: getfem_mesh_fem.h:444
bgeot::geotrans_interpolation_context
the geotrans_interpolation_context structure is passed as the argument of geometric transformation in...
Definition: bgeot_geometric_trans.h:411
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
getfem_partial_mesh_fem.h
a subclass of getfem::mesh_fem which allows to eliminate a number of dof of the original mesh_fem.
getfem::mesh_fem::extension_matrix
const EXTENSION_MATRIX & extension_matrix() const
Return the extension matrix corresponding to reduction applied (RE=I).
Definition: getfem_mesh_fem.h:201
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::mesh_fem::nb_basic_dof
virtual size_type nb_basic_dof() const
Return the total number of basic degrees of freedom (before the optional reduction).
Definition: getfem_mesh_fem.h:556
bgeot::geotrans_interpolation_context::set_ii
void set_ii(size_type ii__)
change the current point (assuming a geotrans_precomp_ is used)
Definition: bgeot_geometric_trans.h:463
getfem::select_dofs_from_im
dal::bit_vector select_dofs_from_im(const mesh_fem &mf, const mesh_im &mim, unsigned P=unsigned(-1))
Return a selection of dof who contribute significantly to the mass-matrix that would be computed with...
Definition: getfem_partial_mesh_fem.cc:116
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::im_none
pintegration_method im_none(void)
return IM_NONE
Definition: getfem_integration.cc:1293
getfem::classical_approx_im
pintegration_method classical_approx_im(bgeot::pgeometric_trans pgt, dim_type degree)
try to find an approximate integration method for the geometric transformation pgt which is able to i...
Definition: getfem_integration.cc:1280
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
getfem::mesh_fem::ind_basic_dof_of_element
virtual ind_dof_ct ind_basic_dof_of_element(size_type cv) const
Give an array of the dof numbers a of convex.
Definition: getfem_mesh_fem.h:450
getfem::mesh_im::int_method_of_element
virtual pintegration_method int_method_of_element(size_type cv) const
return the integration method associated with an element (in no integration is associated,...
Definition: getfem_mesh_im.h:117
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::reduction_matrix
const REDUCTION_MATRIX & reduction_matrix() const
Return the reduction matrix applied to the dofs.
Definition: getfem_mesh_fem.h:198
getfem::mesh_fem::nb_dof
virtual size_type nb_dof() const
Return the total number of degrees of freedom.
Definition: getfem_mesh_fem.h:562