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
refhypercube.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 __refhypercube_H
30 #define __refhypercube_H 1
31 
32 namespace Feel
33 {
34 template<uint16_type Dim, uint16_type Order, uint16_type RDim, typename T>
35 class Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T>
36  :
37 public Hypercube<Dim, Order, RDim>
38 {
39 public:
40 
41 
45 
46  typedef Hypercube<Dim, Order, RDim> super;
47 
48  static const uint16_type nDim = super::nDim;
49  static const uint16_type nOrder = super::nOrder;
50  static const uint16_type nRealDim = super::nRealDim;
51 
52  static const uint16_type topological_dimension = super::topological_dimension;
53  static const uint16_type real_dimension = super::real_dimension;
54 
55  static const size_type Shape = super::Shape;
56  static const size_type Geometry = super::Geometry;
57 
58  typedef T value_type;
59  typedef Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T> self_type;
60 
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;
65 
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;
69 
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;
76 
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;
82 
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;
86 
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;
92 
96 
97  Reference()
98  :
99  super(),
100  M_id( 0 ),
101  M_vertices( nDim, numVertices ),
102  M_points( nDim, numPoints ),
103  M_normals( numNormals ),
104  M_barycenter( nDim ),
105  M_barycenterfaces( nDim, numTopologicalFaces ),
106  M_meas( 0 )
107  {
108  if ( nDim == 1 )
109  {
110  M_vertices( 0, 0 ) = -1.0;
111  M_vertices( 0, 1 ) = 1.0;
112 
113  M_points = make_line_points();
114  }
115 
116  if ( nDim == 2 )
117  {
118  M_vertices( 0, 0 ) = -1.0;
119  M_vertices( 1, 0 ) = -1.0;
120 
121  M_vertices( 0, 1 ) = 1.0;
122  M_vertices( 1, 1 ) = -1.0;
123 
124  M_vertices( 0, 2 ) = 1.0;
125  M_vertices( 1, 2 ) = 1.0;
126 
127  M_vertices( 0, 3 ) = -1.0;
128  M_vertices( 1, 3 ) = 1.0;
129 
130  M_points = make_quad_points();
131  }
132 
133  if ( nDim == 3 )
134  {
135  // z = -1
136  M_vertices( 0, 0 ) = -1.0;
137  M_vertices( 1, 0 ) = -1.0;
138  M_vertices( 2, 0 ) = -1.0;
139 
140  M_vertices( 0, 1 ) = 1.0;
141  M_vertices( 1, 1 ) = -1.0;
142  M_vertices( 2, 1 ) = -1.0;
143 
144  M_vertices( 0, 2 ) = 1.0;
145  M_vertices( 1, 2 ) = 1.0;
146  M_vertices( 2, 2 ) = -1.0;
147 
148  M_vertices( 0, 3 ) = -1.0;
149  M_vertices( 1, 3 ) = 1.0;
150  M_vertices( 2, 3 ) = -1.0;
151 
152  // z = 1
153  M_vertices( 0, 4 ) = -1.0;
154  M_vertices( 1, 4 ) = -1.0;
155  M_vertices( 2, 4 ) = 1.0;
156 
157  M_vertices( 0, 5 ) = 1.0;
158  M_vertices( 1, 5 ) = -1.0;
159  M_vertices( 2, 5 ) = 1.0;
160 
161  M_vertices( 0, 6 ) = 1.0;
162  M_vertices( 1, 6 ) = 1.0;
163  M_vertices( 2, 6 ) = 1.0;
164 
165  M_vertices( 0, 7 ) = -1.0;
166  M_vertices( 1, 7 ) = 1.0;
167  M_vertices( 2, 7 ) = 1.0;
168 
169  M_points = make_hexa_points();
170  }
171 
172  //std::cout << "P = " << M_points << "\n";
173  make_normals();
174  computeMeasure();
175  }
176 
177  Reference( element_type const& e, uint16_type __f )
178  :
179  super(),
180  M_id( __f ),
181  M_vertices( nRealDim, numVertices ),
182  M_points( nRealDim, numPoints ),
183  M_normals( numNormals ),
184  M_barycenter( nDim ),
185  M_barycenterfaces( nDim, numTopologicalFaces ),
186  M_meas( 0 )
187  {
188  if ( __f >= element_type::numTopologicalFaces )
189  {
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() );
194  }
195 
196  for ( int i = 0; i < numVertices; ++i )
197  {
198  if ( real_dimension == 3 )
199  ublas::column( M_vertices, i ) = e.vertex( element_type::f2p( __f, i ) );
200 
201  else
202  ublas::column( M_vertices, i ) = e.vertex( element_type::e2p( __f, i ) );
203  }
204 
205  for ( int i = 0; i < numPoints; ++i )
206  {
207  if ( real_dimension == 3 )
208  ublas::column( M_points, i ) = e.point( element_type::f2p( __f, i ) );
209 
210  else
211  ublas::column( M_points, i ) = e.point( element_type::e2p( __f, i ) );
212  }
213 
214  M_points = M_vertices;
215  make_normals();
216  computeMeasure();
217  }
218 
219  Reference( Reference const & r )
220  :
221  super( r ),
222  M_id( r.M_id ),
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 ),
228  M_meas( r.M_meas )
229  {
230 
231  }
232 
233  ~Reference() {}
234 
236 
240 
241  Reference& operator=( Reference const& r )
242  {
243  if ( this != &r )
244  {
245  M_id = r.M_id;
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;
251  M_meas = r.M_meas;
252  }
253 
254  return *this;
255  }
256 
258 
262 
263  uint16_type topologicalDimension() const
264  {
265  return topological_dimension;
266  }
267  uint16_type dimension() const
268  {
269  return real_dimension;
270  }
271 
272  uint16_type nVertices() const
273  {
274  return numVertices;
275  }
276  uint16_type nPoints() const
277  {
278  return numPoints;
279  }
280  uint16_type nEdges() const
281  {
282  return numEdges;
283  }
284  uint16_type nFaces() const
285  {
286  return numFaces;
287  }
288 
289  points_type const& vertices() const
290  {
291  return M_vertices;
292  }
293 
294  ublas::matrix_column<points_type const> vertex( uint16_type __i ) const
295  {
296  return ublas::column( M_vertices, __i );
297  }
298 
299  ublas::matrix_column<points_type const> edgeVertex( uint16_type __e, uint16_type __p ) const
300  {
301  return ublas::column( M_vertices, edge_to_point_t::e2p( __e,__p ) );
302  }
303 
304  points_type const& points() const
305  {
306  return M_points;
307  }
308 
309  ublas::matrix_column<points_type const> point( uint16_type __i ) const
310  {
311  return ublas::column( M_points, __i );
312  }
313 
317  double measure() const
318  {
319  return M_meas;
320  }
321 
325  matrix_node_type faceVertices( uint16_type f ) const
326  {
327  const int d[3] = { 1,2,4 };
328  matrix_node_type v( nDim, d[nDim-1] );
329  // there is exactely nDim vertices on each face on a d-simplex
330 
331  for ( int p = 0; p < d[nDim]; ++p )
332  {
333  switch ( nDim )
334  {
335  case 1:
336  case 3:
337  ublas::column( v, p ) = ublas::column( M_vertices, face_to_point_t::f2p( f,p ) );
338  break;
339 
340  case 2:
341  ublas::column( v, p ) = ublas::column( M_vertices, face_to_point_t::e2p( f,p ) );
342  break;
343  }
344  }
345 
346  return v;
347  }
348 
352  node_type barycenter() const
353  {
354  return M_barycenter;
355  }
356 
360  points_type barycenterFaces() const
361  {
362  return M_barycenterfaces;
363  }
364 
368  ublas::matrix_column<matrix_node_type const> faceBarycenter( uint16_type f ) const
369  {
370  return ublas::column( M_barycenterfaces, f );
371  }
372 
379  normals_type const& normals() const
380  {
381  return M_normals;
382  }
383 
391  node_type const& normal( uint16_type __n ) const
392  {
393  return M_normals[__n];
394  }
395 
402  normal_const_iterator beginNormal() const
403  {
404  return M_normals.begin();
405  }
406 
413  normal_const_iterator endNormal() const
414  {
415  return M_normals.end();
416  }
417 
418  topological_face_type topologicalFace( uint16_type __f ) const
419  {
420  topological_face_type ref( *this, __f );
421  return ref;
422  }
423 
424  points_type const& G() const
425  {
426  return M_points;
427  }
428 
429  size_type id() const
430  {
431  return 0;
432  }
433  size_type marker() const
434  {
435  return 0;
436  }
437  flag_type marker2() const
438  {
439  return 0;
440  }
441  flag_type marker3() const
442  {
443  return 0;
444  }
445  uint16_type ad_first() const
446  {
447  return -1;
448  }
449  uint16_type ad_second() const
450  {
451  return -1;
452  }
453  uint16_type pos_first() const
454  {
455  return -1;
456  }
457  uint16_type pos_second() const
458  {
459  return -1;
460  }
461  permutation_type permutation( uint16_type /*f*/ ) const
462  {
463  return permutation_type();
464  }
465  double h() const
466  {
467  // FIXME: should be computed once for all in constructor
468  double __max = 0.0;
469 
470  for ( int __e = 0; __e < numEdges; ++ __e )
471  {
472  double __len = ublas::norm_2( edgeVertex( __e, 1 ) - edgeVertex( __e, 0 ) );
473  __max = ( __max < __len )?__len:__max;
474  }
475 
476  return __max;
477  }
478 
479 
480  double hFace( int /*__f*/ ) const
481  {
482 #if 0
483 
484  // FIXME: should be computed once for all in constructor
485  if ( nDim == 1 )
486  return 0.0;
487 
488  else
489  return face( __f ).h();
490 
491 #else
492  return 0.0;
493 #endif
494  }
495 
496  EntityRange<self_type> entityRange( uint16_type d ) const
497  {
498  return EntityRange<self_type>( d );
499  }
500 
502 
506 
508 
512 
517  boost::tuple<bool, value_type>
518  isIn( typename node<value_type>::type const& pt ) const
519  {
520  return isIn( pt, mpl::int_<nDim>() );
521  }
522 
523  boost::tuple<bool, value_type>
524  isIn( typename node<value_type>::type const& pt, mpl::int_<1> ) const
525  {
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 ) );
529 
530  // compute the measure(times 2) of the segments, triangles or
531  // tetras formed by the point pt and the vertices of the
532  // simplex. The simplices are formed such that when the
533  // point is in the simplex the measure are all
534  // positive. In the case the point is outside the simplex
535  // the measure is negative
536  double meas_times = 1e30;
537 
538  for ( int n = 0; n < numVertices; ++n )
539  {
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;
545  }
546 
547  return boost::make_tuple( meas_times>=-1e10, meas_times );
548  }
549  boost::tuple<bool, value_type>
550  isIn( typename node<value_type>::type const& pt, mpl::int_<2> ) const
551  {
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 ) );
555 
556  // compute the measure(times 2) of the segments, triangles or
557  // tetras formed by the point pt and the vertices of the
558  // simplex. The simplices are formed such that when the
559  // point is in the simplex the measure are all
560  // positive. In the case the point is outside the simplex
561  // the measure is negative
562  double meas_times = 1e30;
563 
564  for ( int n = 0; n < numVertices; ++n )
565  {
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;
572  }
573 
574  return boost::make_tuple( meas_times>=-1e10, meas_times );
575  }
576  boost::tuple<bool, value_type>
577  isIn( typename node<value_type>::type const& pt, mpl::int_<3> ) const
578  {
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 ) );
582 
583  // compute the measure(times 2) of the segments, triangles or
584  // tetras formed by the point pt and the vertices of the
585  // simplex. The simplices are formed such that when the
586  // point is in the simplex the measure are all
587  // positive. In the case the point is outside the simplex
588  // the measure is negative
589  double meas_times = 1e30;
590 
591  for ( int n = 0; n < numVertices; ++n )
592  {
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;
600  }
601 
602  return boost::make_tuple( meas_times>=-1e10, meas_times );
603  }
604 
605  points_type makePoints( uint16_type topo_dim, uint16_type __id, int interior = 1 ) const
606  {
607  // vertices
608  if ( topo_dim == 0 )
609  {
610 
611  //std::cerr << "coucpu 2 " << topo_dim << " " << topological_dimension << " " << __id << "\n";
612  points_type G( M_vertices.size1(), 1 );
613  ublas::column( G, 0 ) = ublas::column( M_vertices, __id );
614  return G;
615  }
616 
617  // interior points of the convex
618  else if ( topo_dim == topological_dimension )
619  {
620  //std::cerr << "coucpu 2 " << topo_dim << " " << topological_dimension << " " << __id << "\n";
621  if ( __id == 0 )
622  return makeLattice<Shape>( interior );
623 
624  throw std::logic_error( "cannot make those points" );
625  return points_type();
626  }
627 
628  // all the other points
629  else
630  {
631  points_type G;
632 
633  if ( topo_dim == 1 )
634  {
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() );
639 
640  for ( size_type i = 0; i < G.size2(); ++i )
641  ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
642 
643  return Gret;
644  }
645 
646  else if ( topo_dim == 2 )
647  {
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() );
652 
653  for ( size_type i = 0; i < G.size2(); ++i )
654  ublas::column( Gret, i ) = p_to_e( ublas::column( G, i ) );
655 
656  return Gret;
657  }
658  }
659 
660  return points_type();
661  }
662 
663  template<size_type shape>
664  points_type
665  makeLattice( uint16_type interior = 0 ) const
666  {
667  if ( nOrder > 0 )
668  {
669  if ( shape == SHAPE_LINE )
670  return make_line_points( interior );
671 
672  else if ( shape == SHAPE_QUAD )
673  return make_quad_points( interior );
674 
675  else if ( shape == SHAPE_HEXA )
676  return make_hexa_points( interior );
677  }
678 
679  else if ( nOrder == 0 )
680  return glas::average( M_vertices );
681  }
683 
684 private:
685 
686  int n_line_points( int interior = 0 ) const
687  {
688  return std::max( 0, int( Order )+1-2*interior );
689  }
690  int n_quad_points( int interior = 0 ) const
691  {
692  if ( interior == 1 )
693  return std::max( 0, ( int( Order )+1-2*interior )*( int( Order )+1-2*interior ) );
694 
695  return ( Order+1 )*( Order+1 );
696  }
697  int n_hexa_points( int interior = 0 ) const
698  {
699  if ( interior == 1 )
700  return std::max( 0, ( int( Order )+1-2*interior )*( int( Order )+1-2*interior )*( int( Order )+1-2*interior ) );
701 
702  return ( Order+1 )*( Order+1 )*( Order+1 );
703  }
704 
705 
706  points_type
707  make_line_points( int interior = 0 ) const
708  {
709  if ( nOrder > 0 )
710  {
711  ublas::vector<node_type> h ( 1 );
712  h( 0 ) = vertex( 1 ) - vertex( 0 );
713 
714  points_type p( 1, n_line_points( interior ) );
715 
716  for ( int i = interior, indp = 0; i < int( Order )+1-interior; ++i, ++indp )
717  {
718  ublas::column( p, indp ) = vertex( 0 ) + ( h( 0 ) * value_type( i ) )/value_type( Order );
719  }
720 
721  return p;
722  }
723 
724  else
725  return glas::average( M_vertices );
726  }
727 
728  points_type
729  make_quad_points( int interior = 0 ) const
730  {
731  if ( nOrder > 0 )
732  {
733  ublas::vector<node_type> h ( 2 );
734  h( 0 ) = vertex( 1 ) - vertex( 0 );
735  h( 1 ) = vertex( 3 ) - vertex( 0 );
736  //DVLOG(2) << "h = " << h << "\n";
737  //DVLOG(2) << "n quad pts = " << n_quad_points( interior ) << "\n";
738  points_type G( 2, n_quad_points( interior ) );
739 
740  for ( int i = interior, p = 0; i < int( Order )+1-interior; ++i )
741  {
742  for ( int j = interior; j < int( Order ) + 1 -interior; ++j, ++p )
743  {
744  ublas::column( G, p ) = vertex( 0 ) + ( value_type( i ) * h( 1 ) +
745  value_type( j ) * h( 0 ) )/ value_type( Order );
746  }
747  }
748 
749  return G;
750  }
751 
752  else
753  return glas::average( M_vertices );
754  }
755 
756  points_type
757  make_hexa_points( int interior = 0 ) const
758  {
759  if ( nOrder > 0 )
760  {
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 ) );
766 
767  //DVLOG(2) << "n hexa pts = " << n_hexa_points( interior ) << "\n";
768  for ( int i = interior, p = 0; i < int( Order )+1-interior; ++i )
769  {
770  for ( int j = interior; j < int( Order ) + 1 - interior; ++j )
771  {
772  for ( int k = interior; k < int( Order ) + 1 - interior; ++k, ++p )
773  {
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 );
777 
778  }
779  }
780  }
781 
782  //std::cout << "G = " << G << "\n";
783  return G;
784  }
785 
786  else
787  return glas::average( M_vertices );
788  }
789 
790  template<size_type shape>
791  struct pt_to_edge
792  {
793  pt_to_edge( std::vector<uint16_type> vert_ids )
794  :
795  h( 1.0 ),
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 ] ) ),
800  diff( v-u )
801  {
802  h = 1.0/( b[0]-a[0] );
803  }
804  node_type
805  operator()( node_type const& x ) const
806  {
807  return u + h * ( x[ 0 ] - a[ 0 ] ) * diff;
808  }
809  value_type h;
810  node_type a, b;
811  node_type u, v, diff;
812  };
813 
814  //
815  // pt_to_face hexa
816  //
817  template<size_type shape>
818  struct pt_to_face
819  {
820  pt_to_face( std::vector<uint16_type> vert_ids )
821  :
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 ] ) ),
825  diff( 2 )
826  {
827  diff[0] = v-u;
828  diff[1] = w-u;
829  }
830  node_type
831  operator()( node_type const& x ) const
832  {
833  return u + 0.5*( x[ 0 ]+1.0 ) * diff[ 0 ] + 0.5*( x[ 1 ]+1.0 ) * diff[ 1 ];
834  }
835  node_type u, v, w;
836  ublas::vector<node_type> diff;
837  };
838 
839  template<size_type shape>
840  struct pt_to_element
841  {
842  pt_to_element() {}
843  pt_to_element( std::vector<uint16_type> const& ) {}
844  node_type operator()( node_type const& x ) const
845  {
846  return x;
847  }
848  };
849 
850  template<size_type shape,uint16_type topo_dim>
851  struct pt_to_entity
852  {
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> > >
858  >::type // 2
859  >::type::type _type;
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;
862 
863  pt_to_entity( uint16_type entity_id )
864  :
865  mapping( typename mpl::at<list_v, mpl::int_<topo_dim> >::type().entity( topo_dim, entity_id ) )
866  {}
867 
868  node_type operator()( node_type const& x ) const
869  {
870  return mapping( x );
871  }
872  mapping_type mapping;
873  };
874 
875  void make_normals()
876  {
877  M_normals.resize( numNormals );
878 
879  for ( int n = 0; n < numNormals; ++n )
880  {
881  M_normals[n].resize( nDim );
882  M_normals[n].clear();
883  }
884 
885  if ( nDim == 1 )
886  {
887  M_normals[0][0] = -1;
888  M_normals[1][0] = 1;
889  }
890 
891  if ( nDim == 2 )
892  {
893  M_normals[0][1] = -1;
894  M_normals[1][0] = 1;
895  M_normals[2][1] = 1;
896  M_normals[3][0] = -1;
897  }
898 
899  if ( nDim == 3 )
900  {
901  M_normals[0][2] = -1;
902  M_normals[1][1] = -1;
903  M_normals[2][0] = 1;
904  M_normals[3][1] = 1;
905  M_normals[4][0] = -1;
906  M_normals[5][2] = 1;
907  }
908  }
909 
911  void computeBarycenters();
912 
914  void computeMeasure();
915 private:
916 
917  uint16_type M_id;
918 
919  points_type M_vertices;
920 
921  points_type M_points;
922 
923  normals_type M_normals;
924 
925 
926  node_type M_barycenter;
927 
928  points_type M_barycenterfaces;
929 
930  value_type M_meas;
931 };
932 
933 template<uint16_type Dim, uint16_type Order, uint16_type RDim, typename T>
934 const uint16_type Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T>::nbPtsPerVertex;
935 template<uint16_type Dim, uint16_type Order, uint16_type RDim, typename T>
936 const uint16_type Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T>::nbPtsPerEdge;
937 template<uint16_type Dim, uint16_type Order, uint16_type RDim, typename T>
938 const uint16_type Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T>::nbPtsPerFace;
939 template<uint16_type Dim, uint16_type Order, uint16_type RDim, typename T>
940 const uint16_type Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T>::numGeometricFaces;
941 
942 
943 
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> {};
946 
947 template<uint16_type Dim, uint16_type Order, uint16_type RDim, typename T>
948 void
949 Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T>::computeBarycenters()
950 {
951  M_barycenter = ublas::column( glas::average( M_vertices ), 0 );
952 
953  for ( int f = 0; f < numTopologicalFaces; ++f )
954  {
955  //std::cout << "face " << f << " vertices " << faceVertices( f ) << "\n";
956  ublas::column( M_barycenterfaces, f ) = ublas::column( glas::average( faceVertices( f ) ), 0 );
957  }
958 }
959 
960 template<uint16_type Dim, uint16_type Order, uint16_type RDim, typename T>
961 void
962 Reference<Hypercube<Dim, Order, RDim>, Dim, Order, RDim, T>::computeMeasure()
963 {
964  if ( nDim == nRealDim )
965  {
966 
967  typename matrix_node<value_type>::type M( nDim,nDim );
968 
969  //double factor = 1;
970  switch ( nDim )
971  {
972  case 1:
973  M_meas = 2;
974  break;
975 
976  case 2:
982  ublas::column( M, 0 ) = this->vertex( 2 )-this->vertex( 0 );
983  ublas::column( M, 1 ) = this->vertex( 3 )-this->vertex( 1 );
984  //factor = 2;
985  M_meas = 4;
986  break;
987 
988  case 3:
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 );
994  M_meas = 8;
995  break;
996  }
997 
998  //M_meas = math::abs( details::det( M, mpl::int_<nDim>() ) );
999  }
1000 
1001  else
1002  {
1003 #if 0
1004 
1010  typename matrix_node<value_type>::type M( nRealDim,nRealDim );
1011  value_type factor( 1 );
1012 
1013  switch ( nDim )
1014  {
1015  case 1:
1016  M_meas = ublas::norm2( this->vertex( 1 )-this->vertex( 0 ) );
1017  break;
1018 
1019  case 2:
1020  ublas::column( M, 0 ) = this->vertex( 0 )-this->vertex( 1 );
1021  ublas::column( M, 1 ) = this->vertex( 1 )-this->vertex( 2 );
1022  factor = 2;
1023  break;
1024  }
1025 
1026 #endif
1027  }
1028 }
1029 
1030 }
1031 #endif /* __refhypercube_H */

Generated on Fri Oct 25 2013 14:24:24 for Feel++ by doxygen 1.8.4