Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
discontinuousinterfaces.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2009-01-23
7 
8  Copyright (C) 2009 Universite Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __DiscontinuousInterfaces_H
30 #define __DiscontinuousInterfaces_H 1
31 
33 
34 namespace Feel
35 {
43 template<typename A0>
44 class DiscontinuousInterfaces : public detail::continuity_base
45 {
46 public:
47 
48 
52 
53  static const bool is_continuous = false;
54  static const bool is_discontinuous_locally = true;
55  static const bool is_discontinuous_totally = false;
56 
57  static const uint16_type n_discontinuities = fusion::result_of::size<A0>::type::value;
58 
59 
60 
61 
63 
67  typedef A0 discontinuity_markers_type;
68 
70 
74 
77  :
78  M_d_faces()
79  {}
80 
83  :
84  M_d_faces( d.M_d_faces )
85  {}
86 
89  {}
90 
92 
96 
99  {
100  if ( this != &d )
101  {
102  M_d_faces = d.M_d_faces;
103  }
104 
105  return *this;
106  }
107 
109 
113 
114  discontinuity_markers_type const& discontinuityMarkers() const
115  {
116  return M_d_faces;
117  }
118 
119 
121 
125 
126 
128 
132 
133 
135 
136 
137 
138  template<typename MeshType, typename DofType>
139  class apply
140  {
141  public:
142  typedef size_type result_type;
143  typedef MeshType mesh_type;
144  typedef DofType dof_type;
145  typedef typename dof_type::fe_type fe_type;
146 
147  apply( MeshType& M, DofType& D )
148  :
149  M_mesh( M ),
150  M_dof( D )
151  {}
152  template<typename T>
153 #if BOOST_VERSION < 104200
154  result_type operator()( const T& t, const size_type& start )
155 #else
156  result_type operator()( const size_type& start, const T& t )
157 #endif
158  {
159  boost::tuple<size_type,size_type,size_type> disc = boost::make_tuple( mpl::at<T,mpl::int_<0> >::type::value,
160  mpl::at<T,mpl::int_<1> >::type::value,
161  mpl::at<T,mpl::int_<2> >::type::value );
162 
163  return build( disc, start );
164  }
165  private:
166 
167  // -loop on faces marked as discontinuous
168  // - verify that the face is connected to two elements
169  // - for each element associated to the face
170  // - extract the id of the element
171  // - loop on dof associated to the face in the current element
172  // - insert dof by "just" increment next_free_dof, be careful that
173  // map_gdof must be updated to ensure that these will _never_ be revisited (don't use insertDof)
174  // - return the next_free_dof
175  size_type
176  build( boost::tuple<size_type,size_type,size_type> const& marker, size_type start )
177  {
178  size_type next_free_dof = start;
179  size_type n_dof = 0;
180 
181  typedef typename mesh_type::element_type element_type;
182  //typedef typename mesh_type::marker_element_iterator element_iterator;
183  typedef typename mesh_type::element_const_iterator element_const_iterator;
184 
185  element_const_iterator fit, fen;
186  //boost::tie( fit, fen ) = M_mesh.elementsWithMarker( boost::get<1>( marker ), M_mesh.rank() );
187  boost::tie( fit, fen ) = M_mesh.elementsRange();
188  DVLOG(2) << "[DiscontinuousInterfaces::build] n_elements = " << std::distance( fit, fen )
189  << " with marker " << boost::get<1>( marker ) << "\n";
190 #if 0
191 
192  while ( fit != fen )
193  {
194  if ( fit->marker().value() != boost::get<1>( marker ) )
195  {
196  ++fit;
197  continue;
198  }
199 
200  DVLOG(2) << "found element with marker " << fit->marker().value() << "\n";
201  typename element_type::face_const_iterator it, en;
202  boost::tie( it, en ) = fit->faces();
203 
204 
205  for ( ; it != en; ++it )
206  {
207  DVLOG(2) << "face with marker " << ( *it )->marker().value() << "\n";
208  //if ( (*it)->marker().value() == boost::get<0>( marker ) )
209  {
210  DVLOG(2) << "------------------------------------------------------------\n";
211  DVLOG(2) << "face " << ( *it )->id()
212  << " marker = " << boost::get<0>( marker )
213  << " elt marker 0 = " << boost::get<1>( marker )
214  << " elt marker 1 = " << boost::get<2>( marker ) << "\n";
215  DVLOG(2) << "element marker " << fit->marker() << "\n";
216 
217  addVertexDof( *fit, *( *it ), next_free_dof, 0, mpl::bool_<fe_type::nDofPerVertex>() );
218  addEdgeDof( *fit, *( *it ), next_free_dof, 0, mpl::bool_<fe_type::nDofPerEdge>(), mpl::int_<mesh_type::nDim>() );
219 
220 
221  }
222  }
223 
224  ++fit;
225  }
226 
227 #else
228  boost::tie( fit, fen ) = M_mesh.elementsRange();
229 
230  while ( fit != fen )
231  {
232  if ( fit->marker().value() == ( int )boost::get<1>( marker ) )
233  {
234  M_dof.addDofFromElement( *fit, next_free_dof, 0 );
235  }
236 
237  ++fit;
238  }
239 
240 #endif
241  //n_dof = next_free_dof-start;
242  n_dof = next_free_dof;
243 
244  DVLOG(2) << "[DiscontinuousInterfaces::build] n_dof = " << n_dof << "\n";
245 
246  //boost::tie( fit, fen ) = M_mesh.elementsWithMarker( boost::get<2>( marker ), M_mesh.rank() );
247  boost::tie( fit, fen ) = M_mesh.elementsRange();
248  DVLOG(2) << "[DiscontinuousInterfaces::build] n_elements = " << std::distance( fit, fen )
249  << " with marker " << boost::get<2>( marker ) << "\n";
250 #if 0
251 
252  while ( fit != fen )
253  {
254  if ( fit->marker().value() != boost::get<2>( marker ) )
255  {
256  ++fit;
257  continue;
258  }
259 
260  DVLOG(2) << "found element with marker " << fit->marker().value() << "\n";
261  typename element_type::face_const_iterator it, en;
262  boost::tie( it, en ) = fit->faces();
263 
264 
265  for ( ; it != en; ++it )
266  {
267  DVLOG(2) << "face with marker " << ( *it )->marker().value() << "\n";
268  //if ( (*it)->marker().value() == boost::get<0>( marker ) )
269  {
270  DVLOG(2) << "------------------------------------------------------------\n";
271  DVLOG(2) << "face " << ( *it )->id()
272  << " marker = " << boost::get<0>( marker )
273  << " elt marker 0 = " << boost::get<1>( marker )
274  << " elt marker 1 = " << boost::get<2>( marker ) << "\n";
275  DVLOG(2) << "element marker " << fit->marker() << "\n";
276 
277  addVertexDof( *fit, *( *it ), next_free_dof, n_dof, mpl::bool_<fe_type::nDofPerVertex>() );
278  addEdgeDof( *fit, *( *it ), next_free_dof, n_dof, mpl::bool_<fe_type::nDofPerEdge>(), mpl::int_<mesh_type::nDim>() );
279 
280 
281  }
282 
283  }
284 
285  ++fit;
286  }
287 
288 #else
289 
290  typedef typename DofType::dof_map_type dof_map_type;
291  dof_map_type m1( M_dof.mapGDof() );
292  M_dof.clearMapGDof();
293  boost::tie( fit, fen ) = M_mesh.elementsRange();
294 
295  while ( fit != fen )
296  {
297  if ( fit->marker().value() == boost::get<2>( marker ) )
298  {
299  M_dof.addDofFromElement( *fit, next_free_dof );
300  }
301 
302  ++fit;
303  }
304 
305 #if 0
306  dof_map_type m2( M_dof.mapGDof() );
307  M_dof.clearMapGDof();
308  std::merge( m1.begin(), m1.end(), m2.begin(), m2.end(), M_dof.mapGDof().begin() );
309 #else
310  dof_map_type m2( M_dof.mapGDof() );
311  M_dof.clearMapGDof();
312  typename dof_map_type::iterator it = m1.begin();
313  typename dof_map_type::iterator en = m1.end();
314 
315  for ( ; it != en; ++it )
316  {
317  M_dof.mapGDof().insert( *it );
318  }
319 
320  it = m2.begin();
321  en = m2.end();
322 
323  for ( ; it != en; ++it )
324  {
325  M_dof.mapGDof().insert( *it );
326  }
327 
328  DVLOG(2) << "size dictionnary = " << M_dof.mapGDof().size() << " next_free_dof = " << next_free_dof+n_dof << "\n";
329 #endif
330 #endif
331 
332  return next_free_dof;
333  }
334  template<typename element_type, typename face_type>
335  void
336  addVertexDof( element_type const& elt, face_type const& face, size_type& next_free_dof, size_type shift, mpl::bool_<false> )
337  {
338 
339  }
340  template<typename element_type, typename face_type>
341  void
342  addVertexDof( element_type const& elt, face_type const& face, size_type& next_free_dof, size_type shift, mpl::bool_<true> )
343  {
344  uint16_type iFaEl = invalid_uint16_type_value;
345 
346  if ( face.ad_first() == elt.id() )
347  iFaEl = face.pos_first();
348 
349  else if ( face.ad_second() == elt.id() )
350  iFaEl = face.pos_second();
351 
352  for ( uint16_type iVeFa = 0; iVeFa < face_type::numVertices; ++iVeFa )
353  {
354  // local vertex number (in element)
355  uint16_type iVeEl = element_type::fToP( iFaEl, iVeFa );
356  Feel::detail::ignore_unused_variable_warning( iVeEl );
357 
358  FEELPP_ASSERT( iVeEl != invalid_uint16_type_value ).error( "invalid local dof" );
359 
360  // Loop number of Dof per vertex
361  for ( uint16_type l = 0; l < fe_type::nDofPerVertex; ++l )
362  {
363  uint16_type lid = iVeEl * fe_type::nDofPerVertex + l;
364  const size_type gDof = ( elt.point( iVeEl ).id() ) * fe_type::nDofPerVertex + l;
365 
366  DVLOG(2) << "add vertex discontinuous dof " << next_free_dof << " in element " << elt.id() << " lid = " << lid << "\n";
367  bool inserted = M_dof.insertDof( elt.id(), lid, iVeEl, boost::make_tuple( 0, 0, gDof ), 0, next_free_dof, 1, false, shift );
368 
369  if ( shift )
370  {
371  FEELPP_ASSERT( inserted == false )( elt.id() )
372  ( lid )( gDof )( next_free_dof ).error( "should have inserted unique dof" );
373  }
374 
375  DVLOG(2) << "vertex discontinuous dof inserted : " << inserted << "\n";
376 
377  DVLOG(2) << "added vertex discontinuous dof " << elt.id() << ", "
378  << lid << ", "
379  << boost::get<0>( M_dof.localToGlobal( elt.id(), lid, 0 ) ) << "\n";
380  }
381 
382  }
383 
384 
385  }
386  template<typename element_type, typename face_type>
387  void
388  addEdgeDof( element_type const& elt, face_type const& face, size_type& next_free_dof, size_type shift, mpl::bool_<false>, mpl::int_<1> )
389  {}
390  template<typename element_type, typename face_type>
391  void
392  addEdgeDof( element_type const& elt, face_type const& face, size_type& next_free_dof, size_type shift, mpl::bool_<true>, mpl::int_<1> )
393  {}
394 
395  template<typename element_type, typename face_type>
396  void
397  addEdgeDof( element_type const& elt, face_type const& face, size_type& next_free_dof, size_type shift, mpl::bool_<false>, mpl::int_<2> )
398  {}
399  template<typename element_type, typename face_type>
400  void
401  addEdgeDof( element_type const& elt, face_type const& face, size_type& next_free_dof, size_type shift, mpl::bool_<true>, mpl::int_<2> )
402  {
403  uint16_type iFaEl = invalid_uint16_type_value;
404 
405  if ( face.ad_first() == elt.id() )
406  iFaEl = face.pos_first();
407 
408  else if ( face.ad_second() == elt.id() )
409  iFaEl = face.pos_second();
410 
411  // Loop number of Dof per edge
412  for ( uint16_type l = 0; l < fe_type::nDofPerEdge; ++l )
413  {
414  uint16_type lid = element_type::numVertices*fe_type::nDofPerVertex + iFaEl * fe_type::nDofPerEdge + l;
415  const size_type gDof = ( elt.edge( iFaEl ).id() ) * fe_type::nDofPerEdge + l;
416 
417  DVLOG(2) << "add edge discontinuous dof " << next_free_dof << " in element " << elt.id() << " lid = " << lid << "\n";
418  bool inserted = M_dof.insertDof( elt.id(), lid, iFaEl, boost::make_tuple( 1, 0, gDof ), 0, next_free_dof, 1, false, shift );
419  DVLOG(2) << "edge discontinuous dof inserted (1 or 0) : " << inserted << "\n";
420 
421  DVLOG(2) << "added edge discontinuous dof "
422  << elt.id() << ", "
423  << lid << ", "
424  << boost::get<0>( M_dof.localToGlobal( elt.id(), lid, 0 ) ) << "\n";
425 
426  }
427 
428  }
429  private:
430  MeshType& M_mesh;
431  DofType& M_dof;
432  };
433 
434 private:
435  discontinuity_markers_type M_d_faces;
436 };
437 } // Feel
438 #endif /* __DiscontinuousInterfaces_H */

Generated on Sun Oct 20 2013 08:24:56 for Feel++ by doxygen 1.8.4