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
refsimplex.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: 2010-05-06
7 
8  Copyright (C) 2010 Université 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 2.1 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 __refsimplex_H
30 #define __refsimplex_H 1
31 
32 namespace Feel
33 {
34 template<uint16_type Dim, uint16_type Order, uint16_type RDim, typename T>
35 class Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T>
36  :
37 public Simplex<Dim, Order, RDim>
38 {
39 public:
40  typedef Simplex<Dim, Order, RDim> super;
41 
45 
46  static const uint16_type nDim = super::nDim;
47  static const uint16_type nOrder = super::nOrder;
48  static const uint16_type nRealDim = super::nRealDim;
49 
50  static const uint16_type topological_dimension = super::topological_dimension;
51  static const uint16_type real_dimension = super::real_dimension;
52 
53  static const size_type Shape = super::Shape;
54  static const size_type Geometry = super::Geometry;
55 
56  typedef T value_type;
57  typedef Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T> self_type;
58 
59  typedef typename mpl::if_<boost::is_same<typename super::element_type, boost::none_t>,
60  mpl::identity<boost::none_t>,
61  mpl::identity<Reference<typename super::element_type, nDim+1, nOrder, nRealDim, T> > >::type::type element_type;
62  typedef Reference<typename super::topological_face_type,
63  super::topological_face_type::nDim, nOrder, nRealDim, T> topological_face_type;
64 
65  typedef typename super::edge_to_point_t edge_to_point_t;
66  typedef typename super::face_to_point_t face_to_point_t;
67  typedef typename super::face_to_edge_t face_to_edge_t;
68 
69  static const uint16_type numVertices = super::numVertices;
70  static const uint16_type numFaces = super::numFaces;
71  static const uint16_type numGeometricFaces = super::numGeometricFaces;
72  static const uint16_type numTopologicalFaces = super::numTopologicalFaces;
73  static const uint16_type numEdges = super::numEdges;
74  static const uint16_type numNormals = super::numNormals;
75 
76  static const uint16_type numPoints = super::numPoints;
77  static const uint16_type nbPtsPerVertex = super::nbPtsPerVertex;
78  static const uint16_type nbPtsPerEdge = super::nbPtsPerEdge;
79  static const uint16_type nbPtsPerFace = super::nbPtsPerFace;
80  static const uint16_type nbPtsPerVolume = super::nbPtsPerVolume;
81 
82  typedef typename node<value_type>::type node_type;
83  typedef typename matrix_node<value_type>::type points_type;
84  typedef points_type matrix_node_type;
85 
86  typedef node_type normal_type;
87  typedef ublas::vector<normal_type> normals_type;
88  typedef typename normals_type::const_iterator normal_const_iterator;
89 
90  typedef typename super::permutation_type permutation_type;
92 
96 
97  Reference()
98  :
99  super(),
100  M_id( 0 ),
101  M_vertices( nRealDim, numVertices ),
102  M_points( nRealDim, numPoints ),
103  M_normals( numNormals ),
104  M_tangents( numNormals ),
105  M_barycenter( nRealDim ),
106  M_barycenterfaces( nRealDim, numTopologicalFaces ),
107  M_meas( 0 )
108  {
109  M_vertices *= 0;
110  M_points *= 0;
111 
112  if ( nDim == 1 )
113  {
114  M_vertices( 0, 0 ) = -1.0;
115  M_vertices( 0, 1 ) = 1.0;
116 
117  M_points = make_line_points();
118  }
119 
120  if ( nDim == 2 )
121  {
122  M_vertices( 0, 0 ) = -1.0;
123  M_vertices( 1, 0 ) = -1.0;
124  M_vertices( 0, 1 ) = 1.0;
125  M_vertices( 1, 1 ) = -1.0;
126  M_vertices( 0, 2 ) = -1.0;
127  M_vertices( 1, 2 ) = 1.0;
128 
129  M_points = make_triangle_points();
130  }
131 
132  if ( nDim == 3 )
133  {
134  M_vertices( 0, 0 ) = -1.0;
135  M_vertices( 1, 0 ) = -1.0;
136  M_vertices( 2, 0 ) = -1.0;
137  M_vertices( 0, 1 ) = 1.0;
138  M_vertices( 1, 1 ) = -1.0;
139  M_vertices( 2, 1 ) = -1.0;
140  M_vertices( 0, 2 ) = -1.0;
141  M_vertices( 1, 2 ) = 1.0;
142  M_vertices( 2, 2 ) = -1.0;
143  M_vertices( 0, 3 ) = -1.0;
144  M_vertices( 1, 3 ) = -1.0;
145  M_vertices( 2, 3 ) = 1.0;
146 
147  M_points = make_tetrahedron_points();
148  }
149 
150  //std::cout << "P = " << M_points << "\n";
151  make_normals();
152 
153  if ( nDim == 2 )
154  make_tangents();
155 
156  computeBarycenters();
157  computeMeasure();
158  }
159 
160  Reference( element_type const& e, uint16_type __f )
161  :
162  super(),
163  M_id( __f ),
164  M_vertices( nRealDim, numVertices ),
165  M_points( nRealDim, numPoints ),
166  M_normals( numNormals ),
167  M_tangents( numNormals ),
168  M_barycenter( nRealDim ),
169  M_barycenterfaces( nRealDim, numTopologicalFaces ),
170  M_meas( 0 )
171  {
172  if ( __f >= element_type::numTopologicalFaces )
173  {
174  std::ostringstream str;
175  str << "invalid face number " << __f << "\n"
176  << "must be 0 <= f < " << element_type::numTopologicalFaces << "\n";
177  throw std::invalid_argument( str.str() );
178  }
179 
180  for ( int i = 0; i < numVertices; ++i )
181  {
182  if ( real_dimension == 3 )
183  ublas::column( M_vertices, i ) = e.vertex( element_type::f2p( __f, i ) );
184 
185  else
186  ublas::column( M_vertices, i ) = e.vertex( element_type::e2p( __f, i ) );
187  }
188 
189  for ( int i = 0; i < numPoints; ++i )
190  {
191  if ( real_dimension == 3 )
192  ublas::column( M_points, i ) = e.point( element_type::f2p( __f, i ) );
193 
194  else
195  ublas::column( M_points, i ) = e.point( element_type::e2p( __f, i ) );
196  }
197 
198  make_normals();
199  computeBarycenters();
200  computeMeasure();
201  }
202 
203  Reference( Reference const & r )
204  :
205  super( r ),
206  M_id( r.M_id ),
207  M_vertices( r.M_vertices ),
208  M_points( r.M_points ),
209  M_normals( r.M_normals ),
210  M_tangents( r.M_tangents ),
211  M_barycenter( r.M_barycenter ),
212  M_barycenterfaces( r.M_barycenterfaces ),
213  M_meas( r.M_meas )
214  {
215 
216  }
217 
218  ~Reference() {}
219 
221 
225 
226  Reference& operator=( Reference const& r )
227  {
228  if ( this != &r )
229  {
230  M_id = r.M_id;
231  M_vertices = r.M_vertices;
232  M_points = r.M_points;
233  M_normals = r.M_normals;
234  M_tangents = r.M_tangents;
235  M_barycenter = r.M_barycenter;
236  M_barycenterfaces = r.M_barycenterfaces;
237  M_meas = r.M_meas;
238  }
239 
240  return *this;
241  }
242 
244 
248 
249  uint16_type topologicalDimension() const
250  {
251  return topological_dimension;
252  }
253  uint16_type dimension() const
254  {
255  return real_dimension;
256  }
257 
258  uint16_type nVertices() const
259  {
260  return numVertices;
261  }
262  uint16_type nPoints() const
263  {
264  return numPoints;
265  }
266  uint16_type nEdges() const
267  {
268  return numEdges;
269  }
270  uint16_type nFaces() const
271  {
272  return numFaces;
273  }
274 
275  points_type const& vertices() const
276  {
277  return M_vertices;
278  }
279 
280  ublas::matrix_column<points_type const> vertex( uint16_type __i ) const
281  {
282  return ublas::column( M_vertices, __i );
283  }
284 
285  ublas::matrix_column<points_type const> edgeVertex( uint16_type __e, uint16_type __p ) const
286  {
287  return ublas::column( M_vertices, edge_to_point_t::e2p( __e,__p ) );
288  }
289 
293  ublas::matrix_column<points_type const> faceVertex( uint16_type f, uint16_type p ) const
294  {
295  return ublas::column( M_vertices, face_to_point_t::f2p( f,p ) );
296  }
297 
301  matrix_node_type faceVertices( uint16_type f ) const
302  {
303  matrix_node_type v( nRealDim, nDim );
304 
305  // there is exactely nDim vertices on each face on a d-simplex
306  for ( int p = 0; p < nDim; ++p )
307  {
308  switch ( nDim )
309  {
310  case 1:
311  case 3:
312  ublas::column( v, p ) = ublas::column( M_vertices, face_to_point_t::f2p( f,p ) );
313  break;
314 
315  case 2:
316  ublas::column( v, p ) = ublas::column( M_vertices, face_to_point_t::e2p( f,p ) );
317  break;
318  }
319  }
320 
321  return v;
322  }
323 
324  points_type const& points() const
325  {
326  return M_points;
327  }
328 
329  ublas::matrix_column<points_type const> point( uint16_type __i ) const
330  {
331  return ublas::column( M_points, __i );
332  }
333 
337  node_type barycenter() const
338  {
339  return M_barycenter;
340  }
341 
345  points_type barycenterFaces() const
346  {
347  return M_barycenterfaces;
348  }
349 
353  ublas::matrix_column<matrix_node_type const> faceBarycenter( uint16_type f ) const
354  {
355  return ublas::column( M_barycenterfaces, f );
356  }
357 
364  normals_type const& normals() const
365  {
366  return M_normals;
367  }
368 
376  node_type const& normal( uint16_type __n ) const
377  {
378  return M_normals[__n];
379  }
380 
388  node_type const& tangent( uint16_type __n ) const
389  {
390  return M_tangents[__n];
391  }
392 
399  normal_const_iterator beginNormal() const
400  {
401  return M_normals.begin();
402  }
403 
410  normal_const_iterator endNormal() const
411  {
412  return M_normals.end();
413  }
414 
415  topological_face_type topologicalFace( uint16_type __f ) const
416  {
417  topological_face_type ref( *this, __f );
418  return ref;
419  }
420 
421  points_type const& G() const
422  {
423  return M_points;
424  }
425 
429  double measure() const
430  {
431  return M_meas;
432  }
433 
434  size_type id() const
435  {
436  return 0;
437  }
438  flag_type marker() const
439  {
440  return 0;
441  }
442  flag_type marker2() const
443  {
444  return 0;
445  }
446  flag_type marker3() const
447  {
448  return 0;
449  }
450  uint16_type ad_first() const
451  {
452  return -1;
453  }
454  uint16_type ad_second() const
455  {
456  return -1;
457  }
458  uint16_type pos_first() const
459  {
460  return -1;
461  }
462  uint16_type pos_second() const
463  {
464  return -1;
465  }
466  permutation_type permutation( uint16_type /*f*/ ) const
467  {
468  return permutation_type();
469  }
470  double h() const
471  {
472  // FIXME: should be computed once for all in constructor
473  double __max = 0.0;
474 
475  for ( int __e = 0; __e < numEdges; ++ __e )
476  {
477  double __len = ublas::norm_2( edgeVertex( __e, 1 ) - edgeVertex( __e, 0 ) );
478  __max = ( __max < __len )?__len:__max;
479  }
480 
481  return __max;
482  }
483  double h( int e ) const
484  {
485  return ublas::norm_2( edgeVertex( e, 1 ) - edgeVertex( e, 0 ) );
486  }
487 
488 
489  double hFace( int /*__f*/ ) const
490  {
491 #if 0
492 
493  // FIXME: should be computed once for all in constructor
494  if ( nDim == 1 )
495  return 0.0;
496 
497  else
498  return face( __f ).h();
499 
500 #else
501  return 0.0;
502 #endif
503  }
504 
505  EntityRange<self_type> entityRange( uint16_type d ) const
506  {
507  return EntityRange<self_type>( d );
508  }
509 
511 
515 
517 
521 
522  points_type makePoints( uint16_type topo_dim, uint16_type __id, int interior = 1 ) const
523  {
524  // vertices
525  if ( topo_dim == 0 )
526  {
527 
528  //std::cerr << "coucpu 2 " << topo_dim << " " << topological_dimension << " " << __id << "\n";
529  points_type G( M_vertices.size1(), 1 );
530  ublas::column( G, 0 ) = ublas::column( M_vertices, __id );
531  return G;
532  }
533 
534  // interior points of the convex
535  else if ( topo_dim == topological_dimension )
536  {
537  //std::cerr << "coucpu 2 " << topo_dim << " " << topological_dimension << " " << __id << "\n";
538  if ( __id == 0 )
539  return this->template makeLattice<Shape>( interior );
540 
541  throw std::logic_error( "cannot make those points" );
542  return points_type();
543  }
544 
545  // all the other points
546  else
547  {
548  points_type G;
549 
550  if ( topo_dim == 1 )
551  {
552  Reference<Simplex<1, Order, 1>, 1, Order, 1, T> refline;
553  G = refline.template makeLattice<SHAPE_LINE>( interior );
554  pt_to_entity<Shape,1> p_to_e( __id );
555  points_type Gret( nRealDim, G.size2() );
556 
557  for ( size_type i = 0; i < G.size2(); ++i )
558  ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
559 
560  return Gret;
561  }
562 
563  else if ( topo_dim == 2 )
564  {
565  Reference<Simplex<2, Order, 2>, 2, Order, 2, T> refface;
566  G = refface.template makeLattice<SHAPE_TRIANGLE>( interior );
567  pt_to_entity<Shape,2> p_to_e( __id );
568  points_type Gret( nRealDim, G.size2() );
569 
570  for ( size_type i = 0; i < G.size2(); ++i )
571  ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
572 
573  return Gret;
574  }
575  }
576 
577  return points_type();
578  }
579 
580  template<size_type shape>
581  points_type
582  makeLattice( uint16_type interior = 0 ) const
583  {
584  if ( nOrder > 0 )
585  {
586  if ( shape == SHAPE_LINE )
587  return make_line_points( interior );
588 
589  else if ( shape == SHAPE_TRIANGLE )
590  return make_triangle_points( interior );
591 
592  else if ( shape == SHAPE_TETRA )
593  return make_tetrahedron_points( interior );
594  }
595 
596  else if ( nOrder == 0 )
597  return glas::average( M_vertices );
598  }
599 
652  boost::tuple<bool, value_type>
653  isIn( typename node<value_type>::type const& pt ) const
654  {
655  //std::cout << "[isIn] pt =" << pt << "\n";
656  ublas::vector<double> D( nDim + 1, 0 );
657  typename matrix_node<value_type>::type M( nDim+1,nDim+1 );
658  typename matrix_node<value_type>::type P( nDim,nDim+1 );
659  typename matrix_node<value_type>::type P0( nDim,nDim+1 );
660  std::fill( M.data().begin(), M.data().end(), value_type( 1.0 ) );
661 
662  for ( int p = 0; p < numVertices; ++p )
663  {
664  ublas::column( P0, p ) = this->vertex( p );
665  }
666 
667  ublas::subrange( M, 0, nDim, 0, nDim+1 ) = P0;
668  //VLOG(1) << "M=" << M << "\n";
669 
670 
671 
672  // compute the measure(times 2) of the segments, triangles or
673  // tetras formed by the point pt and the vertices of the
674  // simplex. The simplices are formed such that when the
675  // point is in the simplex the measure are all
676  // positive. In the case the point is outside the simplex
677  // the measure is negative
678  double meas_times = details::det( M, mpl::int_<nDim+1>() );
679 
680  for ( int n = 0; n < numVertices; ++n )
681  {
682  ublas::noalias( P ) = P0;
683  ublas::column( P, n ) = pt;
684  ublas::subrange( M, 0, nDim, 0, nDim+1 ) = P;
685  //VLOG(1) << "M=" << M << "\n";
686 
687  // multiply by -1 such that the volume of the
688  // reference is > 0
689  D( n )= meas_times*details::det( M, mpl::int_<nDim+1>() );
690 
691  //std::cout.setf( std::ios::scientific );
692  //std::cout.precision( 10 );
693  //std::cout << "D " << n << "=" << D(n) << "\n";
694  //VLOG(1) << "sign " << n << "=" << sign << "\n";
695  //meas_times = ( meas_times < res )?meas_times:res;
696  }
697 
698  double dmin = *std::min_element( D.begin(), D.end() );
699  ublas::vector<double>::const_iterator Dit = std::find_if( D.begin(), D.end(), lambda::_1 < -5e-7 );
700 
701  //VLOG(1) << "meas=" << meas_times << "\n";
702  //return meas_times;
703  // sign must be > 0 if the point is inside the simplex
704  return boost::make_tuple( Dit == D.end(), dmin );
705  }
707 
708 private:
709 
710  int n_line_points( int interior = 0 ) const
711  {
712  return std::max( 0, int( Order )+1-2*interior );
713  }
714  int n_triangle_points( int interior = 0 ) const
715  {
716  if ( interior == 1 )
717  return std::max( 0, ( int( Order )+1-2*interior )*( int( Order )-2*interior )/2 );
718 
719  return ( Order+1 )*( Order+2 )/2;
720  }
721  int n_tetrahedron_points( int interior = 0 ) const
722  {
723  if ( interior == 1 )
724  return std::max( 0, ( int( Order )+1-2*interior )*( int( Order )-2*interior )*( int( Order )-1-2*interior )/6 );
725 
726  return ( Order+1 )*( Order+2 )*( Order+3 )/6;
727  }
728 
729  points_type
730  make_line_points( int interior = 0 ) const
731  {
732  if ( nOrder > 0 )
733  {
734  ublas::vector<node_type> h ( 1 );
735  h( 0 ) = vertex( 1 ) - vertex( 0 );
736 
737  points_type p( nRealDim, n_line_points( interior ) );
738 
739  for ( int i = interior, indp = 0; i < int( Order )+1-interior; ++i, ++indp )
740  {
741  ublas::column( p, indp ) = vertex( 0 ) + ( h( 0 ) * value_type( i ) )/value_type( Order );
742  }
743 
744  return p;
745  }
746 
747  else
748  return glas::average( M_vertices );
749 
750  }
751 
752  points_type
753  make_triangle_points( int interior = 0 ) const
754  {
755  if ( nOrder > 0 )
756  {
757  ublas::vector<node_type> h ( 2 );
758  h *= 0;
759  h( 0 ) = vertex( 1 ) - vertex( 0 );
760  h( 1 ) = vertex( 2 ) - vertex( 0 );
761  //std::cout << "h = " << h << "\n";
762  //DVLOG(2) << "n triangle pts = " << n_triangle_points( interior ) << "\n";
763  points_type G( nRealDim, n_triangle_points( interior ) );
764 
765  for ( int i = interior, p = 0; i < int( Order )+1-interior; ++i )
766  {
767  for ( int j = interior; j < int( Order ) + 1 - i-interior; ++j, ++p )
768  {
769  ublas::column( G, p ) = vertex( 0 ) + ( value_type( j ) * h( 0 )+
770  value_type( i ) * h( 1 ) )/ value_type( Order );
771  }
772  }
773 
774  return G;
775  }
776 
777  else
778  return glas::average( M_vertices );
779  }
780 
781  points_type
782  make_tetrahedron_points( int interior = 0 ) const
783  {
784  if ( nOrder > 0 )
785  {
786  ublas::vector<node_type> h ( 3 );
787  h( 0 ) = vertex( 1 ) - vertex( 0 );
788  h( 1 ) = vertex( 2 ) - vertex( 0 );
789  h( 2 ) = vertex( 3 ) - vertex( 0 );
790  points_type G( 3, n_tetrahedron_points( interior ) );
791 
792  //DVLOG(2) << "n tetra pts = " << n_tetrahedron_points( interior ) << "\n";
793  for ( int i = interior, p = 0; i < int( Order )+1-interior; ++i )
794  {
795  for ( int j = interior; j < int( Order ) + 1 - i - interior; ++j )
796  {
797  for ( int k = interior; k < int( Order ) + 1 - i - j - interior; ++k, ++p )
798  {
799  ublas::column( G, p ) = vertex( 0 ) + ( value_type( i ) * h( 2 ) +
800  value_type( j ) * h( 1 ) +
801  value_type( k ) * h( 0 ) ) / value_type( Order );
802 
803  }
804  }
805  }
806 
807  //std::cout << "G = " << G << "\n";
808  return G;
809  }
810 
811  else
812  return glas::average( M_vertices );
813  }
814 
815  template<size_type shape>
816  struct pt_to_edge
817  {
818  pt_to_edge( std::vector<uint16_type> vert_ids )
819  :
820  h( 1.0 ),
821  a( Entity<SHAPE_LINE, value_type>().vertex( 0 ) ),
822  b( Entity<SHAPE_LINE, value_type>().vertex( 1 ) ),
823  u( Entity<shape, value_type>().vertex( vert_ids[ 0 ] ) ),
824  v( Entity<shape, value_type>().vertex( vert_ids[ 1 ] ) ),
825  diff( v-u )
826  {
827  h = 1.0/( b[0]-a[0] );
828  }
829  node_type
830  operator()( node_type const& x ) const
831  {
832  return u + h * ( x[ 0 ] - a[ 0 ] ) * diff;
833  }
834  value_type h;
835  node_type a, b;
836  node_type u, v, diff;
837  };
838 
839  //
840  // pt_to_face tetrahedron
841  //
842  template<size_type shape>
843  struct pt_to_face
844  {
845  pt_to_face( std::vector<uint16_type> vert_ids )
846  :
847  u( Entity<shape, value_type>().vertex( vert_ids[ 0 ] ) ),
848  v( Entity<shape, value_type>().vertex( vert_ids[ 1 ] ) ),
849  w( Entity<shape, value_type>().vertex( vert_ids[ 2 ] ) ),
850  diff( 2 )
851  {
852  diff[0] = v-u;
853  diff[1] = w-u;
854  }
855  node_type
856  operator()( node_type const& x ) const
857  {
858  return u + 0.5*( x[ 0 ]+1.0 ) * diff[ 0 ] + 0.5*( x[ 1 ]+1.0 ) * diff[ 1 ];
859  }
860  node_type u, v, w;
861  ublas::vector<node_type> diff;
862  };
863 
864  template<size_type shape>
865  struct pt_to_element
866  {
867  pt_to_element() {}
868  pt_to_element( std::vector<uint16_type> const& ) {}
869  node_type operator()( node_type const& x ) const
870  {
871  return x;
872  }
873  };
874 
875  template<size_type shape,uint16_type topo_dim>
876  struct pt_to_entity
877  {
878  typedef typename mpl::if_<mpl::equal_to<mpl::size_t<shape>, mpl::size_t<SHAPE_LINE> >,
879  mpl::identity<mpl::vector<boost::none_t,pt_to_edge<shape>,pt_to_edge<shape> > >,
880  typename mpl::if_<mpl::equal_to<mpl::size_t<shape>, mpl::size_t<SHAPE_TRIANGLE> >,
881  mpl::identity<mpl::vector<boost::none_t, pt_to_edge<shape>, pt_to_element<shape> > >,
882  mpl::identity<mpl::vector<boost::none_t, pt_to_edge<shape>, pt_to_face<shape>, pt_to_element<shape> > >
883  >::type // 2
884  >::type::type _type;
885  typedef typename mpl::at<_type, mpl::int_<topo_dim> >::type mapping_type;
886  typedef mpl::vector<boost::none_t, edge_to_point_t, face_to_point_t> list_v;
887 
888  pt_to_entity( uint16_type entity_id )
889  :
890  mapping( typename mpl::at<list_v, mpl::int_<topo_dim> >::type().entity( topo_dim, entity_id ) )
891  {}
892 
893  node_type operator()( node_type const& x ) const
894  {
895  return mapping( x );
896  }
897  mapping_type mapping;
898  };
899 
900  void make_tangents()
901  {
902  uint16_type reindex1[][4] = { { 1, 0, 0, 0},
903  { 0, 1, 2, 0},
904  { 0, 1, 2, 3}
905  };
906 
907  for ( uint16_type __n = 0; __n < numNormals; ++__n )
908  {
909  const int ind_normal = reindex1[nDim-1][__n];
910  M_tangents[ind_normal].resize( nDim );
911  }
912 
913  M_tangents[0][0]=-1/math::sqrt( 2. );
914  M_tangents[0][1]= 1/math::sqrt( 2. );
915  M_tangents[1][0]= 0;
916  M_tangents[1][1]= -1;
917  M_tangents[2][0]= 1;
918  M_tangents[2][1]= 0;
919  }
920  void make_normals()
921  {
922  /*uint16_type reindex[][4] = { { 1, 0, 0, 0},
923  { 1, 2, 0, 0},
924  { 2, 3, 1, 0} };
925  */
926  uint16_type reindex1[][4] = { { 1, 0, 0, 0},
927  { 0, 1, 2, 0},
928  { 0, 1, 2, 3}
929  };
930 
931 
932 #if 0
933  node_type zero = ublas::zero_vector<value_type>( nDim );
934  uint16_type d = nDim;
935  std::for_each( M_normals.begin(), M_normals.end(),
936  ( lambda::bind( &node_type::resize, lambda::_1,
937  lambda::constant( d ), false ),
938  lambda::_1 = lambda::constant( zero ) ) );
939 #endif
940 
941  for ( uint16_type __n = 0; __n < numNormals; ++__n )
942  {
943  //const uint16_type ind_normal = reindex[nDim-1][__n];
944  const int ind_normal = reindex1[nDim-1][__n];
945  M_normals[ind_normal].resize( nDim );
946  M_normals[ind_normal].clear();// = ublas::zero_vector<value_type>( nDim );
947 
948  if ( __n > 0 )
949  {
950  M_normals[ind_normal][__n-1] = -1;
951  }
952 
953  else
954  {
955  typedef typename mpl::if_<mpl::equal_to<mpl::int_<nDim>, mpl::int_<0> >,
956  mpl::int_<1>,
957  mpl::int_<nDim> >::type denom;
958  M_normals[ind_normal] = ublas::scalar_vector<value_type>( nDim, math::sqrt( value_type( 1.0 )/value_type( denom::value ) ) );
959  }
960 
961  //DVLOG(2) << "normal[" << ind_normal << "]=" << M_normals[ind_normal] << "\n";
962  }
963  }
964 
965  void computeBarycenters();
966  void computeMeasure();
967 private:
968 
969  uint16_type M_id;
970 
971  points_type M_vertices;
972 
973  points_type M_points;
974 
975  normals_type M_normals;
976  normals_type M_tangents;
977 
978  node_type M_barycenter;
979 
980  points_type M_barycenterfaces;
981 
982  value_type M_meas;
983 };
984 
985 template<uint16_type Dim, uint16_type Order, uint16_type RDim, typename T>
986 const uint16_type Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T>::nbPtsPerVertex;
987 template<uint16_type Dim, uint16_type Order, uint16_type RDim, typename T>
988 const uint16_type Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T>::nbPtsPerEdge;
989 template<uint16_type Dim, uint16_type Order, uint16_type RDim, typename T>
990 const uint16_type Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T>::nbPtsPerFace;
991 template<uint16_type Dim, uint16_type Order, uint16_type RDim, typename T>
992 const uint16_type Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T>::numGeometricFaces;
993 
994 
995 
996 template<typename T> class Entity<SHAPE_LINE, T>: public Reference<Simplex<1, 1, 1>,1,1, 1, T> {};
997 template<typename T> class Entity<SHAPE_TRIANGLE, T>: public Reference<Simplex<2, 1, 2>,2,1, 2, T> {};
998 template<typename T> class Entity<SHAPE_TETRA, T>: public Reference<Simplex<3, 1, 3>,3,1, 3, T> {};
999 
1000 
1001 template<uint16_type Dim, uint16_type Order, uint16_type RDim, typename T>
1002 void
1003 Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T>::computeBarycenters()
1004 {
1005  M_barycenter = ublas::column( glas::average( M_vertices ), 0 );
1006 
1007  for ( int f = 0; f < numTopologicalFaces; ++f )
1008  {
1009  //std::cout << "face " << f << " vertices " << faceVertices( f ) << "\n";
1010  ublas::column( M_barycenterfaces, f ) = ublas::column( glas::average( faceVertices( f ) ), 0 );
1011  }
1012 }
1013 template<uint16_type Dim, uint16_type Order, uint16_type RDim, typename T>
1014 void
1015 Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T>::computeMeasure()
1016 {
1017  if ( nDim == nRealDim )
1018  {
1019 
1020  typename matrix_node<value_type>::type M( nDim,nDim );
1021  value_type factor( 1 );
1022 
1023  switch ( nDim )
1024  {
1025  case 1:
1026  ublas::column( M, 0 ) = this->vertex( 1 )-this->vertex( 0 );
1027  factor = 1;
1028  break;
1029 
1030  case 2:
1031  ublas::column( M, 0 ) = this->vertex( 0 )-this->vertex( 1 );
1032  ublas::column( M, 1 ) = this->vertex( 1 )-this->vertex( 2 );
1033  factor = 2;
1034  break;
1035 
1036  case 3:
1042  ublas::column( M, 0 ) = this->vertex( 0 )-this->vertex( 1 );
1043  ublas::column( M, 1 ) = this->vertex( 1 )-this->vertex( 2 );
1044  ublas::column( M, 2 ) = this->vertex( 2 )-this->vertex( 3 );
1045  factor = 6;
1046  break;
1047  }
1048 
1049  M_meas = math::abs( details::det( M, mpl::int_<nDim>() ) )/factor;
1050  }
1051 
1052  else
1053  {
1054 #if 0
1055 
1061  typename matrix_node<value_type>::type M( nRealDim,nRealDim );
1062  value_type factor( 1 );
1063 
1064  switch ( nDim )
1065  {
1066  case 1:
1067  M_meas = ublas::norm2( this->vertex( 1 )-this->vertex( 0 ) );
1068  break;
1069 
1070  case 2:
1071  ublas::column( M, 0 ) = this->vertex( 0 )-this->vertex( 1 );
1072  ublas::column( M, 1 ) = this->vertex( 1 )-this->vertex( 2 );
1073  factor = 2;
1074  break;
1075  }
1076 
1077 #endif
1078  }
1079 }
1080 
1081 }
1082 #endif /* __refsimplex_H */

Generated on Sun Oct 20 2013 08:25:04 for Feel++ by doxygen 1.8.4