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
meshmover.hpp
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2007-03-31
7 
8  Copyright (C) 2007 Universit�� Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __MeshMover_H
30 #define __MeshMover_H 1
31 
32 
33 namespace Feel
34 {
41 template<typename MeshType>
42 class MeshMover
43 {
44 public:
45 
46 
50 
51  typedef MeshType mesh_type;
52  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
53 
54  typedef typename mesh_type::value_type value_type;
55 
56  typedef typename mesh_type::gm_type gm_type;
57  typedef boost::shared_ptr<gm_type> gm_ptrtype;
58 
59  typedef typename mesh_type::element_iterator element_iterator;
60  typedef typename mesh_type::element_const_iterator element_const_iterator;
62 
66 
67  MeshMover()
68  {}
69 
70  MeshMover( mesh_ptrtype const& mesh )
71  :
72  M_mesh( mesh )
73  {}
74 
75  MeshMover( MeshMover const & )
76  {}
77 
78  virtual ~MeshMover()
79  {}
80 
82 
86 
87 
89 
93 
94 
96 
100 
101 
103 
107 
108  template<typename DisplType>
109  value_type tryit( DisplType const& u );
110 
111  template<typename DisplType>
112  //boost::tuple<mesh_ptrtype,value_type> apply( mesh_ptrtype const& imesh, DisplType const& u );
113  void apply( mesh_ptrtype& imesh, DisplType const& u );
114 
116 private:
117  mesh_ptrtype M_mesh;
118 };
119 
120 
121 
122 template<typename MeshType>
123 template<typename DisplType>
124 typename MeshMover<MeshType>::value_type
125 MeshMover<MeshType>::tryit( DisplType const& u )
126 {
127  // check that the displacement function is vectorial
128  //BOOST_MPL_ASSERT_MSG( DisplType::is_vectorial, INVALID_DISPLACEMENT_TYPE, DisplType );
129 
130 
131 }
132 template<typename MeshType>
133 template<typename DisplType>
134 //boost::tuple<typename MeshMover<MeshType>::mesh_ptrtype,
135 // typename MeshMover<MeshType>::value_type>
136 void
137 MeshMover<MeshType>::apply( mesh_ptrtype& imesh, DisplType const& u )
138 {
139  // check that the displacement function is vectorial
140  //BOOST_MPL_ASSERT_MSG( DisplType::is_vectorial, INVALID_DISPLACEMENT_TYPE, DisplType );
141 
142  DVLOG(2) << "[Dof::generateDofPoints] generating dof coordinates\n";
143  typedef typename mesh_type::element_type element_type;
144  typedef typename gm_type::template Context<vm::POINT, element_type> gm_context_type;
145  typedef boost::shared_ptr<gm_context_type> gm_context_ptrtype;
146 
147  typedef typename DisplType::functionspace_type::fe_type fe_type;
148  typedef typename fe_type::template Context<vm::POINT, fe_type, gm_type, element_type> fecontext_type;
149 
150 
151  gm_ptrtype gm( new gm_type );
152  fe_type fe;
153 
154  //
155  // Precompute some data in the reference element for
156  // geometric mapping and reference finite element
157  //
158  typename gm_type::precompute_ptrtype __geopc( new typename gm_type::precompute_type( gm, gm->points() ) );
159 
160  //const uint16_type ndofv = fe_type::nDof;
161  //mesh_ptrtype omesh( new mesh_type );
162  //*omesh = *imesh;
163 
164  element_iterator it_elt = imesh->beginElementWithProcessId( imesh->worldComm().localRank() );
165  element_iterator en_elt = imesh->endElementWithProcessId( imesh->worldComm().localRank() );
166  typedef typename DisplType::pc_type pc_type;
167  typedef boost::shared_ptr<pc_type> pc_ptrtype;
168  pc_ptrtype __pc( new pc_type( u.functionSpace()->fe(), gm->points() ) );
169  gm_context_ptrtype __c( new gm_context_type( gm, *it_elt, __geopc ) );
170 
171  typedef typename mesh_type::element_type geoelement_type;
172  typedef typename fe_type::template Context<0, fe_type, gm_type, geoelement_type,gm_context_type::context> fectx_type;
173  typedef boost::shared_ptr<fectx_type> fectx_ptrtype;
174  fectx_ptrtype __ctx( new fectx_type( u.functionSpace()->fe(),
175  __c,
176  __pc ) );
177  typedef typename fectx_type::id_type m_type;
178  typedef boost::multi_array<m_type,1> array_type;
179  array_type uvalues( u.idExtents( *__ctx ) );
180  std::fill( uvalues.data(), uvalues.data()+uvalues.num_elements(), m_type::Zero() );
181  u.id( *__ctx, uvalues );
182 
183  //const uint16_type ndofv = fe_type::nDof;
184 
185 
186  std::map<int,bool> points_done;
187 
188  //if ( !u.areGlobalValuesUpdated() )
189  u.updateGlobalValues();
190 
191 
192  uint16_type nptsperelem = gm->points().size2();
193  ublas::vector<value_type> val( fe_type::nComponents );
194 
195  for ( ; it_elt != en_elt; ++it_elt )
196  {
197  __c->update( *it_elt );
198  __ctx->update( __c );
199  std::fill( uvalues.data(), uvalues.data()+uvalues.num_elements(), m_type::Zero() );
200  u.id( *__ctx, uvalues );
201 
202  for ( uint16_type l =0; l < nptsperelem; ++l )
203  {
204  for ( uint16_type comp = 0; comp < fe_type::nComponents; ++comp )
205  {
206  val[ comp ] = uvalues[l]( comp,0 );
207  }
208 
209  if ( points_done.find( it_elt->point( l ).id() ) == points_done.end() )
210  {
211  //std::cout << "Pt: " << thedof << "Elem " << it_elt->id() << " G=" << it_elt->G() << "\n";
212  imesh->elements().modify( it_elt,
213  lambda::bind( &element_type::applyDisplacement,
214  lambda::_1,
215  l,
216  val ) );
217  points_done[it_elt->point( l ).id()] = true;
218  //std::cout << "Pt: " << thedof << " Moved Elem " << it_elt->id() << " G=" << it_elt->G() << "\n";
219  }
220 
221  else
222  {
223  imesh->elements().modify( it_elt,
224  lambda::bind( &element_type::applyDisplacementG,
225  lambda::_1,
226  l,
227  val ) );
228  }
229  }
230  }
231 
232  imesh->gm()->initCache( imesh.get() );
233  imesh->gm1()->initCache( imesh.get() );
234  // notify observers that the mesh has changed
235  imesh->meshChanged( MESH_CHANGES_POINTS_COORDINATES );
236  //return boost::make_tuple( omesh, 1.0 );
237 
238  imesh->tool_localization()->reset();
239 }
240 
241 template<typename MeshType, typename DisplType>
242 boost::shared_ptr<MeshType>
243 meshMove( boost::shared_ptr<MeshType>& m, DisplType const& u )
244 {
245  MeshMover<MeshType> mover( m );
246  mover.apply( m, u );
247  return m;
248 }
249 } // Feel
250 #endif /* __MeshMover_H */
Expr< Val< typename mpl::if_< boost::is_arithmetic< ExprT1 >, mpl::identity< Cst< ExprT1 > >, mpl::identity< ExprT1 > >::type::type > > val(ExprT1 const &__e1)
precompute expression tensor
Definition: val.hpp:304
Move mesh according to a given map.
Definition: meshmover.hpp:42

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