21 #ifndef __MeshHighOrder
22 #define __MeshHighOrder 1
36 template<
class Convex >
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;
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;
49 typedef typename mpl::if_< mpl::bool_< is_simplex >,
Simplex<2, Order>, Hypercube<2, Order> >::type new_convex_type;
52 typedef boost::shared_ptr<new_mesh_type> new_mesh_ptrtype;
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;
59 typedef typename element_type::edge_permutation_type edge_permutation_type;
60 typedef typename element_type::node_type node_type;
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;
66 typedef PointSetEquiSpaced<Hypercube<1,1>, Order,
double> oned_pointset_type;
67 typedef PointSetEquiSpaced<convex_type, Order, double> interior_pointset_type;
69 typedef typename oned_pointset_type::points_type node_points_type;
70 typedef typename interior_pointset_type::points_type points_type;
77 std::string projectionStrategy =
"axis" );
85 void updateP1Mesh ( mesh_ptrtype
const& mesh );
87 new_mesh_ptrtype getMesh()
const;
89 template<
typename elem_type >
90 void generateMesh( std::vector<flag_type>
const& flags, std::vector<elem_type>
const& polyBoundary )
95 oned_pointset_type pts_space;
96 node_points_type span = pts_space.pointsBySubEntity( 1,0 );
99 std::vector<bool> vertexAdd( old_mesh->numVertices() );
102 std::vector<bool> edgeAdd( old_mesh->numEdges() );
105 std::vector<uint16_type> edgesMap ( old_mesh->numEdges() );
107 std::fill( vertexAdd.begin(), vertexAdd.end(), 0 );
108 std::fill( edgeAdd.begin(), edgeAdd.end(), 0 );
109 std::fill( edgesMap.begin(), edgesMap.end(), 0 );
111 uint16_type nodesCount = old_mesh->numVertices();
114 for ( element_const_iterator elt = old_mesh->beginElement();
115 elt != old_mesh->endElement(); ++elt )
117 new_element_type new_element;
118 new_element.setId( elt->id() );
119 new_element.setOnBoundary( elt->isOnBoundary() );
122 #if !defined ( NDEBUG )
123 DVLOG(2) <<
"Add vertices...\n";
125 this->addVertices( *elt, new_element, new_mesh, vertexAdd );
129 #if !defined ( NDEBUG )
130 DVLOG(2) <<
"Add edges...\n";
133 for ( uint16_type i = 0; i< element_type::numEdges; ++i )
135 edge_type old_edge = elt->edge( i );
137 bool hasEdge = edgeAdd[ old_edge.id() ];
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() );
148 for ( uint16_type j = 0; j<2; ++j )
150 uint16_type localId = j;
152 if ( !swapEdge( elt->id(), i ) )
153 new_edge.setPoint( j, new_mesh->point( old_edge.point( localId ).id() ) );
158 new_edge.setPoint( j, new_mesh->point( old_edge.point( localId ).id() ) );
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();
168 edgesMap[ new_edge.id() ] = nodesCount;
172 node_type
node( Dim );
174 for ( uint16_type j=2; j<=Order; ++j )
176 double t = span( 0,j-2 );
178 if ( swapEdge( elt->id(), i ) )
179 node = this->affineEdge( old_edge.point( 1 ).node(), old_edge.point( 0 ).node(), t );
182 node = this->affineEdge( old_edge.point( 0 ).node(), old_edge.point( 1 ).node(), t );
185 if ( old_edge.isOnBoundary() )
187 this->shiftCoordinates( node, old_edge.point( 0 ).node(), old_edge.point( 1 ).node(),
188 old_edge.marker(), flags, polyBoundary );
191 point_type new_point( nodesCount, node, old_edge.isOnBoundary() );
192 new_point.marker().assign( old_edge.marker().value() );
197 new_mesh->addPoint( new_point );
199 new_edge.setPoint( j, new_mesh->point( nodesCount ) );
201 #if !defined ( NDEBUG )
202 DVLOG(2) <<
"[AddPointToMesh] Point " << j
203 <<
": id=" << new_mesh->point( nodesCount ).id()
204 <<
"; " << new_mesh->point( nodesCount ).node();
211 if ( !swapEdge( elt->id(), i ) )
212 new_point.setId( edgesMap[ new_edge.id() ] + j - 2 );
215 new_point.setId( edgesMap[ new_edge.id() ] + Order - j );
218 new_element.setPoint( element_type::numVertices + ( Order-1 )*i +j-2, new_mesh->point( new_point.id() ) );
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";
230 new_mesh->addFace( new_edge );
232 edgeAdd[ new_edge.id() ] = 1;
235 #if !defined ( NDEBUG )
236 DVLOG(2) <<
"[AddToMesh] Edge of id " << new_edge.id() <<
" has been added\n";
242 if ( Order > ( 1 + is_simplex ) )
244 interior_pointset_type pts_interior( 1 );
246 points_type pts = pts_interior.points();
248 this->GordonHall( *elt, pts, flags, polyBoundary );
251 for ( uint16_type i=0; i < pts.size2(); ++i )
254 node_type
node ( Dim );
255 node[0] = pts( 0,i );
256 node[1] = pts( 1,i );
258 point_type new_point( nodesCount, node,
false );
259 new_point.marker().assign( 0 );
261 new_mesh->addPoint( new_point );
263 new_element.setPoint( element_type::numVertices + ( Order-1 )*element_type::numEdges +i,
264 new_mesh->point( new_point.id() ) );
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";
278 new_mesh->addElement( new_element );
280 #if !defined ( NDEBUG )
281 DVLOG(2) <<
"[AddToMesh] Element of id " << new_element.id() <<
" has been added\n";
282 DVLOG(2) <<
"-------------------------------------------------------\n\n";
287 new_mesh->setNumVertices( old_mesh->numVertices() );
289 #if !defined ( NDEBUG )
290 DVLOG(2) <<
"Number of elements in the new mesh: " << new_mesh->numElements() <<
"\n";
293 new_mesh->components().set( MESH_CHECK | MESH_RENUMBER | MESH_UPDATE_EDGES | MESH_UPDATE_FACES );
294 new_mesh->updateForUse();
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 )
302 points_type final_pts = 0*pts;
304 node_type node1 ( Dim );
305 node_type node2 ( Dim );
306 ublas::vector<double> ones( ublas::scalar_vector<double>( pts.size2(), 1.0 ) );
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 );
313 for ( uint16_type i = 0; i< element_type::numEdges; ++i )
315 std::vector<flag_type>::const_iterator result;
316 result = find( flags.begin(), flags.end(), elt.edge( i ).marker().value() );
318 uint16_type pos = distance( flags.begin(), result );
320 node1 = elt.edge( i ).point( 0 ).node();
321 node2 = elt.edge( i ).point( 1 ).node();
323 if ( this->swapEdge( elt.id(), i ) )
324 std::swap( node1,node2 );
326 bool isAffine = ( result == flags.end() );
328 this->applyDeformation( i, node1, node2,
329 xi, eta, ones, polyBoundary[pos], isAffine,
330 final_pts, mpl::bool_<is_simplex>() );
339 mesh_ptrtype old_mesh;
341 new_mesh_ptrtype new_mesh;
343 std::string strategy;
345 std::map< size_type, std::map< size_type, bool > > M_swap_edges;
348 node_type affineEdge( node_type
const& node1, node_type
const& node2,
double const& x )
const;
350 void createSwapEdgesMap ( mpl::bool_<true> );
352 void createSwapEdgesMap ( mpl::bool_<false> );
354 void updatePts( ublas::vector<double>
const& x, points_type
const& pts, points_type& final_pts )
const;
356 void addVertices( element_type
const& elt, new_element_type& new_element,
357 new_mesh_ptrtype& new_mesh, std::vector<bool>& vertexAdd )
const;
362 template<
typename elem_type >
368 typedef typename matrix_node<double>::type matrix_node_t_type;
369 typedef typename elem_type::functionspace_type::node_type oned_node_type;
371 ProjectPoints( elem_type& _p, node_type& _node, node_type
const& node0, node_type
const& node1 )
374 m( ( node1[1] - node0[1] )/( node1[0] - node0[0] ) ),
378 using namespace Feel::vf;
385 double operator()(
const node_t_type& x )
const
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] );
395 void operator()(
const node_t_type& x, node_t_type& gr )
const
397 typename elem_type::grad_type
gradient( p.grad( x ) );
399 std::cout <<
"Gradient " << g_v1_x <<
"\n";
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;
423 template<
typename elem_type >
424 void shiftCoordinates( node_type&
node,
425 node_type
const& node0,
426 node_type
const& node1,
428 std::vector<flag_type>
const& flags,
429 std::vector<elem_type>
const& polys )
const
431 std::vector<flag_type>::const_iterator result;
432 result = find( flags.begin(), flags.end(), marker.value() );
434 if ( result != flags.end() )
437 uint16_type pos = distance( flags.begin(), result );
439 typedef typename elem_type::functionspace_type::node_type oned_node_type;
441 oned_node_type pt( 1 );
445 elem_type p = polys[pos];
447 if ( strategy ==
"axis" )
450 #if !defined ( NDEBUG )
451 DVLOG(2) <<
"Point (" << node[0] <<
"," << node[1]
452 <<
") moves to (" << node[0] <<
"," << p( pt )( 0,0,0 )
456 node[1] = p( pt )( 0,0,0 );
461 if ( strategy ==
"normal" )
464 ProjectPoints<elem_type>
Projector( p, node, node0, node1 );
467 iter->setMaximumNumberOfIterations( 50 );
468 iter->setRelativePrecision( 1e-8 );
470 bfgs( Projector, Projector, pt, 10, *iter );
472 std::cout <<
"Result: " << pt <<
"\n";
475 double n0 = node0[0];
476 double n1 = node1[0];
481 DVLOG(2) <<
"Initial interval: [" << n0 <<
"," << n1 <<
"]\n";
482 std::pair<double, double> interval = std::make_pair( n0,n1 );
484 double m = ( node1[1] - node0[1] )/( node1[0] - node0[0] );
485 DVLOG(2) <<
"Slope: " << m <<
"\n";
488 uint16_type iter = 0;
490 while ( error > 1e-14 && iter<100 )
492 double midpoint = ( interval.first+interval.second )/2;
495 double feval = p( pt )( 0,0,0 ) - node[1] + ( 1/m )*( pt[0] - node[0] );
497 pt[0] = interval.first;
498 double feval0 = p( pt )( 0,0,0 ) - node[1] + ( 1/m )*( pt[0] - node[0] );
500 if ( feval*feval0 > 0 )
501 interval.first = midpoint;
504 interval.second = midpoint;
506 error = math::abs( feval );
508 DVLOG(2) <<
"Abs(residual) = " << error <<
"\n";
513 DVLOG(2) <<
"Finished in " << iter <<
" iterations...\n";
515 DVLOG(2) <<
"Point (" << node[0] <<
"," << node[1]
516 <<
") moves to (" << pt[0] <<
"," << p( pt )( 0,0,0 )
519 node[1] = p( pt )( 0,0,0 );
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
535 ublas::vector<double> one( ublas::scalar_vector<double>( x.size(), 1.0 ) );
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 ) );
546 typedef typename elem_type::functionspace_type::node_type oned_node_type;
547 oned_node_type pt( 1 );
549 for ( uint16_type i = 0; i < x.size(); i++ )
551 pt[0] = 0.5*( ( b-a )*x( i ) + b+a );
553 pts( 1,i ) = p( pt )( 0,0,0 );
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
565 points_type nodes3 = points_type( Dim, final_pts.size2() );
572 deformEdge( node1, node2, p, -xi, nodes3, isAffine );
575 updatePts( ones + 0.5*( xi+eta ), nodes3, final_pts );
577 deformEdge( node1, node2, p, -( ones+xi+eta ), nodes3, isAffine );
578 updatePts( - 0.5*( ones + xi ), nodes3, final_pts );
584 deformEdge( node1, node2, p, -eta, nodes3, isAffine );
585 updatePts( 0.5*( ones-xi ), nodes3, final_pts );
587 deformEdge( node1, node2, p, xi, nodes3, isAffine );
588 updatePts( - 0.5*( ones + eta ), nodes3, final_pts );
590 deformEdge( node1, node2, p, ones, nodes3, isAffine );
591 updatePts( 0.5*( eta+xi ), nodes3, final_pts );
597 deformEdge( node1, node2, p, xi, nodes3, isAffine );
598 updatePts( 0.5*( ones-eta ), nodes3, final_pts );
600 deformEdge( node1, node2, p, -eta, nodes3, isAffine );
601 updatePts( - 0.5*( ones + xi ), nodes3, final_pts );
603 deformEdge( node1, node2, p, ones, nodes3, isAffine );
604 updatePts( 0.5*( ones+xi ), nodes3, final_pts );
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
618 points_type nodes3 = points_type( Dim, final_pts.size2() );
625 deformEdge( node1, node2, p, xi, nodes3, isAffine );
628 updatePts( 0.5*( ones-eta ), nodes3, final_pts );
631 deformEdge( node1, node2, p, ones, nodes3, isAffine );
634 updatePts( -0.25*ublas::element_prod( ones+xi, ones - eta ), nodes3, final_pts );
640 deformEdge( node1, node2, p, eta, nodes3, isAffine );
643 updatePts( 0.5*( ones+xi ), nodes3, final_pts );
646 deformEdge( node1, node2, p, ones, nodes3, isAffine );
649 updatePts( -0.25*ublas::element_prod( ones+xi, ones + eta ), nodes3, final_pts );
655 deformEdge( node1, node2, p, -xi, nodes3, isAffine );
658 updatePts( 0.5*( ones+eta ), nodes3, final_pts );
661 deformEdge( node1, node2, p, ones, nodes3, isAffine );
664 updatePts( -0.25*ublas::element_prod( ones-xi, ones + eta ), nodes3, final_pts );
670 deformEdge( node1, node2, p, -eta, nodes3, isAffine );
673 updatePts( 0.5*( ones-xi ), nodes3, final_pts );
676 deformEdge( node1, node2, p, ones, nodes3, isAffine );
679 updatePts( -0.25*ublas::element_prod( ones-xi, ones - eta ), nodes3, final_pts );
687 #if !defined(FEELPP_INSTANTIATION_MODE)
688 # include <feel/feeldiscr/meshhighorderimpl.hpp>
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
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