27 #include <boost/numeric/ublas/storage.hpp>
30 #include <feel/feelmesh/geo0d.hpp>
31 #include <feel/feelpoly/geomap.hpp>
32 #include <feel/feelmesh/marker.hpp>
39 template<
int Dim,
int Order,
int RealDim,
template<u
int16_type,u
int16_type,u
int16_type>
class Entity,
typename T>
struct GT_Lagrange;
48 template <
typename GeoShape>
52 uint16_type operate( uint16_type
const & point )
54 return ( point < GeoShape::numVertices ) ?
55 GeoShape::numVertices - point :
56 GeoShape::numPoints - point + GeoShape::numVertices;
69 template <uint16_type Dim,
72 typename POINTTYPE = Geo0D<Dim, T> >
84 typedef GEOSHAPE GeoShape;
85 typedef POINTTYPE PointType;
87 typedef PointType point_type;
88 typedef typename super::face_type face_type;
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;
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;
108 static const uint16_type nDim = super::nDim;
109 static const uint16_type nOrder = super::nOrder;
110 static const uint16_type nRealDim = super::nRealDim;
112 template<
int GmOrder>
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;
121 typedef typename GetGm<nOrder>::type gm_type;
122 typedef typename GetGm<nOrder>::ptrtype gm_ptrtype;
124 typedef typename GetGm<1>::type gm1_type;
125 typedef typename GetGm<1>::ptrtype gm1_ptrtype;
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>,
132 mpl::identity<vertex_permutation_type>,
133 typename mpl::if_<mpl::equal_to<mpl::int_<nDim>,
135 mpl::identity<edge_permutation_type>,
136 mpl::identity<face_permutation_type> >::type>::type::type permutation_type;
143 M_points( numPoints ),
145 M_G( nRealDim, numPoints ),
146 M_barycenter( nRealDim ),
147 M_barycenterfaces( nRealDim, numTopologicalFaces ),
149 M_h_face( numTopologicalFaces, 1 ),
150 M_h_edge( numLocalEdges, 1 ),
152 M_measurefaces( numTopologicalFaces ),
153 M_normals( nRealDim, numTopologicalFaces ),
154 M_has_points( false ),
173 M_points( numPoints ),
175 M_G( nRealDim, numPoints ),
176 M_barycenter( nRealDim ),
177 M_barycenterfaces( nRealDim, numTopologicalFaces ),
179 M_h_face( numTopologicalFaces, 1 ),
180 M_h_edge( numLocalEdges, 1 ),
182 M_measurefaces( numTopologicalFaces ),
183 M_normals( nRealDim, numTopologicalFaces ),
184 M_has_points( false ),
197 M_points( numPoints ),
199 M_G( nRealDim, numPoints ),
200 M_barycenter( e.M_barycenter ),
201 M_barycenterfaces( e.M_barycenterfaces ),
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 ),
210 M_marker1( e.M_marker1 ),
211 M_marker2( e.M_marker2 ),
212 M_marker3( e.M_marker3 ),
218 for ( uint16_type i = 0; i < numLocalPoints; ++i )
219 M_points[ i ] = e.M_points[ i ];
250 void setMesh(
MeshBase const* m )
const
256 gm_ptrtype
gm()
const
296 super::operator=( G );
298 for ( uint16_type i = 0; i < numLocalPoints; ++i )
299 M_points[ i ] = G.M_points[ i ];
304 M_barycenter = G.M_barycenter;
305 M_barycenterfaces = G.M_barycenterfaces;
307 M_h_face = G.M_h_face;
308 M_h_edge = G.M_h_edge;
310 M_has_points = G.M_has_points;
312 M_neighbors = G.M_neighbors;
314 M_marker1 = G.M_marker1;
315 M_marker2 = G.M_marker2;
316 M_marker3 = G.M_marker3;
351 std::pair<size_type,size_type>
const&
neighbor( uint16_type n )
const
353 return M_neighbors[n];
361 M_neighbors[n] = std::make_pair( neigh_id, proc_id );
364 bool isNeighbor( self_type
const& G )
const
366 for ( uint16_type i = 0; i< this->
nNeighbors() ; ++i )
367 if ( this->
neighbor( i ).first==G.id() )
return true;
385 return ublas::column( M_barycenterfaces, f );
393 return M_barycenterfaces;
401 return permutation_type();
410 return *(
static_cast<POINTTYPE *
>( M_points[ i ] ) );
420 return *(
static_cast<POINTTYPE *
>( M_points[ i ] ) );
428 return M_points[ i ];
437 return M_points[ i ];
442 PointType
const & facePoint ( uint16_type __f, uint16_type
const __i )
const
444 return M_face_points[__f][__i];
459 return *(
static_cast<POINTTYPE *
>( M_points[ detail::ReversePoint<GEOSHAPE>::operate( i ) ] ) );
473 return *(
static_cast<POINTTYPE *
>( M_points[ detail::ReversePoint<GEOSHAPE>::operate( i ) ] ) );
481 void setPoint( uint16_type
const i, point_type
const & p );
491 std::ostream &
showMe(
bool verbose =
false, std::ostream & c = std::cout )
const;
502 void swapPoints(
const uint16_type & pt1,
const uint16_type & pt2 );
523 matrix_node_type
const&
G()
const
539 return ublas::subrange( M_G, 0, nRealDim, 0, numVertices );
550 matrix_node_type &
G()
555 point_iterator beginPoint()
557 return M_points.begin();
559 point_const_iterator beginPoint()
const
561 return M_points.begin();
563 point_iterator endPoint()
565 return M_points.end();
567 point_const_iterator endPoint()
const
569 return M_points.end();
595 double hEdge( uint16_type f )
const
602 static uint16_type
fToP( uint16_type
const _localFace, uint16_type
const _point )
621 return M_measurefaces[f];
629 return M_measurefaces;
643 ublas::matrix_column<matrix_node_type const>
normal( uint16_type f )
const
645 return ublas::column( M_normals, f );
656 static uint16_type
fToP( uint16_type
const _localFace, uint16_type
const _point )
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 );
677 return super::nbOppositePointsPerFace;
687 return super::faceToOppositePoint( _localFace, _point );
698 ublas::matrix<T> orientation_matrix ( nRealDim,nRealDim );
700 for (
int i = 0; i < nRealDim ; ++i )
702 ublas::row( orientation_matrix, i ) = ( ublas::column( this->
G(), i+1 ) -
703 ublas::column( this->
G(), 0 ) );
710 return ( sgn > 0 ) ? 1 : 0;
713 void applyDisplacement(
int i, ublas::vector<double>
const& u )
715 ublas::column( M_G, i ) += u;
716 ( *M_points[ i ] ) += u;
718 void applyDisplacementG(
int i, ublas::vector<double>
const& u )
720 ublas::column( M_G, i ) += u;
730 M_marker1.assign( tags[0] );
732 if ( tags.size() > 1 )
733 M_marker2.assign( tags[1] );
735 if ( tags.size() > 2 )
742 std::vector<int> p( tags[2]-1 );
744 for (
size_type i = 0; i < p.size(); ++i )
754 Marker1
const& marker()
const
762 void setMarker( flag_type v )
764 return M_marker1.assign( v );
767 Marker2
const& marker2()
const
775 void setMarker2( flag_type v )
777 return M_marker2.assign( v );
780 Marker3
const& marker3()
const
788 void setMarker3( flag_type v )
790 return M_marker3.assign( v );
796 return M_pneighbors.size();
806 M_meas_pneighbors = meas;
811 return M_meas_pneighbors;
815 void updateWithPc(
typename gm_type::precompute_ptrtype
const& pc,
typename gm_type::faces_precompute_type & pcf );
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> );
825 friend class boost::serialization::access;
826 template<
class Archive>
827 void serialize( Archive & ar,
const unsigned int version )
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";
834 DVLOG(2) <<
" - G...\n";
836 DVLOG(2) <<
" - marker1...\n";
838 DVLOG(2) <<
" - marker1: " << M_marker1.value() <<
"...\n";
839 DVLOG(2) <<
" - marker2...\n";
841 DVLOG(2) <<
" - marker2: " << M_marker2.value() <<
"...\n";
842 DVLOG(2) <<
" - marker3...\n";
844 DVLOG(2) <<
" - marker3: " << M_marker3.value() <<
"...\n";
849 std::vector<point_type*> M_points;
852 std::vector<std::vector<point_type*> > M_face_points;
855 matrix_node_type M_G;
856 node_type M_barycenter;
857 matrix_node_type M_barycenterfaces;
860 std::vector<double> M_h_face;
861 std::vector<double> M_h_edge;
864 std::vector<double> M_measurefaces;
865 matrix_node_type M_normals;
873 std::vector<std::pair<size_type,size_type> > M_neighbors;
875 std::set<size_type> M_pneighbors;
877 value_type M_meas_pneighbors;
884 mutable MeshBase
const* M_mesh;
885 mutable gm_ptrtype M_gm;
886 mutable gm1_ptrtype M_gm1;
889 template <u
int16_type Dim,
typename GEOSHAPE,
typename T,
typename POINTTYPE>
890 const uint16_type GeoND<Dim,GEOSHAPE, T, POINTTYPE>::numLocalPoints;
892 template <u
int16_type Dim,
typename GEOSHAPE,
typename T,
typename POINTTYPE>
893 const uint16_type GeoND<Dim,GEOSHAPE, T, POINTTYPE>::numLocalVertices;
895 template <u
int16_type Dim,
typename GEOSHAPE,
typename T,
typename POINTTYPE>
900 M_points[ i ] =
const_cast<point_type *
>( &p );
902 FEELPP_ASSERT( const_cast<point_type *>( &p ) != 0 ).error(
"invalid Geo0D<>" );
903 ublas::column( M_G, i ) = M_points[i]->node();
908 template <u
int16_type Dim,
typename GEOSHAPE,
typename T,
typename POINTTYPE>
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";
919 for (
int i = 0; i < numVertices; i++ )
921 out <<
"POINT id = " << i << std::endl;
922 out << point( i ).showMe( verbose, out );
925 out <<
"----- END OF GeoND data ---" << std::endl << std::endl;
929 template <u
int16_type Dim,
typename GEOSHAPE,
typename T,
typename POINTTYPE>
932 point_type * tmp( M_points[ pt1 ] );
933 M_points[ pt1 ] = M_points[ pt2 ];
934 M_points[ pt2 ] = tmp;
937 ublas::column( M_G, pt1 ).swap( ublas::column( M_G, pt2 ) );
941 template <u
int16_type Dim,
typename GEOSHAPE,
typename T,
typename POINTTYPE>
944 point_type * tmp[ numPoints ];
946 for (
unsigned int i = 0; i < numPoints; ++i )
948 tmp[ i ] = M_points[ i ];
951 for (
unsigned int i = 0; i < numPoints; ++i )
953 M_points[ i ] = tmp[ otn[ i ] ];
954 ublas::column( M_G, i ) = M_points[i]->node();
958 template <u
int16_type Dim,
typename GEOSHAPE,
typename T,
typename POINTTYPE>
962 if ( !M_gm.use_count() )
963 M_gm = gm_ptrtype(
new gm_type );
965 if ( !M_gm1.use_count() )
966 M_gm1 = gm1_ptrtype(
new gm1_type );
968 auto pc = M_gm->preCompute( M_gm, M_gm->referenceConvex().vertices() );
969 auto pcf = M_gm->preComputeOnFaces( M_gm, M_gm->referenceConvex().barycenterFaces() );
971 updateWithPc( pc, pcf );
974 template <u
int16_type Dim,
typename GEOSHAPE,
typename T,
typename POINTTYPE>
976 GeoND<Dim,GEOSHAPE, T, POINTTYPE>::updateWithPc(
typename gm_type::precompute_ptrtype
const& pc,
977 typename gm_type::faces_precompute_type& pcf )
981 for ( uint16_type __e = 0; __e < numLocalEdges; ++__e )
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];
989 auto M = glas::average( M_G );
990 M_barycenter = ublas::column( M, 0 );
991 M_pneighbors.clear();
993 for ( uint16_type __p = 0; __p < numPoints; ++__p )
995 std::copy( M_points[__p]->
elements().begin(),
997 std::inserter( M_pneighbors, M_pneighbors.begin() ) );
1000 auto ctx = M_gm->template context<vm::JACOBIAN>( *
this, pc );
1002 double w = ( nDim == 3 )?4./3.:2;
1003 M_measure = w*ctx->J( 0 );
1005 updatep( pcf,
typename mpl::equal_to<mpl::int_<nDim>, mpl::int_<nRealDim> >::type() );
1007 template <u
int16_type Dim,
typename GEOSHAPE,
typename T,
typename POINTTYPE>
1009 GeoND<Dim,GEOSHAPE, T, POINTTYPE>::updatep(
typename gm_type::faces_precompute_type& pcf, mpl::bool_<true> )
1019 int nEdges = GEOSHAPE::topological_face_type::numEdges;
1021 for ( uint16_type __f = 0; __f < numTopologicalFaces; ++__f )
1025 for ( uint16_type e = 0; e < nEdges; ++e )
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 );
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 );
1044 M_h_face[__f] = ( M_h_face[__f] > __l )?M_h_face[__f]:__l;
1050 auto ctx = M_gm->template context<vm::POINT|vm::NORMAL|vm::KB|vm::JACOBIAN>(
1056 std::vector<double> f2( numTopologicalFaces, 2 );
1057 std::vector<double> f3( numTopologicalFaces, 2 );
1059 if ( GEOSHAPE::is_simplex )
1061 f2[0] = 2.82842712474619;
1062 f3[0] = 3.464101615137754;
1065 for (
int f = 0; f < numTopologicalFaces; ++f )
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 );
1072 ublas::column( M_barycenterfaces, f ) = ublas::column( glas::average( this->face( f ).G() ) );
1074 double w = ( nDim == 3 )?f3[f]:( ( nDim==2 )?f2[f]:1 );
1075 M_measurefaces[f] = w*ctx->J( 0 )*ctx->normalNorm( 0 );
1079 template <u
int16_type Dim,
typename GEOSHAPE,
typename T,
typename POINTTYPE>
1081 GeoND<Dim,GEOSHAPE, T, POINTTYPE>::updatep(
typename gm_type::faces_precompute_type& pcf, mpl::bool_<false> )
1085 template <u
int16_type Dim,
typename GEOSHAPE,
typename T,
typename POINTTYPE>
1088 operator<<( DebugStream& __os, GeoND<Dim,GEOSHAPE, T, POINTTYPE>
const& __n )
1090 if ( __os.doPrint() )
1092 std::ostringstream __str;
1094 __str << __n.showMe(
true, __str );
1096 __os << __str.str() <<
"\n";
1102 template <u
int16_type Dim,
typename GEOSHAPE,
typename T,
typename POINTTYPE>
1105 operator<<( NdebugStream& __os, GeoND<Dim,GEOSHAPE, T, POINTTYPE>
const& __n )
1111 template <u
int16_type Dim,
typename GEOSHAPE,
typename T,
typename POINTTYPE>
1114 operator<<( std::ostream& __os, GeoND<Dim,GEOSHAPE, T, POINTTYPE>
const& __n )
1116 return __n.showMe(
true, __os );
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
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
matrix_node_type const & G() const
Definition: geond.hpp:523
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