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
geond.hpp
Go to the documentation of this file.
1 /*
2  This file is part of the Feel library
3 
4  Copyright (C) 2001,2002,2003,2004 EPFL, INRIA and Politechnico di Milano
5  Copyright (C) 2008 Université Joseph Fourier (Grenoble I)
6 
7  This library is free software; you can redistribute it and/or
8  modify it under the terms of the GNU Lesser General Public
9  License as published by the Free Software Foundation; either
10  version 3.0 of the License, or (at your option) any later version.
11 
12  This library is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15  Lesser General Public License for more details.
16 
17  You should have received a copy of the GNU Lesser General Public
18  License along with this library; if not, write to the Free Software
19  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20 */
24 #ifndef _GEOND_HH_
25 #define _GEOND_HH_
26 
27 #include <boost/numeric/ublas/storage.hpp>
28 
30 #include <feel/feelmesh/geo0d.hpp>
31 #include <feel/feelpoly/geomap.hpp>
32 #include <feel/feelmesh/marker.hpp>
34 
35 namespace Feel
36 {
37 class MeshBase;
38 
39 template<int Dim, int Order, int RealDim, template<uint16_type,uint16_type,uint16_type> class Entity, typename T> struct GT_Lagrange;
40 
42 namespace detail
43 {
48 template <typename GeoShape>
49 class ReversePoint
50 {
51 public:
52  uint16_type operate( uint16_type const & point )
53  {
54  return ( point < GeoShape::numVertices ) ?
55  GeoShape::numVertices - point :
56  GeoShape::numPoints - point + GeoShape::numVertices;
57  }
58 
59 };
60 }
62 
63 
69 template <uint16_type Dim,
70  typename GEOSHAPE,
71  typename T = double,
72  typename POINTTYPE = Geo0D<Dim, T> >
73 class GeoND
74  :
75 public GeoEntity<GEOSHAPE>
76 {
77  typedef GeoEntity<GEOSHAPE> super;
78 public:
79 
80  typedef T value_type;
82  typedef self_type element_type;
83 
84  typedef GEOSHAPE GeoShape;
85  typedef POINTTYPE PointType;
86 
87  typedef PointType point_type;
88  typedef typename super::face_type face_type;
89 
90  static const size_type Shape = super::Shape;
91  static const uint16_type numPoints = super::numPoints;
92  static const uint16_type numVertices = super::numVertices;
93  static const uint16_type numLocalPoints = super::numPoints;
94  static const uint16_type numLocalEdges = super::numEdges;
95  static const uint16_type numLocalVertices = super::numVertices;
96  static const int numFaces = super::numFaces;
97  static const int numEdges = super::numEdges;
98  static const int numTopologicalFaces = super::numTopologicalFaces;
99  static const uint16_type numNeighbors = super::numTopologicalFaces;
100 
101 
102  typedef typename ublas::bounded_array<point_type*, numPoints>::iterator point_iterator;
103  typedef typename ublas::bounded_array<point_type*, numPoints>::const_iterator point_const_iterator;
104 
105  typedef typename matrix_node<value_type>::type matrix_node_type;
106  typedef typename node<value_type>::type node_type;
107 
108  static const uint16_type nDim = super::nDim;
109  static const uint16_type nOrder = super::nOrder;
110  static const uint16_type nRealDim = super::nRealDim;
111 
112  template<int GmOrder>
113  struct GetGm
114  {
115 
116  typedef typename mpl::if_<mpl::bool_<GeoShape::is_hypercube>,
117  mpl::identity<GT_Lagrange<nDim, GmOrder, nRealDim, Hypercube, T> >,
118  mpl::identity<GT_Lagrange<nDim, GmOrder, nRealDim, Simplex, T> > >::type::type type;
119  typedef boost::shared_ptr<type> ptrtype;
120  };
121  typedef typename GetGm<nOrder>::type gm_type;
122  typedef typename GetGm<nOrder>::ptrtype gm_ptrtype;
123 
124  typedef typename GetGm<1>::type gm1_type;
125  typedef typename GetGm<1>::ptrtype gm1_ptrtype;
126 
127  typedef typename super::vertex_permutation_type vertex_permutation_type;
128  typedef typename super::edge_permutation_type edge_permutation_type;
129  typedef typename super::face_permutation_type face_permutation_type;
130  typedef typename mpl::if_<mpl::equal_to<mpl::int_<nDim>,
131  mpl::int_<1> >,
132  mpl::identity<vertex_permutation_type>,
133  typename mpl::if_<mpl::equal_to<mpl::int_<nDim>,
134  mpl::int_<2> >,
135  mpl::identity<edge_permutation_type>,
136  mpl::identity<face_permutation_type> >::type>::type::type permutation_type;
141  :
142  super( 0 ),
143  M_points( numPoints ),
144  //M_face_points( numTopologicalFaces ),
145  M_G( nRealDim, numPoints ),
146  M_barycenter( nRealDim ),
147  M_barycenterfaces( nRealDim, numTopologicalFaces ),
148  M_h( 1 ),
149  M_h_face( numTopologicalFaces, 1 ),
150  M_h_edge( numLocalEdges, 1 ),
151  M_measure( 1 ),
152  M_measurefaces( numTopologicalFaces ),
153  M_normals( nRealDim, numTopologicalFaces ),
154  M_has_points( false ),
155  M_neighbors( numNeighbors, std::make_pair( invalid_size_type_value, 0 ) ),
156  M_marker1(),
157  M_marker2(),
158  M_marker3(),
159  M_gm(),
160  M_gm1()
161  {
162  }
163 
170  explicit GeoND( size_type id )
171  :
172  super( id ),
173  M_points( numPoints ),
174  //M_face_points( numTopologicalFaces ),
175  M_G( nRealDim, numPoints ),
176  M_barycenter( nRealDim ),
177  M_barycenterfaces( nRealDim, numTopologicalFaces ),
178  M_h( 1 ),
179  M_h_face( numTopologicalFaces, 1 ),
180  M_h_edge( numLocalEdges, 1 ),
181  M_measure( 1 ),
182  M_measurefaces( numTopologicalFaces ),
183  M_normals( nRealDim, numTopologicalFaces ),
184  M_has_points( false ),
185  M_neighbors( numNeighbors, std::make_pair( invalid_size_type_value, 0 ) ),
186  M_marker1(),
187  M_marker2(),
188  M_marker3(),
189  M_gm(),
190  M_gm1()
191  {
192  }
193 
194  GeoND( GeoND const& e )
195  :
196  super( e ),
197  M_points( numPoints ),
198  //M_face_points( e.M_face_points ),
199  M_G( nRealDim, numPoints ),
200  M_barycenter( e.M_barycenter ),
201  M_barycenterfaces( e.M_barycenterfaces ),
202  M_h( e.M_h ),
203  M_h_face( e.M_h_face ),
204  M_h_edge( e.M_h_edge ),
205  M_measure( e.M_measure ),
206  M_measurefaces( numTopologicalFaces ),
207  M_normals( e.M_normals ),
208  M_has_points( false ),
209  M_neighbors( numNeighbors, std::make_pair( invalid_size_type_value, 0 ) ),
210  M_marker1( e.M_marker1 ),
211  M_marker2( e.M_marker2 ),
212  M_marker3( e.M_marker3 ),
213  M_gm(),
214  M_gm1()
215  {
216  M_G = e.M_G;
217 
218  for ( uint16_type i = 0; i < numLocalPoints; ++i )
219  M_points[ i ] = e.M_points[ i ];
220  }
221 
226  {
227 
228  }
229 
230 #if 0
231 
234  void setMeshAndGm( MeshBase const* m, gm_ptrtype const& gm ) const
235  {
236  M_mesh = m;
237  M_gm = gm;
238  }
239 #endif
240 
243  void setMeshAndGm( MeshBase const* m, gm_ptrtype const& gm, gm1_ptrtype const& gm1 ) const
244  {
245  M_mesh = m;
246  M_gm = gm;
247  M_gm1 = gm1;
248  }
249 
250  void setMesh( MeshBase const* m ) const
251  {
252  M_mesh = m;
253  }
254 
256  gm_ptrtype gm() const
257  {
258  return M_gm;
259  }
260 
262  gm1_ptrtype gm1() const
263  {
264  return M_gm1;
265  }
266 
270  MeshBase const* mesh() const
271  {
272  return M_mesh;
273  }
274 
279  bool hasPoints() const
280  {
281  return M_has_points;
282  }
283 
284 #if 0
285 
292  GeoND & operator=( GeoND const & G )
293  {
294  if ( this != &G )
295  {
296  super::operator=( G );
297 
298  for ( uint16_type i = 0; i < numLocalPoints; ++i )
299  M_points[ i ] = G.M_points[ i ];
300 
301  //M_face_points = G.M_face_points;
302  M_G = G.M_G;
303 
304  M_barycenter = G.M_barycenter;
305  M_barycenterfaces = G.M_barycenterfaces;
306  M_h = G.M_h;
307  M_h_face = G.M_h_face;
308  M_h_edge = G.M_h_edge;
309 
310  M_has_points = G.M_has_points;
311 
312  M_neighbors = G.M_neighbors;
313 
314  M_marker1 = G.M_marker1;
315  M_marker2 = G.M_marker2;
316  M_marker3 = G.M_marker3;
317 
318  M_gm = G.M_gm;
319  M_gm1 = G.M_gm1;
320  }
321 
322  return *this;
323  }
324 #endif
325 
328  uint16_type nPoints() const
329  {
330  return numPoints;
331  }
332 
340  uint16_type nNeighbors() const
341  {
342  return numNeighbors;
343  }
344 
351  std::pair<size_type,size_type> const& neighbor( uint16_type n ) const
352  {
353  return M_neighbors[n];
354  }
355 
359  void setNeighbor( uint16_type n, size_type neigh_id, size_type proc_id )
360  {
361  M_neighbors[n] = std::make_pair( neigh_id, proc_id );
362  }
363 
364  bool isNeighbor( self_type const& G ) const
365  {
366  for ( uint16_type i = 0; i< this->nNeighbors() ; ++i )
367  if ( this->neighbor( i ).first==G.id() ) return true;
368 
369  return false;
370  }
371 
375  node_type barycenter() const
376  {
377  return M_barycenter;
378  }
379 
383  node_type faceBarycenter( uint16_type f ) const
384  {
385  return ublas::column( M_barycenterfaces, f );
386  }
387 
391  matrix_node_type faceBarycenters() const
392  {
393  return M_barycenterfaces;
394  }
395 
399  permutation_type permutation( uint16_type /*f*/ ) const
400  {
401  return permutation_type();
402  }
403 
408  PointType & point( uint16_type i )
409  {
410  return *( static_cast<POINTTYPE *>( M_points[ i ] ) );
411  }
412 
413 
418  PointType const & point ( uint16_type i ) const
419  {
420  return *( static_cast<POINTTYPE *>( M_points[ i ] ) );
421  }
426  PointType* pointPtr( uint16_type i )
427  {
428  return M_points[ i ];
429  }
430 
435  PointType const* pointPtr ( uint16_type i ) const
436  {
437  return M_points[ i ];
438  }
439 
442  PointType const & facePoint ( uint16_type __f, uint16_type const __i ) const
443  {
444  return M_face_points[__f][__i];
445  }
446 
447 
448 
457  PointType & reversepoint( uint16_type const i )
458  {
459  return *( static_cast<POINTTYPE *>( M_points[ detail::ReversePoint<GEOSHAPE>::operate( i ) ] ) );
460  }
461 
462 
471  PointType const & reversepoint ( uint16_type const i ) const
472  {
473  return *( static_cast<POINTTYPE *>( M_points[ detail::ReversePoint<GEOSHAPE>::operate( i ) ] ) );
474  }
475 
476 
481  void setPoint( uint16_type const i, point_type const & p );
482 
491  std::ostream & showMe( bool verbose = false, std::ostream & c = std::cout ) const;
492 
502  void swapPoints( const uint16_type & pt1, const uint16_type & pt2 );
503 
513  void exchangePoints( const uint16_type otn[ numPoints ] );
514 
523  matrix_node_type const& G() const
524  {
525  return M_G;
526  }
527 
536  //matrix_node_type vertices() const { return ublas::subrange( M_G, 0, nRealDim, 0, numVertices ); }
537  matrix_node_type vertices() const
538  {
539  return ublas::subrange( M_G, 0, nRealDim, 0, numVertices );
540  }
541 
550  matrix_node_type & G()
551  {
552  return M_G;
553  }
554 
555  point_iterator beginPoint()
556  {
557  return M_points.begin();
558  }
559  point_const_iterator beginPoint() const
560  {
561  return M_points.begin();
562  }
563  point_iterator endPoint()
564  {
565  return M_points.end();
566  }
567  point_const_iterator endPoint() const
568  {
569  return M_points.end();
570  }
571 
578  double h() const
579  {
580  return M_h;
581  }
582 
590  double hFace( uint16_type f ) const
591  {
592  return M_h_face[f];
593  }
594 
595  double hEdge( uint16_type f ) const
596  {
597  return M_h_edge[f];
598  }
599 
600  struct tt
601  {
602  static uint16_type fToP( uint16_type const _localFace, uint16_type const _point )
603  {
604  return super::eToP( _localFace, _point );
605  }
606  };
607 
611  double measure() const
612  {
613  return M_measure;
614  }
615 
619  double faceMeasure( uint16_type f ) const
620  {
621  return M_measurefaces[f];
622  }
623 
627  std::vector<double> const& faceMeasures() const
628  {
629  return M_measurefaces;
630  }
631 
635  matrix_node_type const& normals() const
636  {
637  return M_normals;
638  }
639 
643  ublas::matrix_column<matrix_node_type const> normal( uint16_type f ) const
644  {
645  return ublas::column( M_normals, f );
646  }
647 
656  static uint16_type fToP( uint16_type const _localFace, uint16_type const _point )
657  {
658 #if 1
659  typedef typename mpl::if_<mpl::not_equal_to<mpl::int_<super::nDim>, mpl::int_<2> >,
660  mpl::identity<super>,
661  mpl::identity<tt> >::type the_type;
662  return the_type::type::fToP( _localFace, _point );
663 #else
664  return super::fToP( _localFace, _point );
665 #endif
666 
667  }
668 
669 
675  uint16_type nOppositePointsPerFace() const
676  {
677  return super::nbOppositePointsPerFace;
678  }
679 
685  uint16_type faceToOppositePoint( uint16_type const _localFace, uint16_type const _point ) const
686  {
687  return super::faceToOppositePoint( _localFace, _point );
688  }
689 
695  {
696 
697  // Calculate vectors originating from vertex zero
698  ublas::matrix<T> orientation_matrix ( nRealDim,nRealDim );
699 
700  for ( int i = 0; i < nRealDim ; ++i )
701  {
702  ublas::row( orientation_matrix, i ) = ( ublas::column( this->G(), i+1 ) -
703  ublas::column( this->G(), 0 ) );
704 
705  }
706 
707  LU< ublas::matrix<T> > lu( orientation_matrix );
708  T sgn=lu.det();
709 
710  return ( sgn > 0 ) ? 1 : 0;
711  }
712 
713  void applyDisplacement( int i, ublas::vector<double> const& u )
714  {
715  ublas::column( M_G, i ) += u;
716  ( *M_points[ i ] ) += u;
717  }
718  void applyDisplacementG( int i, ublas::vector<double> const& u )
719  {
720  ublas::column( M_G, i ) += u;
721  }
728  void setTags( std::vector<int> const& tags )
729  {
730  M_marker1.assign( tags[0] );
731 
732  if ( tags.size() > 1 )
733  M_marker2.assign( tags[1] );
734 
735  if ( tags.size() > 2 )
736  {
737  this->setNumberOfPartitions( tags[2] );
738  this->setProcessId( tags[3] );
739 
740  if ( tags[2] > 1 )
741  {
742  std::vector<int> p( tags[2]-1 );
743 
744  for ( size_type i = 0; i < p.size(); ++i )
745  {
746  p[i] = tags[4+i];
747  }
748 
749  this->setNeighborPartitionIds( p );
750  }
751 
752  }
753  }
754  Marker1 const& marker() const
755  {
756  return M_marker1;
757  }
758  Marker1& marker()
759  {
760  return M_marker1;
761  }
762  void setMarker( flag_type v )
763  {
764  return M_marker1.assign( v );
765  }
766 
767  Marker2 const& marker2() const
768  {
769  return M_marker2;
770  }
771  Marker2& marker2()
772  {
773  return M_marker2;
774  }
775  void setMarker2( flag_type v )
776  {
777  return M_marker2.assign( v );
778  }
779 
780  Marker3 const& marker3() const
781  {
782  return M_marker3;
783  }
784  Marker3& marker3()
785  {
786  return M_marker3;
787  }
788  void setMarker3( flag_type v )
789  {
790  return M_marker3.assign( v );
791  }
792 
795  {
796  return M_pneighbors.size();
797  }
799  std::set<size_type> const& pointElementNeighborIds() const
800  {
801  return M_pneighbors;
802  }
804  void setMeasurePointElementNeighbors( value_type meas )
805  {
806  M_meas_pneighbors = meas;
807  }
809  value_type measurePointElementNeighbors() const
810  {
811  return M_meas_pneighbors;
812  }
813 
814  void update();
815  void updateWithPc( typename gm_type::precompute_ptrtype const& pc, typename gm_type::faces_precompute_type & pcf );
816 private:
817 
818  void updatep( typename gm_type::faces_precompute_type & pcf, mpl::bool_<true> );
819  void updatep( typename gm_type::faces_precompute_type & pcf, mpl::bool_<false> );
820 
821 private:
822 
823 private:
824 
825  friend class boost::serialization::access;
826  template<class Archive>
827  void serialize( Archive & ar, const unsigned int version )
828  {
829  DVLOG(2) << "Serializing GeoND...\n";
830  DVLOG(2) << " - base class...\n";
831  ar & boost::serialization::base_object<super>( *this );
832  DVLOG(2) << " - points...\n";
833  ar & M_points;
834  DVLOG(2) << " - G...\n";
835  ar & M_G;
836  DVLOG(2) << " - marker1...\n";
837  ar & M_marker1;
838  DVLOG(2) << " - marker1: " << M_marker1.value() << "...\n";
839  DVLOG(2) << " - marker2...\n";
840  ar & M_marker2;
841  DVLOG(2) << " - marker2: " << M_marker2.value() << "...\n";
842  DVLOG(2) << " - marker3...\n";
843  ar & M_marker3;
844  DVLOG(2) << " - marker3: " << M_marker3.value() << "...\n";
845  }
846 
847 private:
849  std::vector<point_type*> M_points;
850 
852  std::vector<std::vector<point_type*> > M_face_points;
853 
855  matrix_node_type M_G;
856  node_type M_barycenter;
857  matrix_node_type M_barycenterfaces;
858 
859  double M_h;
860  std::vector<double> M_h_face;
861  std::vector<double> M_h_edge;
862 
863  double M_measure;
864  std::vector<double> M_measurefaces;
865  matrix_node_type M_normals;
866 
867 
868  bool M_has_points;
869 
873  std::vector<std::pair<size_type,size_type> > M_neighbors;
875  std::set<size_type> M_pneighbors;
877  value_type M_meas_pneighbors;
878 
879  Marker1 M_marker1;
880  Marker2 M_marker2;
881  Marker3 M_marker3;
882 
883  // mesh to which the geond element belongs to
884  mutable MeshBase const* M_mesh;
885  mutable gm_ptrtype M_gm;
886  mutable gm1_ptrtype M_gm1;
887 };
888 
889 template <uint16_type Dim, typename GEOSHAPE, typename T, typename POINTTYPE>
890 const uint16_type GeoND<Dim,GEOSHAPE, T, POINTTYPE>::numLocalPoints;
891 
892 template <uint16_type Dim, typename GEOSHAPE, typename T, typename POINTTYPE>
893 const uint16_type GeoND<Dim,GEOSHAPE, T, POINTTYPE>::numLocalVertices;
894 
895 template <uint16_type Dim, typename GEOSHAPE, typename T, typename POINTTYPE>
896 inline
897 void
898 GeoND<Dim,GEOSHAPE, T, POINTTYPE>::setPoint( uint16_type const i, point_type const & p )
899 {
900  M_points[ i ] = const_cast<point_type *>( &p );
901  //VLOG(1) << "[setPoint] üpdate point index " << i << " with "<< M_points[i]->id() << "\n";
902  FEELPP_ASSERT( const_cast<point_type *>( &p ) != 0 ).error( "invalid Geo0D<>" );
903  ublas::column( M_G, i ) = M_points[i]->node();
904  M_has_points = true;
905 }
906 
907 
908 template <uint16_type Dim, typename GEOSHAPE, typename T, typename POINTTYPE>
909 std::ostream &
910 GeoND<Dim,GEOSHAPE, T, POINTTYPE>::showMe( bool verbose, std::ostream & out ) const
911 {
912  out << "----- BEGIN OF GeoND data ---" << std::endl << std::endl;
913  out << " GeoND object of shape " << Shape << std::endl;
914  out << " Number of Vertices = " << numVertices << std::endl;
915  out << " Number of Points = " << numPoints << std::endl;
916  out << " id = " << this->id() << std::endl;
917  out << " G = " << M_G << "\n";
918 
919  for ( int i = 0; i < numVertices; i++ )
920  {
921  out << "POINT id = " << i << std::endl;
922  out << point( i ).showMe( verbose, out );
923  }
924 
925  out << "----- END OF GeoND data ---" << std::endl << std::endl;
926  return out;
927 }
928 
929 template <uint16_type Dim, typename GEOSHAPE, typename T, typename POINTTYPE>
930 void GeoND<Dim,GEOSHAPE, T, POINTTYPE>::swapPoints( const uint16_type & pt1, const uint16_type & pt2 )
931 {
932  point_type * tmp( M_points[ pt1 ] );
933  M_points[ pt1 ] = M_points[ pt2 ];
934  M_points[ pt2 ] = tmp;
935 
936  // swap also the entries in G
937  ublas::column( M_G, pt1 ).swap( ublas::column( M_G, pt2 ) );
938 }
939 
940 
941 template <uint16_type Dim, typename GEOSHAPE, typename T, typename POINTTYPE>
942 void GeoND<Dim,GEOSHAPE, T, POINTTYPE>::exchangePoints( const uint16_type otn[ numPoints ] )
943 {
944  point_type * tmp[ numPoints ];
945 
946  for ( unsigned int i = 0; i < numPoints; ++i )
947  {
948  tmp[ i ] = M_points[ i ];
949  }
950 
951  for ( unsigned int i = 0; i < numPoints; ++i )
952  {
953  M_points[ i ] = tmp[ otn[ i ] ];
954  ublas::column( M_G, i ) = M_points[i]->node();
955  }
956 }
957 
958 template <uint16_type Dim, typename GEOSHAPE, typename T, typename POINTTYPE>
959 void
961 {
962  if ( !M_gm.use_count() )
963  M_gm = gm_ptrtype( new gm_type );
964 
965  if ( !M_gm1.use_count() )
966  M_gm1 = gm1_ptrtype( new gm1_type );
967 
968  auto pc = M_gm->preCompute( M_gm, M_gm->referenceConvex().vertices() );
969  auto pcf = M_gm->preComputeOnFaces( M_gm, M_gm->referenceConvex().barycenterFaces() );
970 
971  updateWithPc( pc, pcf );
972 }
973 
974 template <uint16_type Dim, typename GEOSHAPE, typename T, typename POINTTYPE>
975 void
976 GeoND<Dim,GEOSHAPE, T, POINTTYPE>::updateWithPc( typename gm_type::precompute_ptrtype const& pc,
977  typename gm_type::faces_precompute_type& pcf )
978 {
979  M_h = 0;
980 
981  for ( uint16_type __e = 0; __e < numLocalEdges; ++__e )
982  {
983  node_type const& __x1 = this->point( this->eToP( __e, 0 ) ).node();
984  node_type const& __x2 = this->point( this->eToP( __e, 1 ) ).node();
985  M_h_edge[__e] = ublas::norm_2( __x1-__x2 );
986  M_h = ( M_h > M_h_edge[__e] )?M_h:M_h_edge[__e];
987  }
988 
989  auto M = glas::average( M_G );
990  M_barycenter = ublas::column( M, 0 );
991  M_pneighbors.clear();
992 
993  for ( uint16_type __p = 0; __p < numPoints; ++__p )
994  {
995  std::copy( M_points[__p]->elements().begin(),
996  M_points[__p]->elements().end(),
997  std::inserter( M_pneighbors, M_pneighbors.begin() ) );
998  }
999 
1000  auto ctx = M_gm->template context<vm::JACOBIAN>( *this, pc );
1001  //M_gm->preCompute( M_gm, M_gm->referenceConvex().vertices() ) );
1002  double w = ( nDim == 3 )?4./3.:2;
1003  M_measure = w*ctx->J( 0 );
1004 
1005  updatep( pcf, typename mpl::equal_to<mpl::int_<nDim>, mpl::int_<nRealDim> >::type() );
1006 }
1007 template <uint16_type Dim, typename GEOSHAPE, typename T, typename POINTTYPE>
1008 void
1009 GeoND<Dim,GEOSHAPE, T, POINTTYPE>::updatep( typename gm_type::faces_precompute_type& pcf, mpl::bool_<true> )
1010 {
1011  if ( nRealDim==1 )
1012  {
1013  M_h_face[0] = 1;
1014  M_h_face[1] = 1;
1015  }
1016 
1017  else
1018  {
1019  int nEdges = GEOSHAPE::topological_face_type::numEdges;
1020 
1021  for ( uint16_type __f = 0; __f < numTopologicalFaces; ++__f )
1022  {
1023  M_h_face[__f] = 0;
1024 
1025  for ( uint16_type e = 0; e < nEdges; ++e )
1026  {
1027  double __l = 0;
1028 
1029  if ( Dim == 2 )
1030  {
1031  node_type const& __x1 = this->point( this->eToP( this->f2e( __f, __f ), 0 ) ).node();
1032  node_type const& __x2 = this->point( this->eToP( this->f2e( __f, __f ), 1 ) ).node();
1033  __l = ublas::norm_2( __x1-__x2 );
1034  }
1035 
1036  else
1037  {
1038  node_type const& __x1 = this->point( this->eToP( this->f2e( __f, e ), 0 ) ).node();
1039  node_type const& __x2 = this->point( this->eToP( this->f2e( __f, e ), 1 ) ).node();
1040  __l = ublas::norm_2( __x1-__x2 );
1041  }
1042 
1043  //std::cout << "face " << __f << " edge " << e << " edge " << this->f2e( __f, __f ) << " length " << __l << std::endl;
1044  M_h_face[__f] = ( M_h_face[__f] > __l )?M_h_face[__f]:__l;
1045  }
1046  }
1047  }
1048 
1049  //auto pc = M_gm->preComputeOnFaces( M_gm, M_gm->referenceConvex().barycenterFaces() );
1050  auto ctx = M_gm->template context<vm::POINT|vm::NORMAL|vm::KB|vm::JACOBIAN>(
1051  *this,
1052  pcf,
1053  0 );
1054 
1055  // jacobian of transformation
1056  std::vector<double> f2( numTopologicalFaces, 2 );
1057  std::vector<double> f3( numTopologicalFaces, 2 );
1058 
1059  if ( GEOSHAPE::is_simplex )
1060  {
1061  f2[0] = 2.82842712474619;
1062  f3[0] = 3.464101615137754;
1063  }
1064 
1065  for ( int f = 0; f < numTopologicalFaces; ++f )
1066  {
1067  ctx->update( *this, f );
1068  ublas::column( M_normals, f ) = ctx->unitNormal( 0 );
1069 #if 1 // doesn't work (vincent)
1070  ublas::column( M_barycenterfaces, f ) = ctx->xReal( 0 );
1071 #else
1072  ublas::column( M_barycenterfaces, f ) = ublas::column( glas::average( this->face( f ).G() ) );
1073 #endif
1074  double w = ( nDim == 3 )?f3[f]:( ( nDim==2 )?f2[f]:1 );
1075  M_measurefaces[f] = w*ctx->J( 0 )*ctx->normalNorm( 0 );
1076  }
1077 }
1078 
1079 template <uint16_type Dim, typename GEOSHAPE, typename T, typename POINTTYPE>
1080 void
1081 GeoND<Dim,GEOSHAPE, T, POINTTYPE>::updatep( typename gm_type::faces_precompute_type& pcf, mpl::bool_<false> )
1082 {
1083 }
1084 
1085 template <uint16_type Dim, typename GEOSHAPE, typename T, typename POINTTYPE>
1086 inline
1087 DebugStream&
1088 operator<<( DebugStream& __os, GeoND<Dim,GEOSHAPE, T, POINTTYPE> const& __n )
1089 {
1090  if ( __os.doPrint() )
1091  {
1092  std::ostringstream __str;
1093 
1094  __str << __n.showMe( true, __str );
1095 
1096  __os << __str.str() << "\n";
1097  }
1098 
1099  return __os;
1100 }
1101 
1102 template <uint16_type Dim, typename GEOSHAPE, typename T, typename POINTTYPE>
1103 inline
1104 NdebugStream&
1105 operator<<( NdebugStream& __os, GeoND<Dim,GEOSHAPE, T, POINTTYPE> const& __n )
1106 {
1107  return __os;
1108 }
1109 
1110 
1111 template <uint16_type Dim, typename GEOSHAPE, typename T, typename POINTTYPE>
1112 inline
1113 std::ostream&
1114 operator<<( std::ostream& __os, GeoND<Dim,GEOSHAPE, T, POINTTYPE> const& __n )
1115 {
1116  return __n.showMe( true, __os );
1117 }
1118 
1119 } // Feel
1120 #endif
void setMeasurePointElementNeighbors(value_type meas)
set the measure of point element neighbors
Definition: geond.hpp:804
double h() const
Definition: geond.hpp:578
node_type faceBarycenter(uint16_type f) const
Definition: geond.hpp:383
elements_type const & elements() const
Definition: elements.hpp:355
~GeoND()
Definition: geond.hpp:225
PointType & reversepoint(uint16_type const i)
Definition: geond.hpp:457
GeoND()
Definition: geond.hpp:140
ublas::matrix_column< matrix_node_type const > normal(uint16_type f) const
Definition: geond.hpp:643
void setNumberOfPartitions(uint16_type np)
Definition: geoentity.hpp:611
base mesh class
Definition: meshbase.hpp:67
PointType * pointPtr(uint16_type i)
Definition: geond.hpp:426
bool hasPoints() const
Definition: geond.hpp:279
value_type measurePointElementNeighbors() const
Definition: geond.hpp:809
node_type barycenter() const
Definition: geond.hpp:375
static uint16_type fToP(uint16_type const __localFace, uint16_type const __point)
Definition: geoentity.hpp:650
Definition: glas.hpp:157
uint16_type nNeighbors() const
Definition: geond.hpp:340
double hFace(uint16_type f) const
Definition: geond.hpp:590
std::vector< double > const & faceMeasures() const
Definition: geond.hpp:627
matrix_node_type const & normals() const
Definition: geond.hpp:635
MeshBase const * mesh() const
Definition: geond.hpp:270
uint16_type nPoints() const
Definition: geond.hpp:328
void setMeshAndGm(MeshBase const *m, gm_ptrtype const &gm, gm1_ptrtype const &gm1) const
Definition: geond.hpp:243
matrix_node_type faceBarycenters() const
Definition: geond.hpp:391
const size_type invalid_size_type_value
Definition: feelcore/feel.hpp:360
Definition: lu.hpp:219
matrix_node_type const & G() const
Definition: geond.hpp:523
Definition: glas.hpp:125
matrix_node_type vertices() const
Definition: geond.hpp:537
uint16_type nOppositePointsPerFace() const
Definition: geond.hpp:675
size_type numberOfPointElementNeighbors() const
Definition: geond.hpp:794
void setNeighborPartitionIds(std::vector< int > const &npids)
Definition: geoentity.hpp:628
value_type det()
Definition: lu.hpp:420
PointType const & reversepoint(uint16_type const i) const
Definition: geond.hpp:471
void setPoint(uint16_type const i, point_type const &p)
Definition: geond.hpp:898
Base class for Multi-dimensional basis Geometrical Entities.
Definition: geond.hpp:73
double measure() const
Definition: geond.hpp:611
std::pair< size_type, size_type > const & neighbor(uint16_type n) const
Definition: geond.hpp:351
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
GeoND(size_type id)
Definition: geond.hpp:170
PointType const * pointPtr(uint16_type i) const
Definition: geond.hpp:435
void setProcessId(uint16_type pid)
Definition: geoentity.hpp:460
std::set< size_type > const & pointElementNeighborIds() const
Definition: geond.hpp:799
bool isAnticlockwiseOriented() const
Definition: geond.hpp:694
static uint16_type fToP(uint16_type const _localFace, uint16_type const _point)
Definition: geond.hpp:656
PointType & point(uint16_type i)
Definition: geond.hpp:408
uint16_type faceToOppositePoint(uint16_type const _localFace, uint16_type const _point) const
Definition: geond.hpp:685
double faceMeasure(uint16_type f) const
Definition: geond.hpp:619
matrix_node_type & G()
Definition: geond.hpp:550
Definition: geoelement.hpp:483
void setTags(std::vector< int > const &tags)
Definition: geond.hpp:728
static uint16_type eToP(uint16_type const __localEdge, uint16_type const __point)
Definition: geoentity.hpp:642
gm1_ptrtype gm1() const
return the geometric mapping if a mesh was set
Definition: geond.hpp:262
base class for all geometric entities
Definition: geoentity.hpp:47
void setNeighbor(uint16_type n, size_type neigh_id, size_type proc_id)
Definition: geond.hpp:359
std::ostream & showMe(bool verbose=false, std::ostream &c=std::cout) const
Definition: geond.hpp:910
void exchangePoints(const uint16_type otn[numPoints])
Definition: geond.hpp:942
gm_ptrtype gm() const
return the geometric mapping if a mesh was set
Definition: geond.hpp:256
PointType const & point(uint16_type i) const
Definition: geond.hpp:418
permutation_type permutation(uint16_type) const
Definition: geond.hpp:399
void swapPoints(const uint16_type &pt1, const uint16_type &pt2)
Definition: geond.hpp:930

Generated on Sun Dec 22 2013 13:11:04 for Feel++ by doxygen 1.8.5