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
meshhighorder.hpp
1 /*
2  This file is part of the Feel library
3 
4  Copyright (C) 2007,2008 EPFL
5 
6  This library is free software; you can redistribute it and/or
7  modify it under the terms of the GNU Lesser General Public
8  License as published by the Free Software Foundation; either
9  version 3.0 of the License, or (at your option) any later version.
10 
11  This library is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public
17  License along with this library; if not, write to the Free Software
18  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20 
21 #ifndef __MeshHighOrder
22 #define __MeshHighOrder 1
23 
24 
25 #include <feel/feeldiscr/mesh.hpp>
27 #include <feel/feelvf/vf.hpp>
28 
29 
30 namespace Feel
31 {
36 template< class Convex >
38 {
39  static const bool is_simplex = Convex::is_simplex;
40  static const uint16_type Dim = Convex::nDim;
41  static const uint16_type Order = Convex::nOrder;
42 
43 
44  typedef typename mpl::if_< mpl::bool_< is_simplex >, Simplex<Dim, 1>, Hypercube<Dim, 1> >::type convex_type;
46  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
47 
48 
49  typedef typename mpl::if_< mpl::bool_< is_simplex >, Simplex<2, Order>, Hypercube<2, Order> >::type new_convex_type;
51 
52  typedef boost::shared_ptr<new_mesh_type> new_mesh_ptrtype;
53 
54  typedef typename mesh_type::element_type element_type;
55  typedef typename mesh_type::element_const_iterator element_const_iterator;
56  typedef typename mesh_type::edge_type edge_type;
57  typedef typename mesh_type::point_type point_type;
58 
59  typedef typename element_type::edge_permutation_type edge_permutation_type;
60  typedef typename element_type::node_type node_type;
61 
62  typedef typename new_mesh_type::edge_type new_edge_type;
63  typedef typename new_mesh_type::element_type new_element_type;
64  typedef typename new_mesh_type::element_const_iterator new_element_const_iterator;
65 
66  typedef PointSetEquiSpaced<Hypercube<1,1>, Order, double> oned_pointset_type;
67  typedef PointSetEquiSpaced<convex_type, Order, double> interior_pointset_type;
68 
69  typedef typename oned_pointset_type::points_type node_points_type;
70  typedef typename interior_pointset_type::points_type points_type;
71 
72 public:
73 
75 
76  MeshHighOrder( mesh_ptrtype& mesh,
77  std::string projectionStrategy = "axis" );
78 
79  MeshHighOrder( MeshHighOrder const& tc );
80 
81  ~MeshHighOrder();
82 
83  void clearMesh();
84 
85  void updateP1Mesh ( mesh_ptrtype const& mesh );
86 
87  new_mesh_ptrtype getMesh() const;
88 
89  template< typename elem_type >
90  void generateMesh( std::vector<flag_type> const& flags, std::vector<elem_type> const& polyBoundary )
91  {
92  new_mesh->clear();
93 
94  //Create pointset
95  oned_pointset_type pts_space;
96  node_points_type span = pts_space.pointsBySubEntity( 1,0 );
97 
98  // stores a boolean indicating if the vertex is in the new mesh
99  std::vector<bool> vertexAdd( old_mesh->numVertices() );
100 
101  // stores a boolean indicating if the edge is in the new mesh
102  std::vector<bool> edgeAdd( old_mesh->numEdges() );
103 
104  // connectivity map for edges
105  std::vector<uint16_type> edgesMap ( old_mesh->numEdges() );
106 
107  std::fill( vertexAdd.begin(), vertexAdd.end(), 0 );
108  std::fill( edgeAdd.begin(), edgeAdd.end(), 0 );
109  std::fill( edgesMap.begin(), edgesMap.end(), 0 );
110 
111  uint16_type nodesCount = old_mesh->numVertices();
112 
113  // for in the elements
114  for ( element_const_iterator elt = old_mesh->beginElement();
115  elt != old_mesh->endElement(); ++elt )
116  {
117  new_element_type new_element;
118  new_element.setId( elt->id() );
119  new_element.setOnBoundary( elt->isOnBoundary() );
120 
121  // add vertices to mesh
122 #if !defined ( NDEBUG )
123  DVLOG(2) << "Add vertices...\n";
124 #endif
125  this->addVertices( *elt, new_element, new_mesh, vertexAdd );
126 
127 
128  // add edges and nodes in edges to mesh
129 #if !defined ( NDEBUG )
130  DVLOG(2) << "Add edges...\n";
131 #endif
132 
133  for ( uint16_type i = 0; i< element_type::numEdges; ++i )
134  {
135  edge_type old_edge = elt->edge( i );
136 
137  bool hasEdge = edgeAdd[ old_edge.id() ];
138 
139  new_edge_type new_edge;
140  new_edge.setId( old_edge.id() );
141  new_edge.marker().assign( old_edge.marker().value() );
142  new_edge.setOnBoundary( old_edge.isOnBoundary() );
143 
144  // if the edge is not found in the list, create it and add the points in the boundary
145  if ( hasEdge == 0 )
146  {
147  // add points to edge:
148  for ( uint16_type j = 0; j<2; ++j )
149  {
150  uint16_type localId = j;
151 
152  if ( !swapEdge( elt->id(), i ) )
153  new_edge.setPoint( j, new_mesh->point( old_edge.point( localId ).id() ) );
154 
155  else
156  {
157  localId = 1-j;
158  new_edge.setPoint( j, new_mesh->point( old_edge.point( localId ).id() ) );
159  }
160 
161 #if !defined ( NDEBUG )
162  DVLOG(2) << "[AddPointToEdge] Point localId=" << j
163  << ": globalToMeshId=" << new_mesh->point( old_edge.point( localId ).id() ).id()
164  << "; " << new_mesh->point( old_edge.point( localId ).id() ).node();
165 #endif
166  }
167 
168  edgesMap[ new_edge.id() ] = nodesCount;
169  }
170 
171 
172  node_type node( Dim );
173 
174  for ( uint16_type j=2; j<=Order; ++j )
175  {
176  double t = span( 0,j-2 );
177 
178  if ( swapEdge( elt->id(), i ) )
179  node = this->affineEdge( old_edge.point( 1 ).node(), old_edge.point( 0 ).node(), t );
180 
181  else
182  node = this->affineEdge( old_edge.point( 0 ).node(), old_edge.point( 1 ).node(), t );
183 
184  //strategies for curved edges:
185  if ( old_edge.isOnBoundary() )
186  {
187  this->shiftCoordinates( node, old_edge.point( 0 ).node(), old_edge.point( 1 ).node(),
188  old_edge.marker(), flags, polyBoundary );
189  }
190 
191  point_type new_point( nodesCount, node, old_edge.isOnBoundary() );
192  new_point.marker().assign( old_edge.marker().value() );
193  //uint16_type localId = j;
194 
195  if ( hasEdge == 0 )
196  {
197  new_mesh->addPoint( new_point );
198 
199  new_edge.setPoint( j, new_mesh->point( nodesCount ) );
200 
201 #if !defined ( NDEBUG )
202  DVLOG(2) << "[AddPointToMesh] Point " << j
203  << ": id=" << new_mesh->point( nodesCount ).id()
204  << "; " << new_mesh->point( nodesCount ).node();
205 #endif
206  nodesCount++;
207  }
208 
209  else
210  {
211  if ( !swapEdge( elt->id(), i ) )
212  new_point.setId( edgesMap[ new_edge.id() ] + j - 2 );
213 
214  else
215  new_point.setId( edgesMap[ new_edge.id() ] + Order - j );
216  }
217 
218  new_element.setPoint( element_type::numVertices + ( Order-1 )*i +j-2, new_mesh->point( new_point.id() ) );
219 
220 #if !defined ( NDEBUG )
221  DVLOG(2) << "[AddPointToElement] Add point with local id "
222  << element_type::numVertices + ( Order-1 )*i +j-2
223  << " with global id " << new_mesh->point( new_point.id() ).id()
224  << " to element " << new_element.id() << "\n";
225 #endif
226  }
227 
228  if ( hasEdge == 0 )
229  {
230  new_mesh->addFace( new_edge );
231 
232  edgeAdd[ new_edge.id() ] = 1;
233  }
234 
235 #if !defined ( NDEBUG )
236  DVLOG(2) << "[AddToMesh] Edge of id " << new_edge.id() << " has been added\n";
237  DVLOG(2) << "\n";
238 #endif
239  } // end for in edges
240 
241  // for in the points in the interior of the element
242  if ( Order > ( 1 + is_simplex ) )
243  {
244  interior_pointset_type pts_interior( 1 );
245 
246  points_type pts = pts_interior.points();
247 
248  this->GordonHall( *elt, pts, flags, polyBoundary );
249 
250  // add interior points given by the GordonHall transformation to element and mesh
251  for ( uint16_type i=0; i < pts.size2(); ++i )
252  {
253  // add point to mesh
254  node_type node ( Dim );
255  node[0] = pts( 0,i );
256  node[1] = pts( 1,i );
257 
258  point_type new_point( nodesCount, node, false );
259  new_point.marker().assign( 0 );
260 
261  new_mesh->addPoint( new_point );
262 
263  new_element.setPoint( element_type::numVertices + ( Order-1 )*element_type::numEdges +i,
264  new_mesh->point( new_point.id() ) );
265 
266 #if !defined ( NDEBUG )
267  DVLOG(2) << "Added Point "
268  << element_type::numVertices + ( Order-1 )*element_type::numEdges +i
269  << ": id=" << new_mesh->point( new_point.id() ).id() << "; "
270  << new_mesh->point( new_point.id() ).node() << "\n";
271 #endif
272 
273  nodesCount++;
274  }
275 
276  }
277 
278  new_mesh->addElement( new_element );
279 
280 #if !defined ( NDEBUG )
281  DVLOG(2) << "[AddToMesh] Element of id " << new_element.id() << " has been added\n";
282  DVLOG(2) << "-------------------------------------------------------\n\n";
283 #endif
284 
285  } // end for in elements
286 
287  new_mesh->setNumVertices( old_mesh->numVertices() );
288 
289 #if !defined ( NDEBUG )
290  DVLOG(2) << "Number of elements in the new mesh: " << new_mesh->numElements() << "\n";
291 #endif
292 
293  new_mesh->components().set( MESH_CHECK | MESH_RENUMBER | MESH_UPDATE_EDGES | MESH_UPDATE_FACES );
294  new_mesh->updateForUse();
295  }
296 
297 
298  template< typename elem_type >
299  void GordonHall( element_type const& elt, points_type& pts,
300  std::vector<flag_type> const& flags, std::vector<elem_type> const& polyBoundary )
301  {
302  points_type final_pts = 0*pts;
303 
304  node_type node1 ( Dim );
305  node_type node2 ( Dim );
306  ublas::vector<double> ones( ublas::scalar_vector<double>( pts.size2(), 1.0 ) );
307 
308  ublas::vector<double> xi( pts.size2(), 0.0 );
309  ublas::vector<double> eta( pts.size2(), 0.0 );
310  xi = ublas::row( pts, 0 );
311  eta = ublas::row( pts, 1 );
312 
313  for ( uint16_type i = 0; i< element_type::numEdges; ++i )
314  {
315  std::vector<flag_type>::const_iterator result;
316  result = find( flags.begin(), flags.end(), elt.edge( i ).marker().value() );
317 
318  uint16_type pos = distance( flags.begin(), result );
319 
320  node1 = elt.edge( i ).point( 0 ).node();
321  node2 = elt.edge( i ).point( 1 ).node();
322 
323  if ( this->swapEdge( elt.id(), i ) )
324  std::swap( node1,node2 );
325 
326  bool isAffine = ( result == flags.end() );
327 
328  this->applyDeformation( i, node1, node2,
329  xi, eta, ones, polyBoundary[pos], isAffine,
330  final_pts, mpl::bool_<is_simplex>() );
331  }
332 
333  pts = final_pts;
334  }
335 
336 
337 private:
338 
339  mesh_ptrtype old_mesh;
340 
341  new_mesh_ptrtype new_mesh;
342 
343  std::string strategy;
344 
345  std::map< size_type, std::map< size_type, bool > > M_swap_edges;
346 
347  // affine transformation from [-1,1] to the edge defined by node1 and node2
348  node_type affineEdge( node_type const& node1, node_type const& node2, double const& x ) const;
349 
350  void createSwapEdgesMap ( mpl::bool_<true> );
351 
352  void createSwapEdgesMap ( mpl::bool_<false> );
353 
354  void updatePts( ublas::vector<double> const& x, points_type const& pts, points_type& final_pts ) const;
355 
356  void addVertices( element_type const& elt, new_element_type& new_element,
357  new_mesh_ptrtype& new_mesh, std::vector<bool>& vertexAdd ) const;
358 
359  bool swapEdge ( size_type eltId, size_type edgeId );
360 
361 #if 0
362  template< typename elem_type >
363  class ProjectPoints
364  {
365  public:
366 
367  typedef typename node<double>::type node_t_type;
368  typedef typename matrix_node<double>::type matrix_node_t_type;
369  typedef typename elem_type::functionspace_type::node_type oned_node_type;
370 
371  ProjectPoints( elem_type& _p, node_type& _node, node_type const& node0, node_type const& node1 )
372  :
373  p( _p ),
374  m( ( node1[1] - node0[1] )/( node1[0] - node0[0] ) ),
375  node( _node )
376 
377  {
378  using namespace Feel::vf;
379  //interpolate(p.functionSpace(), gradv(p), der );
380  //der = project( p.functionSpace(), p.functionSpace()->mesh(), gradv(p) );
381  //std::cout << der << "\n";
382 
383  }
384 
385  double operator()( const node_t_type& x ) const
386  {
387  //oned_node_type pt(1);
388  //pt[0] = x[0];
389 
390  //double feval = p(pt)(0,0,0) - node[1] + (1/m)*(pt[0] - node[0]);
391  std::cout << "Eval f( " << x[0] << " ) = " << p( x )( 0,0,0 ) - node[1] + ( 1/m )*( x[0] - node[0] ) << "\n";
392  return p( x )( 0,0,0 ) - node[1] + ( 1/m )*( x[0] - node[0] );
393  }
394 
395  void operator()( const node_t_type& x, node_t_type& gr ) const
396  {
397  typename elem_type::grad_type gradient( p.grad( x ) );
398  double g_v1_x = gradient( 0,0,0 );
399  std::cout << "Gradient " << g_v1_x << "\n";
400 
401  /*oned_node_type pt(1);
402  pt[0] = x[0];*/
403  std::cout << "Eval der( " << x[0] << " ) = " << -x[0]/math::sqrt( 1-x[0]*x[0] ) + 1/m << "\n";
404  gr[0] = -x[0]/math::sqrt( 1-x[0]*x[0] ) + 1/m;
405  }
406 
407  private:
408 
409  elem_type p, der;
410  double m;
411  node_type node;
412 
413  };
414 
415 #endif
416 
417 
418 
419 
420 
421 
422 
423  template< typename elem_type >
424  void shiftCoordinates( node_type& node,
425  node_type const& node0,
426  node_type const& node1,
427  Marker1& marker,
428  std::vector<flag_type> const& flags,
429  std::vector<elem_type> const& polys ) const
430  {
431  std::vector<flag_type>::const_iterator result;
432  result = find( flags.begin(), flags.end(), marker.value() );
433 
434  if ( result != flags.end() )
435  {
436  // determine the position of the flag found in vector flags
437  uint16_type pos = distance( flags.begin(), result );
438 
439  typedef typename elem_type::functionspace_type::node_type oned_node_type;
440 
441  oned_node_type pt( 1 );
442 
443  pt[0] = node[0];
444 
445  elem_type p = polys[pos];
446 
447  if ( strategy == "axis" )
448  {
449 
450 #if !defined ( NDEBUG )
451  DVLOG(2) << "Point (" << node[0] << "," << node[1]
452  << ") moves to (" << node[0] << "," << p( pt )( 0,0,0 )
453  << ")" << "\n";
454 #endif
455 
456  node[1] = p( pt )( 0,0,0 );
457  }
458 
459  else
460  {
461  if ( strategy == "normal" )
462  {
463 #if 0
464  ProjectPoints<elem_type> Projector( p, node, node0, node1 );
465 
466  iteration_ptrtype iter( Iteration<double>::New() );
467  iter->setMaximumNumberOfIterations( 50 );
468  iter->setRelativePrecision( 1e-8 );
469 
470  bfgs( Projector, Projector, pt, 10, *iter );
471 
472  std::cout << "Result: " << pt << "\n";
473 #endif
474 
475  double n0 = node0[0];
476  double n1 = node1[0];
477 
478  if ( n0 > n1 )
479  std::swap( n0,n1 );
480 
481  DVLOG(2) << "Initial interval: [" << n0 << "," << n1 << "]\n";
482  std::pair<double, double> interval = std::make_pair( n0,n1 );
483 
484  double m = ( node1[1] - node0[1] )/( node1[0] - node0[0] );
485  DVLOG(2) << "Slope: " << m << "\n";
486 
487  double error = 1;
488  uint16_type iter = 0;
489 
490  while ( error > 1e-14 && iter<100 )
491  {
492  double midpoint = ( interval.first+interval.second )/2;
493 
494  pt[0] = midpoint;
495  double feval = p( pt )( 0,0,0 ) - node[1] + ( 1/m )*( pt[0] - node[0] );
496 
497  pt[0] = interval.first;
498  double feval0 = p( pt )( 0,0,0 ) - node[1] + ( 1/m )*( pt[0] - node[0] );
499 
500  if ( feval*feval0 > 0 )
501  interval.first = midpoint;
502 
503  else
504  interval.second = midpoint;
505 
506  error = math::abs( feval );
507 
508  DVLOG(2) << "Abs(residual) = " << error << "\n";
509 
510  ++iter;
511  }
512 
513  DVLOG(2) << "Finished in " << iter << " iterations...\n";
514 
515  DVLOG(2) << "Point (" << node[0] << "," << node[1]
516  << ") moves to (" << pt[0] << "," << p( pt )( 0,0,0 )
517  << ")" << "\n";
518 
519  node[1] = p( pt )( 0,0,0 );
520  node[0] = pt[0];
521 
522  }
523  }
524  }
525  }
526 
527 
528  template< typename elem_type >
529  void deformEdge( node_type const& node1, node_type const& node2,
530  elem_type const& p, ublas::vector<double> const& x,
531  points_type& pts, bool const& isAffine = 1 ) const
532  {
533  if ( isAffine )
534  {
535  ublas::vector<double> one( ublas::scalar_vector<double>( x.size(), 1.0 ) );
536 
537  ublas::row( pts, 0 ) = 0.5*( node1[0]*( one - x ) + node2[0]*( one+x ) );
538  ublas::row( pts, 1 ) = 0.5*( node1[1]*( one - x ) + node2[1]*( one+x ) );
539  }
540 
541  else
542  {
543  double a = node1[0];
544  double b = node2[0];
545 
546  typedef typename elem_type::functionspace_type::node_type oned_node_type;
547  oned_node_type pt( 1 );
548 
549  for ( uint16_type i = 0; i < x.size(); i++ )
550  {
551  pt[0] = 0.5*( ( b-a )*x( i ) + b+a );
552  pts( 0,i ) = pt[0];
553  pts( 1,i ) = p( pt )( 0,0,0 );
554  }
555  }
556  }
557 
558  template< typename elem_type >
559  void applyDeformation( uint16_type i, node_type const& node1, node_type const& node2,
560  ublas::vector<double> const& xi, ublas::vector<double> const& eta,
561  ublas::vector<double> const& ones,
562  elem_type const& p, bool const& isAffine,
563  points_type& final_pts, mpl::bool_<true> ) const
564  {
565  points_type nodes3 = points_type( Dim, final_pts.size2() );
566 
567  switch ( i )
568  {
569  case 0:
570  {
571  // calculates nodes3 = \pi_0 ( -\xi )
572  deformEdge( node1, node2, p, -xi, nodes3, isAffine );
573 
574  // calculates ( 1 + 0.5*(\xi+\eta) )*nodes3
575  updatePts( ones + 0.5*( xi+eta ), nodes3, final_pts );
576 
577  deformEdge( node1, node2, p, -( ones+xi+eta ), nodes3, isAffine );
578  updatePts( - 0.5*( ones + xi ), nodes3, final_pts );
579  break;
580  }
581 
582  case 1:
583  {
584  deformEdge( node1, node2, p, -eta, nodes3, isAffine );
585  updatePts( 0.5*( ones-xi ), nodes3, final_pts );
586 
587  deformEdge( node1, node2, p, xi, nodes3, isAffine );
588  updatePts( - 0.5*( ones + eta ), nodes3, final_pts );
589 
590  deformEdge( node1, node2, p, ones, nodes3, isAffine );
591  updatePts( 0.5*( eta+xi ), nodes3, final_pts );
592  break;
593  }
594 
595  default:
596  {
597  deformEdge( node1, node2, p, xi, nodes3, isAffine );
598  updatePts( 0.5*( ones-eta ), nodes3, final_pts );
599 
600  deformEdge( node1, node2, p, -eta, nodes3, isAffine );
601  updatePts( - 0.5*( ones + xi ), nodes3, final_pts );
602 
603  deformEdge( node1, node2, p, ones, nodes3, isAffine );
604  updatePts( 0.5*( ones+xi ), nodes3, final_pts );
605  }
606  }
607  }
608 
609 
610 
611  template< typename elem_type >
612  void applyDeformation( uint16_type i, node_type const& node1, node_type const& node2,
613  ublas::vector<double> const& xi, ublas::vector<double> const& eta,
614  ublas::vector<double> const& ones,
615  elem_type const& p, bool const& isAffine,
616  points_type& final_pts, mpl::bool_<false> ) const
617  {
618  points_type nodes3 = points_type( Dim, final_pts.size2() );
619 
620  switch ( i )
621  {
622  case 0:
623  {
624  // calculates nodes3 = \pi_0 ( xi )
625  deformEdge( node1, node2, p, xi, nodes3, isAffine );
626 
627  // calculates ( 0.5*(1-eta) )*nodes3
628  updatePts( 0.5*( ones-eta ), nodes3, final_pts );
629 
630  // calculates nodes3 = \pi_0 ( 1 )
631  deformEdge( node1, node2, p, ones, nodes3, isAffine );
632 
633  // calculates ( -0.25*( (1+xi)*(1-eta)) )*nodes3
634  updatePts( -0.25*ublas::element_prod( ones+xi, ones - eta ), nodes3, final_pts );
635  }
636 
637  case 1:
638  {
639  // calculates nodes3 = \pi_1 ( eta )
640  deformEdge( node1, node2, p, eta, nodes3, isAffine );
641 
642  // calculates ( 0.5*(1+xi) )*nodes3
643  updatePts( 0.5*( ones+xi ), nodes3, final_pts );
644 
645  // calculates nodes3 = \pi_1 ( 1 )
646  deformEdge( node1, node2, p, ones, nodes3, isAffine );
647 
648  // calculates ( -0.25*( (1+xi)*(1+eta)) )*nodes3
649  updatePts( -0.25*ublas::element_prod( ones+xi, ones + eta ), nodes3, final_pts );
650  }
651 
652  case 2:
653  {
654  // calculates nodes3 = \pi_2 ( -xi )
655  deformEdge( node1, node2, p, -xi, nodes3, isAffine );
656 
657  // calculates ( 0.5*(1+eta) )*nodes3
658  updatePts( 0.5*( ones+eta ), nodes3, final_pts );
659 
660  // calculates nodes3 = \pi_2 ( 1 )
661  deformEdge( node1, node2, p, ones, nodes3, isAffine );
662 
663  // calculates ( -0.25*( (1-xi)*(1+eta)) )*nodes3
664  updatePts( -0.25*ublas::element_prod( ones-xi, ones + eta ), nodes3, final_pts );
665  }
666 
667  default:
668  {
669  // calculates nodes3 = \pi_3 ( -eta )
670  deformEdge( node1, node2, p, -eta, nodes3, isAffine );
671 
672  // calculates ( 0.5*(1-xi) )*nodes3
673  updatePts( 0.5*( ones-xi ), nodes3, final_pts );
674 
675  // calculates nodes3 = \pi_3 ( 1 )
676  deformEdge( node1, node2, p, ones, nodes3, isAffine );
677 
678  // calculates ( -0.25*( (1-xi)*(1-eta)) )*nodes3
679  updatePts( -0.25*ublas::element_prod( ones-xi, ones - eta ), nodes3, final_pts );
680  }
681  }
682  }
683 
684 };
685 }
686 
687 #if !defined(FEELPP_INSTANTIATION_MODE)
688 # include <feel/feeldiscr/meshhighorderimpl.hpp>
689 #endif //
690 
691 #endif // __MeshHighOrder
simplex of dimension Dim
Definition: simplex.hpp:73
brief description
Definition: iteration.hpp:57
Polynomial< Poly, Vectorial > gradient(Polynomial< Poly, Scalar > const &p)
compute the gradient of a scalar polynomial p
Definition: operations.hpp:133
Definition: glas.hpp:125
unifying mesh class
Definition: createsubmesh.hpp:41
Projection made easy.
Definition: projector.hpp:38
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
Definition: meshhighorder.hpp:37

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