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
meshhighorderimpl.hpp
1 /*
2  This file is part of the Feel library
3 
4  Copyright (C) 2007,2008 University of Coimbra
5  Copyright (C) 2010 Université de Grenoble 1
6 
7  This library is free software; you can redistribute it and/or
8  modify it under the terms of the GNU Lesser General Public
9  License as published by the Free Software Foundation; either
10  version 3.0 of the License, or (at your option) any later version.
11 
12  This library is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15  Lesser General Public License for more details.
16 
17  You should have received a copy of the GNU Lesser General Public
18  License along with this library; if not, write to the Free Software
19  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20 */
21 #include <boost/preprocessor/comparison/greater_equal.hpp>
22 #include <feel/feeldiscr/meshhighorder.hpp>
23 
24 namespace Feel
25 {
26 template < class Convex >
27 MeshHighOrder<Convex>::MeshHighOrder( mesh_ptrtype& mesh,
28  std::string projectionStrategy )
29  :
30  old_mesh( mesh ),
31  new_mesh( new new_mesh_type ),
32  strategy ( projectionStrategy )
33 {
34  this->createSwapEdgesMap( mpl::bool_< is_simplex >() );
35 }
36 
37 template < class Convex >
38 MeshHighOrder<Convex>::MeshHighOrder( MeshHighOrder const& tc )
39  :
40  old_mesh( tc.old_mesh ),
41  new_mesh( tc.new_mesh ),
42  strategy( tc.strategy ),
43  M_swap_edges( tc.M_swap_edges )
44 {
45 }
46 
47 template < class Convex >
48 MeshHighOrder<Convex>::~MeshHighOrder()
49 {
50 }
51 
52 
53 template < class Convex >
54 void
55 MeshHighOrder<Convex>::clearMesh()
56 {
57  new_mesh->clear();
58 }
59 
60 template < class Convex >
61 void
62 MeshHighOrder<Convex>::updateP1Mesh ( mesh_ptrtype const& mesh )
63 {
64  old_mesh = mesh;
65 }
66 
67 template < class Convex >
68 typename MeshHighOrder<Convex>::new_mesh_ptrtype
69 MeshHighOrder<Convex>::getMesh() const
70 {
71  return new_mesh;
72 }
73 
74 
75 // affine transformation from [-1,1] to the edge defined by node1 and node2
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
79 {
80  return ( 0.5*( ( 1-x )*node1 + ( x+1 )*node2 ) );
81 }
82 
83 
84 template < class Convex >
85 void
86 MeshHighOrder<Convex>::createSwapEdgesMap ( mpl::bool_<true> )
87 {
88  // Construct a map that returns if for a given element and an edge,
89  // the edge points should be inverted (for triangles only)
90  for ( element_const_iterator elt = old_mesh->beginElement();
91  elt != old_mesh->endElement(); ++elt )
92  {
93  //check if point(0) is elt->point(1)
94  // if they are not the same, points need to be swapped
95  if ( elt->edge( 0 ).point( 0 ).id() != elt->point( 1 ).id() )
96  M_swap_edges[elt->id()][0] = 1;
97 
98  //check if point(0) is elt->point(2)
99  // if they are not the same, points need to be swapped
100  if ( elt->edge( 1 ).point( 0 ).id() != elt->point( 2 ).id() )
101  M_swap_edges[elt->id()][1] = 1;
102 
103  //check if point(0) is elt->point(0)
104  // if they are not the same, points need to be swapped
105  if ( elt->edge( 2 ).point( 0 ).id() != elt->point( 0 ).id() )
106  M_swap_edges[elt->id()][2] = 1;
107  }
108 }
109 
110 template < class Convex >
111 void
112 MeshHighOrder<Convex>::createSwapEdgesMap ( mpl::bool_<false> )
113 {
114  // Construct a map that returns if for a given element and an edge,
115  // the edge points should be inverted (for triangles only)
116  for ( element_const_iterator elt = old_mesh->beginElement();
117  elt != old_mesh->endElement(); ++elt )
118  {
119 #if 0
120 
121  //check if point(0) is elt->point(1)
122  // if they are not the same, points need to be swapped
123  if ( elt->edge( 0 ).point( 0 ).id() != elt->point( 1 ).id() )
124  M_swap_edges[elt->id()][0] = 1;
125 
126  //check if point(0) is elt->point(2)
127  // if they are not the same, points need to be swapped
128  if ( elt->edge( 1 ).point( 0 ).id() != elt->point( 2 ).id() )
129  M_swap_edges[elt->id()][1] = 1;
130 
131  //check if point(0) is elt->point(0)
132  // if they are not the same, points need to be swapped
133  if ( elt->edge( 2 ).point( 0 ).id() != elt->point( 0 ).id() )
134  M_swap_edges[elt->id()][2] = 1;
135 
136 #endif
137  }
138 }
139 
140 
141 
142 template < class Convex >
143 bool
144 MeshHighOrder<Convex>::swapEdge ( size_type eltId, size_type edgeId )
145 {
146  return M_swap_edges[eltId][edgeId];
147 }
148 
149 template < class Convex >
150 void
151 MeshHighOrder<Convex>::updatePts( ublas::vector<double> const& x, points_type const& pts,
152  points_type& final_pts ) const
153 {
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 );
156 }
157 
158 
159 template < class Convex >
160 void
161 MeshHighOrder<Convex>::addVertices( element_type const& elt, new_element_type& new_element,
162  new_mesh_ptrtype& new_mesh, std::vector<bool>& vertexAdd ) const
163 {
164  for ( uint16_type i = 0; i< element_type::numVertices; ++i )
165  {
166  point_type old_point = elt.point( i );
167 
168  point_type new_point( old_point.id(), old_point.node(), old_point.isOnBoundary() );
169  new_point.marker().assign( old_point.marker().value() );
170 
171  // if the created point has not been added to the mesh then add it
172  if ( vertexAdd[ old_point.id() ] == 0 )
173  {
174  new_mesh->addPoint( new_point );
175 
176 #if !defined ( NDEBUG )
177  DVLOG(2) << "[AddPointToMesh] Vertex of id: " << new_point.id() << " has coordinates "
178  << new_point.node();
179 #endif
180 
181  vertexAdd[ new_point.id() ] = 1;
182  }
183 
184  // add point to the element it belongs
185  new_element.setPoint( i, new_mesh->point( new_point.id() ) );
186 
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";
191 #endif
192 
193  }
194 }
195 
196 }
197 
198 
199 
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319

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