29 #ifndef __refsimplex_H
30 #define __refsimplex_H 1
34 template<u
int16_type Dim, u
int16_type Order, u
int16_type RDim,
typename T>
35 class Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T>
37 public Simplex<Dim, Order, RDim>
40 typedef Simplex<Dim, Order, RDim> super;
46 static const uint16_type nDim = super::nDim;
47 static const uint16_type nOrder = super::nOrder;
48 static const uint16_type nRealDim = super::nRealDim;
50 static const uint16_type topological_dimension = super::topological_dimension;
51 static const uint16_type real_dimension = super::real_dimension;
53 static const size_type Shape = super::Shape;
54 static const size_type Geometry = super::Geometry;
57 typedef Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T> self_type;
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;
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;
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;
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;
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;
86 typedef node_type normal_type;
87 typedef ublas::vector<normal_type> normals_type;
88 typedef typename normals_type::const_iterator normal_const_iterator;
90 typedef typename super::permutation_type permutation_type;
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 ),
114 M_vertices( 0, 0 ) = -1.0;
115 M_vertices( 0, 1 ) = 1.0;
117 M_points = make_line_points();
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;
129 M_points = make_triangle_points();
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;
147 M_points = make_tetrahedron_points();
156 computeBarycenters();
160 Reference( element_type
const& e, uint16_type __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 ),
172 if ( __f >= element_type::numTopologicalFaces )
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() );
180 for (
int i = 0; i < numVertices; ++i )
182 if ( real_dimension == 3 )
183 ublas::column( M_vertices, i ) = e.vertex( element_type::f2p( __f, i ) );
186 ublas::column( M_vertices, i ) = e.vertex( element_type::e2p( __f, i ) );
189 for (
int i = 0; i < numPoints; ++i )
191 if ( real_dimension == 3 )
192 ublas::column( M_points, i ) = e.point( element_type::f2p( __f, i ) );
195 ublas::column( M_points, i ) = e.point( element_type::e2p( __f, i ) );
199 computeBarycenters();
203 Reference( Reference
const & r )
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 ),
226 Reference&
operator=( Reference
const& r )
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;
249 uint16_type topologicalDimension()
const
251 return topological_dimension;
253 uint16_type dimension()
const
255 return real_dimension;
258 uint16_type nVertices()
const
262 uint16_type nPoints()
const
266 uint16_type nEdges()
const
270 uint16_type nFaces()
const
275 points_type
const& vertices()
const
280 ublas::matrix_column<points_type const> vertex( uint16_type __i )
const
282 return ublas::column( M_vertices, __i );
285 ublas::matrix_column<points_type const> edgeVertex( uint16_type __e, uint16_type __p )
const
287 return ublas::column( M_vertices, edge_to_point_t::e2p( __e,__p ) );
293 ublas::matrix_column<points_type const> faceVertex( uint16_type f, uint16_type p )
const
295 return ublas::column( M_vertices, face_to_point_t::f2p( f,p ) );
301 matrix_node_type faceVertices( uint16_type f )
const
303 matrix_node_type v( nRealDim, nDim );
306 for (
int p = 0; p < nDim; ++p )
312 ublas::column( v, p ) = ublas::column( M_vertices, face_to_point_t::f2p( f,p ) );
316 ublas::column( v, p ) = ublas::column( M_vertices, face_to_point_t::e2p( f,p ) );
324 points_type
const&
points()
const
329 ublas::matrix_column<points_type const> point( uint16_type __i )
const
331 return ublas::column( M_points, __i );
337 node_type barycenter()
const
345 points_type barycenterFaces()
const
347 return M_barycenterfaces;
353 ublas::matrix_column<matrix_node_type const> faceBarycenter( uint16_type f )
const
355 return ublas::column( M_barycenterfaces, f );
364 normals_type
const& normals()
const
376 node_type
const& normal( uint16_type __n )
const
378 return M_normals[__n];
388 node_type
const& tangent( uint16_type __n )
const
390 return M_tangents[__n];
399 normal_const_iterator beginNormal()
const
401 return M_normals.begin();
410 normal_const_iterator endNormal()
const
412 return M_normals.end();
415 topological_face_type topologicalFace( uint16_type __f )
const
417 topological_face_type ref( *
this, __f );
421 points_type
const& G()
const
429 double measure()
const
438 flag_type marker()
const
442 flag_type marker2()
const
446 flag_type marker3()
const
450 uint16_type ad_first()
const
454 uint16_type ad_second()
const
458 uint16_type pos_first()
const
462 uint16_type pos_second()
const
466 permutation_type permutation( uint16_type )
const
468 return permutation_type();
475 for (
int __e = 0; __e < numEdges; ++ __e )
477 double __len = ublas::norm_2( edgeVertex( __e, 1 ) - edgeVertex( __e, 0 ) );
478 __max = ( __max < __len )?__len:__max;
483 double h(
int e )
const
485 return ublas::norm_2( edgeVertex( e, 1 ) - edgeVertex( e, 0 ) );
489 double hFace(
int )
const
498 return face( __f ).h();
505 EntityRange<self_type> entityRange( uint16_type d )
const
507 return EntityRange<self_type>( d );
522 points_type makePoints( uint16_type topo_dim, uint16_type __id,
int interior = 1 )
const
529 points_type G( M_vertices.size1(), 1 );
530 ublas::column( G, 0 ) = ublas::column( M_vertices, __id );
535 else if ( topo_dim == topological_dimension )
539 return this->
template makeLattice<Shape>( interior );
541 throw std::logic_error(
"cannot make those points" );
542 return points_type();
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() );
557 for (
size_type i = 0; i < G.size2(); ++i )
558 ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
563 else if ( topo_dim == 2 )
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() );
570 for (
size_type i = 0; i < G.size2(); ++i )
571 ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
577 return points_type();
580 template<
size_type shape>
582 makeLattice( uint16_type interior = 0 )
const
586 if ( shape == SHAPE_LINE )
587 return make_line_points( interior );
589 else if ( shape == SHAPE_TRIANGLE )
590 return make_triangle_points( interior );
592 else if ( shape == SHAPE_TETRA )
593 return make_tetrahedron_points( interior );
596 else if ( nOrder == 0 )
597 return glas::average( M_vertices );
652 boost::tuple<bool, value_type>
653 isIn(
typename node<value_type>::type
const& pt )
const
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 ) );
662 for (
int p = 0; p < numVertices; ++p )
664 ublas::column( P0, p ) = this->vertex( p );
667 ublas::subrange( M, 0, nDim, 0, nDim+1 ) = P0;
678 double meas_times = details::det( M, mpl::int_<nDim+1>() );
680 for (
int n = 0; n < numVertices; ++n )
682 ublas::noalias( P ) = P0;
683 ublas::column( P, n ) = pt;
684 ublas::subrange( M, 0, nDim, 0, nDim+1 ) = P;
689 D( n )= meas_times*details::det( M, mpl::int_<nDim+1>() );
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 );
704 return boost::make_tuple( Dit == D.end(), dmin );
710 int n_line_points(
int interior = 0 )
const
712 return std::max( 0,
int( Order )+1-2*interior );
714 int n_triangle_points(
int interior = 0 )
const
717 return std::max( 0, (
int( Order )+1-2*interior )*(
int( Order )-2*interior )/2 );
719 return ( Order+1 )*( Order+2 )/2;
721 int n_tetrahedron_points(
int interior = 0 )
const
724 return std::max( 0, (
int( Order )+1-2*interior )*(
int( Order )-2*interior )*(
int( Order )-1-2*interior )/6 );
726 return ( Order+1 )*( Order+2 )*( Order+3 )/6;
730 make_line_points(
int interior = 0 )
const
734 ublas::vector<node_type> h ( 1 );
735 h( 0 ) = vertex( 1 ) - vertex( 0 );
737 points_type p( nRealDim, n_line_points( interior ) );
739 for (
int i = interior, indp = 0; i < int( Order )+1-interior; ++i, ++indp )
741 ublas::column( p, indp ) = vertex( 0 ) + ( h( 0 ) * value_type( i ) )/value_type( Order );
748 return glas::average( M_vertices );
753 make_triangle_points(
int interior = 0 )
const
757 ublas::vector<node_type> h ( 2 );
759 h( 0 ) = vertex( 1 ) - vertex( 0 );
760 h( 1 ) = vertex( 2 ) - vertex( 0 );
763 points_type G( nRealDim, n_triangle_points( interior ) );
765 for (
int i = interior, p = 0; i < int( Order )+1-interior; ++i )
767 for (
int j = interior; j < int( Order ) + 1 - i-interior; ++j, ++p )
769 ublas::column( G, p ) = vertex( 0 ) + ( value_type( j ) * h( 0 )+
770 value_type( i ) * h( 1 ) )/ value_type( Order );
778 return glas::average( M_vertices );
782 make_tetrahedron_points(
int interior = 0 )
const
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 ) );
793 for (
int i = interior, p = 0; i < int( Order )+1-interior; ++i )
795 for (
int j = interior; j < int( Order ) + 1 - i - interior; ++j )
797 for (
int k = interior; k < int( Order ) + 1 - i - j - interior; ++k, ++p )
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 );
812 return glas::average( M_vertices );
815 template<
size_type shape>
818 pt_to_edge( std::vector<uint16_type> vert_ids )
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 ] ) ),
827 h = 1.0/( b[0]-a[0] );
830 operator()( node_type
const& x )
const
832 return u + h * ( x[ 0 ] - a[ 0 ] ) * diff;
836 node_type u, v, diff;
842 template<
size_type shape>
845 pt_to_face( std::vector<uint16_type> vert_ids )
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 ] ) ),
856 operator()( node_type
const& x )
const
858 return u + 0.5*( x[ 0 ]+1.0 ) * diff[ 0 ] + 0.5*( x[ 1 ]+1.0 ) * diff[ 1 ];
861 ublas::vector<node_type> diff;
864 template<
size_type shape>
868 pt_to_element( std::vector<uint16_type>
const& ) {}
869 node_type operator()( node_type
const& x )
const
875 template<
size_type shape,u
int16_type topo_dim>
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> > >
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;
888 pt_to_entity( uint16_type entity_id )
890 mapping( typename mpl::at<list_v, mpl::int_<topo_dim> >::type().entity( topo_dim, entity_id ) )
893 node_type operator()( node_type
const& x )
const
897 mapping_type mapping;
902 uint16_type reindex1[][4] = { { 1, 0, 0, 0},
907 for ( uint16_type __n = 0; __n < numNormals; ++__n )
909 const int ind_normal = reindex1[nDim-1][__n];
910 M_tangents[ind_normal].resize( nDim );
913 M_tangents[0][0]=-1/math::sqrt( 2. );
914 M_tangents[0][1]= 1/math::sqrt( 2. );
916 M_tangents[1][1]= -1;
926 uint16_type reindex1[][4] = { { 1, 0, 0, 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 ) ) );
941 for ( uint16_type __n = 0; __n < numNormals; ++__n )
944 const int ind_normal = reindex1[nDim-1][__n];
945 M_normals[ind_normal].resize( nDim );
946 M_normals[ind_normal].clear();
950 M_normals[ind_normal][__n-1] = -1;
955 typedef typename mpl::if_<mpl::equal_to<mpl::int_<nDim>, mpl::int_<0> >,
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 ) ) );
965 void computeBarycenters();
966 void computeMeasure();
971 points_type M_vertices;
973 points_type M_points;
975 normals_type M_normals;
976 normals_type M_tangents;
978 node_type M_barycenter;
980 points_type M_barycenterfaces;
985 template<u
int16_type Dim, u
int16_type Order, u
int16_type RDim,
typename T>
986 const uint16_type Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T>::nbPtsPerVertex;
987 template<u
int16_type Dim, u
int16_type Order, u
int16_type RDim,
typename T>
988 const uint16_type Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T>::nbPtsPerEdge;
989 template<u
int16_type Dim, u
int16_type Order, u
int16_type RDim,
typename T>
990 const uint16_type Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T>::nbPtsPerFace;
991 template<u
int16_type Dim, u
int16_type Order, u
int16_type RDim,
typename T>
992 const uint16_type Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T>::numGeometricFaces;
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> {};
1001 template<u
int16_type Dim, u
int16_type Order, u
int16_type RDim,
typename T>
1003 Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T>::computeBarycenters()
1005 M_barycenter = ublas::column( glas::average( M_vertices ), 0 );
1007 for (
int f = 0; f < numTopologicalFaces; ++f )
1010 ublas::column( M_barycenterfaces, f ) = ublas::column( glas::average( faceVertices( f ) ), 0 );
1013 template<u
int16_type Dim, u
int16_type Order, u
int16_type RDim,
typename T>
1015 Reference<Simplex<Dim, Order, RDim>, Dim, Order, RDim, T>::computeMeasure()
1017 if ( nDim == nRealDim )
1020 typename matrix_node<value_type>::type M( nDim,nDim );
1021 value_type factor( 1 );
1026 ublas::column( M, 0 ) = this->vertex( 1 )-this->vertex( 0 );
1031 ublas::column( M, 0 ) = this->vertex( 0 )-this->vertex( 1 );
1032 ublas::column( M, 1 ) = this->vertex( 1 )-this->vertex( 2 );
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 );
1049 M_meas = math::abs( details::det( M, mpl::int_<nDim>() ) )/factor;
1061 typename matrix_node<value_type>::type M( nRealDim,nRealDim );
1062 value_type factor( 1 );
1067 M_meas = ublas::norm2( this->vertex( 1 )-this->vertex( 0 ) );
1071 ublas::column( M, 0 ) = this->vertex( 0 )-this->vertex( 1 );
1072 ublas::column( M, 1 ) = this->vertex( 1 )-this->vertex( 2 );
boost::tuple< mpl::size_t< MESH_POINTS >, typename MeshTraits< MeshType >::point_const_iterator, typename MeshTraits< MeshType >::point_const_iterator > points(MeshType const &mesh)
Definition: filters.hpp:1296
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
Elements & operator=(Elements const &e)
Definition: elements.hpp:335