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
geoelement.hpp
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  Copyright (C) 2001,2002,2003,2004 EPFL, INRIA and Politechnico di Milano
6  Copyright (C) 2006-2010 Université Joseph Fourier (Grenoble I)
7 
8  This library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU Lesser General Public
10  License as published by the Free Software Foundation; either
11  version 3.0 of the License, or (at your option) any later version.
12 
13  This library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  Lesser General Public License for more details.
17 
18  You should have received a copy of the GNU Lesser General Public
19  License along with this library; if not, write to the Free Software
20  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21 */
28 #ifndef _GEOELEMENT_HH_
29 #define _GEOELEMENT_HH_
30 
31 #include <boost/optional.hpp>
32 #include <boost/tuple/tuple.hpp>
33 
34 #include <feel/feelcore/feel.hpp>
35 #include <feel/feelalg/matrix.hpp>
36 #include <feel/feelmesh/marker.hpp>
37 #include <feel/feelmesh/geond.hpp>
38 
39 #include <feel/feelalg/lu.hpp>
40 
41 namespace Feel
42 {
43 
44 class SubFaceOfNone
45 {
46 public:
47  static const uint16_type nDim = 1;
48  template<typename ET>
49  struct Element
50  {
51  typedef ET type;
52  };
53  typedef boost::tuple<size_type, size_type, uint16_type> element_connectivity_type;
54 
55  boost::none_t element( uint16_type /* e */ ) const
56  {
57  return boost::none_t();
58  }
59 
60  SubFaceOfNone() {}
61  SubFaceOfNone( SubFaceOfNone const& ) {}
62 
63  virtual ~SubFaceOfNone() {}
64 
65  template<typename SFO>
66  SubFaceOfNone( SFO const& /*sf*/ )
67  {
68  }
69  bool
70  isGhostFace( size_type p ) const
71  {
72  return false;
73  }
74  bool
75  isInterProcessDomain( size_type /*p*/ ) const
76  {
77  return false;
78  }
79 private:
80  friend class boost::serialization::access;
81  template<class Archive>
82  void serialize( Archive & ar, const unsigned int version )
83  {
84  }
85 };
86 
87 template<typename ElementType>
88 class SubFaceOf
89 {
90 public:
91  static const uint16_type nDim = ElementType::nDim;
92  template<typename ET>
93  struct Element
94  {
95  typedef ElementType type;
96  };
97  typedef ElementType entity_type;
98  typedef boost::tuple<ElementType const*, size_type, uint16_type, size_type> element_connectivity_type;
99 
100  SubFaceOf()
101  :
104  {}
105 
106  SubFaceOf( element_connectivity_type const& connect0 )
107  :
108  M_element0( connect0 ),
110  {}
111  SubFaceOf( element_connectivity_type const& connect0,
112  element_connectivity_type const& connect1 )
113  :
114  M_element0( connect0 ),
115  M_element1( connect1 )
116  {}
117 
118  SubFaceOf( SubFaceOf const& sf )
119  :
120  M_element0( sf.M_element0 ),
121  M_element1( sf.M_element1 )
122  {
123  }
124  SubFaceOf( SubFaceOfNone const& /*sf*/ )
125  :
128  {
129  }
130  virtual ~SubFaceOf() {}
131 
132  SubFaceOf& operator=( SubFaceOf const& sf )
133  {
134  if ( this != &sf )
135  {
136  M_element0 = sf.M_element0;
137  M_element1 = sf.M_element1;
138 
139  }
140 
141  return *this;
142  }
143  SubFaceOf& operator=( SubFaceOfNone const& /*sf*/ )
144  {
145  return *this;
146  }
147  entity_type const& element( uint16_type e ) const
148  {
149  if ( e == 0 )
150  return *boost::get<0>( M_element0 );
151 
152  else
153  return *boost::get<0>( M_element1 );
154  }
155 
156  entity_type const& element0() const
157  {
158  return *boost::get<0>( M_element0 );
159  }
160  entity_type const& element1() const
161  {
162  return *boost::get<0>( M_element1 );
163  }
164 
165  size_type ad_first() const
166  {
167  return boost::get<1>( M_element0 );
168  }
169  uint16_type pos_first() const
170  {
171  return boost::get<2>( M_element0 );
172  }
173  size_type proc_first() const
174  {
175  return boost::get<3>( M_element0 );
176  }
177 
178  size_type ad_second() const
179  {
180  return boost::get<1>( M_element1 );
181  }
182  uint16_type pos_second() const
183  {
184  return boost::get<2>( M_element1 );
185  }
186  size_type proc_second() const
187  {
188  return boost::get<3>( M_element1 );
189  }
190 
191  element_connectivity_type const& connection0() const
192  {
193  return M_element0;
194  }
195  element_connectivity_type const& connection1() const
196  {
197  return M_element1;
198  }
199 
200  void setConnection( uint16_type f, element_connectivity_type const& connect )
201  {
202  if ( f == 0 )
203  M_element0 = connect;
204 
205  else
206  M_element1 = connect;
207 
208  }
209 
210  void setConnection0( element_connectivity_type const& connect )
211  {
212  M_element0 = connect;
213  }
214  void setConnection1( element_connectivity_type const& connect )
215  {
216  M_element1 = connect;
217  }
218 
219  bool isConnected() const { return isConnectedTo0() && isConnectedTo1(); }
220 
221  bool isConnectedTo0() const
222  {
223  return ( boost::get<1>( M_element0 ) != invalid_size_type_value &&
224  boost::get<2>( M_element0 ) != invalid_uint16_type_value &&
225  boost::get<3>( M_element0 ) != invalid_size_type_value );
226  }
227  bool isConnectedTo1() const
228  {
229  return ( boost::get<1>( M_element1 ) != invalid_size_type_value &&
230  boost::get<2>( M_element1 ) != invalid_uint16_type_value &&
231  boost::get<3>( M_element1 ) != invalid_size_type_value );
232  }
233 
234  bool
235  isGhostFace( size_type p ) const
236  {
237  return ( ( boost::get<3>( M_element1 ) != invalid_size_type_value ) &&
238  ( ( ( boost::get<3>( M_element0 ) == p ) && ( boost::get<3>( M_element1 ) < p ) ) ||
239  ( ( boost::get<3>( M_element0 ) < p ) && ( boost::get<3>( M_element1 ) == p ) ) ) );
240  }
241 
242  bool
243  isInterProcessDomain( size_type p ) const
244  {
245  return ( ( boost::get<3>( M_element1 ) != invalid_size_type_value ) &&
246  ( ( boost::get<3>( M_element0 ) == p ) || ( boost::get<3>( M_element1 ) == p ) ) &&
247  ( boost::get<3>( M_element0 ) != boost::get<3>( M_element1 ) ) );
248  }
249  bool
250  isIntraProcessDomain( size_type p ) const
251  {
252  return ( ( boost::get<3>( M_element0 ) == p ) &&
253  ( boost::get<3>( M_element1 ) == p ) );
254  }
255 
256  void disconnect0()
257  {
258  M_element0 = boost::make_tuple( ( ElementType const* )0,
262  }
263 
264  void disconnect1()
265  {
266  M_element1 = boost::make_tuple( ( ElementType const* )0,
270  }
271 
272  void disconnect()
273  {
274  disconnect0();
275  disconnect1();
276  }
277 
278  void disconnect( ElementType const& elem )
279  {
280  if(boost::get<0>( M_element0 ) == boost::addressof(elem))
281  {
282  DVLOG(2) << "connecting 1 to 0 and disconnecting 1..\n";
283  M_element0 = M_element1;
284  disconnect1();
285  }
286  else
287  {
288  DVLOG(2) << "disconnecting 1..\n";
289  disconnect1();
290  }
291  }
292 
293 private:
294  friend class boost::serialization::access;
295  template<class Archive>
296  void serialize( Archive & ar, const unsigned int version )
297  {
298 #if 1
299  //ar & M_element0.template get<0>();
300  ar & M_element0.template get<1>();
301  ar & M_element0.template get<2>();
302  ar & M_element0.template get<3>();
303  M_element0.template get<0>() = 0;
304 
305 
306  //ar & M_element1.template get<0>();
307  ar & M_element1.template get<1>();
308  ar & M_element1.template get<2>();
309  ar & M_element1.template get<3>();
310  M_element1.template get<0>() = 0;
311 #endif
312  }
313 private:
314 
315  element_connectivity_type M_element0;
316  element_connectivity_type M_element1;
317 
318 };
319 #if 0
320 template<typename ElementType>
321 class SubFaceOfMany
322 {
323 public:
324  static const uint16_type nDim = ElementType::nDim;
325  static const uint16_type nRealDim = ElementType::nRealDim;
326  template<typename ET>
327  struct Element
328  {
329  typedef ElementType type;
330  };
331  typedef ElementType entity_type;
332  typedef boost::tuple<ElementType const*, size_type, uint16_type, size_type> element_connectivity_type;
333 
334  SubFaceOfMany()
335  :
336  M_elements()
337  {}
338 
339  SubFaceOfMany( SubFaceOfMany const& sf )
340  :
341  M_elements( sf.M_elements )
342  {
343  }
344  SubFaceOfMany( SubFaceOfNone const& /*sf*/ )
345  :
346  M_elements()
347  {
348  }
349  virtual ~SubFaceOfMany() {}
350 
351  SubFaceOfMany& operator=( SubFaceOfMany const& sf )
352  {
353  if ( this != &sf )
354  {
355  M_elements = sf.M_elements;
356  }
357 
358  return *this;
359  }
360  SubFaceOfMany& operator=( SubFaceOfNone const& /*sf*/ )
361  {
362  return *this;
363  }
364 
365  entity_type const& element0() const
366  {
367  return *boost::get<0>( *M_elements.begin() );
368  }
369  entity_type const& element1() const
370  {
371  return *boost::get<0>( *boost::next(M_elements.begin()) );
372  }
373 
374  void setConnection( element_connectivity_type const& connect )
375  {
376  this->insert( connect );
377  }
378 
379  size_type ad_first() const
380  {
381  return boost::get<1>( *M_elements.begin() );
382  }
383  uint16_type pos_first() const
384  {
385  return boost::get<2>( *M_elements.begin() );
386  }
387  size_type proc_first() const
388  {
389  return boost::get<3>( *M_elements.begin() );
390  }
391 
392  size_type ad_second() const
393  {
394  return boost::get<1>( *boost::next(M_elements.begin()) );
395  }
396  uint16_type pos_second() const
397  {
398  return boost::get<2>( *boost::next(M_elements.begin()) );
399  }
400  size_type proc_second() const
401  {
402  return boost::get<3>( *boost::next(M_elements.begin()) );
403  }
404 
405 
406  void setConnection0( element_connectivity_type const& connect )
407  {
408  M_elements.insert( connect );
409  }
410  void setConnection1( element_connectivity_type const& connect )
411  {
412  M_elements.insert( connect );
413  }
414 
415  element_connectivity_type const& connection0() const
416  {
417  return *M_elements.begin();
418  }
419  element_connectivity_type const& connection1() const
420  {
421  return *boost::next(M_elements.begin());
422  }
423 
424  bool isConnectedTo0() const
425  {
426  return ( boost::get<1>( *M_elements.begin() ) != invalid_size_type_value &&
427  boost::get<2>( *M_elements.begin() ) != invalid_uint16_type_value &&
428  boost::get<3>( *M_elements.begin() ) != invalid_size_type_value );
429  }
430  bool isConnectedTo1() const
431  {
432  return ( boost::get<1>( *boost::next(M_elements.begin()) ) != invalid_size_type_value &&
433  boost::get<2>( *boost::next(M_elements.begin()) ) != invalid_uint16_type_value &&
434  boost::get<3>( *boost::next(M_elements.begin()) ) != invalid_size_type_value );
435  }
436  bool
437  isInterProcessDomain( size_type p ) const
438  {
439  return false;
440  }
441  bool
442  isIntraProcessDomain( size_type p ) const
443  {
444  return true;
445  }
446 
447  entity_type const& element( uint16_type e ) const
448  {
449  return *boost::get<0>( *M_elements.begin() );
450  }
451  void disconnect()
452  {
453  M_elements.clear();
454  }
455 
456 private:
457  friend class boost::serialization::access;
458  template<class Archive>
459  void serialize( Archive & ar, const unsigned int version )
460  {
461  }
462 private:
463 
464  std::set<element_connectivity_type> M_elements;
465 };
466 #else
467 //template<typename ElementType> class SubFaceOfMany: public SubFaceOf<ElementType> {};
468 #define SubFaceOfMany SubFaceOf
469 
470 #endif
471 
472 // *********** Geometrical Elements *****************
476 
480 template <uint16_type Dim,
481  typename SubFace = SubFaceOfNone,
482  typename T = double>
484  :
485  //public GeoND<Dim, GEOSHAPE, T, GeoElement0D<Dim, SubFaceOfNone, T> >,
486 public Geo0D<Dim,T>,
487 public SubFace
488 {
489 public:
490 
491  static const uint16_type nDim = 0;
492  static const uint16_type nRealDim = Dim;
493 
494  typedef Geo0D<Dim,T> geo0d_type;
495  typedef typename geo0d_type::node_type node_type;
496 
497  typedef geo0d_type super;
498  typedef SubFace super2;
499 
501  typedef typename mpl::if_<mpl::equal_to<mpl::int_<SubFace::nDim>, mpl::int_<0> >, mpl::identity<self_type>, mpl::identity<typename SubFace::template Element<self_type>::type> >::type::type element_type;
502  typedef self_type point_type;
503 
504  static const uint16_type numLocalVertices = super::numVertices;
505 
506  GeoElement0D()
507  :
508  super(),
509  super2(),
510  M_facept()
511  {}
512 
513 
515  GeoElement0D( size_type id, bool boundary = false )
516  :
517  super( id, boundary ),
518  super2(),
519  M_facept()
520  {}
521 
522  GeoElement0D( size_type id, node_type const& n, bool boundary = false )
523  :
524  super( id, n, boundary ),
525  super2(),
526  M_facept()
527  {}
528 
531  GeoElement0D( size_type id, Real x, Real y, Real z, bool boundary = false )
532  :
533  super( id, x, y, z, boundary ),
534  super2(),
535  M_facept()
536  {}
537 
538  template<typename SF>
540  :
541  super( g ),
542  super2( g ),
543  M_facept( g.M_facept )
544  {
545  }
546 
547  ~GeoElement0D()
548  {}
549 
550  template<typename SF>
551  GeoElement0D & operator = ( GeoElement0D<Dim,SF,T> const & g )
552  {
553  super::operator=( g );
554  super2::operator=( g );
555  M_facept = g.M_facept;
556  return *this;
557  }
558 
559  //void setMesh( MeshBase const* m ) { super::setMesh( m ); }
560  MeshBase const* mesh() const
561  {
562  return super::mesh();
563  }
564 
568  size_type id() const
569  {
570  return super::id();
571  }
572 
576  uint16_type processId() const
577  {
578  return super::processId();
579  }
580 
581  bool isGhostFace() const
582  {
583  return super2::isGhostFace( super::processId() );
584  }
585 
589  bool isInterProcessDomain() const
590  {
591  return super2::isInterProcessDomain( super::processId() );
592  }
593 
597  bool isOnBoundary() const
598  {
599  return super::isOnBoundary();
600  }
601 
605  uint16_type boundaryEntityDimension() const
606  {
607  return 0;
608  }
609 
613  bool isGhostCell() const
614  {
615  return super::isGhostCell();
616  }
617 
621  geo0d_type const& point( uint16_type /*i*/ ) const
622  {
623  return M_facept;
624  }
625 
629  //void setPoint( uint16_type /*i*/, GeoElement0D<Dim,SubFaceOfNone,T> const& e ) { M_facept = e; }
630  void setPoint( uint16_type /*i*/, geo0d_type const& e )
631  {
632  M_facept = e;
633  }
634 
638  Marker1 const& marker() const
639  {
640  return super::marker();
641  }
645  Marker1& marker()
646  {
647  return super::marker();
648  }
652  Marker2 const& marker2() const
653  {
654  return super::marker2();
655  }
659  Marker2& marker2()
660  {
661  return super::marker2();
662  }
666  Marker3 const& marker3() const
667  {
668  return super::marker3();
669  }
673  Marker3& marker3()
674  {
675  return super::marker3();
676  }
677 
678 //private:
679  geo0d_type M_facept;
680 
681 private:
682 
683  friend class boost::serialization::access;
684  template<class Archive>
685  void serialize( Archive & ar, const unsigned int version )
686  {
687  ar & boost::serialization::base_object<super>( *this );
688  ar & boost::serialization::base_object<super2>( *this );
689  }
690 
691 
692 };
693 
694 
695 
696 
704 template<uint16_type Dim,
705  typename GEOSHAPE,
706  typename SubFace = SubFaceOfNone,
707  typename T = double>
709  :
710 public GeoND<Dim, GEOSHAPE, T, GeoElement0D<Dim, SubFaceOfNone, T> >,
711 public SubFace
712 {
713 public:
714 
716  typedef SubFace super2;
717 
718  static const uint16_type nDim = super::nDim;
719  static const uint16_type nOrder = super::nOrder;
720  static const uint16_type nRealDim = super::nRealDim;
721 
722  static const bool condition = ( Dim==nRealDim );
723  BOOST_MPL_ASSERT_MSG( ( condition ), INVALID_ELEMENT_REAL_DIMENSION, ( mpl::int_<Dim>, mpl::int_<nRealDim>, GEOSHAPE ) );
724 
725  typedef GEOSHAPE GeoShape;
727  //typedef typename SubFace::template Element<self_type>::type element_type;
728  typedef self_type element_type;
729  typedef typename mpl::if_<mpl::equal_to<mpl::int_<nRealDim>,mpl::int_<1> >,
730  mpl::identity<GeoElement0D<Dim, SubFaceOf<self_type>, T> >,
731  mpl::identity<GeoElement0D<Dim, SubFaceOfMany<self_type>, T> > >::type::type point_type;
732  typedef point_type GeoBElement;
733 
734  static const uint16_type numLocalVertices = super::numLocalVertices;
735  static const uint16_type numLocalEdges = super::numEdges;
736  static const uint16_type numLocalFaces = super::numLocalVertices;
737 
738 
739  typedef typename super::node_type node_type;
740  typedef typename super::vertex_permutation_type vertex_permutation_type;
741  typedef typename super::edge_permutation_type edge_permutation_type;
742  typedef typename super::face_permutation_type face_permutation_type;
743  typedef typename super::permutation_type permutation_type;
744 
749  explicit GeoElement1D( size_type id = 0 )
750  :
751  super( id ),
752  super2(),
753  M_vertices( numLocalVertices, 0 ),
754  M_vertex_permutation( numLocalVertices )
755  {
756  }
761  :
762  super( g ),
763  super2( g ),
764  M_vertices( g.M_vertices ),
765  M_vertex_permutation( g.M_vertex_permutation )
766  {}
767 
772  {}
773 
778  {
779  if ( this != &g )
780  {
781  super::operator=( g );
782  super2::operator=( g );
783  M_vertices = g.M_vertices;
784  M_vertex_permutation = g.M_vertex_permutation;
785  }
786 
787  return *this;
788  }
789 
790 
791  //void setMesh( MeshBase const* m ) { super::setMesh( m ); }
792  MeshBase const* mesh() const
793  {
794  return super::mesh();
795  }
796 
797 
801  size_type id() const
802  {
803  return super::id();
804  }
805 
806  bool isGhostFace() const
807  {
808  return super2::isGhostFace( super::processId() );
809  }
810 
811 
815  bool isInterProcessDomain() const
816  {
817  return super2::isInterProcessDomain( super::processId() );
818  }
819 
823  bool isOnBoundary() const
824  {
825  return super::isOnBoundary();
826  }
827 
831  uint16_type boundaryEntityDimension() const
832  {
833  return 0;
834  }
835 
839  bool isGhostCell() const
840  {
841  return super::isGhostCell();
842  }
843 
847  uint16_type processId() const
848  {
849  return super::processId();
850  }
851 
852  void setMap( uint8_type k_1, uint8_type k_2 )
853  {
854  M_map[k_1] = k_2;
855  }
856 
857  uint8_type map( uint8_type k_1 ) const
858  {
859  return M_map[ k_1 ];
860  }
861 
862  Marker1 const& marker() const
863  {
864  return super::marker();
865  }
866  Marker1& marker()
867  {
868  return super::marker();
869  }
870  Marker2 const& marker2() const
871  {
872  return super::marker2();
873  }
874  Marker3 const& marker3() const
875  {
876  return super::marker3();
877  }
878 
882  void setFace( uint16_type const i, point_type const & p )
883  {
884  FEELPP_ASSERT( i < numLocalVertices )( i ).error( "invalid local point index" );
885  M_vertices[i] = const_cast<point_type*>( boost::addressof( p ) );
886  }
887 
888  edge_permutation_type permutation( uint16_type /*i*/ ) const
889  {
890  return edge_permutation_type();
891  }
892 
893  point_type const& edge( uint16_type i ) const
894  {
895  return *M_vertices[i];
896  }
897  point_type const& face( uint16_type i ) const
898  {
899  return *M_vertices[i];
900  }
901  point_type const* facePtr( uint16_type i ) const
902  {
903  FEELPP_ASSERT( i < numLocalVertices )( this->id() )( i ).error( "invalid local vertex index" );
904  return M_vertices[i];
905  }
906  point_type* facePtr( uint16_type i )
907  {
908  FEELPP_ASSERT( i < numLocalVertices )( this->id() )( i ).error( "invalid local vertex index" );
909  return M_vertices[i];
910  }
911 
912  typedef typename ublas::bounded_array<point_type*, numLocalVertices>::iterator face_iterator;
913  typedef typename ublas::bounded_array<point_type*, numLocalVertices>::const_iterator face_const_iterator;
914 
918  std::pair<face_iterator,face_iterator>
920  {
921  return std::make_pair( M_vertices.begin(), M_vertices.end() );
922  }
923 
927  std::pair<face_const_iterator,face_const_iterator>
928  faces() const
929  {
930  return std::make_pair( M_vertices.begin(), M_vertices.end() );
931  }
932 
936  vertex_permutation_type facePermutation( uint16_type i ) const
937  {
938  FEELPP_ASSERT( i < numLocalVertices )( i )( numLocalVertices ).error( "invalid local vertex index" );
939  return M_vertex_permutation[i];
940  }
941 
942 private:
943 
944  friend class boost::serialization::access;
945  template<class Archive>
946  void serialize( Archive & ar, const unsigned int version )
947  {
948  DVLOG(2) << "Serializing Geoelement1D id: " << this->id() << "...\n";
949  ar & boost::serialization::base_object<super>( *this );
950  ar & boost::serialization::base_object<super2>( *this );
951  }
952 
953 
954 
955 private:
956 
957  std::vector<uint8_type> M_map;
958  ublas::bounded_array<point_type*, numLocalVertices> M_vertices;
959  ublas::bounded_array<vertex_permutation_type, numLocalVertices> M_vertex_permutation;
960 
961 };
962 template<uint16_type Dim,
963  typename GEOSHAPE,
964  typename SubFace,
965  typename T>
966 const uint16_type GeoElement1D<Dim,GEOSHAPE,SubFace,T>::numLocalVertices;
967 
968 template<uint16_type Dim,
969  typename GEOSHAPE,
970  typename SubFace,
971  typename T>
972 const uint16_type GeoElement1D<Dim,GEOSHAPE,SubFace,T>::nDim;
973 
981 template<uint16_type Dim,
982  typename GEOSHAPE,
983  typename SubFace = SubFaceOfNone,
984  typename T = double>
986  :
987 public GeoND<Dim, GEOSHAPE, T, GeoElement0D<Dim, SubFaceOfNone, T> >,
988 public SubFace
989 {
990 public:
991 
992 
994  typedef SubFace super2;
995 
996  static const uint16_type nDim = super::nDim;
997  static const uint16_type nOrder = super::nOrder;
998  static const uint16_type nRealDim = super::nRealDim;
999 
1000  static const bool condition = ( Dim==nRealDim );
1001  BOOST_MPL_ASSERT_MSG( ( condition ), INVALID_ELEMENT_REAL_DIMENSION, ( mpl::int_<Dim>, mpl::int_<nRealDim>, GEOSHAPE ) );
1002 
1004  static const uint16_type numLocalEdges = super::numEdges;
1005  static const uint16_type numLocalFaces = super::numFaces;
1006 
1007  typedef GEOSHAPE GeoShape;
1008  typedef typename super::face_type entity_face_type;
1010  //typedef typename SubFace::template Element<self_type>::type element_type;
1011  typedef self_type element_type;
1012  typedef typename mpl::if_<mpl::equal_to<mpl::int_<nRealDim>,mpl::int_<2> >,
1013  mpl::identity<GeoElement1D<Dim, entity_face_type, SubFaceOf<self_type>, T> >,
1014  mpl::identity<GeoElement1D<Dim, entity_face_type, SubFaceOfMany<self_type>, T> > >::type::type edge_type;
1015  //typedef GeoElement1D<Dim, entity_face_type, SubFaceOf<self_type>, T > edge_type;
1017 #if 0
1018  BOOST_MPL_ASSERT_MSG( ( boost::is_same<point_type,typename edge_type::point_type>::value ),
1019  INCOMPATIBLE_POINT_TYPE,
1020  ( point_type, typename edge_type::point_type, edge_type, element_type, self_type ) );
1021  BOOST_STATIC_ASSERT( ( boost::is_same<point_type,typename edge_type::point_type>::value ) );
1022 #endif // 0
1023  typedef typename super::node_type node_type;
1024 
1025  typedef typename super::vertex_permutation_type vertex_permutation_type;
1026  typedef typename super::edge_permutation_type edge_permutation_type;
1027  typedef typename super::edge_permutation_type permutation_type;
1028  typedef typename super::face_permutation_type face_permutation_type;
1029  typedef typename super2::element_connectivity_type element_connectivity_type;
1030 
1035  explicit GeoElement2D( size_type id = 0 )
1036  :
1037  super( id ),
1038  super2(),
1039  M_edges( numLocalEdges, nullptr ),
1040  M_edge_permutation( numLocalEdges, edge_permutation_type( edge_permutation_type::IDENTITY ) )
1041  {
1042  }
1043 
1048  :
1049  super( g ),
1050  super2( g ),
1051  M_edges( g.M_edges ),
1052  M_edge_permutation( g.M_edge_permutation )
1053  {}
1054 
1059  {}
1060 
1065  {
1066  if ( this != &g )
1067  {
1068  super::operator=( g );
1069  super2::operator=( g );
1070  M_edges = g.M_edges;
1071  M_edge_permutation = g.M_edge_permutation;
1072  }
1073 
1074  return *this;
1075  }
1076 
1077  //void setMesh( MeshBase const* m ) { super::setMesh( m ); }
1078  MeshBase const* mesh() const
1079  {
1080  return super::mesh();
1081  }
1085  size_type id() const
1086  {
1087  return super::id();
1088  }
1089 
1090  Marker1 const& marker() const
1091  {
1092  return super::marker();
1093  }
1094  Marker1& marker()
1095  {
1096  return super::marker();
1097  }
1098  Marker2 const& marker2() const
1099  {
1100  return super::marker2();
1101  }
1102  Marker3 const& marker3() const
1103  {
1104  return super::marker3();
1105  }
1106 
1107  bool isGhostFace() const
1108  {
1109  return super2::isGhostFace( super::processId() );
1110  }
1111 
1112 
1117  {
1118  return super2::isInterProcessDomain( super::processId() );
1119  }
1120 
1124  bool isOnBoundary() const
1125  {
1126  return super::isOnBoundary();
1127  }
1128 
1132  uint16_type boundaryEntityDimension() const
1133  {
1135  }
1136 
1140  bool isGhostCell() const
1141  {
1142  return super::isGhostCell();
1143  }
1144 
1148  uint16_type processId() const
1149  {
1150  return super::processId();
1151  }
1152 
1153 
1154 
1158  edge_type const& edge( uint16_type i ) const
1159  {
1160  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1161  DCHECK( M_edges[i] != nullptr ) << "invalid edge (null pointer) for edge local id " << i << " in element " << this->id();
1162  return boost::cref( *M_edges[i] );
1163  }
1164 
1168  edge_type& edge( uint16_type i )
1169  {
1170  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1171  DCHECK( M_edges[i] != nullptr ) << "invalid edge (null pointer) for edge local id " << i << " in element " << this->id();
1172  return boost::ref( *M_edges[i] );
1173  }
1174 
1175  edge_type & face( uint16_type i )
1176  {
1177  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1178  DCHECK( M_edges[i] != nullptr ) << "invalid edge (null pointer) for edge local id " << i << " in element " << this->id();
1179 
1180  return boost::ref( *M_edges[i] );
1181  }
1182 
1186  edge_type const& face( uint16_type i ) const
1187  {
1188  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1189  DCHECK( M_edges[i] != nullptr ) << "invalid edge (null pointer) for edge local id " << i << " in element " << this->id();
1190 
1191  return boost::cref( *M_edges[i] );
1192  }
1193 
1194  edge_type const* facePtr( uint16_type i ) const
1195  {
1196  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1197  DCHECK( M_edges[i] != nullptr ) << "invalid edge (null pointer) for edge local id " << i << " in element " << this->id();
1198 
1199  return M_edges[i];
1200  }
1201 
1206  void setFace( uint16_type const i, edge_type const & p )
1207  {
1208  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1209 
1210  M_edges[i] = const_cast<edge_type*>( boost::addressof( p ) );
1211  }
1212 
1216  edge_permutation_type edgePermutation( uint16_type i ) const
1217  {
1218  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1219 
1220  return M_edge_permutation[i];
1221  }
1225  edge_permutation_type facePermutation( uint16_type i ) const
1226  {
1227  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1228 
1229  return M_edge_permutation[i];
1230  }
1231 
1235  edge_permutation_type permutation( uint16_type i ) const
1236  {
1237  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges << " in element id " << this->id();
1238  return M_edge_permutation[i];
1239  }
1240 
1245  void setEdge( uint16_type i, edge_type const & p )
1246  {
1247  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1248  M_edges[i] = const_cast<edge_type*>( boost::addressof( p ) );
1249  }
1250 
1251  void setEdgePermutation( uint16_type i, edge_permutation_type o )
1252  {
1253  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1254 
1255  M_edge_permutation[i] = o;
1256  }
1257 
1258  typedef typename std::vector<edge_type*>::iterator face_iterator;
1259  typedef typename std::vector<edge_type*>::const_iterator face_const_iterator;
1260 
1264  std::pair<face_iterator,face_iterator>
1266  {
1267  return std::make_pair( M_edges.begin(), M_edges.end() );
1268  }
1269 
1273  std::pair<face_const_iterator,face_const_iterator>
1274  faces() const
1275  {
1276  return std::make_pair( M_edges.begin(), M_edges.end() );
1277  }
1278 
1279  void disconnectSubEntities()
1280  {
1281  for(unsigned int i = 0; i<numLocalEdges;++i)
1282  {
1283  M_edges[i]->disconnect(*this);
1284  }
1285  }
1286 
1287 private:
1288 
1289  friend class boost::serialization::access;
1290  template<class Archive>
1291  void serialize( Archive & ar, const unsigned int version )
1292  {
1293  DVLOG(2) << "Serializing Geoelement2D id: " << this->id() << "...\n";
1294  ar & boost::serialization::base_object<super>( *this );
1295  ar & boost::serialization::base_object<super2>( *this );
1296  ar & M_edges;
1297  }
1298 
1299 
1300 private:
1301 
1302  std::vector<edge_type*> M_edges;
1303  std::vector<edge_permutation_type> M_edge_permutation;
1304 };
1305 
1306 /*-------------------------------------------------------------------------
1307  GeoElement2D
1308  --------------------------------------------------------------------------*/
1309 template <uint16_type Dim, typename GEOSHAPE, typename SFO, typename T>
1310 const uint16_type GeoElement2D<Dim, GEOSHAPE, SFO, T>::numLocalEdges;
1311 template <uint16_type Dim, typename GEOSHAPE, typename SFO, typename T>
1312 const uint16_type GeoElement2D<Dim, GEOSHAPE, SFO, T>::nDim;
1313 
1314 
1320 template<uint16_type Dim,
1321  typename GEOSHAPE,
1322  typename T = double>
1324  :
1325 public GeoND<Dim, GEOSHAPE, T, GeoElement0D<Dim, SubFaceOfNone, T> >,
1326 public SubFaceOfNone
1327 {
1328 public:
1329 
1330  static const uint16_type nDim = Dim;
1331 
1333  typedef SubFaceOfNone super2;
1334 
1335  typedef GEOSHAPE GeoShape;
1336 
1337  typedef typename super::face_type entity_face_type;
1338 
1340  typedef self_type element_type;
1344 
1345  typedef typename super::node_type node_type;
1346 
1347  typedef typename super::vertex_permutation_type vertex_permutation_type;
1348  typedef typename super::edge_permutation_type edge_permutation_type;
1349  typedef typename super::face_permutation_type face_permutation_type;
1350  typedef typename super::face_permutation_type permutation_type;
1351 
1353  static const uint16_type numLocalVertices = super::numVertices;
1355  static const uint16_type numLocalFaces = super::numFaces;
1357  static const uint16_type numLocalEdges = super::numEdges;
1358 
1364  explicit GeoElement3D( size_type id = 0 )
1365  :
1366  super( id ),
1367  super2(),
1368  M_edges( numLocalEdges, nullptr ),
1369  M_faces( numLocalFaces ),
1370  M_edge_permutation( numLocalEdges, edge_permutation_type( edge_permutation_type::IDENTITY ) ),
1371  M_face_permutation( numLocalFaces )
1372  {
1373  std::fill( M_faces.begin(), M_faces.end(), ( face_type* )0 );
1374 
1375  std::fill( M_face_permutation.begin(), M_face_permutation.end(), face_permutation_type( face_permutation_type::IDENTITY ) );
1376  }
1377 
1382  :
1383  super( g ),
1384  super2( g ),
1385  M_edges( g.M_edges ),
1386  M_faces( g.M_faces ),
1387  M_edge_permutation( g.M_edge_permutation ),
1388  M_face_permutation( g.M_face_permutation )
1389  {}
1390 
1395  {}
1396 
1401  {
1402  if ( this != &g )
1403  {
1404  super::operator=( g );
1405  M_edges = g.M_edges;
1406  M_faces = g.M_faces;
1407  M_edge_permutation = g.M_edge_permutation;
1408  M_face_permutation = g.M_face_permutation;
1409  }
1410 
1411  return *this;
1412  }
1413 
1414  //void setMesh( MeshBase const* m ) { super::setMesh( m ); }
1415  MeshBase const* mesh() const
1416  {
1417  return super::mesh();
1418  }
1422  size_type id() const
1423  {
1424  return super::id();
1425  }
1426 
1427  Marker1 const& marker() const
1428  {
1429  return super::marker();
1430  }
1431  Marker1& marker()
1432  {
1433  return super::marker();
1434  }
1435  Marker2 const& marker2() const
1436  {
1437  return super::marker2();
1438  }
1439  Marker3 const& marker3() const
1440  {
1441  return super::marker3();
1442  }
1443 
1444  bool isGhostFace() const
1445  {
1446  return super2::isGhostFace( super::processId() );
1447  }
1448 
1449 
1454  {
1455  return super2::isInterProcessDomain( super::processId() );
1456  }
1457 
1461  bool isOnBoundary() const
1462  {
1463  return super::isOnBoundary();
1464  }
1465 
1469  uint16_type boundaryEntityDimension() const
1470  {
1472  }
1473 
1477  bool isGhostCell() const
1478  {
1479  return super::isGhostCell();
1480  }
1481 
1485  uint16_type processId() const
1486  {
1487  return super::processId();
1488  }
1489 
1490  size_type ad_first() const
1491  {
1492  return invalid_size_type_value;
1493  }
1494  uint16_type pos_first() const
1495  {
1497  }
1498  size_type ad_second() const
1499  {
1500  return invalid_size_type_value;
1501  }
1502  uint16_type pos_second() const
1503  {
1505  }
1506 
1507  edge_type const& edge( uint16_type i ) const
1508  {
1509  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1510  DCHECK( M_edges[i] != nullptr ) << "invalid edge (null pointer) for edge local id " << i << " in element " << this->id();
1511 
1512  return *M_edges[i];
1513  }
1514 
1515  edge_type& edge( uint16_type i )
1516  {
1517  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1518  DCHECK( M_edges[i] != nullptr ) << "invalid edge (null pointer) for edge local id " << i << " in element " << this->id();
1519 
1520  return *M_edges[i];
1521  }
1522 
1523  edge_type const* edgePtr( uint16_type i ) const
1524  {
1525  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1526  DCHECK( M_edges[i] != nullptr ) << "invalid edge (null pointer) for edge local id " << i << " in element " << this->id();
1527 
1528  return M_edges[i];
1529  }
1530 
1531  edge_permutation_type edgePermutation( uint16_type i ) const
1532  {
1533  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1534  DCHECK( M_edges[i] != nullptr ) << "invalid edge (null pointer) for edge local id " << i << " in element " << this->id();
1535 
1536 
1537  return M_edge_permutation[i];
1538  }
1539 
1543  void setEdge( uint16_type const i, edge_type const & p )
1544  {
1545  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1546  DCHECK( boost::addressof( p ) ) << "invalid edge (null pointer) for edge local id " << i << " in element " << this->id();
1547  M_edges[i] = const_cast<edge_type*>( boost::addressof( p ) );
1548  }
1549 
1550  void setEdgePermutation( uint16_type i, edge_permutation_type o )
1551  {
1552  DCHECK( i < numLocalEdges ) << "invalid local edge index " << i << " should be less than " << numLocalEdges ;
1553 
1554  M_edge_permutation[i] = o;
1555  }
1556 
1557  face_type const& face( uint16_type i ) const
1558  {
1559  FEELPP_ASSERT( i < numLocalFaces )( this->id() )( i ).error( "invalid local edge index" );
1560  FEELPP_ASSERT( M_faces[i] )( this->id() )( i ).error( "invalid edge (null pointer)" );
1561  return *M_faces[i];
1562  }
1563 
1564  face_type& face( uint16_type i )
1565  {
1566  FEELPP_ASSERT( i < numLocalFaces )( this->id() )( i ).error( "invalid local edge index" );
1567  FEELPP_ASSERT( M_faces[i] )( this->id() )( i ).error( "invalid edge (null pointer)" );
1568  return *M_faces[i];
1569  }
1570 
1571  face_type const* facePtr( uint16_type i ) const
1572  {
1573  FEELPP_ASSERT( i < numLocalFaces )( this->id() )( i ).error( "invalid local edge index" );
1574  //FEELPP_ASSERT( M_faces[i] )( i ).error( "invalid edge (null pointer)" );
1575  return M_faces[i];
1576  }
1577 
1578  face_permutation_type facePermutation( uint16_type i ) const
1579  {
1580  FEELPP_ASSERT( i < numLocalFaces )( this->id() )( i ).error( "invalid local face index" );
1581  FEELPP_ASSERT( M_faces[i] )( this->id() )( i ).error( "invalid face (null pointer)" );
1582  return M_face_permutation[i];
1583  }
1584  face_permutation_type permutation( uint16_type i ) const
1585  {
1586  FEELPP_ASSERT( i < numLocalFaces )( this->id() )( i ).error( "invalid local face index" );
1587  FEELPP_ASSERT( M_faces[i] )( this->id() )( i ).error( "invalid face (null pointer)" );
1588  return M_face_permutation[i];
1589  }
1590 
1594  void setFace( uint16_type const i, face_type const & p )
1595  {
1596  M_faces[i] = const_cast<face_type*>( boost::addressof( p ) );
1597  }
1598 
1599  void setFacePermutation( uint16_type i, face_permutation_type o )
1600  {
1601  FEELPP_ASSERT( i < numLocalFaces )( this->id() )( i ).error( "invalid local face index" );
1602  M_face_permutation[i] = o;
1603  }
1604 
1605  typedef typename ublas::bounded_array<face_type*, numLocalFaces>::iterator face_iterator;
1606  typedef typename ublas::bounded_array<face_type*, numLocalFaces>::const_iterator face_const_iterator;
1607 
1611  std::pair<face_iterator,face_iterator>
1613  {
1614  return std::make_pair( M_faces.begin(), M_faces.end() );
1615  }
1616 
1620  std::pair<face_const_iterator,face_const_iterator>
1621  faces() const
1622  {
1623  return std::make_pair( M_faces.begin(), M_faces.end() );
1624  }
1625 private:
1626 
1627  friend class boost::serialization::access;
1628  template<class Archive>
1629  void serialize( Archive & ar, const unsigned int version )
1630  {
1631  ar & boost::serialization::base_object<super>( *this );
1632  }
1633 
1634 
1635 private:
1636 
1637  ublas::bounded_array<edge_type*, numLocalEdges> M_edges;
1638  ublas::bounded_array<face_type*, numLocalFaces> M_faces;
1639 
1640  ublas::bounded_array<edge_permutation_type, numLocalEdges> M_edge_permutation;
1641  ublas::bounded_array<face_permutation_type, numLocalFaces> M_face_permutation;
1642 };
1648 /*-------------------------------------------------------------------------
1649  GeoElement3D
1650  --------------------------------------------------------------------------*/
1651 template <uint16_type Dim, typename GEOSHAPE, typename T>
1652 const uint16_type GeoElement3D<Dim, GEOSHAPE, T>::numLocalVertices;
1653 template <uint16_type Dim, typename GEOSHAPE, typename T>
1654 const uint16_type GeoElement3D<Dim, GEOSHAPE, T>::numLocalFaces;
1655 template <uint16_type Dim, typename GEOSHAPE, typename T>
1656 const uint16_type GeoElement3D<Dim, GEOSHAPE, T>::numLocalEdges;
1657 template <uint16_type Dim, typename GEOSHAPE, typename T>
1658 const uint16_type GeoElement3D<Dim, GEOSHAPE, T>::nDim;
1659 
1660 } // Feel
1661 #endif

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