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
geomap.hpp
1 /*
2  This file is part of the Feel library
3 
4  Copyright (C) 2001,2002,2003,2004 EPFL, INRIA and Politechnico di Milano
5  Copyright (C) 2005,2006 EPFL
6  Copyright (C) 2006-2010 Université Joseph Fourier Grenoble 1
7 
8  This library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU Lesser General Public
10  License as published by the Free Software Foundation; either
11  version 3.0 of the License, or (at your option) any later version.
12 
13  This library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  Lesser General Public License for more details.
17 
18  You should have received a copy of the GNU Lesser General Public
19  License along with this library; if not, write to the Free Software
20  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21 */
22 #ifndef _GEOMAP_H
23 #define _GEOMAP_H
24 
25 #include <boost/version.hpp>
26 #if (BOOST_VERSION >= 103400)
27 #include <boost/none.hpp>
28 #else
29 #include <boost/none_t.hpp>
30 #endif /* BOOST_VERSION >= 103400 */
31 
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>
36 
37 //#include <boost/numeric/bindings/traits/traits.hpp>
38 //#include <boost/numeric/bindings/traits/ublas_vector.hpp>
39 //#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
40 //#include <boost/numeric/bindings/atlas/cblas3.hpp>
41 
42 #include <feel/feelcore/feel.hpp>
43 
44 #include <feel/feelalg/svd.hpp>
46 #include <feel/feelalg/lu.hpp>
48 
49 #include <feel/feelpoly/context.hpp>
52 #include <feel/feelpoly/fekete.hpp>
53 #include <feel/feelmesh/marker.hpp>
54 
55 
56 
57 namespace Feel
58 {
62 enum class GeomapStrategyType
63 {
64  GEOMAP_OPT = 0,
65  GEOMAP_O1 = 1,
66  GEOMAP_HO = 2
67 };
68 
69 //namespace blas = boost::numeric::bindings::blas;
70 //namespace traits = boost::numeric::bindings::traits;
71 
72 struct GeomapInverse
73 {
74  //GeomapInverse
75 };
76 
77 template<uint16_type O,
78  template<uint16_type Dim> class PolySetType,
79  typename ContinuityType,
80  template<class, uint16_type, class> class Pts,
81  uint16_type TheTag > class Lagrange;
82 
83 
93 template<uint16_type Dim,
94  uint16_type Order,
95  uint16_type RealDim,
96  typename T = double,
97  template<uint16_type, uint16_type, uint16_type> class Entity = Simplex,
98  template<uint16_type, template<uint16_type RDim> class PolySetType, typename ContinuityType,
99  template<class, uint16_type, class> class Pts, uint16_type > class PP = Lagrange>
100 class GeoMap
101  :
102  public PP<Order,Scalar, Continuous,PointSetEquiSpaced, 0>::template apply<Dim,RealDim/*Dim*/, T, Entity<Dim,Order,/*RealDim*/Dim> >::result_type,
103  public boost::enable_shared_from_this<GeoMap<Dim, Order, RealDim, T, Entity, PP > >
104  //public PP<Order,Scalar, PointSetFekete>::template apply<Dim, T, Entity<Dim,Order,Dim> >::result_type
105  {
106  //typedef typename PP<Order, Scalar, PointSetFekete>::template apply<Dim, T, Entity<Dim,Order,Dim> >::result_type super;
107  typedef typename PP<Order, Scalar, Continuous, PointSetEquiSpaced, 0>::template apply<Dim, RealDim/*Dim*/, T, Entity<Dim,Order,/*RealDim*/Dim> >::result_type super;
108 
109  static const uint16_type nRealDimCheck2d = mpl::if_< mpl::less_equal<mpl::int_<2>,mpl::int_<RealDim> >,
110  mpl::int_<RealDim>,
111  mpl::int_<Dim> >::type::value;
112 
113 
114  typedef mpl::vector<boost::none_t, boost::none_t, boost::none_t,
116 
117  typedef mpl::vector<boost::none_t, boost::none_t,
119  GeoMap<2, Order, nRealDimCheck2d/*RealDim*/, T, Entity, PP> > geomap_faces_t;
120 
121  typedef mpl::vector<GeoMap<1, Order, RealDim,T, Entity, PP>,
122  GeoMap<2, Order, nRealDimCheck2d/*RealDim*/,T, Entity, PP>,
124  boost::none_t> geomap_elements_t;
125  public:
126 
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;
132 
133  typedef typename super::value_type value_type;
134 
135  typedef typename super::PreCompute precompute_type;
136  typedef boost::shared_ptr<precompute_type> precompute_ptrtype;
137 
138  typedef typename super::convex_type convex_type;
139  typedef Reference<convex_type,nDim,Order,nDim/*nRealDim*/> reference_convex_type;
140 
143  typedef boost::shared_ptr<geometric_mapping_type> geometric_mapping_ptrtype;
144 
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;
148 
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;
151 
152  template<int N>
153  struct face_gm
154  {
155  typedef typename mpl::at<geomap_faces_t, mpl::int_<N> >::type type;
156  };
157 
158 
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;
161 
162  template<int N>
163  struct edge_gm
164  {
165  typedef typename mpl::at<geomap_edges_t, mpl::int_<N> >::type type;
166  };
167 
168 
169  typedef typename node<value_type>::type node_t_type;
170  typedef typename matrix_node<value_type>::type matrix_node_t_type;
171 
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;
178 
179  typedef typename ublas::vector<value_type> vector_type;
180  typedef typename ublas::matrix<value_type> matrix_type;
181 
182 
183 
186  :
187  super(),
188  M_is_cached( false ),
189  _elementMap( ),
190  _boundaryMap(),
191  M_g_linear( nNodes, nDim ),
192  M_refconvex()
193 {
194  if ( trans == fem::LINEAR )
195  {
196  //M_g_linear.resize( nNodes, nDim );
197  matrix_node_t_type __dummy_pts( ublas::zero_matrix<value_type>( nDim, 1 ) );
198 
199  ublas::vector<ublas::matrix<value_type> > m = super::derivate( __dummy_pts );
200 
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" );
203 
204  FEELPP_ASSERT( M_g_linear.size1() == m( 0 ).size1() )( M_g_linear.size1() )( m( 0 ).size1() ).error( "invalid number of DOF" );
205 
206  //std::cout << "nNodes= " << nNodes << "\n"
207  //<< "nDim= " << nDim << "\n";
208  //std::cout << "M_g_linear = " << M_g_linear << "\n"
209  //<< "m(0) = " << m( 0 ) << "\n";
210 
211  for ( uint16_type i = 0; i < nNodes; ++i )
212  {
213  for ( uint16_type n = 0; n < nDim; ++n )
214  {
215  //std::cout << "m(n)= " << m( n ) << "\n";
216  M_g_linear( i, n ) = m( n )( i, 0 );
217  }
218  }
219 
220 #if 0
221 
222  for ( uint16_type i = 0; i < nNodes; ++i )
223  {
224  for ( uint16_type n = 0; n < nDim; ++n )
225  {
226  M_g_linear( i, n ) = this->dPhi( i, n, __dummy_pt );
227  }
228  }
229 
230 #endif // 0
231  }
232 
233 }
235 GeoMap( element_gm_ptrtype const& e, face_gm_ptrtype const& f )
236  :
237  super(),
238  M_is_cached( false ),
239  _elementMap( e ),
240  _boundaryMap( f ),
241  M_g_linear( nNodes, nDim ),
242  M_refconvex()
243 {
244  if ( trans == fem::LINEAR )
245  {
246  //M_g_linear.resize( nNodes, nDim );
247  node_t_type __dummy_pt( nDim );
248 
249  matrix_node_t_type __dummy_pts( ublas::zero_matrix<value_type>( nDim, 1 ) );
250 
251  //std::cout << "geomap::derivate<> pts=" << __dummy_pts << "\n";
252  //std::cout << "geomap::derivate<> m=" << super::derivate( __dummy_pts ) << "\n";
253  ublas::vector<ublas::matrix<value_type> > m = super::derivate( __dummy_pts );
254  //std::cout << "nNodes= " << nNodes << "\n"
255  //<< "nDim= " << nDim << "\n";
256  //std::cout << "M_g_linear = " << M_g_linear << "\n"
257  //<< "m(0) = " << m( 0 ) << "\n";
258 
259  for ( uint16_type i = 0; i < nNodes; ++i )
260  {
261  for ( uint16_type n = 0; n < nDim; ++n )
262  {
263  //std::cout << "m(n)= " << m( n ) << "\n";
264  M_g_linear( i, n ) = m( n )( i, 0 );
265  }
266  }
267 
268 #if 0
269 
270  for ( uint16_type i = 0; i < nNodes; ++i )
271  {
272  for ( uint16_type n = 0; n < nDim; ++n )
273  {
274  M_g_linear( i, n ) = this->dPhi( i, n, __dummy_pt );
275  }
276  }
277 
278 #endif // 0
279  }
280 
281 }
286 {}
287 
291 uint16_type dim() const
292 {
293  return nDim;
294 }
295 
296 uint16_type realDim() const
297 {
298  return nRealDim;
299 }
300 
304 bool isLinear() const
305 {
306  return trans == fem::LINEAR;
307 }
308 
309 
313 const element_gm_ptrtype elementMap() const
314 {
315  return _elementMap;
316 }
317 
321 const face_gm_ptrtype& boundaryMap() const
322 {
323  return _boundaryMap;
324 }
325 
330 ublas::vector<value_type> refNode( uint16_type i ) const
331 {
332  return ublas::column( this->points(), i );
333 }
334 
339 {
340  return M_refconvex;
341 }
342 
347 boost::tuple<bool,value_type> isIn( typename node<value_type>::type const& pt ) const
348 {
349  return M_refconvex.isIn( pt );
350 }
351 
352 
357 node_t_type transform( const node_t_type &__ref_p,
358  matrix_node_t_type const& __G ) const
359 {
360  namespace lambda = boost::lambda;
361 
362  typename node<value_type>::type __real_p( __G.size1() );
363  __real_p.clear();
364 
365  for ( uint16_type __i = 0; __i < nNodes; ++__i )
366  {
367  //value_type __phi_at_pt = super::phi( __i, __ref_p );
368  value_type __phi_at_pt = super::evaluate( __i, __ref_p )( 0 );
369  __real_p.plus_assign( __phi_at_pt*ublas::column( __G, __i ) );
370  }
371 
372  return __real_p;
373 }
374 
379 node_t_type transform( uint16_type __idref,
380  matrix_node_t_type const& __G,
381  precompute_type const* __pc ) const
382 {
383  node_t_type __real_p( __G.size1() );
384  __real_p.clear();
385 
386  for ( uint16_type __i = 0; __i < nNodes; ++__i )
387  {
388  // evaluate transformation at point pt
389  value_type __phi_at_pt = __pc->phi( __idref, __i );
390  __real_p.plus_assign( __phi_at_pt*ublas::column( __G, __i ) );
391  }
392 
393  return __real_p;
394 }
395 
399 void transform( matrix_node_t_type const& G,
400  precompute_type const* pc,
401  matrix_type & x ) const
402 {
403 #if 0
404  blas::gemm( traits::NO_TRANSPOSE, traits::NO_TRANSPOSE,
405  1.0, G, pc->phi(),
406  0.0, x );
407 #else
408  ublas::axpy_prod( G, pc->phi(), x, true );
409 #endif
410 }
411 
412 
421 void gradient( const node_t_type& __pt,
422  matrix_type& __g ) const
423 {
424  namespace lambda = boost::lambda;
425 
426  if ( trans == fem::LINEAR )
427  {
428  __g = M_g_linear;
429  }
430 
431  else
432  {
433  FEELPP_ASSERT( __pt.size() == dim() )( __pt.size() )( dim() ).error( "invalid dimension" );
434 
435  matrix_node_t_type __pts( nDim, 1 );
436  ublas::column( __pts, 0 ) = __pt;
437 
438  ublas::vector<ublas::matrix<value_type> > m = super::derivate( __pts );
439 
440  for ( uint16_type n = 0; n < nDim; ++n )
441  {
442  ublas::column( __g, n ) = ublas::column( m( n ), 0 );
443  }
444  }
445 }
446 
455 void gradient( uint16_type __idref,
456  matrix_type& __g,
457  precompute_type const* __pc ) const
458 {
459  FEELPP_ASSERT( __pc )( __idref ).error( "a PreCompute must be set first before using this function" );
460 
461  if ( trans == fem::LINEAR )
462  {
463  __g = M_g_linear;
464  }
465 
466  else
467  {
468  FEELPP_ASSERT( __pc->dim() == dim() )( __pc->dim() )( dim() ).error( "invalid dimension" );
469 
470  for ( size_type i = 0; i < nNodes; ++i )
471  {
472  for ( uint16_type n = 0; n < nDim; ++n )
473  {
474  __g( i, n ) = __pc->grad( i, 0, n, __idref );
475  }
476  }
477  }
478 }
479 
480 
487 value_type radiusEstimate( matrix_node_t_type const& G ) const
488 {
489  size_type N = G.size1();
490  size_type Nref = this->nbNodes();
491  //size_type n = ( this->isLinear()) ? 1 : Nref;
492 
493  // compute K
494  matrix_node_t_type __g( Nref, dim() );
495  // xref is not defined at this point: it serves only as a dummy point
496  // in the case of linear inversion (__g is a constant matrix in this case)
497  this->gradient( this->refNode( 0 ), __g );
498 
499  FEELPP_ASSERT( __g.size1() == G.size2() )( __g.size1() )( G.size2() ).error( "invalid sizes" );
500 
501  matrix_node_t_type K( G.size1(), __g.size2() );
502  ublas::axpy_prod( G, __g, K );
503 
504  SVD<matrix_node_t_type> __svd( K );
505  value_type __max;
506  value_type __min;
507  __svd.conditionNumber( __max, __min );
508  return __max*sqrt( value_type( N ) ) / value_type( N );
509 }
510 
511 
512 bool M_is_cached;
513 bool isCached() const
514 {
515  return M_is_cached;
516 }
517 
518 std::vector<bool> M_cached;
519 std::vector<double> M_J;
520 //boost::multi_array<double,3> M_K;
521 //boost::multi_array<double,3> M_B;
522 std::vector<matrix_type> M_K;
523 std::vector<matrix_type> M_B;
524 
525 template<typename MeshType>
526 void initCache( MeshType const* mesh )
527 {
528  size_type nelts = mesh->numElements();
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 );
532 
533  M_J.resize( nelts );
534  M_K.resize( nelts );
535  M_B.resize( nelts );
536  M_is_cached = true;
537 }
538 bool isCacheValid() const
539 {
540  if ( !( M_cached.size() > 0 &&
541  M_J.size() > 0 &&
542  M_K.size() > 0 &&
543  M_B.size() > 0 ) )
544  {
545  LOG(INFO) << "invalid cache\n";
546  return false;
547  }
548 
549  return true;
550 }
551 bool cached( int e ) const
552 {
553  FEELPP_ASSERT( this->isCacheValid() )( e ).error( "invalid cache" );
554  return M_cached[e];
555 }
556 void setCached( int e, bool v )
557 {
558  FEELPP_ASSERT( this->isCacheValid() )( e ).error( "invalid cache" );
559  M_cached[e] = v;
560 }
561 double J( int e ) const
562 {
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";
565  return M_J[e];
566 }
567 matrix_type const& B( int e ) const
568 {
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";
571  return M_B[e];
572 }
573 matrix_type const& K( int e ) const
574 {
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";
577  return M_K[e];
578 }
579 void addJ( int e, double v )
580 {
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";
583  M_J[e] = v;
584 }
585 void addK( int e, matrix_type const& K )
586 {
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() );
590 
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 );
594 }
595 void addB( int e, matrix_type const& B )
596 {
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() );
600 
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 );
604 }
605 
612 template<size_type context_v, typename ElementType>
613 class Context
614 {
615  public:
616  static const size_type contextv = context_v;
617  static const size_type context = context_v;
618  // reference space dimension
619  static const uint16_type PDim = ElementType::nDim;
620  // real space dimension
621  static const uint16_type NDim = ElementType::nRealDim;
622  static const uint16_type nDim = NDim;
623  // type of transformation (linear or not)
624  static const fem::transformation_type trans = geometric_mapping_type::trans;
625  static const bool is_linear = ( trans == fem::LINEAR );
626 
627  static const bool condition = ( ( PDim==NDim )||( ( NDim>=1 )&&( PDim==NDim-1 ) ) );
628  //BOOST_MPL_ASSERT_MSG( condition, INVALID_DIM, (mpl::int_<NDim>, mpl::int_<PDim>, ElementType ) );
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> > >,
633  //typename mpl::if_<mpl::equal_to<mpl::int_<PDim>, mpl::int_<NDim-1> >,
634  // mpl::identity<typename GeoMap<Dim, Order, T, Entity, PP >::template face_gm<NDim>::type>,
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> > >,
638  //typename mpl::if_<mpl::equal_to<mpl::int_<PDim>, mpl::int_<NDim-1> >,
639  // mpl::identity<typename GeoMap<Dim, Order, T, Entity, PP >::template face_gm<NDim>::type>,
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;
643 
644  typedef typename gm_type::value_type value_type;
645 
646  typedef typename gm_type::precompute_ptrtype precompute_ptrtype;
647 
649  typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
650 
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;
657 
658  typedef ElementType element_type;
659  typedef typename element_type::permutation_type permutation_type;
660 
661  Context( gm_ptrtype __gm,
662  element_type const& __e,
663  precompute_ptrtype const& __pc )
664  :
665  M_gm( __gm ),
666  M_element( boost::addressof( __e ) ),
667  M_pc( __pc ),
668  M_pc_faces(),
669  M_npoints( M_pc->nPoints() ),
670 
671  //M_xref( PDim ),
672  //M_xreal( NDim ),
673  //M_x0( NDim ),
674  M_J( 0 ),
675  M_G( ( gm_type::nNodes == element_type::numVertices ) ?__e.vertices() : __e.G() ),
676  M_n( M_gm->referenceConvex().normals() ),
677  M_n_real( NDim ),
678  M_u_n_real( NDim ),
679  M_n_norm( 0 ),
680  M_t_real( NDim ),
681  M_xrefq( PDim, nPoints() ),
682  M_xrealq( NDim, nPoints() ),
683  M_nrealq( NDim, nPoints() ),
684  M_unrealq( NDim, nPoints() ),
685  M_nnormq( nPoints() ),
686 
687  M_g( M_G.size2(), PDim ),
688  M_K( NDim, PDim ),
689  M_CS( PDim, PDim ),
690  M_CSi( PDim, PDim ),
691  M_B( NDim, PDim ),
692  M_B3( boost::extents[NDim][NDim][PDim][PDim] ),
693  M_id( __e.id() ),
694  M_e_marker( __e.marker() ),
695  M_e_marker2( __e.marker2() ),
696  M_e_marker3( __e.marker3() ),
697  M_elem_id_1( invalid_size_type_value ),// __e.ad_first() ),
698  M_pos_in_elem_id_1( invalid_uint16_type_value ), //__e.pos_first() ),
699  M_elem_id_2( invalid_size_type_value ), //__e.ad_second() ),
700  M_pos_in_elem_id_2( invalid_uint16_type_value ), //__e.pos_second() ),
701  M_face_id( invalid_uint16_type_value ),
702  M_h( __e.h() ),
703  M_h_face( 0 ),
704  M_meas( __e.measure() ),
705  M_measface( 0 ),
706  M_Jt(),
707  M_Bt(),
708  M_perm( )
709 {
710 
711  if ( is_linear )
712  {
713  M_gm->gradient( node_t_type(), M_g_linear );
714  }
715 
716  else
717  {
718  M_Jt.resize( nPoints() );
719  M_Bt.resize( nPoints() );
720  }
721 
722  update( __e );
723 }
724 
725 Context( gm_ptrtype __gm,
726  element_type const& __e,
727  std::vector<std::map<permutation_type, precompute_ptrtype> > & __pc,
728  uint16_type __f )
729  :
730  M_gm( __gm ),
731  M_element( boost::addressof( __e ) ),
732  M_pc(),
733  M_pc_faces( __pc ),
734  M_npoints( __pc[__f][__e.permutation( __f )]->nPoints() ),
735 
736  //M_xref( PDim ),
737  //M_xreal( NDim ),
738  //M_x0( NDim ),
739  M_J( 0 ),
740  M_G( ( gm_type::nNodes == element_type::numVertices ) ?__e.vertices() : __e.G() ),
741  M_n( M_gm->referenceConvex().normals() ),
742  M_n_real( NDim ),
743  M_u_n_real( NDim ),
744  M_n_norm( 0 ),
745  M_t_real( NDim ),
746  M_xrefq( PDim, nPoints() ),
747  M_xrealq( NDim, nPoints() ),
748  M_nrealq( NDim, nPoints() ),
749  M_unrealq( NDim, nPoints() ),
750  M_nnormq( nPoints() ),
751 
752  M_g( M_G.size2(), PDim ),
753  M_K( NDim, PDim ),
754  M_CS( PDim, PDim ),
755  M_CSi( PDim, PDim ),
756  M_B( NDim, PDim ),
757  M_B3( boost::extents[NDim][NDim][PDim][PDim] ),
758  M_id( __e.id() ),
759  M_e_marker( __e.marker() ),
760  M_e_marker2( __e.marker2() ),
761  M_e_marker3( __e.marker3() ),
762  M_elem_id_1( invalid_size_type_value ),// __e.ad_first() ),
763  M_pos_in_elem_id_1( invalid_uint16_type_value ), //__e.pos_first() ),
764  M_elem_id_2( invalid_size_type_value ), //__e.ad_second() ),
765  M_pos_in_elem_id_2( invalid_uint16_type_value ), //__e.pos_second() ),
766  M_face_id( __f ),
767  M_h( __e.h() ),
768  M_h_face( 0 ),
769  M_meas( __e.measure() ),
770  M_measface( __e.faceMeasure( __f ) ),
771  M_Jt(),
772  M_Bt(),
773  M_perm( )
774 {
775 
776  if ( is_linear )
777  {
778  M_gm->gradient( node_t_type(), M_g_linear );
779  }
780 
781  else
782  {
783  M_Jt.resize( nPoints() );
784  M_Bt.resize( nPoints() );
785  }
786 
787  update( __e, __f );
788 }
789 Context( gmc_ptrtype& p )
790  :
791  M_gm( p->M_gm ),
792  M_element( p->M_element ),
793  M_pc( p->M_pc ),
794  M_pc_faces( p->M_pc_faces ),
795  M_npoints( M_pc->nPoints() ),
796  //M_xref( PDim ),
797  //M_xreal( NDim ),
798  //M_x0( NDim ),
799  M_J( p->M_J ),
800  M_G( ( gm_type::nNodes == element_type::numVertices ) ?M_element->vertices() : M_element->G() ),
801  M_n( p->M_n ),
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 ),
811  M_g( p->M_g ),
812  M_K( p->M_K ),
813  M_CS( p->M_CS ),
814  M_CSi( p->M_CSi ),
815  M_B( p->M_B ),
816  M_B3( p->M_B3 ),
817  M_id( p->M_id ),
818  M_e_marker( p->M_e_marker ),
819  M_e_marker2( p->M_e_marker2 ),
820  M_e_marker3( p->M_e_marker3 ),
821  M_elem_id_1( invalid_size_type_value ),// M_element.ad_first() ),
822  M_pos_in_elem_id_1( invalid_uint16_type_value ), //M_element.pos_first() ),
823  M_elem_id_2( invalid_size_type_value ), //M_element.ad_second() ),
824  M_pos_in_elem_id_2( invalid_uint16_type_value ), //M_element.pos_second() ),
825  M_face_id( invalid_uint16_type_value ),
826  M_h( p->M_h ),
827  M_h_face( p->M_h_face ),
828  M_meas( p->M_meas ),
829  M_measface( p->M_measface ),
830  M_Jt(),
831  M_Bt(),
832  M_perm( p->M_perm )
833 {
834  if ( is_linear )
835  {
836  M_gm->gradient( node_t_type(), M_g_linear );
837  }
838 
839  else
840  {
841  M_Jt.resize( nPoints() );
842  M_Bt.resize( nPoints() );
843  }
844 
845  update( *M_element );
846 }
847 
851 gmc_ptrtype clone()
852 {
853  return gmc_ptrtype( new gmc_type( *this ) );
854 }
855 
868 void update( element_type const& __e, uint16_type __f )
869 {
870  //M_element_c = boost::shared_ptr<element_type const>(&__e);
871  M_element = boost::addressof( __e );
872  M_face_id = __f;
873 
874  M_perm = __e.permutation( M_face_id );
875 
876  M_h_face = __e.hFace( M_face_id );
877  //M_h_edge = __e.hEdge( M_face_id );
878 
879  M_pc = M_pc_faces[__f][M_perm];
880  //M_G = __e.G();
881  M_G = ( gm_type::nNodes == element_type::numVertices ) ?__e.vertices() : __e.G();
882  M_id = __e.id();
883  M_e_marker = __e.marker();
884  M_e_marker2 = __e.marker2();
885  M_e_marker3 = __e.marker3();
886  M_h = __e.h();
887  M_meas = __e.measure();
888  M_measface = __e.faceMeasure( __f );
889  M_xrefq = M_pc->nodes();
890 
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" );
893 
894  if ( vm::has_point<context>::value )
895  {
896 
897  //ublas::axpy_prod( M_G, pc->phi(), M_xrealq, true );
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();
902 
903  for ( uint16_type i = 0; i < size1; ++i )
904  for ( uint16_type j = 0; j < size2; ++j )
905  {
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 );
908  }
909  }
910 
911  if ( vm::has_jacobian<context>::value )
912  {
913  updateJKBN( mpl::bool_<is_linear>() );
914  }
915 
916 }
917 
918 void update( element_type const& __e, uint16_type __f, permutation_type __perm )
919 {
920  //M_element_c = boost::shared_ptr<element_type const>(&__e);
921  M_element = boost::addressof( __e );
922  M_face_id = __f;
923 
924  M_perm = __perm;
925 
926  M_h_face = __e.hFace( M_face_id );
927  //M_h_edge = __e.hEdge( M_face_id );
928 
929  M_pc = M_pc_faces[__f][M_perm];
930  //M_G = __e.G();
931  M_G = ( gm_type::nNodes == element_type::numVertices ) ?__e.vertices() : __e.G();
932  M_id = __e.id();
933  M_e_marker = __e.marker();
934  M_e_marker2 = __e.marker2();
935  M_e_marker3 = __e.marker3();
936  M_h = __e.h();
937  M_meas = __e.measure();
938  M_measface = __e.faceMeasure( __f );
939  M_xrefq = M_pc->nodes();
940 
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" );
943 
944  if ( vm::has_point<context>::value )
945  {
946 
947  //ublas::axpy_prod( M_G, pc->phi(), M_xrealq, true );
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();
952 
953  for ( uint16_type i = 0; i < size1; ++i )
954  for ( uint16_type j = 0; j < size2; ++j )
955  {
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 );
958  }
959  }
960 
961  if ( vm::has_jacobian<context>::value )
962  {
963  updateJKBN( mpl::bool_<is_linear>() );
964  }
965 
966 }
967 
968 void update( element_type const& __e,
969  precompute_ptrtype const& __pc )
970 {
971  M_pc = __pc;
972 
973 
974  if ( M_npoints != M_pc->nPoints() )
975  {
976  M_npoints = M_pc.get()->nPoints();
977 
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() );
983 
984  if ( is_linear )
985  {
986  M_gm->gradient( node_t_type(), M_g_linear );
987  }
988 
989  else
990  {
991  M_Jt.resize( nPoints() );
992  M_Bt.resize( nPoints() );
993  }
994  }
995 
996  update( __e );
997 }
998 
1011 void update( element_type const& __e )
1012 {
1013  M_G = ( gm_type::nNodes == element_type::numVertices ) ?__e.vertices() : __e.G();
1014  //M_G = __e.G();
1015  M_g.resize( M_G.size2(), PDim );
1016  //M_element_c = boost::shared_ptr<element_type const>(&__e);
1017  M_element = boost::addressof( __e );
1018  M_id = __e.id();
1019  M_e_marker = __e.marker();
1020  M_e_marker2 = __e.marker2();
1021  M_e_marker3 = __e.marker3();
1022  M_face_id = invalid_uint16_type_value;
1023  M_h = __e.h();
1024  M_meas = __e.measure();
1025  M_xrefq = M_pc->nodes();
1026 
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" );
1029 
1030  if ( vm::has_point<context>::value )
1031  {
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();
1036 
1037  for ( uint16_type i = 0; i < size1; ++i )
1038  for ( uint16_type j = 0; j < size2; ++j )
1039  {
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 );
1042  }
1043  }
1044 
1045  if ( vm::has_jacobian<context>::value )
1046  {
1047  updateJKBN( mpl::bool_<is_linear>() );
1048  }
1049 
1050 }
1051 
1052 ~Context()
1053 {
1054 }
1055 
1059 
1063 gm_ptrtype const& geometricMapping() const
1064 {
1065  return M_gm;
1066 }
1067 
1071 uint16_type N() const
1072 {
1073  return NDim;
1074 }
1075 
1079 uint16_type P() const
1080 {
1081  return PDim;
1082 }
1083 
1087 element_type const& element() const
1088 {
1089  return *M_element;
1090 }
1091 
1092 element_type const& element_c() const
1093 {
1094  return *M_element;
1095 }
1096 
1100 uint16_type nPoints() const
1101 {
1102  return M_npoints;
1103 }
1104 
1108 matrix_node_t_type const& xRefs() const
1109 {
1110  return M_xrefq;
1111 }
1112 
1116 ublas::matrix_column<matrix_node_t_type const> xRef( int q ) const
1117 {
1118  return ublas::column( M_xrefq, q );
1119 }
1120 
1124 //matrix_node_t_type const& xReal() const
1125 matrix_type const& xReal() const
1126 {
1127  //BOOST_STATIC_ASSERT( vm::has_point<context>::value );
1128  return M_xrealq;
1129 }
1130 
1134 ublas::matrix_column<matrix_type const> xReal( int q ) const
1135 {
1136  // BOOST_STATIC_ASSERT( vm::has_point<context>::value );
1137  return ublas::column( M_xrealq, q );
1138 }
1139 
1146 #if 0
1147 value_type J( int q ) const
1148 {
1149  //BOOST_STATIC_ASSERT( vm::has_jacobian<context>::value );
1150  //return J( q, mpl::bool_<is_linear>() );
1151  if ( is_linear )
1152  return M_J;
1153 
1154  else
1155  return M_Jt[q];
1156 }
1157 #else
1158 value_type J( int q ) const
1159 {
1160  if ( is_linear )
1161  return M_J;
1162 
1163  else
1164  return M_Jt[q];
1165 }
1166 #endif
1167 #if 0
1168 
1171 value_type J( int /*q*/, mpl::bool_<true> ) const
1172 {
1173  return M_J;
1174 }
1175 
1179 value_type J( int q, mpl::bool_<false> ) const
1180 {
1181  return M_Jt[q];
1182 }
1183 #endif
1184 
1188 //matrix_node_t_type const& G() const { return M_G; }
1189 matrix_type const& G() const
1190 {
1191  return M_G;
1192 }
1193 
1200 matrix_type const& B( int i ) const
1201 {
1202  //BOOST_STATIC_ASSERT( vm::has_kb<context>::value );
1203  return B( i, mpl::bool_<is_linear>() );
1204 }
1205 value_type B( int c1, int c2, int q ) const
1206 {
1207  //BOOST_STATIC_ASSERT( vm::has_kb<context>::value );
1208  return B( q, mpl::bool_<is_linear>() )( c1, c2 );
1209 }
1210 
1211 matrix_type const& K( int i ) const
1212 {
1213  return M_K;
1214 }
1215 value_type K( int c1, int c2, int q ) const
1216 {
1217  return M_K( c1, c2 );
1218 }
1222 matrix_type const& B( int /*i*/, mpl::bool_<true> ) const
1223 {
1224  return M_B;
1225 }
1229 matrix_type const& B( int i, mpl::bool_<false> ) const
1230 {
1231  return M_Bt[i];
1232 }
1233 
1241 boost::multi_array<value_type,4> const& B3() const
1242 {
1243  return M_B3;
1244 }
1245 
1249 node_t_type barycenterRef() const
1250 {
1251  node_t_type __barycenter( M_gm->dim() );
1252  __barycenter.clear();
1253 
1254  for ( uint16_type __c = 0; __c < M_gm->refNodes().size(); ++__c )
1255  {
1256  __barycenter += M_gm->refNode( __c );
1257  }
1258 
1259  __barycenter /= M_gm->refNodes().size();
1260  return __barycenter;
1261 }
1262 
1266 node_t_type barycenterReal() const
1267 {
1268  node_t_type __barycenter( M_G.size1() );
1269  __barycenter.clear();
1270 
1271 
1272  for ( uint16_type __c = 0; __c < M_G.size2(); ++__c )
1273  {
1274  __barycenter += ublas::column( M_G, __c );
1275  }
1276 
1277  __barycenter /= M_G.size2();
1278  return __barycenter;
1279 }
1280 
1281 
1286 bool isOnConvexSurface() const
1287 {
1288  if ( trans == fem::LINEAR )
1289  {
1290  // x -x0 - K(0)\bar{x}
1291  return std::abs( ublas::norm_2( xReal()-ublas::column( M_G, 0 )-
1292  ublas::prod( M_K, xRef() ) ) ) < 1e-10;
1293  }
1294 
1295  return false;
1296 }
1297 
1298 node_t_type const& refNormal( int /*q*/ ) const
1299 {
1300  return M_gm->referenceConvex().normal( M_face_id );
1301 }
1302 
1308 value_type normalNorm( int q ) const
1309 {
1310  if ( is_linear )
1311  return M_n_norm;
1312 
1313  else
1314  return M_nnormq[q];
1315 }
1316 
1322 node_t_type const& normal() const
1323 {
1324  //BOOST_STATIC_ASSERT( vm::has_normal<context>::value );
1325  return M_n_real;
1326 }
1327 
1333 ublas::matrix_column<matrix_node_t_type const> normal( int q ) const
1334 {
1335  //BOOST_STATIC_ASSERT( vm::has_normal<context>::value );
1336  if ( is_linear )
1337  return M_n_real;
1338 
1339  else
1340  {
1341  return ublas::column( M_nrealq, q );
1342  }
1343 }
1344 
1350 node_t_type const& unitNormal() const
1351 {
1352  //BOOST_STATIC_ASSERT( vm::has_normal<context>::value );
1353  return M_u_n_real;
1354 }
1355 
1356 //ublas::matrix_column<matrix_node_t_type const> unitNormal( int q ) const
1357 // node_t_type const& unitNormal( int q ) const
1358 node_t_type unitNormal( int q ) const
1359 {
1360  //BOOST_STATIC_ASSERT( vm::has_normal<context>::value );
1361  if ( is_linear )
1362  return M_u_n_real;
1363 
1364  else
1365  return ublas::column( M_unrealq, q );
1366 }
1367 
1368 
1369 value_type const& unitNormal( int n, int q ) const
1370 {
1371  //BOOST_STATIC_ASSERT( vm::has_normal<context>::value );
1372 
1373  if ( is_linear )
1374  return M_u_n_real( n );
1375 
1376  else
1377  return M_unrealq( n, q );
1378 }
1379 
1380 node_t_type const& tangent() const
1381 {
1382  //BOOST_STATIC_ASSERT( vm::has_normal<context>::value );
1383  return M_t_real;
1384 }
1385 //ublas::matrix_column<matrix_node_t_type const> unitNormal( int q ) const
1386 // node_t_type const& unitNormal( int q ) const
1387 node_t_type tangent( int q ) const
1388 {
1389  //BOOST_STATIC_ASSERT( vm::has_normal<context>::value );
1390  if ( is_linear )
1391  return M_t_real;
1392 
1393  //else
1394  //return ublas::column( M_utrealq, q );
1395 }
1396 
1397 
1398 value_type const& unitTangent( int n, int q ) const
1399 {
1400  //BOOST_STATIC_ASSERT( vm::has_normal<context>::value );
1401 #if 0
1402  if ( is_linear )
1403  return M_t_real( n );
1404 
1405  else
1406  return M_utrealq( n, q );
1407 
1408 #endif
1409 }
1410 
1416 size_type id() const
1417 {
1418  return M_id;
1419 }
1420 
1421 /*
1422  * \return the face id
1423  */
1424 uint16_type faceId() const
1425 {
1426  return M_face_id;
1427 }
1428 
1432 bool elementIsAFace() const
1433 {
1434  return M_face_id != invalid_uint16_type_value;
1435 }
1436 
1442 Marker1 marker() const
1443 {
1444  return M_e_marker;
1445 }
1446 
1452 Marker2 marker2() const
1453 {
1454  return M_e_marker2;
1455 }
1456 
1462 Marker2 marker3() const
1463 {
1464  return M_e_marker3;
1465 }
1466 
1473 {
1474  return M_elem_id_1;
1475 }
1476 
1482 uint16_type idIn1() const
1483 {
1484  return M_pos_in_elem_id_1;
1485 }
1486 
1493 {
1494  return M_elem_id_2;
1495 }
1496 
1502 uint16_type idIn2() const
1503 {
1504  return M_pos_in_elem_id_2;
1505 }
1506 
1512 value_type radiusEstimate() const
1513 {
1514  return M_gm->radiusEstimate( M_G );
1515 }
1516 
1523 value_type h() const
1524 {
1525  return M_h;
1526 }
1527 
1533 value_type hFace() const
1534 {
1535  return M_h_face;
1536 }
1537 
1538 /*
1539  * @return the measure of the element
1540  */
1541 value_type meas() const
1542 {
1543  return M_meas;
1544 }
1545 
1546 /*
1547  * @return the measure of the set of elements which share a vertex with \p M_elements including himself
1548  */
1549 value_type measurePointElementNeighbors() const
1550 {
1551  return M_element->measurePointElementNeighbors();
1552 }
1553 
1554 
1555 /*
1556  * @return the measure of the (current) face of the element
1557  */
1558 value_type measFace() const
1559 {
1560  return M_measface;
1561 }
1562 
1563 
1567 permutation_type permutation() const
1568 {
1569  return permutation( mpl::bool_<( nDim>=2 )>() );
1570 }
1571 
1575 precompute_ptrtype const& pc() const
1576 {
1577  return M_pc;
1578 }
1579 
1583 std::vector<std::map<permutation_type, precompute_ptrtype> > const & pcFaces() const
1584 {
1585  return M_pc_faces;
1586 }
1588 
1592 
1596 void setPc( precompute_ptrtype const& __pc )
1597 {
1598  M_pc = __pc;
1599 }
1600 
1601 void setPcFaces( std::vector<std::map<permutation_type, precompute_ptrtype> > const& __pcfaces )
1602 {
1603  M_pc_faces = __pcfaces;
1604 }
1605 
1607 private:
1608 
1612 permutation_type permutation( mpl::bool_<false> ) const
1613 {
1614  return M_perm;
1615 }
1616 
1620 permutation_type permutation( mpl::bool_<true> ) const
1621 {
1622  FEELPP_ASSERT( M_face_id == invalid_uint16_type_value ||
1623  ( M_face_id != invalid_uint16_type_value &&
1624  M_perm != permutation_type( permutation_type::NO_PERMUTATION ) ) )
1625  ( M_face_id ).error( "invalid permutation" );
1626  return M_perm;
1627 }
1628 
1632 void updateJKBN( mpl::bool_<true> )
1633 {
1634 
1635  if ( !M_gm->isCached() ||
1636  ( M_gm->isCached() && M_gm->cached( M_id ) == false ) )
1637  {
1638 #if 0
1639 
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,
1643  0.0, M_K );
1644 
1645  else
1646 #endif
1647  ublas::axpy_prod( M_G, M_g_linear, M_K, true );
1648 
1649  if ( NDim == PDim )
1650  {
1651  M_J = math::abs( det<NDim>( M_K ) );
1652  //if ( vm::has_kb<context>::value )
1653  {
1654 #if 0
1655  inverse<NDim>( M_K, M_CS, M_J );
1656  ublas::noalias( M_B ) = ublas::trans( M_CS );
1657 #else
1658  inverse<NDim>( M_K, M_CS );
1659  ublas::noalias( M_B ) = ublas::trans( M_CS );
1660 #endif
1661  }
1662  }
1663 
1664  else // N != P
1665  {
1666  // CS = K^T K
1667 #if 0
1668  if ( boost::is_arithmetic<value_type>::value )
1669  atlas::gemm( traits::TRANSPOSE, traits::NO_TRANSPOSE,
1670  1.0, M_K, M_K,
1671  0.0, M_CS );
1672 
1673  else
1674 #endif
1675 
1676  ublas::noalias( M_CS ) = ublas::prod( ublas::trans( M_K ), M_K );
1677 
1678  M_J = math::sqrt( math::abs( det<PDim>( M_CS ) ) );
1679  //if ( vm::has_kb<context>::value )
1680  {
1681  inverse<PDim>( M_CS, M_CSi );
1682  // B = K CS
1683 #if 0
1684 
1685  if ( boost::is_arithmetic<value_type>::value )
1686  atlas::gemm( traits::NO_TRANSPOSE, traits::NO_TRANSPOSE,
1687  1.0, M_K, M_CSi,
1688  0.0, M_B );
1689 
1690  else
1691 #endif
1692  ublas::axpy_prod( M_K, M_CSi, M_B, true );
1693  }
1694 
1695  }
1696 
1697  if ( M_gm->isCached() )
1698  {
1699  // cache J, K and B
1700  M_gm->addJ( M_id, M_J );
1701  //if ( vm::has_kb<context>::value )
1702  {
1703  M_gm->addK( M_id, M_K );
1704  M_gm->addB( M_id, M_B );
1705  }
1706  M_gm->setCached( M_id, true );
1707 
1708  //LOG(INFO) << "(add to cache) J[" << M_id << "]=" << M_J << "\n";
1709  //LOG(INFO) << "(add to cache) B[" << M_id << "]=" << M_B << "\n";
1710  }
1711  }
1712 
1713  else
1714  {
1715  M_J = M_gm->J( M_id );
1716  //LOG(INFO) << "(use cache) J[" << M_id << "]=" << M_J << "\n";
1717  //if ( vm::has_kb<context>::value )
1718  {
1719  M_K = M_gm->K( M_id );
1720  M_B = M_gm->B( M_id );
1721  //LOG(INFO) << "(use cache) B[" << M_id << "]=" << M_B << "\n";
1722  }
1723 
1724 
1725  }
1726 
1727  if ( vm::has_hessian<context>::value )
1728  {
1729 
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 );
1735 
1736 #if 0
1737  //B3(k + N_*l, i + P*j) = BB(k, i) * BB(l, j);
1738  std::cout << "element " << this->id() << " B3 = " << "\n";
1739 
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";
1746 
1747 #endif
1748  }
1749 
1750  if ( ( ( NDim != PDim ) || ( vm::has_normal<context>::value ) ) && ( M_face_id != invalid_uint16_type_value ) )
1751  {
1752 #if 0
1753  blas::gemv( traits::NO_TRANSPOSE,
1754  1.0, M_Bt[ q ], M_n[M_face_id],
1755  0.0, M_n_real );
1756 #else
1757  ublas::axpy_prod( M_B,
1758  M_gm->referenceConvex().normal( M_face_id ),
1759  M_n_real,
1760  true );
1761 #endif
1762  M_n_norm = ublas::norm_2( M_n_real );
1763  M_u_n_real = M_n_real/M_n_norm;
1764 
1765  }
1766 
1767  if ( vm::has_tangent<context>::value && ( M_face_id != invalid_uint16_type_value ) )
1768  {
1769  // t = |\hat{e}|*o_K*(K*t_ref)/|e| where o_K is the sign(e*x_K(\hat{e}))
1770  ublas::axpy_prod( M_K,
1771  M_gm->referenceConvex().tangent( M_face_id ),
1772  M_t_real,
1773  true );
1774  double ratio = M_gm->referenceConvex().h( M_face_id )/M_h_face;
1775 
1776  M_t_real *= ratio;
1777  }
1778 
1779 }
1780 
1784 void updateJKBN( mpl::bool_<false> )
1785 {
1786  //std::cout << "nPoints() =" << nPoints() << "\n";
1787  //VLOG(1) << "[geomap] G = "<< M_G << "\n";
1788  //double res = 0;
1789  for ( int q = 0; q < nPoints(); ++q )
1790  {
1791  //std::cout << "q =" << q << "\n";
1792  M_gm->gradient( q, M_g, M_pc.get() );
1793  //VLOG(1) << "[geomap] g[" << q << "] = "<< M_g << "\n";
1794 
1795 #if 0
1796  blas::gemm( M_G, M_g, M_K );
1797 #else
1798  ublas::axpy_prod( M_G, M_g, M_K, true );
1799 #endif
1800 
1801  if ( NDim == PDim )
1802  {
1803  M_J = math::abs( det<NDim>( M_K ) );
1804 
1805  if ( vm::has_kb<context>::value )
1806  {
1807  inverse<NDim>( M_K, M_CS );
1808  //ublas::noalias(M_B) = ublas::trans( M_CS );
1809  M_B = ublas::trans( M_CS );
1810  //std::cout << "========== B[" << q << "]=" << M_B << "\n";
1811  }
1812  }
1813 
1814  else // N != P
1815  {
1816  // CS = K^T K
1817 #if 1
1818  ublas::prod( ublas::trans( M_K ), M_K, M_CS );
1819 #else
1820  blas::gemm( traits::TRANSPOSE, traits::NO_TRANSPOSE,
1821  1.0, M_K, M_K,
1822  0.0, M_CS );
1823 #endif
1824  M_J = math::sqrt( math::abs( det<PDim>( M_CS ) ) );
1825 
1826  if ( vm::has_kb<context>::value )
1827  {
1828  inverse<PDim>( M_CS, M_CSi );
1829  // B = K CS
1830 #if 1
1831  ublas::axpy_prod( M_K, M_CSi, M_B );
1832 #else
1833  blas::gemm( traits::NO_TRANSPOSE, traits::NO_TRANSPOSE,
1834  1.0, M_K, M_CSi,
1835  0.0, M_B );
1836 #endif
1837  }
1838 
1839  }
1840 
1841  //VLOG(1) << "[geomap] J[" << q << "]= "<< M_J << "\n";
1842  //res += M_J;
1843  // store q-th jacobian entry
1844 
1845  M_Jt[q] = M_J;
1846 
1847  if ( vm::has_kb<context>::value )
1848  {
1849  //VLOG(1) << "[geomap] B[" << q << "]= "<< M_B << "\n";
1850  M_Bt[q].resize( M_B.size1(), M_B.size2() );
1851  M_Bt[q] = M_B;
1852  }
1853 
1854 
1855  }
1856 
1857  //VLOG(1) << "[geomap] res(sum J) = " << res << "\n";
1858  if ( ( ( NDim != PDim ) || ( vm::has_normal<context>::value ) ) && ( M_face_id != invalid_uint16_type_value ) )
1859  {
1860  //std::cout << "has normal\n";
1861  for ( int q = 0; q < nPoints(); ++q )
1862  {
1863 #if 0
1864  blas::gemv( traits::NO_TRANSPOSE,
1865  1.0, M_Bt[ q ], M_n[M_face_id],
1866  0.0, M_n_real );
1867 #else
1868 
1869  if ( 0 ) //trans == fem::LINEAR )
1870  {
1871  ublas::axpy_prod( M_B,
1872  M_gm->referenceConvex().normal( M_face_id ),
1873  M_n_real,
1874  true );
1875  }
1876 
1877  else
1878  ublas::axpy_prod( M_Bt[ q ],
1879  M_gm->referenceConvex().normal( M_face_id ),
1880  M_n_real,
1881  true );
1882 
1883 #endif
1884  //std::cout << "[geomap] point " << q << " n_real = " << M_n_real << "\n";
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;
1890 
1891  if ( NDim != PDim )
1892  M_Jt[q] *= M_n_norm;
1893 
1894  if ( vm::has_tangent<context>::value )
1895  {
1896 
1897  }
1898  }
1899 
1900 #if 0
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";
1906 #endif
1907  }
1908 }
1909 
1910 Context();
1911 
1912 private:
1913 
1914 gm_ptrtype M_gm;
1915 
1916 element_type const* M_element;
1917 //boost::shared_ptr<element_type const> M_element_c;
1918 //element_type M_element_c;
1919 
1920 precompute_ptrtype M_pc;
1921 std::vector<std::map<permutation_type, precompute_ptrtype> > M_pc_faces;
1922 uint16_type M_npoints;
1923 
1924 value_type M_J;
1925 
1926 //matrix_node_t_type M_G;
1927 matrix_type M_G;
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;
1933 
1934 //matrix_node_t_type const& M_xrefq;
1935 matrix_node_t_type M_xrefq;
1936 //matrix_node_t_type M_xrealq;
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;
1941 
1942 matrix_type M_g_linear;
1943 matrix_type M_g;
1944 matrix_type M_K;
1945 matrix_type M_CS;
1946 matrix_type M_CSi;
1947 matrix_type M_B;
1948 boost::multi_array<value_type,4> M_B3;
1949 
1950 size_type M_id;
1951 Marker1 M_e_marker;
1952 Marker2 M_e_marker2;
1953 Marker3 M_e_marker3;
1954 size_type M_elem_id_1;
1955 uint16_type M_pos_in_elem_id_1;
1956 size_type M_elem_id_2;
1957 uint16_type M_pos_in_elem_id_2;
1958 
1959 uint16_type M_face_id;
1960 
1961 value_type M_h;
1962 value_type M_h_face;
1963 value_type M_meas;
1964 value_type M_measface;
1965 
1966 vector_type M_Jt;
1967 std::vector<matrix_type> M_Bt;
1968 
1969 permutation_type M_perm;
1970 }; // Context
1971 
1972 
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 )
1976 {
1977  return boost::shared_ptr<Context<context_v,ElementType> >(
1978  new Context<context_v, ElementType>( gm,
1979  e,
1980  pc ) );
1981 }
1982 
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 )
1986 {
1987  return boost::shared_ptr<Context<context_v,ElementType> >(
1988  new Context<context_v, ElementType>( this->shared_from_this(),
1989  e,
1990  pc ) );
1991 }
1992 
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,
1998  uint16_type f )
1999 {
2000  return boost::shared_ptr<Context<context_v,ElementType> >(
2001  new Context<context_v, ElementType>( gm,
2002  e,
2003  pc,
2004  f ) );
2005 }
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,
2010  uint16_type f )
2011 {
2012  return boost::shared_ptr<Context<context_v,ElementType> >(
2013  new Context<context_v, ElementType>( this->shared_from_this(),
2014  e,
2015  pc,
2016  f ) );
2017 }
2018 
2019 
2023 class Inverse
2024 {
2025  public:
2026 
2028  typedef boost::shared_ptr<geometric_mapping_type> geometric_mapping_ptrtype;
2029 
2030  static const fem::transformation_type trans = geometric_mapping_type::trans;
2031 
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;
2036 
2037  template<typename GeoElem>
2038  Inverse( geometric_mapping_ptrtype __gm, GeoElem const& __ge,
2039  WorldComm const& worldComm = Environment::worldComm().subWorldCommSeq() )
2040  :
2041  M_gm( __gm ),
2042  M_xref( __gm->dim() ),
2043  M_xreal( __ge.G().size1() ),
2044  M_is_in( false ),
2045  M_G( __ge.G() ),
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 )
2051  M_nlsolver( SolverNonLinear<double>::build( SOLVERS_PETSC, worldComm ) )
2052 #else
2053  M_nlsolver( SolverNonLinear<double>::build( SOLVERS_GMM, worldComm ) )
2054 #endif
2055 {
2056  FEELPP_ASSERT( M_G.size2() == __gm->nbPoints() )
2057  ( M_G.size2() )( __gm->nbPoints() ).error( "invalid dimensions" );
2058 
2059  if ( M_gm->isLinear() )
2060  {
2061  update();
2062  }
2063  else
2064  {
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 );
2068 #else
2069 
2070 #endif
2071  }
2072 }
2073 
2074 
2075  template<typename GeoElem>
2076  Inverse( geometric_mapping_ptrtype __gm, GeoElem const& __ge, mpl::int_<1> ,
2077  WorldComm const& worldComm = Environment::worldComm().subWorldCommSeq())
2078  :
2079  M_gm( __gm ),
2080  M_xref( __gm->dim() ),
2081  M_xreal( __ge.vertices().size1() ),
2082  M_is_in( false ),
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 )
2089  M_nlsolver( SolverNonLinear<double>::build( SOLVERS_PETSC,worldComm) )
2090 #else
2091  M_nlsolver( SolverNonLinear<double>::build( SOLVERS_GMM,worldComm ) )
2092 #endif
2093 {
2094  FEELPP_ASSERT( M_G.size2() == __gm->nbPoints() )
2095  ( M_G.size2() )( __gm->nbPoints() ).error( "invalid dimensions" );
2096 
2097  if ( M_gm->isLinear() )
2098  {
2099  update();
2100  }
2101  else
2102  {
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 );
2106 #else
2107  //M_nlsolver( SolverNonLinear<double>::build( SOLVERS_GMM ) )
2108 #endif
2109  }
2110 }
2111 
2112 
2113 
2114 
2118 
2122 geometric_mapping_ptrtype const& geometricMapping() const
2123 {
2124  return M_gm;
2125 }
2126 
2130 uint16_type N() const
2131 {
2132  return G().size1();
2133 }
2134 uint16_type P() const
2135 {
2136  return M_gm->dim();
2137 }
2138 
2142 node_t_type const& xRef() const
2143 {
2144  return M_xref;
2145 }
2146 
2150 node_t_type const& xReal() const
2151 {
2152  return M_xreal;
2153 }
2154 
2158 matrix_type const& G() const
2159 {
2160  return M_G;
2161 }
2162 
2166 matrix_type const& K() const
2167 {
2168  return M_K;
2169 }
2170 
2171 value_type J() const
2172 {
2173  return math::abs( det<Dim>( M_K ) );
2174 }
2175 
2180 matrix_type const& B() const
2181 {
2182  return M_B;
2183 }
2184 
2188 node_type barycenterRef() const
2189 {
2190  node_type __barycenter( M_gm->dim() );
2191  __barycenter.clear();
2192 
2193  for ( uint16_type __c = 0; __c < M_gm->referenceConvex().nPoints(); ++__c )
2194  {
2195  __barycenter += M_gm->refNode( __c );
2196  }
2197 
2198  __barycenter /= M_gm->referenceConvex().nPoints();
2199  return __barycenter;
2200 }
2201 
2202 
2206 node_type barycenterReal() const;
2207 
2212 bool isIn() const
2213 {
2214  return M_is_in;
2215 }
2216 
2221 bool isOnConvexSurface() const
2222 {
2223  if ( M_gm->isLinear() )
2224  {
2225  // x -x0 - K(0)\bar{x}
2226  return std::abs( ublas::norm_2( xReal()-ublas::column( M_G, 0 )-
2227  ublas::prod( M_K, xRef() ) ) ) < 1e-10;
2228  }
2229 
2230  return false;
2231 }
2232 
2234 
2238 
2242 void setXReal( node_type const& __xreal )
2243 {
2244  M_xreal = __xreal;
2245 
2246  if ( M_gm->isLinear() )
2247  {
2248  update();
2249  M_is_in = linearInverse();
2250  }
2251 
2252  else
2253  {
2254  M_is_in = nonLinearInversePetsc();
2255 
2256  //bool isin = nonLinearInverse();
2257  }
2258 }
2259 
2260 matrix_node_t_type operator() ( matrix_node_t_type const& real_pts, bool allow_extrapolation = false ) const
2261 {
2262  if ( trans == fem::LINEAR )
2263  return linearInversePoints( real_pts, allow_extrapolation );
2264 
2265  else
2266  return matrix_node_t_type();
2267 }
2268 
2269 
2271 private:
2272 
2273 class GeomapInverseConvex
2274 {
2275  public:
2276 
2277  GeomapInverseConvex( Inverse& __gmi,
2278  const node_t_type& __xr )
2279  :
2280  gmi( __gmi ),
2281  xreal( __xr )
2282 {}
2283 // evaluate
2284 scalar_type operator()( const node_t_type& x ) const
2285 {
2286  node_type r = gmi.geometricMapping()->transform( x, gmi.G() ) - xreal;
2287  return ublas::inner_prod( r, r )/2.;
2288 }
2289 void operator()( const node_t_type& x, node_t_type& gr ) const
2290 {
2291  // gmi.setXreal( xreal );
2292  gmi.M_xref.assign( x );
2293  gmi.update();
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 );
2297 }
2298 private:
2299 
2300 Inverse& gmi;
2301 node_type xreal;
2302 
2303 };
2304 typedef ublas::vector<double> dense_vector_type;
2305 typedef ublas::matrix<double> dense_matrix_type;
2306 
2307 void updateResidual( dense_vector_type const& x, dense_vector_type& r )
2308 {
2309  dense_vector_type y = M_gm->transform( x, M_G );
2310 
2311  if ( N() == P() )
2312  r = y - M_xreal ;
2313 
2314  else
2315  {
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 );
2319  }
2320 
2321 #if 0
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";
2325 
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";
2330 #endif // 0
2331 
2332 }
2333 
2334 void updateJacobian( dense_vector_type const& x, dense_matrix_type& j )
2335 {
2336  M_gm->gradient( x, M_g );
2337 
2338  if ( N() == P() )
2339  ublas::axpy_prod( M_G, M_g, j );
2340 
2341  else
2342  {
2343  ublas::axpy_prod( M_G, M_g, M_K );
2344  ublas::prod( ublas::trans( M_K ), M_K, j );
2345  }
2346 
2347 #if 0
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";
2352 #endif
2353 }
2354 
2367 void update()
2368 {
2369  // xref is not defined at this point: it serves only as a dummy point
2370  // in the case of linear inversion (__g is a constant matrix in this case)
2371  // in the non linear case , xRef() is update in the newton iterations
2372  M_gm->gradient( xRef(), M_g );
2373  DVLOG(2) << "[update] g = " << M_g << "\n";
2374 
2375  checkInvariant();
2376 
2377  ublas::axpy_prod( M_G, M_g, M_K );
2378  DVLOG(2) << "[update] K(0) = " << M_K << "\n";
2379 
2380  // compute B
2381  if ( M_gm->dim() != N() )
2382  {
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 );
2387  }
2388 
2389  else
2390  {
2391  LU<matrix_type> __lu( M_K );
2392  __lu.inverse( M_CS );
2393  M_B = ublas::trans( M_CS );
2394  }
2395 
2396  DVLOG(2) << "[update] B(0) = " << M_B << "\n";
2397 }
2398 
2399 void
2400 checkInvariant() const
2401 {
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" );
2406 }
2407 
2408 
2409 bool linearInverse()
2410 {
2411  checkInvariant();
2412 
2413  size_type N = M_xreal.size();
2414  size_type P = M_xref.size();
2415 
2416  node_type y( M_xreal );
2417 
2418  DVLOG(2) << "y = xreal = " << y << "\n";
2419  //DVLOG(2) << "G(0) = " << node_type( M_x0 << "\n";
2420  y.minus_assign( ublas::column( M_G, 0 ) );
2421  DVLOG(2) << "y - G(0) = " << y << "\n";
2422 
2423  DVLOG(2) << "B(0) = " << M_B << "\n";
2424  DVLOG(2) << "xref = " << ublas::prod( ublas::trans( M_B ), y ) << "\n";
2425 
2426  // xref = B^T * y = B^T * ( x_real - x_0)
2427  M_xref.assign( ublas::prod( ublas::trans( M_B ), y )-ublas::scalar_vector<value_type>( P, 1.0 ) );
2428 
2429  DVLOG(2) << "[GeoMap::Inverse::linearInverse] xref : " << M_xref << "\n";
2430 
2431  bool __isin;
2432  double vmin;
2433  boost::tie( __isin, vmin ) = M_gm->isIn( M_xref );
2434  DVLOG(2) << "[GeoMap::Inverse::linearInverse] isIn : " << __isin << "\n";
2435 
2437  if ( __isin )
2438  {
2439  if ( N == P )
2440  return true;
2441 
2442  else
2443  {
2444  // y = y - K * x_ref
2445  ublas::axpy_prod( M_K, -M_xref, y );
2446 
2447  if ( ublas::norm_2( y ) < 1e-10 )
2448  return true;
2449  }
2450  }
2451 
2452  return false;
2453 }
2454 
2455 matrix_node_t_type linearInversePoints( matrix_node_t_type const& real_pts, bool /*allow_extrapolation*/ = false ) const
2456 {
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 ) );
2461 
2462 }
2463 
2464 /*
2465  *inversion for non-linear geometric transformations
2466  * (Newton on Grad(pgt)(y - pgt(x)) = 0 )
2467  */
2468 bool nonLinearInversePetsc()
2469 {
2470  //LOG(INFO) << "starting new nonlinear inverse\n";
2471  //const double EPS = 1e-10;
2472  const double IN_EPS = 1e-10;
2473  size_type N = M_xreal.size();
2474  size_type P = M_xref.size();
2475 
2476 #if 0
2477  /*
2478  find an initial guess: closest geometric node to M_xreal
2479  */
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 );
2483 
2484  for ( size_type j = 1; j < M_gm->nbPoints(); ++j )
2485  {
2486  scalar_type d2 = ublas::inner_prod( ublas::column( M_G, j ), M_xreal );
2487 
2488  if ( d2 < d )
2489  {
2490  d = d2;
2491  x0 = M_gm->refNode( j );
2492  y = ublas::column( M_G, j );
2493  }
2494  }
2495 
2496  M_xref = x0;
2497  node_type x0 = barycenterRef();
2498 
2499 #else
2500  M_xref = barycenterRef();
2501 #endif
2502 
2503  dense_matrix_type J( P, P );
2504  dense_vector_type R( P );
2505 
2506  updateResidual( M_xref, R );
2507  updateJacobian( M_xref, J );
2508 
2509  // find xref by solving the non linear equation
2510  M_nlsolver->setType( TRUST_REGION );
2511  M_nlsolver->setRelativeResidualTol( 1e-16 );
2512  M_nlsolver->solve( J, M_xref, R, 1e-10, 10 );
2513 
2514 
2515  // compute the location of xref: inside or outside the element
2516  bool __isin;
2517  double vmin;
2518  this->updateResidual( M_xref, R );
2519  boost::tie( __isin, vmin ) = M_gm->isIn( M_xref );
2520 
2521  if ( __isin &&
2522  ( P == N || ublas::norm_2( R ) < IN_EPS ) )
2523  {
2524  //LOG(INFO) << "point " << M_xref << "in IN (" << vmin << ") residual = " << ublas::norm_2(R) << "\n";
2525  return true;
2526  }
2527 
2528  else
2529  {
2530  //LOG(INFO) << "point " << M_xref << "in OUT (" << vmin << ") residual = " << ublas::norm_2(R) << "\n";
2531  }
2532 
2533  //LOG(INFO) << "done in new nonlinear inverse\n";
2534  return false;
2535 
2536 }
2537 /*
2538  *inversion for non-linear geometric transformations
2539  * (Newton on Grad(pgt)(y - pgt(x)) = 0 )
2540  */
2541 bool nonLinearInverse()
2542 {
2543  const double EPS = 1e-10;
2544  const double IN_EPS = 1e-10;
2545  size_type N = M_xreal.size();
2546  size_type P = M_xref.size();
2547 
2548  /*
2549  find an initial guess: closest geometric node to M_xreal
2550  */
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 );
2554 
2555  for ( size_type j = 1; j < M_gm->nbPoints(); ++j )
2556  {
2557  scalar_type d2 = ublas::inner_prod( ublas::column( M_G, j ), M_xreal );
2558 
2559  if ( d2 < d )
2560  {
2561  d = d2;
2562  x0 = M_gm->refNode( j );
2563  y = ublas::column( M_G, j );
2564  }
2565  }
2566 
2567  M_xref = x0;
2568 
2569  node_type vres( N );
2570  node_type rn( M_xreal );
2571  rn.minus_assign( y );
2572 
2573  this->update();
2574 
2575  // vres = K^T rn
2576  ublas::axpy_prod( rn, M_K, vres );
2577  scalar_type res = ublas::norm_2( vres );
2578 
2579  // std::cerr << "DEBUT: res0=" << res << ", X=" << M_xreal << "\nB=" << M_B << ", K=" << M_K << "\n";
2580  unsigned cnt = 50;
2581 
2582  while ( res > EPS/10 && cnt > 0 )
2583  {
2584  node_type xn( P );
2585 
2586  // xn = B^T rn
2587  ublas::axpy_prod( rn, M_B, xn );
2588 
2589  scalar_type newres;
2590 
2591  for ( int16_type i=1; i<=256; i*=2 )
2592  {
2593  M_xref.plus_assign( xn / scalar_type( i ) );
2594  y = M_gm->transform( M_xref, G() );
2595 
2596  rn = M_xreal - y;
2597 
2598  this->update();
2599 
2600  if ( P != N )
2601  {
2602  // vres = K^T rn
2603  ublas::axpy_prod( rn, M_K, vres );
2604  newres = ublas::norm_2( vres );
2605  }
2606 
2607  else
2608  {
2609  newres = ublas::norm_2( rn );
2610  }
2611 
2612  if ( newres < 1.5*res )
2613  break;
2614  }
2615 
2616  res = newres;
2617  // std::cout << "cnt=" << cnt << ", x=" << M_xref << ", res=" << res << "\n";
2618  --cnt;
2619  }
2620 
2621  // std::cerr << " invert_nonlin done\n";
2622  // std::cerr << "cnt=" << cnt << ", P=" << P << ", N=" << N
2623  // << "\nX=" << M_xreal << " Xref=" << M_xref << "\nresidu=" << res << "\n";
2624  //<< ", G=" << G << "\nX=" << M_xreal << " Xref=" << x << "\nresidu=" << res << "\nB=" << B << ", K=" << K << "\n" << "\n-------------------^^^^^^^^\n";
2625  if ( cnt == 0 )
2626  {
2627 #if 0
2628  std::cerr << "BFGS in geotrans_inv_convex!\n";
2629  GeomapInverseConvex b( *this, M_xreal );
2630 
2631  iteration_ptrtype iter( Iteration<double>::New() );
2632  iter->setMaximumNumberOfIterations( 50 );
2633  iter->setRelativePrecision( 1e-8 );
2634 
2635  node_type x( x0 );
2636  bfgs( b,b,x,10,*iter );
2637  rn = M_gm->transform( x,G() ) - M_xreal;
2638 
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)";
2643 
2644 #endif
2645  }
2646 
2647  bool __isin;
2648  double vmin;
2649  boost::tie( __isin, vmin ) = M_gm->isIn( M_xref );
2650 
2651  if ( __isin &&
2652  ( P == N || ublas::norm_2( rn ) < IN_EPS ) )
2653  {
2654 
2655  return true;
2656  }
2657 
2658  else
2659  {
2660  //LOG(INFO) << "point " << M_xref << "in OUT (" << vmin << ")\n";
2661  }
2662 
2663  return false;
2664 }
2665 
2666 
2667 
2668 private:
2669 
2670 geometric_mapping_ptrtype M_gm;
2671 node_type M_xref;
2672 node_type M_xreal;
2673 
2674 bool M_is_in;
2675 
2676 //matrix_type const& M_G;
2677 matrix_type M_G;
2678 matrix_type M_K;
2679 matrix_type M_B;
2680 matrix_type M_CS;
2681 matrix_type M_g;
2682 
2683 boost::shared_ptr<SolverNonLinear<double> > M_nlsolver;
2684 }; // Inverse
2685 
2686 private:
2687 
2688 element_gm_ptrtype _elementMap;
2689 face_gm_ptrtype _boundaryMap;
2690 matrix_type M_g_linear;
2691 
2692 friend class Inverse;
2693 
2694 reference_convex_type M_refconvex;
2695 
2696  };
2697 
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>
2710 
2711 template<int Dim, int Order,int RealDim, template<uint16_type,uint16_type,uint16_type> class Entity = Simplex, typename T = double>
2712 struct GT_Lagrange
2713 {};
2714 
2715 template<int Dim, int Order, int RealDim, typename T = double>
2716 struct GT_QK
2717 {};
2718 
2719 # /* List of dims. */
2720 # define FEELPP_GEOMAP \
2721  BOOST_PP_TUPLE_TO_LIST( \
2722  1, \
2723  ( \
2724  Lagrange \
2725  ) \
2726  ) \
2727 
2728 # /* List of dims. */
2729 # define FEELPP_DIMS \
2730  BOOST_PP_TUPLE_TO_LIST( \
2731  3, \
2732  ( \
2733  1,2,3 \
2734  ) \
2735  ) \
2736 
2737 # /* List of real dims. */
2738 # define FEELPP_REALDIMS \
2739  BOOST_PP_TUPLE_TO_LIST( \
2740  3, \
2741  ( \
2742  1,2,3 \
2743  ) \
2744  ) \
2745 
2746 # /* List of real dims. */
2747 
2748 # define FEELPP_NEWDIMS \
2749  BOOST_PP_TUPLE_TO_LIST( \
2750  6, \
2751  ( \
2752  (1,1),(1,2),(1,3), \
2753  (2,2),(2,3), \
2754  (3,3) \
2755  ) \
2756  ) \
2757 
2758 # /* List of orders. */
2759 # define FEELPP_ORDERS \
2760  BOOST_PP_TUPLE_TO_LIST( \
2761  5, \
2762  ( \
2763  1,2,3,4,5 \
2764  ) \
2765  ) \
2766 
2767 #
2768 # define FEELPP_ENTITY BOOST_PP_TUPLE_TO_LIST( 2, ( Simplex,Hypercube ) )
2769 
2770 #
2771 # define FEELPP_GEN_GT(GEOM,LDIM,LORDER) \
2772  "GT_" #GEOM "(" #LDIM "," #LORDER ")" \
2773 
2774 #
2775 # /* Generates code for all dim and order. */
2776 # define FEELPP_GT_FACTORY_OP(_, GDO) \
2777  FEELPP_GT_FACTORY GDO \
2778 
2779 #
2780 #
2781 # define FEELPP_GT_DIM(T) BOOST_PP_TUPLE_ELEM(2, 0, T) \
2782 
2783 # define FEELPP_GT_REALDIM(T) BOOST_PP_TUPLE_ELEM(2, 1, T) \
2784 
2785 #
2786 #
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> \
2790  : \
2791  public GeoMap<FEELPP_GT_DIM(LDIMS), LORDER, FEELPP_GT_REALDIM(LDIMS), T, ENTITY, GEOM> \
2792  { \
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; \
2799  \
2800  typedef mpl::vector<boost::none_t, \
2801  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; \
2811  \
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; \
2816  template<int N> \
2817  struct face_gm \
2818  { \
2819  typedef typename mpl::at<geomap_faces_t, mpl::int_<N> >::type type; \
2820  }; \
2821  \
2822  BOOST_PP_CAT(GT_,GEOM)() \
2823  : \
2824  super( boost::shared_ptr<element_gm_type>(new element_gm_type()), boost::shared_ptr<face_gm_type>(new face_gm_type() )) \
2825  {} \
2826  }; \
2827 
2828 #if 0
2829 #define FEELPP_GT_FACTORY(GEOM,LDIM,LORDER,LREALDIM,ENTITY) \
2830  FEELPP_GT_FACTORY_BIS(GEOM,LDIM,LORDER,LREALDIM,ENTITY)), \
2831  BOOST_PP_EMPTY ) \
2832 
2833 #endif
2834 
2835 //BOOST_PP_LIST_FOR_EACH_PRODUCT(FEELPP_GT_FACTORY_OP, 5, (FEELPP_GEOMAP, FEELPP_DIMS, FEELPP_ORDERS, FEELPP_REALDIMS, FEELPP_ENTITY))
2836 BOOST_PP_LIST_FOR_EACH_PRODUCT( FEELPP_GT_FACTORY_OP, 4, ( FEELPP_GEOMAP, FEELPP_NEWDIMS, FEELPP_ORDERS, FEELPP_ENTITY ) )
2837 
2838 #undef FEELPP_DIMS
2839 #undef FEELPP_ORDERS
2840 #undef FEELPP_REALDIMS
2841 #undef FEELPP_NEWDIMS
2842 
2843 
2844 template<typename Elem, template<uint16_type,uint16_type,uint16_type> class Entity = Simplex, typename T = double>
2845 class RealToReference
2846 {
2847 public:
2848  static const uint16_type nDim = Elem::nDim;
2849  static const uint16_type nRealDim = Elem::nRealDim;
2850 
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;
2856 
2857  RealToReference( Elem const& elem )
2858  :
2859  M_gm( new gm_type ),
2860  M_igm( M_gm, elem )
2861  {}
2862 
2863  points_type operator()( points_type const& pts ) const
2864  {
2865  return M_igm( pts );
2866  }
2867 
2868  value_type J() const
2869  {
2870  return M_igm.J();
2871  }
2872 
2873 private:
2874 
2875  gm_ptrtype M_gm;
2876  inverse_gm_type M_igm;
2877 
2878 };
2879 
2880 } // Feel
2881 #endif
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
Definition: glas.hpp:157
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
Definition: glas.hpp:125
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

Generated on Sun Dec 22 2013 13:11:03 for Feel++ by doxygen 1.8.5