25 #include <boost/version.hpp>
26 #if (BOOST_VERSION >= 103400)
27 #include <boost/none.hpp>
29 #include <boost/none_t.hpp>
32 #include <boost/mpl/vector.hpp>
33 #include <boost/shared_ptr.hpp>
34 #include <boost/optional.hpp>
35 #include <boost/enable_shared_from_this.hpp>
42 #include <feel/feelcore/feel.hpp>
49 #include <feel/feelpoly/context.hpp>
53 #include <feel/feelmesh/marker.hpp>
62 enum class GeomapStrategyType
77 template<uint16_type O,
78 template<u
int16_type Dim>
class PolySetType,
79 typename ContinuityType,
80 template<
class, u
int16_type,
class>
class Pts,
81 uint16_type TheTag >
class Lagrange;
93 template<uint16_type Dim,
97 template<u
int16_type, u
int16_type, u
int16_type>
class Entity = Simplex,
98 template<u
int16_type,
template<u
int16_type RDim>
class PolySetType,
typename ContinuityType,
99 template<
class, u
int16_type,
class>
class Pts, uint16_type >
class PP = Lagrange>
102 public PP<Order,Scalar, Continuous,PointSetEquiSpaced, 0>::template apply<Dim,RealDim, T, Entity<Dim,Order,Dim> >::result_type,
103 public boost::enable_shared_from_this<GeoMap<Dim, Order, RealDim, T, Entity, PP > >
107 typedef typename PP<Order, Scalar, Continuous, PointSetEquiSpaced, 0>::template apply<Dim, RealDim, T, Entity<Dim,Order,Dim> >::result_type super;
109 static const uint16_type nRealDimCheck2d = mpl::if_< mpl::less_equal<mpl::int_<2>,mpl::int_<RealDim> >,
111 mpl::int_<Dim> >::type::value;
114 typedef mpl::vector<boost::none_t, boost::none_t, boost::none_t,
117 typedef mpl::vector<boost::none_t, boost::none_t,
119 GeoMap<2, Order, nRealDimCheck2d, T, Entity, PP> > geomap_faces_t;
121 typedef mpl::vector<GeoMap<1, Order, RealDim,T, Entity, PP>,
122 GeoMap<2, Order, nRealDimCheck2d,T, Entity, PP>,
124 boost::none_t> geomap_elements_t;
127 static const uint16_type nDim = super::nDim;
128 static const uint16_type nRealDim = super::nRealDim;
129 static const uint16_type nDof = super::nDof;
130 static const uint16_type nNodes = super::nNodes;
131 static const fem::transformation_type trans = super::trans;
133 typedef typename super::value_type value_type;
135 typedef typename super::PreCompute precompute_type;
136 typedef boost::shared_ptr<precompute_type> precompute_ptrtype;
138 typedef typename super::convex_type convex_type;
143 typedef boost::shared_ptr<geometric_mapping_type> geometric_mapping_ptrtype;
145 typedef typename mpl::at<geomap_elements_t, mpl::int_<nDim> >::type element_gm_type;
146 typedef boost::shared_ptr<element_gm_type> element_gm_ptrtype;
147 typedef element_gm_ptrtype gm_ptrtype;
149 typedef typename mpl::at<geomap_faces_t, mpl::int_<nDim> >::type face_gm_type;
150 typedef boost::shared_ptr<face_gm_type> face_gm_ptrtype;
155 typedef typename mpl::at<geomap_faces_t, mpl::int_<N> >::type type;
159 typedef typename mpl::at<geomap_edges_t, mpl::int_<nDim> >::type edge_gm_type;
160 typedef boost::shared_ptr<edge_gm_type> edge_gm_ptrtype;
165 typedef typename mpl::at<geomap_edges_t, mpl::int_<N> >::type type;
172 typedef node_t_type normal_type;
173 typedef ublas::vector<normal_type> normals_type;
174 typedef typename normals_type::const_iterator normal_const_iterator;
175 typedef node_t_type tangent_type;
176 typedef ublas::vector<tangent_type> tangents_type;
177 typedef typename tangents_type::const_iterator tangent_const_iterator;
179 typedef typename ublas::vector<value_type> vector_type;
180 typedef typename ublas::matrix<value_type> matrix_type;
188 M_is_cached( false ),
191 M_g_linear( nNodes, nDim ),
194 if ( trans == fem::LINEAR )
197 matrix_node_t_type __dummy_pts( ublas::zero_matrix<value_type>( nDim, 1 ) );
199 ublas::vector<ublas::matrix<value_type> > m = super::derivate( __dummy_pts );
201 FEELPP_ASSERT( M_g_linear.size2() == m.size() )( M_g_linear.size2() )( m.size() ).error(
"invalid dimension" );
202 FEELPP_ASSERT( m( 0 ).size2() == 1 )( m( 0 ).size2() ).error(
"Invalid number of points" );
204 FEELPP_ASSERT( M_g_linear.size1() == m( 0 ).size1() )( M_g_linear.size1() )( m( 0 ).size1() ).error(
"invalid number of DOF" );
211 for ( uint16_type i = 0; i < nNodes; ++i )
213 for ( uint16_type n = 0; n < nDim; ++n )
216 M_g_linear( i, n ) = m( n )( i, 0 );
222 for ( uint16_type i = 0; i < nNodes; ++i )
224 for ( uint16_type n = 0; n < nDim; ++n )
226 M_g_linear( i, n ) = this->dPhi( i, n, __dummy_pt );
235 GeoMap( element_gm_ptrtype
const& e, face_gm_ptrtype
const& f )
238 M_is_cached( false ),
241 M_g_linear( nNodes, nDim ),
244 if ( trans == fem::LINEAR )
247 node_t_type __dummy_pt( nDim );
249 matrix_node_t_type __dummy_pts( ublas::zero_matrix<value_type>( nDim, 1 ) );
253 ublas::vector<ublas::matrix<value_type> > m = super::derivate( __dummy_pts );
259 for ( uint16_type i = 0; i < nNodes; ++i )
261 for ( uint16_type n = 0; n < nDim; ++n )
264 M_g_linear( i, n ) = m( n )( i, 0 );
270 for ( uint16_type i = 0; i < nNodes; ++i )
272 for ( uint16_type n = 0; n < nDim; ++n )
274 M_g_linear( i, n ) = this->dPhi( i, n, __dummy_pt );
296 uint16_type realDim()
const
306 return trans == fem::LINEAR;
330 ublas::vector<value_type>
refNode( uint16_type i )
const
332 return ublas::column( this->
points(), i );
349 return M_refconvex.isIn( pt );
358 matrix_node_t_type
const& __G )
const
360 namespace lambda = boost::lambda;
365 for ( uint16_type __i = 0; __i < nNodes; ++__i )
368 value_type __phi_at_pt = super::evaluate( __i, __ref_p )( 0 );
369 __real_p.plus_assign( __phi_at_pt*ublas::column( __G, __i ) );
380 matrix_node_t_type
const& __G,
381 precompute_type
const* __pc )
const
383 node_t_type __real_p( __G.size1() );
386 for ( uint16_type __i = 0; __i < nNodes; ++__i )
389 value_type __phi_at_pt = __pc->phi( __idref, __i );
390 __real_p.plus_assign( __phi_at_pt*ublas::column( __G, __i ) );
400 precompute_type
const* pc,
401 matrix_type & x )
const
404 blas::gemm( traits::NO_TRANSPOSE, traits::NO_TRANSPOSE,
408 ublas::axpy_prod( G, pc->phi(), x, true );
422 matrix_type& __g )
const
424 namespace lambda = boost::lambda;
426 if ( trans == fem::LINEAR )
433 FEELPP_ASSERT( __pt.size() ==
dim() )( __pt.size() )(
dim() ).error(
"invalid dimension" );
435 matrix_node_t_type __pts( nDim, 1 );
436 ublas::column( __pts, 0 ) = __pt;
438 ublas::vector<ublas::matrix<value_type> > m = super::derivate( __pts );
440 for ( uint16_type n = 0; n < nDim; ++n )
442 ublas::column( __g, n ) = ublas::column( m( n ), 0 );
457 precompute_type
const* __pc )
const
459 FEELPP_ASSERT( __pc )( __idref ).error(
"a PreCompute must be set first before using this function" );
461 if ( trans == fem::LINEAR )
468 FEELPP_ASSERT( __pc->dim() ==
dim() )( __pc->dim() )(
dim() ).error(
"invalid dimension" );
472 for ( uint16_type n = 0; n < nDim; ++n )
474 __g( i, n ) = __pc->grad( i, 0, n, __idref );
494 matrix_node_t_type __g( Nref,
dim() );
499 FEELPP_ASSERT( __g.size1() == G.size2() )( __g.size1() )( G.size2() ).error(
"invalid sizes" );
501 matrix_node_t_type K( G.size1(), __g.size2() );
502 ublas::axpy_prod( G, __g, K );
508 return __max*sqrt( value_type( N ) ) / value_type( N );
513 bool isCached()
const
518 std::vector<bool> M_cached;
519 std::vector<double> M_J;
522 std::vector<matrix_type> M_K;
523 std::vector<matrix_type> M_B;
525 template<
typename MeshType>
526 void initCache( MeshType
const* mesh )
529 LOG(INFO) <<
"[Geomap] start caching J,K,B for " << nelts <<
" elements\n";
530 M_cached.resize( nelts );
531 std::fill( M_cached.begin(), M_cached.end(), false );
538 bool isCacheValid()
const
540 if ( !( M_cached.size() > 0 &&
545 LOG(INFO) <<
"invalid cache\n";
551 bool cached(
int e )
const
553 FEELPP_ASSERT( this->isCacheValid() )( e ).error(
"invalid cache" );
556 void setCached(
int e,
bool v )
558 FEELPP_ASSERT( this->isCacheValid() )( e ).error(
"invalid cache" );
561 double J(
int e )
const
563 FEELPP_ASSERT( this->isCacheValid() )( e ).error(
"invalid cache" );
564 DCHECK( e >= 0 && e < M_J.size() ) <<
"invalid element id " << e <<
"( " << M_J.size() <<
") for geomap cache, are you using the proper mesh\n";
567 matrix_type
const& B(
int e )
const
569 FEELPP_ASSERT( this->isCacheValid() )( e ).error(
"invalid cache" );
570 DCHECK( e >= 0 && e < M_B.size() ) <<
"invalid element id " << e <<
"( " << M_B.size() <<
") for geomap cache, are you using the proper mesh\n";
573 matrix_type
const& K(
int e )
const
575 FEELPP_ASSERT( this->isCacheValid() )( e ).error(
"invalid cache" );
576 DCHECK( e >= 0 && e < M_K.size() ) <<
"invalid element id " << e <<
"( " << M_K.size() <<
") for geomap cache, are you using the proper mesh\n";
579 void addJ(
int e,
double v )
581 FEELPP_ASSERT( this->isCacheValid() )( e ).error(
"invalid cache" );
582 DCHECK( e >= 0 && e < M_J.size() ) <<
"invalid element id " << e <<
"( " << M_J.size() <<
") for geomap cache, are you using the proper mesh\n";
585 void addK(
int e, matrix_type
const& K )
587 FEELPP_ASSERT( this->isCacheValid() )( e ).error(
"invalid cache" );
588 DCHECK( e >= 0 && e < M_K.size() ) <<
"invalid element id " << e <<
"( " << M_K.size() <<
") for geomap cache, are you using the proper mesh\n";
589 M_K[e].resize( K.size1(), K.size2() );
591 for (
size_type i = 0; i < K.size1(); ++i )
592 for (
size_type j = 0; j < K.size2(); ++j )
593 M_K[e]( i, j ) = K( i, j );
595 void addB(
int e, matrix_type
const& B )
597 FEELPP_ASSERT( this->isCacheValid() )( e ).error(
"invalid cache" );
598 DCHECK( e >= 0 && e < M_B.size() ) <<
"invalid element id " << e <<
"( " << M_B.size() <<
") for geomap cache, are you using the proper mesh\n";
599 M_B[e].resize( B.size1(), B.size2() );
601 for (
size_type i = 0; i < B.size1(); ++i )
602 for (
size_type j = 0; j < B.size2(); ++j )
603 M_B[e]( i, j ) = B( i, j );
612 template<
size_type context_v,
typename ElementType>
616 static const size_type contextv = context_v;
617 static const size_type context = context_v;
619 static const uint16_type PDim = ElementType::nDim;
621 static const uint16_type NDim = ElementType::nRealDim;
622 static const uint16_type nDim = NDim;
624 static const fem::transformation_type trans = geometric_mapping_type::trans;
625 static const bool is_linear = ( trans == fem::LINEAR );
627 static const bool condition = ( ( PDim==NDim )||( ( NDim>=1 )&&( PDim==NDim-1 ) ) );
629 typedef typename mpl::if_<mpl::equal_to<mpl::int_<PDim>, mpl::int_<NDim> >,
630 mpl::identity<GeoMap<Dim, Order, NDim, T, Entity, PP > >,
631 typename mpl::if_<mpl::and_<mpl::greater_equal<mpl::int_<NDim>, mpl::int_<1> >,
632 mpl::equal_to<mpl::int_<PDim>, mpl::int_<NDim-1> > >,
635 mpl::identity<typename GeoMap<NDim, Order, NDim, T, Entity, PP >::face_gm_type>,
636 typename mpl::if_<mpl::and_<mpl::equal_to<mpl::int_<NDim>, mpl::int_<3> >,
637 mpl::equal_to<mpl::int_<PDim>, mpl::int_<NDim-2> > >,
640 mpl::identity<typename GeoMap<NDim, Order, NDim, T, Entity, PP >::edge_gm_type>,
641 mpl::identity<boost::none_t> >::type>::type>::type::type gm_type;
642 typedef boost::shared_ptr<gm_type> gm_ptrtype;
644 typedef typename gm_type::value_type value_type;
646 typedef typename gm_type::precompute_ptrtype precompute_ptrtype;
649 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
651 typedef node_t_type normal_type;
652 typedef ublas::vector<normal_type> normals_type;
653 typedef typename normals_type::const_iterator normal_const_iterator;
654 typedef node_t_type tangent_type;
655 typedef ublas::vector<tangent_type> tangents_type;
656 typedef typename tangents_type::const_iterator tangent_const_iterator;
658 typedef ElementType element_type;
659 typedef typename element_type::permutation_type permutation_type;
662 element_type
const& __e,
663 precompute_ptrtype
const& __pc )
666 M_element( boost::addressof( __e ) ),
669 M_npoints( M_pc->nPoints() ),
675 M_G( ( gm_type::nNodes == element_type::numVertices ) ?__e.vertices() : __e.G() ),
676 M_n( M_gm->referenceConvex().normals() ),
681 M_xrefq( PDim, nPoints() ),
682 M_xrealq( NDim, nPoints() ),
683 M_nrealq( NDim, nPoints() ),
684 M_unrealq( NDim, nPoints() ),
685 M_nnormq( nPoints() ),
687 M_g( M_G.size2(), PDim ),
692 M_B3( boost::extents[NDim][NDim][PDim][PDim] ),
694 M_e_marker( __e.marker() ),
695 M_e_marker2( __e.marker2() ),
696 M_e_marker3( __e.marker3() ),
704 M_meas( __e.measure() ),
713 M_gm->gradient( node_t_type(), M_g_linear );
718 M_Jt.resize( nPoints() );
719 M_Bt.resize( nPoints() );
726 element_type
const& __e,
727 std::vector<std::map<permutation_type, precompute_ptrtype> > & __pc,
731 M_element( boost::addressof( __e ) ),
734 M_npoints( __pc[__f][__e.permutation( __f )]->nPoints() ),
740 M_G( ( gm_type::nNodes == element_type::numVertices ) ?__e.vertices() : __e.G() ),
741 M_n( M_gm->referenceConvex().normals() ),
746 M_xrefq( PDim, nPoints() ),
747 M_xrealq( NDim, nPoints() ),
748 M_nrealq( NDim, nPoints() ),
749 M_unrealq( NDim, nPoints() ),
750 M_nnormq( nPoints() ),
752 M_g( M_G.size2(), PDim ),
757 M_B3( boost::extents[NDim][NDim][PDim][PDim] ),
759 M_e_marker( __e.marker() ),
760 M_e_marker2( __e.marker2() ),
761 M_e_marker3( __e.marker3() ),
769 M_meas( __e.measure() ),
770 M_measface( __e.faceMeasure( __f ) ),
778 M_gm->gradient( node_t_type(), M_g_linear );
783 M_Jt.resize( nPoints() );
784 M_Bt.resize( nPoints() );
792 M_element( p->M_element ),
794 M_pc_faces( p->M_pc_faces ),
795 M_npoints( M_pc->nPoints() ),
800 M_G( ( gm_type::nNodes == element_type::numVertices ) ?M_element->vertices() : M_element->G() ),
802 M_n_real( p->M_n_real ),
803 M_u_n_real( p->M_u_n_real ),
804 M_n_norm( p->M_n_norm ),
805 M_t_real( p->M_t_real ),
806 M_xrefq( p->M_xrefq ),
807 M_xrealq( p->M_xrealq ),
808 M_nrealq( p->M_nrealq ),
809 M_unrealq( p->M_unrealq ),
810 M_nnormq( p->M_nnormq ),
818 M_e_marker( p->M_e_marker ),
819 M_e_marker2( p->M_e_marker2 ),
820 M_e_marker3( p->M_e_marker3 ),
827 M_h_face( p->M_h_face ),
829 M_measface( p->M_measface ),
836 M_gm->gradient( node_t_type(), M_g_linear );
841 M_Jt.resize( nPoints() );
842 M_Bt.resize( nPoints() );
853 return gmc_ptrtype(
new gmc_type( *
this ) );
868 void update( element_type
const& __e, uint16_type __f )
871 M_element = boost::addressof( __e );
874 M_perm = __e.permutation( M_face_id );
876 M_h_face = __e.hFace( M_face_id );
879 M_pc = M_pc_faces[__f][M_perm];
881 M_G = ( gm_type::nNodes == element_type::numVertices ) ?__e.vertices() : __e.G();
883 M_e_marker = __e.marker();
884 M_e_marker2 = __e.marker2();
885 M_e_marker3 = __e.marker3();
887 M_meas = __e.measure();
888 M_measface = __e.faceMeasure( __f );
889 M_xrefq = M_pc->nodes();
891 FEELPP_ASSERT( M_G.size2() == M_gm->nbPoints() )( M_G.size2() )( M_gm->nbPoints() ).error(
"invalid dimensions" );
892 FEELPP_ASSERT( M_pc ).error(
"invalid precompute data structure" );
894 if ( vm::has_point<context>::value )
898 std::fill( M_xrealq.data().begin(), M_xrealq.data().end(), value_type( 0 ) );
899 const uint16_type size1 = M_G.size1();
900 const uint16_type size3 = M_G.size2();
901 const uint16_type size2 = M_pc->nPoints();
903 for ( uint16_type i = 0; i < size1; ++i )
904 for ( uint16_type j = 0; j < size2; ++j )
906 for ( uint16_type k = 0; k < size3; ++k )
907 M_xrealq( i, j ) += M_G( i, k ) * M_pc->phi()[k][j]( 0,0 );
911 if ( vm::has_jacobian<context>::value )
913 updateJKBN( mpl::bool_<is_linear>() );
918 void update( element_type
const& __e, uint16_type __f, permutation_type __perm )
921 M_element = boost::addressof( __e );
926 M_h_face = __e.hFace( M_face_id );
929 M_pc = M_pc_faces[__f][M_perm];
931 M_G = ( gm_type::nNodes == element_type::numVertices ) ?__e.vertices() : __e.G();
933 M_e_marker = __e.marker();
934 M_e_marker2 = __e.marker2();
935 M_e_marker3 = __e.marker3();
937 M_meas = __e.measure();
938 M_measface = __e.faceMeasure( __f );
939 M_xrefq = M_pc->nodes();
941 FEELPP_ASSERT( M_G.size2() == M_gm->nbPoints() )( M_G.size2() )( M_gm->nbPoints() ).error(
"invalid dimensions" );
942 FEELPP_ASSERT( M_pc ).error(
"invalid precompute data structure" );
944 if ( vm::has_point<context>::value )
948 std::fill( M_xrealq.data().begin(), M_xrealq.data().end(), value_type( 0 ) );
949 const uint16_type size1 = M_G.size1();
950 const uint16_type size3 = M_G.size2();
951 const uint16_type size2 = M_pc->nPoints();
953 for ( uint16_type i = 0; i < size1; ++i )
954 for ( uint16_type j = 0; j < size2; ++j )
956 for ( uint16_type k = 0; k < size3; ++k )
957 M_xrealq( i, j ) += M_G( i, k ) * M_pc->phi()[k][j]( 0,0 );
961 if ( vm::has_jacobian<context>::value )
963 updateJKBN( mpl::bool_<is_linear>() );
968 void update( element_type
const& __e,
969 precompute_ptrtype
const& __pc )
974 if ( M_npoints != M_pc->nPoints() )
976 M_npoints = M_pc.get()->nPoints();
978 M_xrefq.resize( PDim, nPoints() );
979 M_xrealq.resize( NDim, nPoints() );
980 M_nrealq.resize( NDim, nPoints() );
981 M_unrealq.resize( NDim, nPoints() );
982 M_nnormq.resize( nPoints() );
986 M_gm->gradient( node_t_type(), M_g_linear );
991 M_Jt.resize( nPoints() );
992 M_Bt.resize( nPoints() );
1013 M_G = ( gm_type::nNodes == element_type::numVertices ) ?__e.vertices() : __e.G();
1015 M_g.resize( M_G.size2(), PDim );
1017 M_element = boost::addressof( __e );
1019 M_e_marker = __e.marker();
1020 M_e_marker2 = __e.marker2();
1021 M_e_marker3 = __e.marker3();
1024 M_meas = __e.measure();
1025 M_xrefq = M_pc->nodes();
1027 FEELPP_ASSERT( M_G.size2() == M_gm->nbPoints() )( M_G.size2() )( M_gm->nbPoints() ).error(
"invalid dimensions" );
1028 FEELPP_ASSERT( M_pc ).error(
"invalid precompute data structure" );
1030 if ( vm::has_point<context>::value )
1032 std::fill( M_xrealq.data().begin(), M_xrealq.data().end(), value_type( 0 ) );
1033 const uint16_type size1 = M_G.size1();
1034 const uint16_type size3 = M_G.size2();
1035 const uint16_type size2 = M_pc->nPoints();
1037 for ( uint16_type i = 0; i < size1; ++i )
1038 for ( uint16_type j = 0; j < size2; ++j )
1040 for ( uint16_type k = 0; k < size3; ++k )
1041 M_xrealq( i, j ) += M_G( i, k ) * M_pc->phi()[k][j]( 0,0 );
1045 if ( vm::has_jacobian<context>::value )
1047 updateJKBN( mpl::bool_<is_linear>() );
1071 uint16_type
N()
const
1079 uint16_type
P()
const
1092 element_type
const& element_c()
const
1100 uint16_type nPoints()
const
1116 ublas::matrix_column<matrix_node_t_type const>
xRef(
int q )
const
1118 return ublas::column( M_xrefq, q );
1134 ublas::matrix_column<matrix_type const>
xReal(
int q )
const
1137 return ublas::column( M_xrealq, q );
1147 value_type
J(
int q )
const
1158 value_type
J(
int q )
const
1171 value_type
J(
int , mpl::bool_<true> )
const
1179 value_type
J(
int q, mpl::bool_<false> )
const
1189 matrix_type
const&
G()
const
1200 matrix_type
const&
B(
int i )
const
1203 return B( i, mpl::bool_<is_linear>() );
1205 value_type
B(
int c1,
int c2,
int q )
const
1208 return B( q, mpl::bool_<is_linear>() )( c1, c2 );
1211 matrix_type
const& K(
int i )
const
1215 value_type K(
int c1,
int c2,
int q )
const
1217 return M_K( c1, c2 );
1222 matrix_type
const&
B(
int , mpl::bool_<true> )
const
1229 matrix_type
const&
B(
int i, mpl::bool_<false> )
const
1241 boost::multi_array<value_type,4>
const&
B3()
const
1251 node_t_type __barycenter( M_gm->dim() );
1252 __barycenter.clear();
1254 for ( uint16_type __c = 0; __c < M_gm->refNodes().size(); ++__c )
1256 __barycenter += M_gm->refNode( __c );
1259 __barycenter /= M_gm->refNodes().size();
1260 return __barycenter;
1268 node_t_type __barycenter( M_G.size1() );
1269 __barycenter.clear();
1272 for ( uint16_type __c = 0; __c < M_G.size2(); ++__c )
1274 __barycenter += ublas::column( M_G, __c );
1277 __barycenter /= M_G.size2();
1278 return __barycenter;
1288 if ( trans == fem::LINEAR )
1291 return std::abs( ublas::norm_2(
xReal()-ublas::column( M_G, 0 )-
1292 ublas::prod( M_K,
xRef() ) ) ) < 1e-10;
1298 node_t_type
const& refNormal(
int )
const
1300 return M_gm->referenceConvex().normal( M_face_id );
1333 ublas::matrix_column<matrix_node_t_type const>
normal(
int q )
const
1341 return ublas::column( M_nrealq, q );
1365 return ublas::column( M_unrealq, q );
1369 value_type
const&
unitNormal(
int n,
int q )
const
1374 return M_u_n_real( n );
1377 return M_unrealq( n, q );
1380 node_t_type
const& tangent()
const
1387 node_t_type tangent(
int q )
const
1398 value_type
const& unitTangent(
int n,
int q )
const
1403 return M_t_real( n );
1406 return M_utrealq( n, q );
1424 uint16_type faceId()
const
1484 return M_pos_in_elem_id_1;
1504 return M_pos_in_elem_id_2;
1514 return M_gm->radiusEstimate( M_G );
1523 value_type
h()
const
1541 value_type meas()
const
1549 value_type measurePointElementNeighbors()
const
1551 return M_element->measurePointElementNeighbors();
1558 value_type measFace()
const
1575 precompute_ptrtype
const&
pc()
const
1583 std::vector<std::map<permutation_type, precompute_ptrtype> >
const &
pcFaces()
const
1596 void setPc( precompute_ptrtype
const& __pc )
1601 void setPcFaces( std::vector<std::map<permutation_type, precompute_ptrtype> >
const& __pcfaces )
1603 M_pc_faces = __pcfaces;
1612 permutation_type
permutation( mpl::bool_<false> )
const
1620 permutation_type
permutation( mpl::bool_<true> )
const
1624 M_perm != permutation_type( permutation_type::NO_PERMUTATION ) ) )
1625 ( M_face_id ).error(
"invalid permutation" );
1632 void updateJKBN( mpl::bool_<true> )
1635 if ( !M_gm->isCached() ||
1636 ( M_gm->isCached() && M_gm->cached( M_id ) == false ) )
1640 if ( boost::is_arithmetic<value_type>::value )
1641 atlas::gemm( traits::NO_TRANSPOSE, traits::NO_TRANSPOSE,
1642 1.0, M_G, M_g_linear,
1647 ublas::axpy_prod( M_G, M_g_linear, M_K,
true );
1651 M_J = math::abs( det<NDim>( M_K ) );
1655 inverse<NDim>( M_K, M_CS, M_J );
1656 ublas::noalias( M_B ) = ublas::trans( M_CS );
1658 inverse<NDim>( M_K, M_CS );
1659 ublas::noalias( M_B ) = ublas::trans( M_CS );
1668 if ( boost::is_arithmetic<value_type>::value )
1669 atlas::gemm( traits::TRANSPOSE, traits::NO_TRANSPOSE,
1676 ublas::noalias( M_CS ) = ublas::prod( ublas::trans( M_K ), M_K );
1678 M_J = math::sqrt( math::abs( det<PDim>( M_CS ) ) );
1681 inverse<PDim>( M_CS, M_CSi );
1685 if ( boost::is_arithmetic<value_type>::value )
1686 atlas::gemm( traits::NO_TRANSPOSE, traits::NO_TRANSPOSE,
1692 ublas::axpy_prod( M_K, M_CSi, M_B,
true );
1697 if ( M_gm->isCached() )
1700 M_gm->addJ( M_id, M_J );
1703 M_gm->addK( M_id, M_K );
1704 M_gm->addB( M_id, M_B );
1706 M_gm->setCached( M_id,
true );
1715 M_J = M_gm->J( M_id );
1719 M_K = M_gm->K( M_id );
1720 M_B = M_gm->B( M_id );
1727 if ( vm::has_hessian<context>::value )
1730 for ( uint16_type k = 0; k < NDim; ++k )
1731 for ( uint16_type l = 0; l < NDim; ++l )
1732 for ( uint16_type i = 0; i < PDim; ++i )
1733 for ( uint16_type j = 0; j < PDim; ++j )
1734 M_B3[k][l][i][j] = M_B( k,i )*M_B( l,j );
1738 std::cout <<
"element " << this->
id() <<
" B3 = " <<
"\n";
1740 for (
int i = 0; i < nDim; ++i )
1741 for (
int j = 0; j < nDim; ++j )
1742 for (
typename boost::multi_array<value_type,4>::index k = 0; k < nDim; ++k )
1743 for (
typename boost::multi_array<value_type,4>::index l = 0; l < nDim; ++l )
1744 std::cout <<
"B3[" << i <<
"][" << j <<
"][" << k <<
"][" << l <<
"]="
1745 << this->
B3()[i][j][k][l] <<
"\n";
1753 blas::gemv( traits::NO_TRANSPOSE,
1754 1.0, M_Bt[ q ], M_n[M_face_id],
1757 ublas::axpy_prod( M_B,
1758 M_gm->referenceConvex().normal( M_face_id ),
1762 M_n_norm = ublas::norm_2( M_n_real );
1763 M_u_n_real = M_n_real/M_n_norm;
1770 ublas::axpy_prod( M_K,
1771 M_gm->referenceConvex().tangent( M_face_id ),
1774 double ratio = M_gm->referenceConvex().h( M_face_id )/M_h_face;
1784 void updateJKBN( mpl::bool_<false> )
1789 for (
int q = 0; q < nPoints(); ++q )
1792 M_gm->gradient( q, M_g, M_pc.get() );
1796 blas::gemm( M_G, M_g, M_K );
1798 ublas::axpy_prod( M_G, M_g, M_K,
true );
1803 M_J = math::abs( det<NDim>( M_K ) );
1805 if ( vm::has_kb<context>::value )
1807 inverse<NDim>( M_K, M_CS );
1809 M_B = ublas::trans( M_CS );
1818 ublas::prod( ublas::trans( M_K ), M_K, M_CS );
1820 blas::gemm( traits::TRANSPOSE, traits::NO_TRANSPOSE,
1824 M_J = math::sqrt( math::abs( det<PDim>( M_CS ) ) );
1826 if ( vm::has_kb<context>::value )
1828 inverse<PDim>( M_CS, M_CSi );
1831 ublas::axpy_prod( M_K, M_CSi, M_B );
1833 blas::gemm( traits::NO_TRANSPOSE, traits::NO_TRANSPOSE,
1847 if ( vm::has_kb<context>::value )
1850 M_Bt[q].resize( M_B.size1(), M_B.size2() );
1861 for (
int q = 0; q < nPoints(); ++q )
1864 blas::gemv( traits::NO_TRANSPOSE,
1865 1.0, M_Bt[ q ], M_n[M_face_id],
1871 ublas::axpy_prod( M_B,
1872 M_gm->referenceConvex().normal( M_face_id ),
1878 ublas::axpy_prod( M_Bt[ q ],
1879 M_gm->referenceConvex().normal( M_face_id ),
1885 M_n_norm = ublas::norm_2( M_n_real );
1886 M_u_n_real = M_n_real/M_n_norm;
1887 ublas::column( M_nrealq, q ) = M_n_real;
1888 ublas::column( M_unrealq, q ) = M_u_n_real;
1889 M_nnormq[q] = M_n_norm;
1892 M_Jt[q] *= M_n_norm;
1894 if ( vm::has_tangent<context>::value )
1901 std::cout <<
"[geomap] face id = " << M_face_id <<
"\n"
1902 <<
"[geomap] ref normal = " << M_gm->referenceConvex().normal( M_face_id ) <<
"\n"
1903 <<
"[geomap] M_n_real = " << M_nrealq <<
"\n"
1904 <<
"[geomap] M_unrealq = " << M_unrealq <<
"\n"
1905 <<
"[geomap] M_nnormq = " << M_nnormq <<
"\n";
1916 element_type
const* M_element;
1920 precompute_ptrtype M_pc;
1921 std::vector<std::map<permutation_type, precompute_ptrtype> > M_pc_faces;
1922 uint16_type M_npoints;
1928 ublas::vector<node_t_type> M_n;
1929 node_t_type M_n_real;
1930 node_t_type M_u_n_real;
1931 value_type M_n_norm;
1932 node_t_type M_t_real;
1935 matrix_node_t_type M_xrefq;
1937 matrix_type M_xrealq;
1938 matrix_node_t_type M_nrealq;
1939 matrix_node_t_type M_unrealq;
1940 ublas::vector<value_type> M_nnormq;
1942 matrix_type M_g_linear;
1948 boost::multi_array<value_type,4> M_B3;
1952 Marker2 M_e_marker2;
1953 Marker3 M_e_marker3;
1955 uint16_type M_pos_in_elem_id_1;
1957 uint16_type M_pos_in_elem_id_2;
1959 uint16_type M_face_id;
1962 value_type M_h_face;
1964 value_type M_measface;
1967 std::vector<matrix_type> M_Bt;
1969 permutation_type M_perm;
1973 template<
size_type context_v,
typename ElementType>
1974 boost::shared_ptr<Context<context_v,ElementType> >
1975 context( geometric_mapping_ptrtype gm, ElementType
const& e, precompute_ptrtype
const& pc )
1977 return boost::shared_ptr<Context<context_v,ElementType> >(
1978 new Context<context_v, ElementType>( gm,
1983 template<
size_type context_v,
typename ElementType>
1984 boost::shared_ptr<Context<context_v,ElementType> >
1985 context( ElementType
const& e, precompute_ptrtype
const& pc )
1987 return boost::shared_ptr<Context<context_v,ElementType> >(
1988 new Context<context_v, ElementType>( this->shared_from_this(),
1993 template<
size_type context_v,
typename ElementType>
1994 boost::shared_ptr<Context<context_v,ElementType> >
1995 context( geometric_mapping_ptrtype gm,
1996 ElementType
const& e,
1997 std::vector<std::map<typename ElementType::permutation_type,precompute_ptrtype> > & pc,
2000 return boost::shared_ptr<Context<context_v,ElementType> >(
2001 new Context<context_v, ElementType>( gm,
2006 template<
size_type context_v,
typename ElementType>
2007 boost::shared_ptr<Context<context_v,ElementType> >
2008 context( ElementType
const& e,
2009 std::vector<std::map<typename ElementType::permutation_type,precompute_ptrtype> > & pc,
2012 return boost::shared_ptr<Context<context_v,ElementType> >(
2013 new Context<context_v, ElementType>( this->shared_from_this(),
2028 typedef boost::shared_ptr<geometric_mapping_type> geometric_mapping_ptrtype;
2030 static const fem::transformation_type trans = geometric_mapping_type::trans;
2032 typedef typename geometric_mapping_type::node_t_type node_t_type;
2033 typedef typename geometric_mapping_type::node_t_type node_type;
2034 typedef typename geometric_mapping_type::matrix_node_t_type matrix_node_t_type;
2035 typedef typename geometric_mapping_type::matrix_node_t_type matrix_node_type;
2037 template<
typename GeoElem>
2038 Inverse( geometric_mapping_ptrtype __gm, GeoElem
const& __ge,
2039 WorldComm
const& worldComm = Environment::worldComm().subWorldCommSeq() )
2042 M_xref( __gm->dim() ),
2043 M_xreal( __ge.G().size1() ),
2046 M_K(
N(), __gm->dim() ),
2047 M_B(
N(), __gm->dim() ),
2048 M_CS( __gm->dim(), __gm->dim() ),
2049 M_g( M_gm->nbPoints(), __gm->dim() ),
2050 #
if defined( FEELPP_HAS_PETSC )
2056 FEELPP_ASSERT( M_G.size2() == __gm->nbPoints() )
2057 ( M_G.size2() )( __gm->nbPoints() ).error(
"invalid dimensions" );
2059 if ( M_gm->isLinear() )
2065 #if defined( FEELPP_HAS_PETSC )
2066 M_nlsolver->dense_residual = boost::bind( &Inverse::updateResidual, boost::ref( *
this ), _1, _2 );
2067 M_nlsolver->dense_jacobian = boost::bind( &Inverse::updateJacobian, boost::ref( *
this ), _1, _2 );
2075 template<
typename GeoElem>
2076 Inverse( geometric_mapping_ptrtype __gm, GeoElem
const& __ge, mpl::int_<1> ,
2077 WorldComm
const& worldComm = Environment::worldComm().subWorldCommSeq())
2080 M_xref( __gm->dim() ),
2081 M_xreal( __ge.vertices().size1() ),
2083 M_G( __ge.vertices() ),
2084 M_K(
N(), __gm->dim() ),
2085 M_B(
N(), __gm->dim() ),
2086 M_CS( __gm->dim(), __gm->dim() ),
2087 M_g( M_gm->nbPoints(), __gm->dim() ),
2088 #
if defined( FEELPP_HAS_PETSC )
2094 FEELPP_ASSERT( M_G.size2() == __gm->nbPoints() )
2095 ( M_G.size2() )( __gm->nbPoints() ).error(
"invalid dimensions" );
2097 if ( M_gm->isLinear() )
2103 #if defined( FEELPP_HAS_PETSC )
2104 M_nlsolver->dense_residual = boost::bind( &Inverse::updateResidual, boost::ref( *
this ), _1, _2 );
2105 M_nlsolver->dense_jacobian = boost::bind( &Inverse::updateJacobian, boost::ref( *
this ), _1, _2 );
2130 uint16_type
N()
const
2134 uint16_type P()
const
2158 matrix_type
const&
G()
const
2166 matrix_type
const&
K()
const
2171 value_type J()
const
2173 return math::abs( det<Dim>( M_K ) );
2180 matrix_type
const&
B()
const
2190 node_type __barycenter( M_gm->dim() );
2191 __barycenter.clear();
2193 for ( uint16_type __c = 0; __c < M_gm->referenceConvex().nPoints(); ++__c )
2195 __barycenter += M_gm->refNode( __c );
2198 __barycenter /= M_gm->referenceConvex().nPoints();
2199 return __barycenter;
2223 if ( M_gm->isLinear() )
2226 return std::abs( ublas::norm_2(
xReal()-ublas::column( M_G, 0 )-
2227 ublas::prod( M_K,
xRef() ) ) ) < 1e-10;
2246 if ( M_gm->isLinear() )
2249 M_is_in = linearInverse();
2254 M_is_in = nonLinearInversePetsc();
2260 matrix_node_t_type operator() ( matrix_node_t_type
const& real_pts,
bool allow_extrapolation =
false )
const
2262 if ( trans == fem::LINEAR )
2263 return linearInversePoints( real_pts, allow_extrapolation );
2266 return matrix_node_t_type();
2273 class GeomapInverseConvex
2277 GeomapInverseConvex( Inverse& __gmi,
2278 const node_t_type& __xr )
2284 scalar_type operator()(
const node_t_type& x )
const
2286 node_type r = gmi.geometricMapping()->transform( x, gmi.G() ) - xreal;
2287 return ublas::inner_prod( r, r )/2.;
2289 void operator()(
const node_t_type& x, node_t_type& gr )
const
2292 gmi.M_xref.assign( x );
2294 node_type r = gmi.geometricMapping()->transform( x, gmi.G() ) - xreal;
2295 gr.resize( x.size() );
2296 ublas::prod( ublas::trans( gmi.K() ), r, gr );
2304 typedef ublas::vector<double> dense_vector_type;
2305 typedef ublas::matrix<double> dense_matrix_type;
2307 void updateResidual( dense_vector_type
const& x, dense_vector_type& r )
2309 dense_vector_type y = M_gm->transform( x, M_G );
2316 M_gm->gradient( x, M_g );
2317 ublas::axpy_prod( M_G, M_g, M_K );
2318 ublas::prod( ublas::trans( M_K ), y-M_xreal, r );
2322 LOG(INFO) <<
"[geomap::residual] begin ------------------------------\n";
2323 LOG(INFO) <<
"[geomap::residual] x =" << x <<
"\n";
2324 LOG(INFO) <<
"[geomap::residual] M_G =" << M_G <<
"\n";
2326 LOG(INFO) <<
"[geomap::residual] y =" << y <<
"\n";
2327 LOG(INFO) <<
"[geomap::residual] xreal =" << M_xreal <<
"\n";
2328 LOG(INFO) <<
"[geomap::residual] r(xreal-y) =" << r <<
"\n";
2329 LOG(INFO) <<
"[geomap::residual] end ------------------------------\n";
2334 void updateJacobian( dense_vector_type
const& x, dense_matrix_type& j )
2336 M_gm->gradient( x, M_g );
2339 ublas::axpy_prod( M_G, M_g, j );
2343 ublas::axpy_prod( M_G, M_g, M_K );
2344 ublas::prod( ublas::trans( M_K ), M_K, j );
2348 LOG(INFO) <<
"[geomap::jacobian] begin ------------------------------\n";
2349 LOG(INFO) <<
"[geomap::jacobian] x =" << x <<
"\n";
2350 LOG(INFO) <<
"[geomap::jacobian] j =" << j <<
"\n";
2351 LOG(INFO) <<
"[geomap::jacobian] end ------------------------------\n";
2372 M_gm->gradient(
xRef(), M_g );
2373 DVLOG(2) <<
"[update] g = " << M_g <<
"\n";
2377 ublas::axpy_prod( M_G, M_g, M_K );
2378 DVLOG(2) <<
"[update] K(0) = " << M_K <<
"\n";
2381 if ( M_gm->dim() !=
N() )
2383 ublas::prod( ublas::trans( M_K ), M_K, M_CS );
2384 LU<matrix_type> __lu( M_CS );
2385 __lu.inverse( M_CS );
2386 ublas::axpy_prod( M_K, M_CS, M_B );
2391 LU<matrix_type> __lu( M_K );
2392 __lu.inverse( M_CS );
2393 M_B = ublas::trans( M_CS );
2396 DVLOG(2) <<
"[update] B(0) = " << M_B <<
"\n";
2400 checkInvariant()
const
2402 FEELPP_ASSERT( M_G.size2() == M_g.size1() )( M_G.size2() )( M_g.size1() ).error(
"G,g invalid dimensions" );
2403 FEELPP_ASSERT( M_G.size1() == M_K.size1() )( M_G.size1() )( M_K.size1() ).error(
"G,K invalid dimensions" );
2404 FEELPP_ASSERT( M_g.size2() == M_K.size2() )( M_g.size2() )( M_K.size2() ).error(
"g,K invalid dimensions" );
2405 FEELPP_ASSERT( M_B.size2() == M_gm->dim() )( M_B.size1() )(
N() ).error(
"B,gm invalid dimensions" );
2409 bool linearInverse()
2416 node_type y( M_xreal );
2418 DVLOG(2) <<
"y = xreal = " << y <<
"\n";
2420 y.minus_assign( ublas::column( M_G, 0 ) );
2421 DVLOG(2) <<
"y - G(0) = " << y <<
"\n";
2423 DVLOG(2) <<
"B(0) = " << M_B <<
"\n";
2424 DVLOG(2) <<
"xref = " << ublas::prod( ublas::trans( M_B ), y ) <<
"\n";
2427 M_xref.assign( ublas::prod( ublas::trans( M_B ), y )-ublas::scalar_vector<value_type>( P, 1.0 ) );
2429 DVLOG(2) <<
"[GeoMap::Inverse::linearInverse] xref : " << M_xref <<
"\n";
2433 boost::tie( __isin, vmin ) = M_gm->isIn( M_xref );
2434 DVLOG(2) <<
"[GeoMap::Inverse::linearInverse] isIn : " << __isin <<
"\n";
2445 ublas::axpy_prod( M_K, -M_xref, y );
2447 if ( ublas::norm_2( y ) < 1e-10 )
2455 matrix_node_t_type linearInversePoints( matrix_node_t_type
const& real_pts,
bool =
false )
const
2457 return ublas::prod( ublas::trans( M_B ),
2458 real_pts-ublas::outer_prod( ublas::column( M_G, 0 ),
2459 ublas::scalar_vector<value_type>( real_pts.size2(), value_type( 1 ) ) ) ) -
2460 ublas::scalar_matrix<value_type>( M_B.size2(), real_pts.size2(), value_type( 1 ) );
2468 bool nonLinearInversePetsc()
2472 const double IN_EPS = 1e-10;
2480 node_type x0 = M_gm->refNode( 0 );
2481 node_type y = ublas::column( M_G, 0 );
2482 scalar_type d = ublas::inner_prod( y, M_xreal );
2484 for (
size_type j = 1; j < M_gm->nbPoints(); ++j )
2486 scalar_type d2 = ublas::inner_prod( ublas::column( M_G, j ), M_xreal );
2491 x0 = M_gm->refNode( j );
2492 y = ublas::column( M_G, j );
2503 dense_matrix_type J( P, P );
2504 dense_vector_type R( P );
2506 updateResidual( M_xref, R );
2507 updateJacobian( M_xref, J );
2510 M_nlsolver->setType( TRUST_REGION );
2511 M_nlsolver->setRelativeResidualTol( 1e-16 );
2512 M_nlsolver->solve( J, M_xref, R, 1e-10, 10 );
2518 this->updateResidual( M_xref, R );
2519 boost::tie( __isin, vmin ) = M_gm->isIn( M_xref );
2522 ( P == N || ublas::norm_2( R ) < IN_EPS ) )
2541 bool nonLinearInverse()
2543 const double EPS = 1e-10;
2544 const double IN_EPS = 1e-10;
2551 node_type x0 = M_gm->refNode( 0 );
2552 node_type y = ublas::column( M_G, 0 );
2553 scalar_type d = ublas::inner_prod( y, M_xreal );
2555 for (
size_type j = 1; j < M_gm->nbPoints(); ++j )
2557 scalar_type d2 = ublas::inner_prod( ublas::column( M_G, j ), M_xreal );
2562 x0 = M_gm->refNode( j );
2563 y = ublas::column( M_G, j );
2569 node_type vres( N );
2570 node_type rn( M_xreal );
2571 rn.minus_assign( y );
2576 ublas::axpy_prod( rn, M_K, vres );
2577 scalar_type res = ublas::norm_2( vres );
2582 while ( res > EPS/10 && cnt > 0 )
2587 ublas::axpy_prod( rn, M_B, xn );
2591 for ( int16_type i=1; i<=256; i*=2 )
2593 M_xref.plus_assign( xn / scalar_type( i ) );
2594 y = M_gm->transform( M_xref,
G() );
2603 ublas::axpy_prod( rn, M_K, vres );
2604 newres = ublas::norm_2( vres );
2609 newres = ublas::norm_2( rn );
2612 if ( newres < 1.5*res )
2628 std::cerr <<
"BFGS in geotrans_inv_convex!\n";
2629 GeomapInverseConvex b( *
this, M_xreal );
2631 iteration_ptrtype iter( Iteration<double>::New() );
2632 iter->setMaximumNumberOfIterations( 50 );
2633 iter->setRelativePrecision( 1e-8 );
2636 bfgs( b,b,x,10,*iter );
2637 rn = M_gm->transform( x,
G() ) - M_xreal;
2639 if ( M_gm->isIn( x ) < IN_EPS &&
2640 N==P && ublas::norm_2( rn ) > IN_EPS )
2641 throw "inversion of non-linear geometric transformation "
2642 "failed ! (too much iterations)";
2649 boost::tie( __isin, vmin ) = M_gm->isIn( M_xref );
2652 ( P == N || ublas::norm_2( rn ) < IN_EPS ) )
2670 geometric_mapping_ptrtype M_gm;
2683 boost::shared_ptr<SolverNonLinear<double> > M_nlsolver;
2688 element_gm_ptrtype _elementMap;
2689 face_gm_ptrtype _boundaryMap;
2690 matrix_type M_g_linear;
2692 friend class Inverse;
2694 reference_convex_type M_refconvex;
2698 #include <boost/preprocessor/comparison/less.hpp>
2699 #include <boost/preprocessor/logical/and.hpp>
2700 #include <boost/preprocessor/control/if.hpp>
2701 #include <boost/preprocessor/list/at.hpp>
2702 #include <boost/preprocessor/list/cat.hpp>
2703 #include <boost/preprocessor/list/for_each_product.hpp>
2704 #include <boost/preprocessor/logical/or.hpp>
2705 #include <boost/preprocessor/tuple/to_list.hpp>
2706 #include <boost/preprocessor/tuple/eat.hpp>
2707 #include <boost/preprocessor/facilities/empty.hpp>
2708 #include <boost/preprocessor/punctuation/comma.hpp>
2709 #include <boost/preprocessor/facilities/identity.hpp>
2711 template<
int Dim,
int Order,
int RealDim,
template<u
int16_type,u
int16_type,u
int16_type>
class Entity = Simplex,
typename T =
double>
2715 template<
int Dim,
int Order,
int RealDim,
typename T =
double>
2720 # define FEELPP_GEOMAP \
2721 BOOST_PP_TUPLE_TO_LIST( \
2729 # define FEELPP_DIMS \
2730 BOOST_PP_TUPLE_TO_LIST( \
2738 # define FEELPP_REALDIMS \
2739 BOOST_PP_TUPLE_TO_LIST( \
2748 # define FEELPP_NEWDIMS \
2749 BOOST_PP_TUPLE_TO_LIST( \
2752 (1,1),(1,2),(1,3), \
2759 # define FEELPP_ORDERS \
2760 BOOST_PP_TUPLE_TO_LIST( \
2768 # define FEELPP_ENTITY BOOST_PP_TUPLE_TO_LIST( 2, ( Simplex,Hypercube ) )
2771 # define FEELPP_GEN_GT(GEOM,LDIM,LORDER) \
2772 "GT_" #GEOM "(" #LDIM "," #LORDER ")" \
2776 # define FEELPP_GT_FACTORY_OP(_, GDO) \
2777 FEELPP_GT_FACTORY GDO \
2781 # define FEELPP_GT_DIM(T) BOOST_PP_TUPLE_ELEM(2, 0, T) \
2783 # define FEELPP_GT_REALDIM(T) BOOST_PP_TUPLE_ELEM(2, 1, T) \
2787 #define FEELPP_GT_FACTORY(GEOM,LDIMS,LORDER,ENTITY) \
2788 template<typename T> \
2789 struct BOOST_PP_CAT(GT_, GEOM)<FEELPP_GT_DIM(LDIMS), LORDER, FEELPP_GT_REALDIM(LDIMS), ENTITY, T> \
2791 public GeoMap<FEELPP_GT_DIM(LDIMS), LORDER, FEELPP_GT_REALDIM(LDIMS), T, ENTITY, GEOM> \
2793 static const uint16_type nDim = FEELPP_GT_DIM(LDIMS); \
2794 static const uint16_type order = LORDER; \
2795 static const uint16_type nRealDim = FEELPP_GT_REALDIM(LDIMS); \
2796 static const uint16_type nRealDimCheck2d = mpl::if_< mpl::less_equal<mpl::int_<2>,mpl::int_<nRealDim> >, \
2797 mpl::int_<nRealDim>, \
2798 mpl::int_<nDim> >::type::value; \
2800 typedef mpl::vector<boost::none_t, \
2802 GeoMap<1, LORDER, nRealDim, T, ENTITY, GEOM>, \
2803 GeoMap<2, LORDER, nRealDimCheck2d, T, ENTITY, GEOM> > geomap_faces_t; \
2804 typedef mpl::vector<GeoMap<1, LORDER, nRealDim, T, ENTITY, GEOM>, \
2805 GeoMap<2, LORDER, nRealDimCheck2d, T, ENTITY, GEOM>, \
2806 GeoMap<3, LORDER, 3, T, ENTITY, GEOM>, \
2807 boost::none_t> geomap_elements_t; \
2808 typedef typename type_traits<T>::value_type value_type; \
2809 typedef GeoMap<nDim, LORDER, nRealDim, T, ENTITY, GEOM> super; \
2810 typedef BOOST_PP_CAT(GT_,GEOM)<nDim-1, LORDER, nRealDim,ENTITY, T> face_geo_type; \
2812 static const uint16_type nDof = super::nDof; \
2813 static const uint16_type nNodes = super::nNodes; \
2814 typedef typename mpl::at<geomap_elements_t, mpl::int_<nDim> >::type element_gm_type; \
2815 typedef typename mpl::at<geomap_faces_t, mpl::int_<nDim> >::type face_gm_type; \
2819 typedef typename mpl::at<geomap_faces_t, mpl::int_<N> >::type type; \
2822 BOOST_PP_CAT(GT_,GEOM)() \
2824 super( boost::shared_ptr<element_gm_type>(new element_gm_type()), boost::shared_ptr<face_gm_type>(new face_gm_type() )) \
2829 #define FEELPP_GT_FACTORY(GEOM,LDIM,LORDER,LREALDIM,ENTITY) \
2830 FEELPP_GT_FACTORY_BIS(GEOM,LDIM,LORDER,LREALDIM,ENTITY)), \
2836 BOOST_PP_LIST_FOR_EACH_PRODUCT( FEELPP_GT_FACTORY_OP, 4, ( FEELPP_GEOMAP, FEELPP_NEWDIMS, FEELPP_ORDERS, FEELPP_ENTITY ) )
2839 #undef FEELPP_ORDERS
2840 #undef FEELPP_REALDIMS
2841 #undef FEELPP_NEWDIMS
2844 template<
typename Elem,
template<u
int16_type,u
int16_type,u
int16_type>
class Entity = Simplex,
typename T =
double>
2845 class RealToReference
2848 static const uint16_type nDim = Elem::nDim;
2849 static const uint16_type nRealDim = Elem::nRealDim;
2851 typedef T value_type;
2852 typedef typename matrix_node<T>::type points_type;
2853 typedef GT_Lagrange<nDim,1, nRealDim, Entity, value_type> gm_type;
2854 typedef boost::shared_ptr<gm_type> gm_ptrtype;
2855 typedef typename gm_type::Inverse inverse_gm_type;
2857 RealToReference( Elem
const& elem )
2859 M_gm( new gm_type ),
2863 points_type operator()( points_type
const& pts )
const
2865 return M_igm( pts );
2868 value_type J()
const
2876 inverse_gm_type M_igm;
const face_gm_ptrtype & boundaryMap() const
Definition: geomap.hpp:321
element_type const & element() const
Definition: geomap.hpp:1087
permutation_type permutation() const
Definition: geomap.hpp:1567
bool isOnConvexSurface() const
Definition: geomap.hpp:1286
bool isIn() const
Definition: geomap.hpp:2212
gm_ptrtype const & geometricMapping() const
Definition: geomap.hpp:1063
ublas::vector< value_type > refNode(uint16_type i) const
Definition: geomap.hpp:330
const uint16_type invalid_uint16_type_value
Definition: feelcore/feel.hpp:338
node_t_type const & xRef() const
Definition: geomap.hpp:2142
Non linear solver base interface.
Definition: solvernonlinear.hpp:62
Structure for the geometrical mapping.
Definition: geomap.hpp:100
size_type id() const
Definition: geomap.hpp:1416
matrix_type const & B() const
Definition: geomap.hpp:2180
size_type id2() const
Definition: geomap.hpp:1492
value_type radiusEstimate(matrix_node_t_type const &G) const
Definition: geomap.hpp:487
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
value_type normalNorm(int q) const
Definition: geomap.hpp:1308
uint16_type idIn2() const
Definition: geomap.hpp:1502
Marker1 marker() const
Definition: geomap.hpp:1442
Marker2 marker2() const
Definition: geomap.hpp:1452
matrix_type const & G() const
Definition: geomap.hpp:2158
node_t_type barycenterRef() const
Definition: geomap.hpp:1249
bool isOnConvexSurface() const
Definition: geomap.hpp:2221
value_type h() const
Definition: geomap.hpp:1523
const size_type invalid_size_type_value
Definition: feelcore/feel.hpp:360
void gradient(const node_t_type &__pt, matrix_type &__g) const
Definition: geomap.hpp:421
uint16_type N() const
Definition: geomap.hpp:1071
uint16_type dim() const
Definition: geomap.hpp:291
value_type radiusEstimate() const
Definition: geomap.hpp:1512
value_type hFace() const
Definition: geomap.hpp:1533
node_t_type barycenterReal() const
Definition: geomap.hpp:1266
ublas::matrix_column< matrix_type const > xReal(int q) const
Definition: geomap.hpp:1134
matrix_node_t_type const & xRefs() const
Definition: geomap.hpp:1108
void update(element_type const &__e, uint16_type __f)
Definition: geomap.hpp:868
node_type barycenterReal() const
gmc_ptrtype clone()
Definition: geomap.hpp:851
matrix_type const & B(int i) const
Definition: geomap.hpp:1200
reference_convex_type const & referenceConvex() const
Definition: geomap.hpp:338
Marker2 marker3() const
Definition: geomap.hpp:1462
precompute_ptrtype const & pc() const
Definition: geomap.hpp:1575
Reference convex.
Definition: refentity.hpp:58
GeoMap()
Definition: geomap.hpp:185
Definition: geomap.hpp:613
size_type id1() const
Definition: geomap.hpp:1472
boost::tuple< bool, value_type > isIn(typename node< value_type >::type const &pt) const
Definition: geomap.hpp:347
uint16_type P() const
Definition: geomap.hpp:1079
node_t_type const & normal() const
Definition: geomap.hpp:1322
Singular Value Decomposition of a rectangular matrix.
Definition: svd.hpp:73
const element_gm_ptrtype elementMap() const
Definition: geomap.hpp:313
node_t_type transform(uint16_type __idref, matrix_node_t_type const &__G, precompute_type const *__pc) const
Definition: geomap.hpp:379
~GeoMap()
Definition: geomap.hpp:285
boost::multi_array< value_type, 4 > const & B3() const
Definition: geomap.hpp:1241
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
node_t_type const & xReal() const
Definition: geomap.hpp:2150
value_type J(int q) const
Definition: geomap.hpp:1158
matrix_type const & K() const
Definition: geomap.hpp:2166
void transform(matrix_node_t_type const &G, precompute_type const *pc, matrix_type &x) const
Definition: geomap.hpp:399
ublas::matrix_column< matrix_node_t_type const > xRef(int q) const
Definition: geomap.hpp:1116
node_type barycenterRef() const
Definition: geomap.hpp:2188
bool isLinear() const
Definition: geomap.hpp:304
geometric_mapping_ptrtype const & geometricMapping() const
Definition: geomap.hpp:2122
node_t_type const & unitNormal() const
Definition: geomap.hpp:1350
matrix_type const & G() const
Definition: geomap.hpp:1189
value_type conditionNumber() const
Definition: svd.hpp:132
Definition: geomap.hpp:2023
bool elementIsAFace() const
Definition: geomap.hpp:1432
void setXReal(node_type const &__xreal)
Definition: geomap.hpp:2242
matrix_type const & xReal() const
Definition: geomap.hpp:1125
uint16_type N() const
Definition: geomap.hpp:2130
void gradient(uint16_type __idref, matrix_type &__g, precompute_type const *__pc) const
Definition: geomap.hpp:455
GeoMap(element_gm_ptrtype const &e, face_gm_ptrtype const &f)
Definition: geomap.hpp:235
node_t_type transform(const node_t_type &__ref_p, matrix_node_t_type const &__G) const
Definition: geomap.hpp:357
uint16_type idIn1() const
Definition: geomap.hpp:1482
void update(element_type const &__e)
Definition: geomap.hpp:1011
ublas::matrix_column< matrix_node_t_type const > normal(int q) const
Definition: geomap.hpp:1333
void setPc(precompute_ptrtype const &__pc)
Definition: geomap.hpp:1596
std::vector< std::map< permutation_type, precompute_ptrtype > > const & pcFaces() const
Definition: geomap.hpp:1583