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
operatorlagrangep1.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: 2008-01-31
7 
8  Copyright (C) 2008-2010 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 __OperatorLagrangeP1_H
30 #define __OperatorLagrangeP1_H 1
31 
32 #include <feel/feeldiscr/functionspace.hpp>
33 #include <feel/feelpoly/fekete.hpp>
36 //#include <feel/feelfilters/exporterquick.hpp>
37 
38 namespace Feel
39 {
40 namespace detailOpLagP1
41 {
46 template<typename SpaceType>
47 struct SpaceToLagrangeP1Space
48 {
49  typedef SpaceType domain_space_type;
50  typedef typename domain_space_type::fe_type::convex_type convex_type;
51  typedef typename domain_space_type::mesh_type domain_mesh_type;
52 
53  typedef typename mpl::if_<mpl::and_<mpl::equal_to<mpl::bool_<convex_type::is_simplex>, mpl::bool_<true> >,
54  mpl::equal_to<mpl::int_<convex_type::nDim>, mpl::int_<1> > >,
55  mpl::identity<Hypercube<convex_type::nDim,convex_type::nOrder,convex_type::nRealDim > >,
56  mpl::identity<convex_type> >::type::type domain_reference_convex_type;
57 
58  template< template<uint16_type Dim> class PsetType>
59  struct SelectConvex
60  {
61  //typedef typename PsetType<convex_type::nDim>::value_type value_type;
62  typedef typename SpaceType::value_type value_type;
63 
64  typedef Lagrange<1,PsetType> type;
65  };
66  typedef typename convex_type::template shape<domain_mesh_type::nDim, 1, domain_mesh_type::nRealDim>::type image_convex_type;
67  typedef Mesh<image_convex_type> image_mesh_type;
68 
69  typedef typename mpl::if_<mpl::bool_<domain_space_type::is_scalar>,
70  mpl::identity<typename SelectConvex<Scalar>::type>,
71  typename mpl::if_<mpl::bool_<domain_space_type::is_vectorial>,
72  mpl::identity<typename SelectConvex<Vectorial>::type>,
73  mpl::identity<typename SelectConvex<Tensor2>::type> >::type>::type::type basis_type;
74 
75 
76  typedef FunctionSpace<image_mesh_type, bases<basis_type> > image_space_type;
77 };
78 
79 }
87 template<typename SpaceType>
89  :
90  //public OperatorInterpolation<SpaceType, typename detail::SpaceToLagrangeP1Space<SpaceType>::image_space_type>
91 public OperatorLinear<SpaceType, typename detailOpLagP1::SpaceToLagrangeP1Space<SpaceType>::image_space_type>
92 {
93  //typedef OperatorInterpolation<SpaceType, typename detail::SpaceToLagrangeP1Space<SpaceType>::image_space_type> super;
95 
96 public:
97 
98 
105  typedef typename super::domain_space_type domain_space_type;
106  typedef typename super::domain_space_ptrtype domain_space_ptrtype;
107  typedef typename domain_space_type::mesh_type domain_mesh_type;
108  typedef typename domain_space_type::mesh_ptrtype domain_mesh_ptrtype;
109  typedef typename domain_space_type::value_type value_type;
110  typedef typename domain_space_type::fe_type domain_fe_type;
111 
112 
113  typedef typename super::backend_ptrtype backend_ptrtype;
117  typedef typename super::dual_image_space_type dual_image_space_type;
118  typedef typename super::dual_image_space_ptrtype dual_image_space_ptrtype;
119  typedef typename dual_image_space_type::mesh_type image_mesh_type;
120  typedef typename dual_image_space_type::mesh_ptrtype image_mesh_ptrtype;
121  typedef typename dual_image_space_type::fe_type image_fe_type;
122 
128  typedef std::vector<std::list<size_type> > el2el_type;
129 
130  typedef typename detailOpLagP1::SpaceToLagrangeP1Space<SpaceType>::convex_type domain_convex_type;
131  typedef typename detailOpLagP1::SpaceToLagrangeP1Space<SpaceType>::domain_reference_convex_type domain_reference_convex_type;
132  typedef typename detailOpLagP1::SpaceToLagrangeP1Space<SpaceType>::image_convex_type image_convex_type;
133  //typedef PointSetWarpBlend<domain_convex_type, domain_space_type::basis_type::nOrder, value_type> pset_type;
134  typedef PointSetFekete<domain_reference_convex_type/*domain_convex_type*/, domain_space_type::basis_type::nOrder, value_type> pset_type;
135  typedef PointSetToMesh<domain_reference_convex_type/*domain_convex_type*/, value_type> p2m_type;
136 
137  typedef typename domain_mesh_type::gm_type gm_type;
138  typedef typename domain_mesh_type::element_type element_type;
139  typedef typename gm_type::template Context<vm::POINT, element_type> gmc_type;
140  typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
141  typedef typename gm_type::precompute_ptrtype gmpc_ptrtype;
142 
143  //typedef typename mpl::at_c<functionspace_vector_type,SpaceIndex>::type functionspace_ptrtype;
144 
145 
147 
151 
156  OperatorLagrangeP1( domain_space_ptrtype const& space,
157  backend_ptrtype const& backend,
158  std::vector<WorldComm> const& worldsComm = std::vector<WorldComm>(1,Environment::worldComm()),
159  std::string pathMeshLagP1=".",
160  std::string prefix="",
161  bool rebuild=true,
162  bool parallelBuild=true );
163 
168  {}
169 
171 
175 
180  {
181  if ( this != &lp1 )
182  {
183  M_el2el = lp1.M_el2el;
184  M_el2pt = lp1.M_el2pt;
185  M_pset = lp1.M_pset;
186  M_gmpc = lp1.M_gmpc;
187  M_p2m = lp1.M_p2m;
188  }
189 
190  return *this;
191  }
192 
193 
195 
199 
200  uint16_type
201  localDof( typename domain_mesh_type::element_const_iterator it,
202  uint16_type localptid ) const
203  {
204  //return localDof( it, localptid, mpl::bool_<is_modal>() );
205  return localDof( it, localptid, mpl::bool_<false>() );
206  }
207 
208  uint16_type
209  localDof( typename domain_mesh_type::element_const_iterator /*it*/,
210  uint16_type localptid,
211  mpl::bool_<false> ) const
212  {
213  return localptid;
214  }
215 
216  uint16_type
217  localDof( typename domain_mesh_type::element_const_iterator it,
218  uint16_type localptid,
219  mpl::bool_<true> ) const;
220 
221 
222  image_mesh_ptrtype
223  mesh()
224  {
225  return M_mesh;
226  }
227 
229 
233 
234 
236 
240 
244  void check() const;
245 
246 
248 
249 private:
251  OperatorLagrangeP1( OperatorLagrangeP1 const & );
252 
253 private:
254  image_mesh_ptrtype M_mesh;
255  el2el_type M_el2el;
256  std::vector<std::vector<size_type> > M_el2pt;
257 
258  pset_type M_pset;
259  gmpc_ptrtype M_gmpc;
260 
261  mutable p2m_type M_p2m;
262 
263 };
264 
265 
266 //
267 // OperatorLagrangeP1
268 //
269 template<typename space_type>
270 OperatorLagrangeP1<space_type>::OperatorLagrangeP1( domain_space_ptrtype const& space,
271  backend_ptrtype const& backend,
272  std::vector<WorldComm> const& worldsComm,
273  std::string pathMeshLagP1,
274  std::string prefix,
275  bool rebuild,
276  bool parallelBuild )
277  :
278  super( space,
279  dual_image_space_ptrtype( new dual_image_space_type( image_mesh_ptrtype( new image_mesh_type ) ) ),
280  backend,
281  false ),
282  M_mesh( new image_mesh_type(worldsComm[0]) ),
283  M_el2el(),
284  M_el2pt(),
285  M_pset( 0 ),
286  M_gmpc( new typename gm_type::precompute_type( this->domainSpace()->mesh()->gm(),
287  M_pset.points() ) ),
288  M_p2m()
289 {
290 
291  std::string nameMeshLagP1base = "meshLagP1-"+space->basisName()+(boost::format("-order%1%")%domain_space_type::basis_type::nOrder).str();
292  std::string nameMeshLagP1 = prefixvm(prefix,nameMeshLagP1base);
293  if ( rebuild )
294  {
295  if (M_mesh->worldComm().globalRank() == M_mesh->worldComm().masterRank() )
296  {
297  // create a triangulation of the current convex
298  // using equispaced points defined on the reference
299  // element
300  //M_p2m.addBoundaryPoints( M_pset.points() );
301  M_p2m.visit( &M_pset );
302 
303  // #if !defined( NDEBUG )
304  // ExporterQuick<image_mesh_type> exp( "vtk", "ensight" );
305  // exp.save( M_p2m.mesh() );
306  // #endif
307 
308  // do not renumber the mesh entities
309  M_p2m.mesh()->components().clear ( MESH_RENUMBER );
310  M_p2m.mesh()->components().clear ( MESH_CHECK );
311  M_p2m.mesh()->updateForUse();
312 
313  M_p2m.mesh()->save( _name=nameMeshLagP1,_path=pathMeshLagP1 );
314  }
315 
316  M_mesh->worldComm().barrier();
317 
318  if (M_mesh->worldComm().globalRank() != M_mesh->worldComm().masterRank() )
319  {
320  M_p2m.mesh()->load( _name=nameMeshLagP1,_path=pathMeshLagP1,
321  _update=MESH_UPDATE_EDGES|MESH_UPDATE_FACES );
322  }
323  }
324  else
325  {
326  M_p2m.mesh()->load( _name=nameMeshLagP1,_path=pathMeshLagP1,
327  _update=MESH_UPDATE_EDGES|MESH_UPDATE_FACES );
328  }
329 
330 
331  VLOG(2) << "[P1 Lagrange] Pointset " << M_pset.points() << "\n";
332 
333  typedef typename image_mesh_type::point_type point_type;
334  typedef typename image_mesh_type::element_type element_type;
335  typedef typename image_mesh_type::face_type face_type;
336 
337 
338 
339  // inherit the table of markersName
340  BOOST_FOREACH( auto itMark, this->domainSpace()->mesh()->markerNames() )
341  {
342  M_mesh->addMarkerName( itMark.first,itMark.second[0],itMark.second[1] );
343  }
344 
345 
346  bool doParallelBuild = (M_mesh->worldComm().size()==1)?false:parallelBuild;
347 
348  //std::vector<bool> pts_done( this->domainSpace()->nLocalDof()/domain_space_type::nComponents, false );
349  //size_type ne = this->domainSpace()->mesh()->numElements() * M_p2m.mesh()->numElements();
350  //std::vector<boost::tuple<size_type, uint16_type, size_type> > dofindices( image_mesh_type::element_type::numVertices*ne, boost::make_tuple( 0,0,0 ) );
351  //std::vector<bool> eltid_done( this->domainSpace()->nLocalDof(), false );
352 
353  //typename domain_mesh_type::element_const_iterator it = this->domainSpace()->mesh()->beginElement();
354  //typename domain_mesh_type::element_const_iterator en = this->domainSpace()->mesh()->endElement();
355  auto it = this->domainSpace()->mesh()->beginElementWithProcessId(this->domainSpace()->mesh()->worldComm().localRank());
356  auto en = this->domainSpace()->mesh()->endElementWithProcessId(this->domainSpace()->mesh()->worldComm().localRank());
357 
358  // memory nodes
359  std::vector<size_type> new_node_numbers( this->domainSpace()->nLocalDof() );
360  std::fill ( new_node_numbers.begin(),new_node_numbers.end(),invalid_size_type_value );
361  // the number of nodes on the new mesh, will be incremented
362  size_type n_new_nodes = 0, n_new_elem = 0, n_new_faces = 0;
363 
364  std::map<size_type,size_type> mapGhostDofIdClusterToProcess;
365  std::map<size_type,std::vector<size_type> > mapActiveEltWhichAreGhostInOtherPartition;
366  // use the geometric transformation to transform
367  // the local equispace mesh to a submesh of the
368  // current element
369  gmc_ptrtype gmc( new gmc_type( this->domainSpace()->mesh()->gm(), *it, M_gmpc ) );
370 
371  for ( size_type elid = 0, pt_image_id = 0; it != en; ++it )
372  {
373  DVLOG(2) << "=========================================\n";
374  DVLOG(2) << "global element " << it->id() << " oriented ok ? : " << it->isAnticlockwiseOriented() << "\n";
375  DVLOG(2) << "global element G=" << it->G() << "\n";
376 
377  gmc->update( *it );
378 
379  // accumulate the local mesh in element *it in the new mesh
380  auto itl = M_p2m.mesh()->beginElement();
381  auto enl = M_p2m.mesh()->endElement();
382  if (doParallelBuild && it->numberOfPartitions() >1)
383  mapActiveEltWhichAreGhostInOtherPartition[it->id()] = std::vector<size_type>(std::distance(itl,enl),invalid_size_type_value);
384 
385  for ( int cptEltp2m=0 ; itl != enl; ++itl, ++elid,++cptEltp2m )
386  {
387  DVLOG(2) << "************************************\n";
388  DVLOG(2) << "local elt = " << elid << "\n";
389  DVLOG(2) << "local element " << itl->id() << " oriented ok ? : " << itl->isAnticlockwiseOriented() << "\n";
390  DVLOG(2) << "local element G=" << itl->G() << "\n";
391 
392  element_type elt;
393  elt.setId( elid );
394  elt.setMarker( it->marker().value() );
395  elt.setMarker2( it->marker2().value() );
396  elt.setMarker3( it->marker3().value() );
397  elt.setProcessIdInPartition( it->pidInPartition() );
398  elt.setNumberOfPartitions( it->numberOfPartitions() );
399  elt.setProcessId(it->processId());
400  elt.setNeighborPartitionIds( it->neighborPartitionIds() );
401 
402  // accumulate the points
403  for ( int p = 0; p < image_mesh_type::element_type::numVertices; ++p )
404  {
405  DVLOG(2) << "local In original element, vertex number " << itl->point( p ).id() << "\n";
406  uint16_type localptid = itl->point( p ).id();
407  uint16_type localptid_dof = localDof( it, localptid );
408  size_type ptid = boost::get<0>( this->domainSpace()->dof()->localToGlobal( it->id(),
409  localptid_dof, 0 ) );
410  FEELPP_ASSERT( ptid < this->domainSpace()->nLocalDof()/domain_space_type::nComponents )( ptid )( this->domainSpace()->nLocalDof()/domain_space_type::nComponents ).warn( "invalid domain dof index" );
411 
412  if (doParallelBuild && !this->domainSpace()->dof()->dofGlobalClusterIsOnProc(this->domainSpace()->dof()->mapGlobalProcessToGlobalCluster(ptid)))
413  mapGhostDofIdClusterToProcess[this->domainSpace()->dof()->mapGlobalProcessToGlobalCluster(ptid)] = ptid;
414 #if 0
415  //if ( pts_done[ptid] == false )
416  {
417  dofindices[pt_image_id] = boost::make_tuple( elid, p, ptid );
418  pt_image_id++;
419  //pts_done[ptid] = pt_image_id;
420  }
421 #endif
422 
423  // add new node if not inserted before
424  if ( new_node_numbers[ptid] == invalid_size_type_value )
425  {
426  new_node_numbers[ptid] = n_new_nodes;
427 
428  point_type __pt( n_new_nodes, boost::get<0>( this->domainSpace()->dof()->dofPoint( ptid ) ) );
429  __pt.setProcessId( it->processId() );
430 
431 
432  DVLOG(2) << "[OperatorLagrangeP1] element id "
433  << elid << "\n";
434  DVLOG(2) << "[OperatorLagrangeP1] local point id "
435  << localptid << " coord " << itl->point( p ).node() << "\n";
436  DVLOG(2) << "[OperatorLagrangeP1] point local id "
437  << localptid << " global id " << ptid << "\n"
438  << " localptid_dof = " << localptid_dof << "\n";
439  DVLOG(2) << "[OperatorLagrangeP1] adding point "
440  << __pt.id() << " : " << __pt.node() << "\n";
441  DVLOG(2) << "[OperatorLagrangeP1] domain point gid "
442  << boost::get<1>( this->domainSpace()->dof()->dofPoint( ptid ) )
443  << " coords: " << boost::get<0>( this->domainSpace()->dof()->dofPoint( ptid ) ) << "\n";
444 
445  DLOG_IF( WARNING, ( ublas::norm_2( boost::get<0>( this->domainSpace()->dof()->dofPoint( ptid ) )-__pt.node() ) > 1e-10 ) )
446  << "inconsistent point coordinates at pt index " << ptid << " in element " << elid << " with "
447  << " dot pt: " << boost::get<0>( this->domainSpace()->dof()->dofPoint( ptid ) )
448  << " mesh pt: " << __pt.node()
449  << " itl->id: " << itl->id()
450  << " local pt id : " << localptid
451  << " local pt id dof : " << localptid_dof << "\n";
452 
453  M_mesh->addPoint( __pt );
454 
455  n_new_nodes++;
456  }
457  //--------------------------------------------------------------------//
458 
459  elt.setPoint( p, M_mesh->point( new_node_numbers[ptid] ) );
460 
461 #if 1
462  FEELPP_ASSERT( ublas::norm_2( boost::get<0>( this->domainSpace()->dof()->dofPoint( ptid ) )-elt.point( p ).node() ) < 1e-10 )
463  ( boost::get<0>( this->domainSpace()->dof()->dofPoint( ptid ) ) )
464  ( p )
465  ( elt.point( p ).node() )
466  ( elid )
467  ( itl->id() )
468  ( localptid )
469  ( localptid_dof )
470  ( ptid ).warn( "[after] inconsistent point coordinates" );
471 
472  FEELPP_ASSERT( ublas::norm_2( elt.point( p ).node()-ublas::column( elt.G(), p ) ) < 1e-10 )
473  ( p )
474  ( elt.point( p ).node() )
475  ( elid )
476  ( elt.G() ).warn( "[after] inconsistent point coordinates" );
477 #endif
478  }
479 
480 #if 0
481  static const int nDim = element_type::nDim;
482  ublas::vector<double> D( nDim + 1, 0 );
483  typename matrix_node<value_type>::type M( nDim+1,nDim+1 );
484  typename matrix_node<value_type>::type P( nDim,nDim+1 );
485  typename matrix_node<value_type>::type P0( nDim,nDim+1 );
486  std::fill( M.data().begin(), M.data().end(), value_type( 1.0 ) );
487 
488  for ( int p = 0; p < element_type::numVertices; ++p )
489  {
490  ublas::column( P0, p ) = elt.vertex( p );
491  }
492 
493  ublas::subrange( M, 0, nDim, 0, nDim+1 ) = P0;
494  double meas_times = details::det( M, mpl::int_<nDim+1>() );
495  FEELPP_ASSERT( meas_times > 0 )( elt.id() )( elt.G() ).warn( "negative area" );
496 #else
497  //FEELPP_ASSERT( elt.isAnticlockwiseOriented() )( elt.id() )( elt.G() ).warn( "negative area" );
498  //FEELPP_ASSERT( ublas::norm_inf( elt.G()- it->G() ) < 1e-10 )( it->G() )( elt.G() ).warn( "global: not same element" );
499  //FEELPP_ASSERT( ublas::norm_inf( elt.G()- itl->G() ) < 1e-10 )( itl->G() )( elt.G() ).warn( "local: not same element" );
500 #endif
501 
502  // set id of element
503  elt.setId ( n_new_elem );
504  // increment the new element counter
505  n_new_elem++;
506  // add element in mesh
507  auto const& theNewElt = M_mesh->addElement ( elt );
508  if (doParallelBuild && it->numberOfPartitions() >1) mapActiveEltWhichAreGhostInOtherPartition[it->id()][cptEltp2m] = theNewElt.id();
509 
510  // Maybe add faces for this element
511  for ( unsigned int s=0; s<itl->numTopologicalFaces; s++ )
512  {
513  if ( !itl->facePtr( s ) ) continue;
514 
515  auto const& theFaceBase = itl->face( s );
516  //size_type global_face_id = theFaceBase.id();
517 
518  face_type new_face;// = theFaceBase;
519  new_face.disconnect();
520  new_face.setOnBoundary( true );
521  new_face.setProcessId( it->processId() );
522  new_face.setNumberOfPartitions( 1 );
523  //new_face.setMarker( theFaceBase.marker().value() );
524  //new_face.setMarker2( theFaceBase.marker2().value() );
525  //new_face.setMarker3( theFaceBase.marker3().value() );
526 
527  for ( uint16_type p = 0; p < theFaceBase.nPoints(); ++p )
528  {
529  uint16_type localptidFace = theFaceBase.point( p ).id();
530  uint16_type localptidFace_dof = localDof( it, localptidFace );
531  const size_type theglobdof = this->domainSpace()->dof()->localToGlobal( it->id(),localptidFace_dof,0 ).template get<0>();
532  new_face.setPoint( p, M_mesh->point( new_node_numbers[theglobdof] ) );
533  }
534  // add it to the list of faces
535  auto addFaceRes = M_mesh->addFace( new_face );
536  }
537 
538  } // for ( int cptEltp2m=0 ; itl != enl; ++itl, ++elid,++cptEltp2m )
539 
540  } // for ( size_type elid = 0, pt_image_id = 0; it != en; ++it )
541 
542 
543 
544  if ( doParallelBuild )
545  {
546  VLOG(2) << "[P1 Lagrange] start parallel build \n";
547 
548  auto const theWorldCommSize = M_mesh->worldComm().size();
549  std::vector<int> nbMsgToSend( theWorldCommSize , 0 );
550  std::vector< std::map<int,size_type> > mapMsg( theWorldCommSize );
551 
552  auto iv = this->domainSpace()->mesh()->beginGhostElement();
553  auto const en = this->domainSpace()->mesh()->endGhostElement();
554  for ( ; iv != en; ++iv )
555  {
556  auto const procGhost = iv->processId();
557  auto const idEltOnGhost = iv->idInOthersPartitions(procGhost);
558  M_mesh->worldComm().localComm().send(procGhost, nbMsgToSend[procGhost], idEltOnGhost);
559  mapMsg[procGhost].insert( std::make_pair( nbMsgToSend[procGhost],iv->id() ) );
560  ++nbMsgToSend[procGhost];
561  }
562 
563  // counter of msg received for each process
564  std::vector<int> nbMsgToRecv;
565  mpi::all_to_all( M_mesh->worldComm().localComm(),
566  nbMsgToSend,
567  nbMsgToRecv );
568 
569  VLOG(2) << "[P1 Lagrange] parallel build : finish first send \n";
570 
571  // recv dof asked and re-send dof in this proc
572  for ( int proc=0; proc<theWorldCommSize; ++proc )
573  {
574  for ( int cpt=0; cpt<nbMsgToRecv[proc]; ++cpt )
575  {
576  //recv
577  size_type idEltRecv;
578  M_mesh->worldComm().localComm().recv( proc, cpt, idEltRecv );
579  auto const& theeltIt = this->domainSpace()->mesh()->elementIterator(idEltRecv);
580 
581  //if (mapActiveEltWhichAreGhostInOtherPartition.find(idEltRecv) == mapActiveEltWhichAreGhostInOtherPartition.end() ) std::cout << "BUG" << std::endl;
582  auto const& idOfNewElt = mapActiveEltWhichAreGhostInOtherPartition[idEltRecv];
583 
584  std::vector<boost::tuple<size_type,ublas::vector<double> > > resultClusterDofsAndNodesToSend(M_p2m.mesh()->numPoints());
585 
586  auto itp = M_p2m.mesh()->beginPoint();
587  auto const enp = M_p2m.mesh()->endPoint();
588  for ( ; itp!=enp ; ++itp )
589  {
590  uint16_type localptid = itp->id();// if (itp->id() >=M_p2m.mesh()->numElements()) std::cout << "aie " << itp->id() << "\n" << std::endl;
591  uint16_type localptid_dof = localDof( theeltIt, localptid );
592  size_type ptid = boost::get<0>( this->domainSpace()->dof()->localToGlobal( theeltIt->id(),localptid_dof, 0 ) );
593  size_type idInProcAsked = this->domainSpace()->dof()->mapGlobalProcessToGlobalCluster(ptid);
594  ublas::vector<double> thedofnode = boost::get<0>( this->domainSpace()->dof()->dofPoint( ptid ) );
595  resultClusterDofsAndNodesToSend[localptid/*itp->id()*/] = boost::make_tuple( idInProcAsked, thedofnode);
596  }
597 
598  auto itl = M_p2m.mesh()->beginElement();
599  auto const enl = M_p2m.mesh()->endElement();
600  // localIdNewElt -> ( localIdPt, globalIdNewElt )
601  std::vector< boost::tuple< std::vector<uint16_type>, size_type > > resultToSendBis(std::distance(itl,enl));
602  for ( size_type elidp2m = 0 ; itl != enl; ++itl, ++elidp2m )
603  {
604  std::vector<uint16_type> vecLocalIdPtToSend(image_mesh_type::element_type::numVertices);
605  // accumulate the points
606  for ( int p = 0; p < image_mesh_type::element_type::numVertices; ++p )
607  {
608  DVLOG(2) << "local In original element, vertex number " << itl->point( p ).id() << "\n";
609  const uint16_type localptid = itl->point( p ).id();
610  vecLocalIdPtToSend[p]=localptid;
611  }
612  resultToSendBis[elidp2m] = boost::make_tuple(vecLocalIdPtToSend,idOfNewElt[elidp2m] );
613  }
614  auto resultToSend = boost::make_tuple(resultToSendBis,resultClusterDofsAndNodesToSend);
615  M_mesh->worldComm().localComm().send( proc, cpt, resultToSend);
616  } // for ( int cpt=0; cpt<nbMsgToRecv[proc]; ++cpt )
617  } // for ( int proc=0; proc<theWorldCommSize; ++proc )
618 
619  VLOG(2) << "[P1 Lagrange] parallel build : finish first recv and resend response \n";
620 
621  std::map<size_type,size_type> new_ghost_node_numbers;
622 
623  // get response to initial request
624  for ( int proc=0; proc<theWorldCommSize; ++proc )
625  {
626  for ( int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
627  {
628  boost::tuple< std::vector< boost::tuple< std::vector<uint16_type>, size_type > >,
629  std::vector<boost::tuple<size_type,ublas::vector<double> > > > resultRecvData;
630  M_mesh->worldComm().localComm().recv( proc, cpt, resultRecvData );
631 
632  auto const& myGhostEltBase = this->domainSpace()->mesh()->element(mapMsg[proc][cpt],proc);
633 
634  auto const& mapLocalToGlobalPointId = resultRecvData.template get<1>();
635  auto resultRecv = resultRecvData.template get<0>();
636 
637  auto itelt = resultRecv.begin();
638  auto const enelt = resultRecv.end();
639  for ( ;itelt!=enelt;++itelt )
640  {
641  element_type elt;
642  elt.setMarker( myGhostEltBase.marker().value() );
643  elt.setMarker2( myGhostEltBase.marker2().value() );
644  elt.setMarker3( myGhostEltBase.marker3().value() );
645  elt.setProcessIdInPartition( myGhostEltBase.pidInPartition() );
646  elt.setNumberOfPartitions(2/*it->numberOfPartitions()*/);
647  elt.setProcessId(proc/*it->processId()*/);
648  elt.setNeighborPartitionIds( std::vector<int>({proc}) );
649  bool findConnection=false;
650 
651  auto itpt = itelt->template get<0>().begin();
652  auto const enpt = itelt->template get<0>().end();
653  for ( int p = 0 ; itpt!=enpt ; ++itpt,++p )
654  {
655  //auto const globclusterdofRecv = itpt->template get<0>();
656  auto const globclusterdofRecv = mapLocalToGlobalPointId[(int)*itpt].template get<0>();
657  size_type theptid = invalid_size_type_value;
658  if ( this->domainSpace()->dof()->dofGlobalClusterIsOnProc(globclusterdofRecv))
659  {
660  theptid = this->domainSpace()->dof()->mapGlobalClusterToGlobalProcess( globclusterdofRecv - this->domainSpace()->dof()->firstDofGlobalCluster() );
661  findConnection = true;
662  }
663  else if (mapGhostDofIdClusterToProcess.find(globclusterdofRecv) != mapGhostDofIdClusterToProcess.end() )
664  {
665  theptid = mapGhostDofIdClusterToProcess[globclusterdofRecv];
666  findConnection = true;
667  }
668 
669  if ( theptid != invalid_size_type_value )
670  elt.setPoint( p, M_mesh->point( new_node_numbers[theptid] ) );
671  else
672  {
673  if (new_ghost_node_numbers.find(globclusterdofRecv) == new_ghost_node_numbers.end() )
674  {
675  new_ghost_node_numbers[globclusterdofRecv] = n_new_nodes;
676 
677  point_type __pt( n_new_nodes, mapLocalToGlobalPointId[(int)*itpt].template get<1>()/*itpt->template get<1>()*/ );
678  __pt.setProcessId( invalid_uint16_type_value );
679 
680  auto const& theNewPt = M_mesh->addPoint( __pt );
681  n_new_nodes++;
682  }
683 
684  elt.setPoint( p, M_mesh->point( new_ghost_node_numbers[globclusterdofRecv] ) );
685  }
686 
687  } // for ( int p = 0 ; itpt!=enpt ; ++itpt,++p )
688 
689  if (!findConnection) continue;
690 
691  // set id of element
692  elt.setId ( n_new_elem );
693 
694  // increment the new element counter
695  n_new_elem++;
696 
697  auto const& theNewElt = M_mesh->addElement ( elt );
698  M_mesh->elements().modify( M_mesh->elementIterator( theNewElt.id(), proc ),
699  Feel::detail::updateIdInOthersPartitions( proc, itelt->template get<1>() ) );
700 
701  } // for ( ;itelt!=enelt;++itelt )
702  } // for ( int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
703  } // for ( int proc=0; proc<theWorldCommSize; ++proc )
704 
705  } // if ( parallelBuild )
706 
707  DVLOG(2) << "[P1 Lagrange] Number of points in mesh: " << M_mesh->numPoints() << "\n";
708 
709  M_mesh->setNumVertices( M_mesh->numPoints() );
710  M_mesh->components().reset();
711  M_mesh->components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
712  M_mesh->updateForUse();
713 
714  this->check();
715 
716  // construct the p1 space and set the operator
717 #if defined(FEELPP_ENABLE_MPI_MODE)
718  auto Xh_image = dual_image_space_type::New(_mesh=M_mesh,_worldscomm=worldsComm);
719 #else
720  auto Xh_image = dual_image_space_type::New(_mesh=M_mesh);
721 #endif
722  this->init( this->domainSpace(),
723  Xh_image,
724  backend,
725  false );
726 
727  this->check();
728 
729  //
730  // Update generates the matrix associated with the interpolation operator
731  // comment out for now
732 #if 0
733  this->update();
734 #endif // 0
735 
736 
737 #if 0
738 
739  for ( size_type i = 0; i < this->domainSpace()->nLocalDof()/domain_space_type::nComponents; ++i )
740  {
741  FEELPP_ASSERT( boost::get<1>( this->domainSpace()->dof()->dofPoint( i ) ) ==
742  boost::get<1>( this->dualImageSpace()->dof()->dofPoint( i ) ) )
743  ( boost::get<0>( this->domainSpace()->dof()->dofPoint( i ) ) )
744  ( boost::get<0>( this->dualImageSpace()->dof()->dofPoint( i ) ) )
745  ( i ).warn( "check inconsistent point id" );
746  FEELPP_ASSERT( ublas::norm_2( boost::get<0>( this->domainSpace()->dof()->dofPoint( i ) )-
747  boost::get<0>( this->dualImageSpace()->dof()->dofPoint( i ) )
748  ) < 1e-10 )
749  ( boost::get<0>( this->domainSpace()->dof()->dofPoint( i ) ) )
750  ( boost::get<0>( this->dualImageSpace()->dof()->dofPoint( i ) ) )
751  ( i ).warn( "check inconsistent point coordinates" );
752  }
753 
754 #endif
755 
756 }
757 template<typename space_type>
758 void
760 {
761 #if !defined(NDEBUG)
762  for ( size_type i = 0; i < this->domainSpace()->nLocalDof()/domain_space_type::nComponents; ++i )
763  {
764  // this test is not good: need to investigate
765 #if 0
766  DLOG_IF( WARNING, (ublas::norm_2( boost::get<0>( this->domainSpace()->dof()->dofPoint( i ) )- M_mesh->point( i ).node() ) > 1e-10 ) )
767  << "inconsistent point coordinates at index " << i
768  << "dofpt : " << boost::get<0>( this->domainSpace()->dof()->dofPoint( i ) )
769  << "mesh pt : " << M_mesh->point( i ).node() << "\n";
770 #endif
771  }
772 #endif
773 
774 }
775 #if 0
776 
777 template<typename space_type>
779 OperatorLagrangeP1<space_type>::operator()( element_type const& u ) const
780 {
781  typename image_space_type::element_type res( this->dualImageSpace(), u.name() );
782  //res.resize( M_mesh->numPoints() );
783  res = ublas::scalar_vector<value_type>( res.size(), value_type( 0 ) );
784  // construct the mesh
785  typename domain_mesh_type::element_const_iterator it = this->domainSpace()->mesh()->beginElement();
786  typename domain_mesh_type::element_const_iterator en = this->domainSpace()->mesh()->endElement();
787 
788  for ( ; it != en; ++it )
789  {
790  gmc_ptrtype gmc( new gmc_type( this->domainSpace()->mesh()->gm(), *it, M_gmpc ) );
791  typename matrix_node<value_type>::type u_at_pts = u.id( *gmc );
792 
793  typename std::list<size_type>::const_iterator ite = M_el2el[it->id()].begin();
794  typename std::list<size_type>::const_iterator ene = M_el2el[it->id()].end();
795 
796  for ( ; ite != ene; ++ite )
797  {
798  typename domain_mesh_type::element_type const& elt = M_mesh->element( *ite );
799 
800  for ( int c = 0; c < nComponents; ++c )
801  for ( int p = 0; p < domain_mesh_type::element_type::numVertices; ++p )
802  {
803  size_type ptid = boost::get<0>( this->dualImageSpace()->dof()->localToGlobal( elt.id(), p, c ) );
804  res( ptid ) = ublas::column( u_at_pts, M_el2pt[*ite][p] )( nComponents*c+c );
805  }
806  }
807  }
808 
809  return res;
810 }
811 #endif
812 
813 template<typename space_type>
814 uint16_type
815 OperatorLagrangeP1<space_type>::localDof( typename domain_mesh_type::element_const_iterator it,
816  uint16_type localptid,
817  mpl::bool_<true> ) const
818 {
819  typedef typename domain_mesh_type::element_type::edge_permutation_type edge_permutation_type;
820  uint16_type localptid_dof = localptid;
821 
822  if ( domain_mesh_type::nDim == 2 &&
823  M_pset.pointToEntity( localptid ).first == 1 &&
824  it->edgePermutation( M_pset.pointToEntity( localptid ).second ) == edge_permutation_type::REVERSE_PERMUTATION )
825  {
826  //size_type edge_id = it->edge( M_pset.pointToEntity( localptid ).second ).id();
827  uint16_type edge_id = M_pset.pointToEntity( localptid ).second;
828 
829  uint16_type l = localptid - ( domain_mesh_type::element_type::numVertices * domain_fe_type::nDofPerVertex +
830  edge_id * domain_fe_type::nDofPerEdge );
831 
832  FEELPP_ASSERT( l != invalid_uint16_type_value && ( l < domain_fe_type::nDofPerEdge ) )
833  ( l )
834  ( domain_fe_type::nDofPerEdge )
835  ( localptid )
836  ( M_pset.pointToEntity( localptid ).second )
837  ( domain_mesh_type::element_type::numVertices * domain_fe_type::nDofPerVertex +
838  edge_id * domain_fe_type::nDofPerEdge ).error( "invalid value" );
839 
840  localptid_dof = ( domain_mesh_type::element_type::numVertices * domain_fe_type::nDofPerVertex +
841  edge_id * domain_fe_type::nDofPerEdge +
842  domain_fe_type::nDofPerEdge - 1 - l );
843 
844  FEELPP_ASSERT( localptid_dof != invalid_uint16_type_value &&
845  ( localptid_dof < domain_fe_type::nLocalDof ) )
846  ( l )
847  ( localptid )
848  ( localptid_dof )
849  ( domain_fe_type::nLocalDof )
850  ( M_pset.pointToEntity( localptid ).second )
851  ( domain_mesh_type::element_type::numVertices * domain_fe_type::nDofPerVertex +
852  edge_id * domain_fe_type::nDofPerEdge ).error( "invalid value" );
853  }
854 
855  return localptid_dof;
856 }
857 
858 //-----------------------------------------------------------------------------------------------------------------//
859 //-----------------------------------------------------------------------------------------------------------------//
860 //-----------------------------------------------------------------------------------------------------------------//
864 template<typename space_type>
865 boost::shared_ptr<OperatorLagrangeP1<space_type> >
866 opLagrangeP1_impl( boost::shared_ptr<space_type> const& Xh,
867  typename OperatorLagrangeP1<space_type>::backend_ptrtype const& backend,
868  std::vector<WorldComm> const& worldsComm,
869  std::string pathMeshLagP1,
870  std::string prefix,
871  bool rebuild,
872  bool parallel
873  )
874 {
875  return boost::shared_ptr<OperatorLagrangeP1<space_type> >( new OperatorLagrangeP1<space_type>( Xh,backend,worldsComm,pathMeshLagP1,prefix,rebuild,parallel ) );
876 }
877 
878 
879 template<typename Args>
880 struct compute_opLagrangeP1_return
881 {
882  typedef typename boost::remove_reference<typename parameter::binding<Args, tag::space>::type>::type::element_type space_type;
883  typedef boost::shared_ptr<OperatorLagrangeP1<space_type> > type;
884 };
885 
886 
887 BOOST_PARAMETER_FUNCTION(
888  ( typename compute_opLagrangeP1_return<Args>::type ), // 1. return type
889  lagrangeP1, // 2. name of the function template
890  tag, // 3. namespace of tag types
891  ( required
892  ( space, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
893  ) // required
894  ( optional
895  ( backend, *, Backend<typename compute_opLagrangeP1_return<Args>::space_type::value_type>::build() )
896  ( worldscomm, *, std::vector<WorldComm>( 1,Environment::worldComm() ) )
897  ( path, *( boost::is_convertible<mpl::_,std::string> ), std::string(".") )
898  ( prefix, *( boost::is_convertible<mpl::_,std::string> ), std::string("") )
899  ( rebuild, *( boost::is_integral<mpl::_> ), 1 )
900  ( parallel, *( boost::is_integral<mpl::_> ), 1 )
901  ) // optionnal
902 )
903 {
904  Feel::detail::ignore_unused_variable_warning( args );
905  return opLagrangeP1_impl(space,backend,worldscomm,path,prefix,rebuild,parallel);
906 }
907 
908 
909 
910 
911 
912 
913 
914 
915 
916 
917 } // Feel
918 #endif /* __OperatorLagrangeP1_H */

Generated on Sun Oct 20 2013 08:25:03 for Feel++ by doxygen 1.8.4