29 #ifndef __refhypercube_H
30 #define __refhypercube_H 1
34 template<u
int16_type Dim, u
int16_type Order, u
int16_type RDim,
typename T>
35 class Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T>
37 public Hypercube<Dim, Order, RDim>
46 typedef Hypercube<Dim, Order, RDim> super;
48 static const uint16_type nDim = super::nDim;
49 static const uint16_type nOrder = super::nOrder;
50 static const uint16_type nRealDim = super::nRealDim;
52 static const uint16_type topological_dimension = super::topological_dimension;
53 static const uint16_type real_dimension = super::real_dimension;
55 static const size_type Shape = super::Shape;
56 static const size_type Geometry = super::Geometry;
59 typedef Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T> self_type;
61 typedef typename mpl::if_<boost::is_same<typename super::element_type, boost::none_t>,
62 mpl::identity<boost::none_t>,
63 mpl::identity<Reference<typename super::element_type, nDim+1, nOrder, nRealDim, T> > >::type::type element_type;
64 typedef Reference<typename super::topological_face_type, super::topological_face_type::nDim, nOrder, nRealDim, T> topological_face_type;
66 typedef typename super::edge_to_point_t edge_to_point_t;
67 typedef typename super::face_to_point_t face_to_point_t;
68 typedef typename super::face_to_edge_t face_to_edge_t;
70 static const uint16_type numVertices = super::numVertices;
71 static const uint16_type numFaces = super::numFaces;
72 static const uint16_type numGeometricFaces = super::numGeometricFaces;
73 static const uint16_type numTopologicalFaces = super::numTopologicalFaces;
74 static const uint16_type numEdges = super::numEdges;
75 static const uint16_type numNormals = super::numNormals;
77 static const uint16_type numPoints = super::numPoints;
78 static const uint16_type nbPtsPerVertex = super::nbPtsPerVertex;
79 static const uint16_type nbPtsPerEdge = super::nbPtsPerEdge;
80 static const uint16_type nbPtsPerFace = super::nbPtsPerFace;
81 static const uint16_type nbPtsPerVolume = super::nbPtsPerVolume;
83 typedef typename node<value_type>::type node_type;
84 typedef typename matrix_node<value_type>::type points_type;
85 typedef points_type matrix_node_type;
87 typedef typename node<value_type>::type normal_type;
88 typedef ublas::vector<normal_type> normals_type;
89 typedef typename normals_type::const_iterator normal_const_iterator;
90 typedef typename super::permutation_type permutation_type;
101 M_vertices( nDim, numVertices ),
102 M_points( nDim, numPoints ),
103 M_normals( numNormals ),
104 M_barycenter( nDim ),
105 M_barycenterfaces( nDim, numTopologicalFaces ),
110 M_vertices( 0, 0 ) = -1.0;
111 M_vertices( 0, 1 ) = 1.0;
113 M_points = make_line_points();
118 M_vertices( 0, 0 ) = -1.0;
119 M_vertices( 1, 0 ) = -1.0;
121 M_vertices( 0, 1 ) = 1.0;
122 M_vertices( 1, 1 ) = -1.0;
124 M_vertices( 0, 2 ) = 1.0;
125 M_vertices( 1, 2 ) = 1.0;
127 M_vertices( 0, 3 ) = -1.0;
128 M_vertices( 1, 3 ) = 1.0;
130 M_points = make_quad_points();
136 M_vertices( 0, 0 ) = -1.0;
137 M_vertices( 1, 0 ) = -1.0;
138 M_vertices( 2, 0 ) = -1.0;
140 M_vertices( 0, 1 ) = 1.0;
141 M_vertices( 1, 1 ) = -1.0;
142 M_vertices( 2, 1 ) = -1.0;
144 M_vertices( 0, 2 ) = 1.0;
145 M_vertices( 1, 2 ) = 1.0;
146 M_vertices( 2, 2 ) = -1.0;
148 M_vertices( 0, 3 ) = -1.0;
149 M_vertices( 1, 3 ) = 1.0;
150 M_vertices( 2, 3 ) = -1.0;
153 M_vertices( 0, 4 ) = -1.0;
154 M_vertices( 1, 4 ) = -1.0;
155 M_vertices( 2, 4 ) = 1.0;
157 M_vertices( 0, 5 ) = 1.0;
158 M_vertices( 1, 5 ) = -1.0;
159 M_vertices( 2, 5 ) = 1.0;
161 M_vertices( 0, 6 ) = 1.0;
162 M_vertices( 1, 6 ) = 1.0;
163 M_vertices( 2, 6 ) = 1.0;
165 M_vertices( 0, 7 ) = -1.0;
166 M_vertices( 1, 7 ) = 1.0;
167 M_vertices( 2, 7 ) = 1.0;
169 M_points = make_hexa_points();
177 Reference( element_type
const& e, uint16_type __f )
181 M_vertices( nRealDim, numVertices ),
182 M_points( nRealDim, numPoints ),
183 M_normals( numNormals ),
184 M_barycenter( nDim ),
185 M_barycenterfaces( nDim, numTopologicalFaces ),
188 if ( __f >= element_type::numTopologicalFaces )
190 std::ostringstream str;
191 str <<
"invalid face number " << __f <<
"\n"
192 <<
"must be 0 <= f < " << element_type::numTopologicalFaces <<
"\n";
193 throw std::invalid_argument( str.str() );
196 for (
int i = 0; i < numVertices; ++i )
198 if ( real_dimension == 3 )
199 ublas::column( M_vertices, i ) = e.vertex( element_type::f2p( __f, i ) );
202 ublas::column( M_vertices, i ) = e.vertex( element_type::e2p( __f, i ) );
205 for (
int i = 0; i < numPoints; ++i )
207 if ( real_dimension == 3 )
208 ublas::column( M_points, i ) = e.point( element_type::f2p( __f, i ) );
211 ublas::column( M_points, i ) = e.point( element_type::e2p( __f, i ) );
214 M_points = M_vertices;
219 Reference( Reference
const & r )
223 M_vertices( r.M_vertices ),
224 M_points( r.M_points ),
225 M_normals( r.M_normals ),
226 M_barycenter( r.M_barycenter ),
227 M_barycenterfaces( r.M_barycenterfaces ),
241 Reference&
operator=( Reference
const& r )
246 M_vertices = r.M_vertices;
247 M_points = r.M_points;
248 M_normals = r.M_normals;
249 M_barycenter = r.M_barycenter;
250 M_barycenterfaces = r.M_barycenterfaces;
263 uint16_type topologicalDimension()
const
265 return topological_dimension;
267 uint16_type dimension()
const
269 return real_dimension;
272 uint16_type nVertices()
const
276 uint16_type nPoints()
const
280 uint16_type nEdges()
const
284 uint16_type nFaces()
const
289 points_type
const& vertices()
const
294 ublas::matrix_column<points_type const> vertex( uint16_type __i )
const
296 return ublas::column( M_vertices, __i );
299 ublas::matrix_column<points_type const> edgeVertex( uint16_type __e, uint16_type __p )
const
301 return ublas::column( M_vertices, edge_to_point_t::e2p( __e,__p ) );
304 points_type
const&
points()
const
309 ublas::matrix_column<points_type const> point( uint16_type __i )
const
311 return ublas::column( M_points, __i );
317 double measure()
const
325 matrix_node_type faceVertices( uint16_type f )
const
327 const int d[3] = { 1,2,4 };
328 matrix_node_type v( nDim, d[nDim-1] );
331 for (
int p = 0; p < d[nDim]; ++p )
337 ublas::column( v, p ) = ublas::column( M_vertices, face_to_point_t::f2p( f,p ) );
341 ublas::column( v, p ) = ublas::column( M_vertices, face_to_point_t::e2p( f,p ) );
352 node_type barycenter()
const
360 points_type barycenterFaces()
const
362 return M_barycenterfaces;
368 ublas::matrix_column<matrix_node_type const> faceBarycenter( uint16_type f )
const
370 return ublas::column( M_barycenterfaces, f );
379 normals_type
const& normals()
const
391 node_type
const& normal( uint16_type __n )
const
393 return M_normals[__n];
402 normal_const_iterator beginNormal()
const
404 return M_normals.begin();
413 normal_const_iterator endNormal()
const
415 return M_normals.end();
418 topological_face_type topologicalFace( uint16_type __f )
const
420 topological_face_type ref( *
this, __f );
424 points_type
const& G()
const
437 flag_type marker2()
const
441 flag_type marker3()
const
445 uint16_type ad_first()
const
449 uint16_type ad_second()
const
453 uint16_type pos_first()
const
457 uint16_type pos_second()
const
461 permutation_type permutation( uint16_type )
const
463 return permutation_type();
470 for (
int __e = 0; __e < numEdges; ++ __e )
472 double __len = ublas::norm_2( edgeVertex( __e, 1 ) - edgeVertex( __e, 0 ) );
473 __max = ( __max < __len )?__len:__max;
480 double hFace(
int )
const
489 return face( __f ).h();
496 EntityRange<self_type> entityRange( uint16_type d )
const
498 return EntityRange<self_type>( d );
517 boost::tuple<bool, value_type>
518 isIn(
typename node<value_type>::type
const& pt )
const
520 return isIn( pt, mpl::int_<nDim>() );
523 boost::tuple<bool, value_type>
524 isIn(
typename node<value_type>::type
const& pt, mpl::int_<1> )
const
526 typename matrix_node<value_type>::type M( nDim+1,nDim+1 );
527 typename matrix_node<value_type>::type P( nDim,nDim+1 );
528 std::fill( M.data().begin(), M.data().end(), value_type( 1.0 ) );
536 double meas_times = 1e30;
538 for (
int n = 0; n < numVertices; ++n )
540 ublas::column( P, 0 ) = pt;
541 ublas::column( P, 1 ) = this->vertex( n );
542 ublas::subrange( M, 0, nDim, 0, nDim+1 ) = P;
543 double res = details::det( M, mpl::int_<nDim+1>() )/2;
544 meas_times = ( meas_times < res )?meas_times:res;
547 return boost::make_tuple( meas_times>=-1e10, meas_times );
549 boost::tuple<bool, value_type>
550 isIn(
typename node<value_type>::type
const& pt, mpl::int_<2> )
const
552 typename matrix_node<value_type>::type M( nDim+1,nDim+1 );
553 typename matrix_node<value_type>::type P( nDim,nDim+1 );
554 std::fill( M.data().begin(), M.data().end(), value_type( 1.0 ) );
562 double meas_times = 1e30;
564 for (
int n = 0; n < numVertices; ++n )
566 ublas::column( P, 0 ) = pt;
567 ublas::column( P, 1 ) = this->vertex( n );
568 ublas::column( P, 2 ) = this->vertex( ( n+1 )%numVertices );
569 ublas::subrange( M, 0, nDim, 0, nDim+1 ) = P;
570 double res = details::det( M, mpl::int_<nDim+1>() )/2;
571 meas_times = ( meas_times < res )?meas_times:res;
574 return boost::make_tuple( meas_times>=-1e10, meas_times );
576 boost::tuple<bool, value_type>
577 isIn(
typename node<value_type>::type
const& pt, mpl::int_<3> )
const
579 typename matrix_node<value_type>::type M( nDim+1,nDim+1 );
580 typename matrix_node<value_type>::type P( nDim,nDim+1 );
581 std::fill( M.data().begin(), M.data().end(), value_type( 1.0 ) );
589 double meas_times = 1e30;
591 for (
int n = 0; n < numVertices; ++n )
593 ublas::column( P, 0 ) = pt;
594 ublas::column( P, 1 ) = this->vertex( n );
595 ublas::column( P, 2 ) = this->vertex( ( n+1 )%numVertices );
596 ublas::column( P, 3 ) = this->vertex( ( n+2 )%numVertices );
597 ublas::subrange( M, 0, nDim, 0, nDim+1 ) = P;
598 double res = details::det( M, mpl::int_<nDim+1>() )/2;
599 meas_times = ( meas_times < res )?meas_times:res;
602 return boost::make_tuple( meas_times>=-1e10, meas_times );
605 points_type makePoints( uint16_type topo_dim, uint16_type __id,
int interior = 1 )
const
612 points_type G( M_vertices.size1(), 1 );
613 ublas::column( G, 0 ) = ublas::column( M_vertices, __id );
618 else if ( topo_dim == topological_dimension )
622 return makeLattice<Shape>( interior );
624 throw std::logic_error(
"cannot make those points" );
625 return points_type();
635 Reference<Hypercube<1, Order, 1>, 1, Order, 1, T> refhyp1;
636 G = refhyp1.template makeLattice<SHAPE_LINE>( interior );
637 pt_to_entity<Shape,1> p_to_e( __id );
638 points_type Gret( nRealDim, G.size2() );
640 for (
size_type i = 0; i < G.size2(); ++i )
641 ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
646 else if ( topo_dim == 2 )
648 Reference<Hypercube<2, Order, 2>, 2, Order, 2, T> refhyp2;
649 G = refhyp2.template makeLattice<SHAPE_QUAD>( interior );
650 pt_to_entity<Shape,2> p_to_e( __id );
651 points_type Gret( nRealDim, G.size2() );
653 for (
size_type i = 0; i < G.size2(); ++i )
654 ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
660 return points_type();
663 template<
size_type shape>
665 makeLattice( uint16_type interior = 0 )
const
669 if ( shape == SHAPE_LINE )
670 return make_line_points( interior );
672 else if ( shape == SHAPE_QUAD )
673 return make_quad_points( interior );
675 else if ( shape == SHAPE_HEXA )
676 return make_hexa_points( interior );
679 else if ( nOrder == 0 )
680 return glas::average( M_vertices );
686 int n_line_points(
int interior = 0 )
const
688 return std::max( 0,
int( Order )+1-2*interior );
690 int n_quad_points(
int interior = 0 )
const
693 return std::max( 0, (
int( Order )+1-2*interior )*(
int( Order )+1-2*interior ) );
695 return ( Order+1 )*( Order+1 );
697 int n_hexa_points(
int interior = 0 )
const
700 return std::max( 0, (
int( Order )+1-2*interior )*(
int( Order )+1-2*interior )*(
int( Order )+1-2*interior ) );
702 return ( Order+1 )*( Order+1 )*( Order+1 );
707 make_line_points(
int interior = 0 )
const
711 ublas::vector<node_type> h ( 1 );
712 h( 0 ) = vertex( 1 ) - vertex( 0 );
714 points_type p( 1, n_line_points( interior ) );
716 for (
int i = interior, indp = 0; i < int( Order )+1-interior; ++i, ++indp )
718 ublas::column( p, indp ) = vertex( 0 ) + ( h( 0 ) * value_type( i ) )/value_type( Order );
725 return glas::average( M_vertices );
729 make_quad_points(
int interior = 0 )
const
733 ublas::vector<node_type> h ( 2 );
734 h( 0 ) = vertex( 1 ) - vertex( 0 );
735 h( 1 ) = vertex( 3 ) - vertex( 0 );
738 points_type G( 2, n_quad_points( interior ) );
740 for (
int i = interior, p = 0; i < int( Order )+1-interior; ++i )
742 for (
int j = interior; j < int( Order ) + 1 -interior; ++j, ++p )
744 ublas::column( G, p ) = vertex( 0 ) + ( value_type( i ) * h( 1 ) +
745 value_type( j ) * h( 0 ) )/ value_type( Order );
753 return glas::average( M_vertices );
757 make_hexa_points(
int interior = 0 )
const
761 ublas::vector<node_type> h ( 3 );
762 h( 0 ) = vertex( 1 ) - vertex( 0 );
763 h( 1 ) = vertex( 3 ) - vertex( 0 );
764 h( 2 ) = vertex( 4 ) - vertex( 0 );
765 points_type G( 3, n_hexa_points( interior ) );
768 for (
int i = interior, p = 0; i < int( Order )+1-interior; ++i )
770 for (
int j = interior; j < int( Order ) + 1 - interior; ++j )
772 for (
int k = interior; k < int( Order ) + 1 - interior; ++k, ++p )
774 ublas::column( G, p ) = vertex( 0 ) + ( value_type( i ) * h( 2 ) +
775 value_type( j ) * h( 1 ) +
776 value_type( k ) * h( 0 ) ) / value_type( Order );
787 return glas::average( M_vertices );
790 template<
size_type shape>
793 pt_to_edge( std::vector<uint16_type> vert_ids )
796 a( Entity<SHAPE_LINE, value_type>().vertex( 0 ) ),
797 b( Entity<SHAPE_LINE, value_type>().vertex( 1 ) ),
798 u( Entity<shape, value_type>().vertex( vert_ids[ 0 ] ) ),
799 v( Entity<shape, value_type>().vertex( vert_ids[ 1 ] ) ),
802 h = 1.0/( b[0]-a[0] );
805 operator()( node_type
const& x )
const
807 return u + h * ( x[ 0 ] - a[ 0 ] ) * diff;
811 node_type u, v, diff;
817 template<
size_type shape>
820 pt_to_face( std::vector<uint16_type> vert_ids )
822 u( Entity<shape, value_type>().vertex( vert_ids[ 0 ] ) ),
823 v( Entity<shape, value_type>().vertex( vert_ids[ 1 ] ) ),
824 w( Entity<shape, value_type>().vertex( vert_ids[ 3 ] ) ),
831 operator()( node_type
const& x )
const
833 return u + 0.5*( x[ 0 ]+1.0 ) * diff[ 0 ] + 0.5*( x[ 1 ]+1.0 ) * diff[ 1 ];
836 ublas::vector<node_type> diff;
839 template<
size_type shape>
843 pt_to_element( std::vector<uint16_type>
const& ) {}
844 node_type operator()( node_type
const& x )
const
850 template<
size_type shape,u
int16_type topo_dim>
853 typedef typename mpl::if_<mpl::equal_to<mpl::size_t<shape>, mpl::size_t<SHAPE_LINE> >,
854 mpl::identity<mpl::vector<boost::none_t,pt_to_edge<shape>,pt_to_edge<shape> > >,
855 typename mpl::if_<mpl::equal_to<mpl::size_t<shape>, mpl::size_t<SHAPE_QUAD> >,
856 mpl::identity<mpl::vector<boost::none_t, pt_to_edge<shape>, pt_to_element<shape> > >,
857 mpl::identity<mpl::vector<boost::none_t, pt_to_edge<shape>, pt_to_face<shape>, pt_to_element<shape> > >
860 typedef typename mpl::at<_type, mpl::int_<topo_dim> >::type mapping_type;
861 typedef mpl::vector<boost::none_t, edge_to_point_t, face_to_point_t> list_v;
863 pt_to_entity( uint16_type entity_id )
865 mapping( typename mpl::at<list_v, mpl::int_<topo_dim> >::type().entity( topo_dim, entity_id ) )
868 node_type operator()( node_type
const& x )
const
872 mapping_type mapping;
877 M_normals.resize( numNormals );
879 for (
int n = 0; n < numNormals; ++n )
881 M_normals[n].resize( nDim );
882 M_normals[n].clear();
887 M_normals[0][0] = -1;
893 M_normals[0][1] = -1;
896 M_normals[3][0] = -1;
901 M_normals[0][2] = -1;
902 M_normals[1][1] = -1;
905 M_normals[4][0] = -1;
911 void computeBarycenters();
914 void computeMeasure();
919 points_type M_vertices;
921 points_type M_points;
923 normals_type M_normals;
926 node_type M_barycenter;
928 points_type M_barycenterfaces;
933 template<u
int16_type Dim, u
int16_type Order, u
int16_type RDim,
typename T>
934 const uint16_type Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T>::nbPtsPerVertex;
935 template<u
int16_type Dim, u
int16_type Order, u
int16_type RDim,
typename T>
936 const uint16_type Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T>::nbPtsPerEdge;
937 template<u
int16_type Dim, u
int16_type Order, u
int16_type RDim,
typename T>
938 const uint16_type Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T>::nbPtsPerFace;
939 template<u
int16_type Dim, u
int16_type Order, u
int16_type RDim,
typename T>
940 const uint16_type Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T>::numGeometricFaces;
944 template<
typename T>
class Entity<SHAPE_QUAD, T>:
public Reference<Hypercube<2, 1, 2>, 2, 1, 2, T> {};
945 template<
typename T>
class Entity<SHAPE_HEXA, T>:
public Reference<Hypercube<3, 1, 3>, 3, 1, 3, T> {};
947 template<u
int16_type Dim, u
int16_type Order, u
int16_type RDim,
typename T>
949 Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T>::computeBarycenters()
951 M_barycenter = ublas::column( glas::average( M_vertices ), 0 );
953 for (
int f = 0; f < numTopologicalFaces; ++f )
956 ublas::column( M_barycenterfaces, f ) = ublas::column( glas::average( faceVertices( f ) ), 0 );
960 template<u
int16_type Dim, u
int16_type Order, u
int16_type RDim,
typename T>
962 Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T>::computeMeasure()
964 if ( nDim == nRealDim )
967 typename matrix_node<value_type>::type M( nDim,nDim );
982 ublas::column( M, 0 ) = this->vertex( 2 )-this->vertex( 0 );
983 ublas::column( M, 1 ) = this->vertex( 3 )-this->vertex( 1 );
991 ublas::column( M, 0 ) = this->vertex( 1 )-this->vertex( 0 );
992 ublas::column( M, 1 ) = this->vertex( 1 )-this->vertex( 2 );
993 ublas::column( M, 2 ) = this->vertex( 2 )-this->vertex( 3 );
1010 typename matrix_node<value_type>::type M( nRealDim,nRealDim );
1011 value_type factor( 1 );
1016 M_meas = ublas::norm2( this->vertex( 1 )-this->vertex( 0 ) );
1020 ublas::column( M, 0 ) = this->vertex( 0 )-this->vertex( 1 );
1021 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