29 #ifndef __OperatorLagrangeP1_H
30 #define __OperatorLagrangeP1_H 1
32 #include <feel/feeldiscr/functionspace.hpp>
40 namespace detailOpLagP1
46 template<
typename SpaceType>
47 struct SpaceToLagrangeP1Space
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;
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;
58 template<
template<u
int16_type Dim>
class PsetType>
62 typedef typename SpaceType::value_type value_type;
64 typedef Lagrange<1,PsetType> type;
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;
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;
76 typedef FunctionSpace<image_mesh_type, bases<basis_type> > image_space_type;
87 template<
typename SpaceType>
91 public OperatorLinear<SpaceType, typename detailOpLagP1::SpaceToLagrangeP1Space<SpaceType>::image_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;
113 typedef typename super::backend_ptrtype backend_ptrtype;
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;
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;
134 typedef PointSetFekete<domain_reference_convex_type, domain_space_type::basis_type::nOrder, value_type>
pset_type;
137 typedef typename domain_mesh_type::gm_type gm_type;
138 typedef typename domain_mesh_type::element_type element_type;
140 typedef boost::shared_ptr<gmc_type> gmc_ptrtype;
141 typedef typename gm_type::precompute_ptrtype gmpc_ptrtype;
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=
"",
162 bool parallelBuild=
true );
183 M_el2el = lp1.M_el2el;
184 M_el2pt = lp1.M_el2pt;
201 localDof(
typename domain_mesh_type::element_const_iterator it,
202 uint16_type localptid )
const
205 return localDof( it, localptid, mpl::bool_<false>() );
209 localDof(
typename domain_mesh_type::element_const_iterator ,
210 uint16_type localptid,
211 mpl::bool_<false> )
const
217 localDof(
typename domain_mesh_type::element_const_iterator it,
218 uint16_type localptid,
219 mpl::bool_<true> )
const;
254 image_mesh_ptrtype M_mesh;
256 std::vector<std::vector<size_type> > M_el2pt;
261 mutable p2m_type M_p2m;
269 template<
typename space_type>
271 backend_ptrtype
const& backend,
272 std::vector<WorldComm>
const& worldsComm,
273 std::string pathMeshLagP1,
282 M_mesh( new image_mesh_type(worldsComm[0]) ),
286 M_gmpc( new typename gm_type::precompute_type( this->domainSpace()->mesh()->gm(),
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);
295 if (M_mesh->worldComm().globalRank() == M_mesh->worldComm().masterRank() )
301 M_p2m.
visit( &M_pset );
309 M_p2m.mesh()->components().clear ( MESH_RENUMBER );
310 M_p2m.mesh()->components().clear ( MESH_CHECK );
311 M_p2m.mesh()->updateForUse();
313 M_p2m.mesh()->save( _name=nameMeshLagP1,_path=pathMeshLagP1 );
316 M_mesh->worldComm().barrier();
318 if (M_mesh->worldComm().globalRank() != M_mesh->worldComm().masterRank() )
320 M_p2m.mesh()->load( _name=nameMeshLagP1,_path=pathMeshLagP1,
321 _update=MESH_UPDATE_EDGES|MESH_UPDATE_FACES );
326 M_p2m.mesh()->load( _name=nameMeshLagP1,_path=pathMeshLagP1,
327 _update=MESH_UPDATE_EDGES|MESH_UPDATE_FACES );
331 VLOG(2) <<
"[P1 Lagrange] Pointset " << M_pset.points() <<
"\n";
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;
340 BOOST_FOREACH(
auto itMark, this->domainSpace()->mesh()->markerNames() )
342 M_mesh->addMarkerName( itMark.first,itMark.second[0],itMark.second[1] );
346 bool doParallelBuild = (M_mesh->worldComm().size()==1)?
false:parallelBuild;
355 auto it = this->domainSpace()->mesh()->beginElementWithProcessId(this->domainSpace()->mesh()->worldComm().localRank());
356 auto en = this->domainSpace()->mesh()->endElementWithProcessId(this->domainSpace()->mesh()->worldComm().localRank());
359 std::vector<size_type> new_node_numbers( this->domainSpace()->nLocalDof() );
362 size_type n_new_nodes = 0, n_new_elem = 0, n_new_faces = 0;
364 std::map<size_type,size_type> mapGhostDofIdClusterToProcess;
365 std::map<size_type,std::vector<size_type> > mapActiveEltWhichAreGhostInOtherPartition;
369 gmc_ptrtype gmc(
new gmc_type( this->domainSpace()->mesh()->gm(), *it, M_gmpc ) );
371 for (
size_type elid = 0, pt_image_id = 0; it != en; ++it )
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";
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);
385 for (
int cptEltp2m=0 ; itl != enl; ++itl, ++elid,++cptEltp2m )
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";
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() );
403 for (
int p = 0; p < image_mesh_type::element_type::numVertices; ++p )
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" );
412 if (doParallelBuild && !this->domainSpace()->dof()->dofGlobalClusterIsOnProc(this->domainSpace()->dof()->mapGlobalProcessToGlobalCluster(ptid)))
413 mapGhostDofIdClusterToProcess[this->domainSpace()->dof()->mapGlobalProcessToGlobalCluster(ptid)] = ptid;
417 dofindices[pt_image_id] = boost::make_tuple( elid, p, ptid );
426 new_node_numbers[ptid] = n_new_nodes;
428 point_type __pt( n_new_nodes, boost::get<0>( this->domainSpace()->dof()->dofPoint( ptid ) ) );
429 __pt.setProcessId( it->processId() );
432 DVLOG(2) <<
"[OperatorLagrangeP1] element id "
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";
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";
453 M_mesh->addPoint( __pt );
459 elt.setPoint( p, M_mesh->point( new_node_numbers[ptid] ) );
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 ) ) )
465 ( elt.point( p ).node() )
470 ( ptid ).warn(
"[after] inconsistent point coordinates" );
472 FEELPP_ASSERT( ublas::norm_2( elt.point( p ).node()-ublas::column( elt.G(), p ) ) < 1e-10 )
474 ( elt.point( p ).node() )
476 ( elt.G() ).warn(
"[after] inconsistent point coordinates" );
481 static const int nDim = element_type::nDim;
482 ublas::vector<double> D( nDim + 1, 0 );
486 std::fill( M.data().begin(), M.data().end(), value_type( 1.0 ) );
488 for (
int p = 0; p < element_type::numVertices; ++p )
490 ublas::column( P0, p ) = elt.vertex( p );
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" );
503 elt.setId ( n_new_elem );
507 auto const& theNewElt = M_mesh->addElement ( elt );
508 if (doParallelBuild && it->numberOfPartitions() >1) mapActiveEltWhichAreGhostInOtherPartition[it->id()][cptEltp2m] = theNewElt.id();
511 for (
unsigned int s=0; s<itl->numTopologicalFaces; s++ )
513 if ( !itl->facePtr( s ) )
continue;
515 auto const& theFaceBase = itl->face( s );
519 new_face.disconnect();
520 new_face.setOnBoundary(
true );
521 new_face.setProcessId( it->processId() );
522 new_face.setNumberOfPartitions( 1 );
527 for ( uint16_type p = 0; p < theFaceBase.nPoints(); ++p )
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] ) );
535 auto addFaceRes = M_mesh->addFace( new_face );
544 if ( doParallelBuild )
546 VLOG(2) <<
"[P1 Lagrange] start parallel build \n";
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 );
552 auto iv = this->domainSpace()->mesh()->beginGhostElement();
553 auto const en = this->domainSpace()->mesh()->endGhostElement();
554 for ( ; iv != en; ++iv )
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];
564 std::vector<int> nbMsgToRecv;
565 mpi::all_to_all( M_mesh->worldComm().localComm(),
569 VLOG(2) <<
"[P1 Lagrange] parallel build : finish first send \n";
572 for (
int proc=0; proc<theWorldCommSize; ++proc )
574 for (
int cpt=0; cpt<nbMsgToRecv[proc]; ++cpt )
578 M_mesh->worldComm().localComm().recv( proc, cpt, idEltRecv );
579 auto const& theeltIt = this->domainSpace()->mesh()->elementIterator(idEltRecv);
582 auto const& idOfNewElt = mapActiveEltWhichAreGhostInOtherPartition[idEltRecv];
584 std::vector<boost::tuple<size_type,ublas::vector<double> > > resultClusterDofsAndNodesToSend(M_p2m.mesh()->numPoints());
586 auto itp = M_p2m.mesh()->beginPoint();
587 auto const enp = M_p2m.mesh()->endPoint();
588 for ( ; itp!=enp ; ++itp )
590 uint16_type localptid = itp->id();
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] = boost::make_tuple( idInProcAsked, thedofnode);
598 auto itl = M_p2m.mesh()->beginElement();
599 auto const enl = M_p2m.mesh()->endElement();
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 )
604 std::vector<uint16_type> vecLocalIdPtToSend(image_mesh_type::element_type::numVertices);
606 for (
int p = 0; p < image_mesh_type::element_type::numVertices; ++p )
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;
612 resultToSendBis[elidp2m] = boost::make_tuple(vecLocalIdPtToSend,idOfNewElt[elidp2m] );
614 auto resultToSend = boost::make_tuple(resultToSendBis,resultClusterDofsAndNodesToSend);
615 M_mesh->worldComm().localComm().send( proc, cpt, resultToSend);
619 VLOG(2) <<
"[P1 Lagrange] parallel build : finish first recv and resend response \n";
621 std::map<size_type,size_type> new_ghost_node_numbers;
624 for (
int proc=0; proc<theWorldCommSize; ++proc )
626 for (
int cpt=0; cpt<nbMsgToSend[proc]; ++cpt )
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 );
632 auto const& myGhostEltBase = this->domainSpace()->mesh()->element(mapMsg[proc][cpt],proc);
634 auto const& mapLocalToGlobalPointId = resultRecvData.template get<1>();
635 auto resultRecv = resultRecvData.template get<0>();
637 auto itelt = resultRecv.begin();
638 auto const enelt = resultRecv.end();
639 for ( ;itelt!=enelt;++itelt )
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);
647 elt.setProcessId(proc);
648 elt.setNeighborPartitionIds( std::vector<int>({proc}) );
649 bool findConnection=
false;
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 )
656 auto const globclusterdofRecv = mapLocalToGlobalPointId[(int)*itpt].
template get<0>();
658 if ( this->domainSpace()->dof()->dofGlobalClusterIsOnProc(globclusterdofRecv))
660 theptid = this->domainSpace()->dof()->mapGlobalClusterToGlobalProcess( globclusterdofRecv - this->domainSpace()->dof()->firstDofGlobalCluster() );
661 findConnection =
true;
663 else if (mapGhostDofIdClusterToProcess.find(globclusterdofRecv) != mapGhostDofIdClusterToProcess.end() )
665 theptid = mapGhostDofIdClusterToProcess[globclusterdofRecv];
666 findConnection =
true;
670 elt.setPoint( p, M_mesh->point( new_node_numbers[theptid] ) );
673 if (new_ghost_node_numbers.find(globclusterdofRecv) == new_ghost_node_numbers.end() )
675 new_ghost_node_numbers[globclusterdofRecv] = n_new_nodes;
677 point_type __pt( n_new_nodes, mapLocalToGlobalPointId[(
int)*itpt].
template get<1>() );
680 auto const& theNewPt = M_mesh->addPoint( __pt );
684 elt.setPoint( p, M_mesh->point( new_ghost_node_numbers[globclusterdofRecv] ) );
689 if (!findConnection)
continue;
692 elt.setId ( n_new_elem );
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>() ) );
707 DVLOG(2) <<
"[P1 Lagrange] Number of points in mesh: " << M_mesh->numPoints() <<
"\n";
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();
717 #if defined(FEELPP_ENABLE_MPI_MODE)
718 auto Xh_image = dual_image_space_type::New(_mesh=M_mesh,_worldscomm=worldsComm);
720 auto Xh_image = dual_image_space_type::New(_mesh=M_mesh);
722 this->init( this->domainSpace(),
739 for (
size_type i = 0; i < this->domainSpace()->nLocalDof()/domain_space_type::nComponents; ++i )
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 ) )
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" );
757 template<
typename space_type>
762 for (
size_type i = 0; i < this->domainSpace()->nLocalDof()/domain_space_type::nComponents; ++i )
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";
777 template<
typename space_type>
781 typename image_space_type::element_type res( this->dualImageSpace(), u.name() );
783 res = ublas::scalar_vector<value_type>( res.size(), value_type( 0 ) );
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();
788 for ( ; it != en; ++it )
790 gmc_ptrtype gmc(
new gmc_type( this->domainSpace()->mesh()->gm(), *it, M_gmpc ) );
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();
796 for ( ; ite != ene; ++ite )
798 typename domain_mesh_type::element_type
const& elt = M_mesh->element( *ite );
800 for (
int c = 0; c < nComponents; ++c )
801 for (
int p = 0; p < domain_mesh_type::element_type::numVertices; ++p )
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 );
813 template<
typename space_type>
815 OperatorLagrangeP1<space_type>::localDof(
typename domain_mesh_type::element_const_iterator it,
816 uint16_type localptid,
817 mpl::bool_<true> )
const
819 typedef typename domain_mesh_type::element_type::edge_permutation_type edge_permutation_type;
820 uint16_type localptid_dof = localptid;
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 )
827 uint16_type edge_id = M_pset.pointToEntity( localptid ).second;
829 uint16_type l = localptid - ( domain_mesh_type::element_type::numVertices * domain_fe_type::nDofPerVertex +
830 edge_id * domain_fe_type::nDofPerEdge );
834 ( domain_fe_type::nDofPerEdge )
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" );
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 );
845 ( localptid_dof < domain_fe_type::nLocalDof ) )
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" );
855 return localptid_dof;
864 template<
typename space_type>
865 boost::shared_ptr<OperatorLagrangeP1<space_type> >
867 typename OperatorLagrangeP1<space_type>::backend_ptrtype
const& backend,
868 std::vector<WorldComm>
const& worldsComm,
869 std::string pathMeshLagP1,
875 return boost::shared_ptr<OperatorLagrangeP1<space_type> >(
new OperatorLagrangeP1<space_type>( Xh,backend,worldsComm,pathMeshLagP1,prefix,rebuild,parallel ) );
879 template<
typename Args>
880 struct compute_opLagrangeP1_return
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;
887 BOOST_PARAMETER_FUNCTION(
888 (
typename compute_opLagrangeP1_return<Args>::type ),
892 ( space, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
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 )
904 Feel::detail::ignore_unused_variable_warning( args );