21 #include <boost/preprocessor/comparison/greater_equal.hpp>
22 #include <feel/feeldiscr/meshhighorder.hpp>
26 template <
class Convex >
27 MeshHighOrder<Convex>::MeshHighOrder( mesh_ptrtype& mesh,
28 std::string projectionStrategy )
31 new_mesh( new new_mesh_type ),
32 strategy ( projectionStrategy )
34 this->createSwapEdgesMap( mpl::bool_< is_simplex >() );
37 template <
class Convex >
38 MeshHighOrder<Convex>::MeshHighOrder( MeshHighOrder
const& tc )
40 old_mesh( tc.old_mesh ),
41 new_mesh( tc.new_mesh ),
42 strategy( tc.strategy ),
43 M_swap_edges( tc.M_swap_edges )
47 template <
class Convex >
48 MeshHighOrder<Convex>::~MeshHighOrder()
53 template <
class Convex >
55 MeshHighOrder<Convex>::clearMesh()
60 template <
class Convex >
62 MeshHighOrder<Convex>::updateP1Mesh ( mesh_ptrtype
const& mesh )
67 template <
class Convex >
68 typename MeshHighOrder<Convex>::new_mesh_ptrtype
69 MeshHighOrder<Convex>::getMesh()
const
76 template <
class Convex >
77 typename MeshHighOrder<Convex>::node_type
78 MeshHighOrder<Convex>::affineEdge( node_type
const& node1, node_type
const& node2,
double const& x )
const
80 return ( 0.5*( ( 1-x )*node1 + ( x+1 )*node2 ) );
84 template <
class Convex >
86 MeshHighOrder<Convex>::createSwapEdgesMap ( mpl::bool_<true> )
90 for ( element_const_iterator elt = old_mesh->beginElement();
91 elt != old_mesh->endElement(); ++elt )
95 if ( elt->edge( 0 ).point( 0 ).id() != elt->point( 1 ).id() )
96 M_swap_edges[elt->id()][0] = 1;
100 if ( elt->edge( 1 ).point( 0 ).id() != elt->point( 2 ).id() )
101 M_swap_edges[elt->id()][1] = 1;
105 if ( elt->edge( 2 ).point( 0 ).id() != elt->point( 0 ).id() )
106 M_swap_edges[elt->id()][2] = 1;
110 template <
class Convex >
112 MeshHighOrder<Convex>::createSwapEdgesMap ( mpl::bool_<false> )
116 for ( element_const_iterator elt = old_mesh->beginElement();
117 elt != old_mesh->endElement(); ++elt )
123 if ( elt->edge( 0 ).point( 0 ).id() != elt->point( 1 ).id() )
124 M_swap_edges[elt->id()][0] = 1;
128 if ( elt->edge( 1 ).point( 0 ).id() != elt->point( 2 ).id() )
129 M_swap_edges[elt->id()][1] = 1;
133 if ( elt->edge( 2 ).point( 0 ).id() != elt->point( 0 ).id() )
134 M_swap_edges[elt->id()][2] = 1;
142 template <
class Convex >
146 return M_swap_edges[eltId][edgeId];
149 template <
class Convex >
151 MeshHighOrder<Convex>::updatePts( ublas::vector<double>
const& x, points_type
const& pts,
152 points_type& final_pts )
const
154 ublas::row( final_pts, 0 ) += element_prod( ublas::row( pts, 0 ), x );
155 ublas::row( final_pts, 1 ) += element_prod( ublas::row( pts, 1 ), x );
159 template <
class Convex >
161 MeshHighOrder<Convex>::addVertices( element_type
const& elt, new_element_type& new_element,
162 new_mesh_ptrtype& new_mesh, std::vector<bool>& vertexAdd )
const
164 for ( uint16_type i = 0; i< element_type::numVertices; ++i )
166 point_type old_point = elt.point( i );
168 point_type new_point( old_point.id(), old_point.node(), old_point.isOnBoundary() );
169 new_point.marker().assign( old_point.marker().value() );
172 if ( vertexAdd[ old_point.id() ] == 0 )
174 new_mesh->addPoint( new_point );
176 #if !defined ( NDEBUG )
177 DVLOG(2) <<
"[AddPointToMesh] Vertex of id: " << new_point.id() <<
" has coordinates "
181 vertexAdd[ new_point.id() ] = 1;
185 new_element.setPoint( i, new_mesh->point( new_point.id() ) );
187 #if !defined ( NDEBUG )
188 DVLOG(2) <<
"[AddVertexElement] Add point: localId=" << i
189 <<
" globalToMeshId=" << new_mesh->point( new_point.id() ).
id()
190 <<
" to element " << new_element.id() <<
"\n";