29 #ifndef _PROJECTOR_HPP_
30 #define _PROJECTOR_HPP_
32 #include <feel/feelcore/parameter.hpp>
38 template<
class DomainSpace,
class DualImageSpace>
class Projector;
42 template<
typename Args>
45 typedef typename vf::detail::clean_type<Args,tag::domainSpace>::type::element_type domain_type;
46 typedef typename vf::detail::clean_type<Args,tag::imageSpace>::type::element_type image_type;
47 typedef boost::shared_ptr<Projector<domain_type,image_type> > return_type;
50 template<
typename Args>
53 typedef typename vf::detail::clean_type<Args,tag::domainSpace>::type::element_type domain_type;
54 typedef boost::shared_ptr<Projector<domain_type,domain_type> > lift_return_type;
65 template<
class DomainSpace,
class DualImageSpace>
66 class Projector :
public OperatorLinear<DomainSpace, DualImageSpace>
68 typedef Projector<DomainSpace,DualImageSpace> super;
77 typedef OperatorLinear<DomainSpace, DualImageSpace> ol_type;
79 typedef typename super::domain_space_type domain_space_type;
80 typedef typename super::dual_image_space_type dual_image_space_type;
81 typedef typename super::domain_space_ptrtype domain_space_ptrtype;
82 typedef typename super::dual_image_space_ptrtype dual_image_space_ptrtype;
83 typedef typename domain_space_type::element_type domain_element_type;
84 typedef typename dual_image_space_type::element_type dual_image_element_type;
86 typedef typename super::backend_type backend_type;
87 typedef typename super::backend_ptrtype backend_ptrtype;
88 typedef typename backend_type::sparse_matrix_type matrix_type;
89 typedef typename backend_type::vector_type vector_type;
90 typedef typename backend_type::vector_ptrtype vector_ptrtype;
91 typedef boost::shared_ptr<matrix_type> matrix_ptrtype;
93 typedef FsFunctionalLinear<DualImageSpace> image_element_type;
100 Projector( domain_space_ptrtype domainSpace,
101 dual_image_space_ptrtype dualImageSpace,
102 backend_ptrtype backend = Backend<double>::build( BACKEND_PETSC ),
103 ProjectorType proj_type=L2,
104 double epsilon = 0.01,
106 DirichletType dirichlet_type = WEAK
109 ol_type( domainSpace, dualImageSpace, backend ),
110 M_backend( backend ),
111 M_epsilon( epsilon ),
113 M_proj_type( proj_type ),
114 M_dir( dirichlet_type )
117 M_matrix = M_backend->newMatrix( _trial=domainSpace, _test=dualImageSpace, _pattern=Pattern::EXTENDED ) ;
127 template<
typename Args,
typename IntEltsDefault>
128 struct integrate_type
130 typedef typename vf::detail::clean_type<Args,tag::expr>::type _expr_type;
131 typedef typename vf::detail::clean2_type<Args,tag::range,IntEltsDefault>::type _range_type;
132 typedef typename vf::detail::clean2_type<Args,tag::quad, _Q< vf::ExpressionOrder<_range_type,_expr_type>::value > >::type _quad_type;
133 typedef typename vf::detail::clean2_type<Args,tag::quad1, _Q< vf::ExpressionOrder<_range_type,_expr_type>::value_1 > >::type _quad1_type;
136 template<
typename Range,
typename Expr>
137 void applyOn( Range range, Expr expr)
139 typedef typename boost::tuples::template element<0, Range>::type idim_type;
140 applyOn( range, expr, mpl::int_<idim_type::value>() );
144 BOOST_PARAMETER_MEMBER_FUNCTION( ( domain_element_type ),
151 ( range, *,
elements( this->dualImageSpace()->mesh() ) )
152 ( quad, *, (
typename integrate_type<Args,decltype(
elements( this->dualImageSpace()->mesh() ) )>::_quad_type() ) )
153 ( quad1, *, (
typename integrate_type<Args,decltype(
elements( this->dualImageSpace()->mesh() ) )>::_quad1_type() ) )
154 ( geomap, *, GeomapStrategyType::GEOMAP_OPT )
160 auto sol = this->domainSpace()->element();
162 ie = M_backend->newVector( this->dualImageSpace() );
164 form1( _test=this->dualImageSpace(), _vector=ie, _init=
true );
166 if ( (M_proj_type != LIFT) )
168 form1( _test=this->dualImageSpace(), _vector=ie ) +=
169 integrate( _range=range, _expr=expr *
id( this->dualImageSpace()->element() ),
170 _quad=quad, _quad1=quad1, _geomap=geomap );
172 else if ( ( M_proj_type == LIFT ) && ( M_dir == WEAK ) )
174 form1( _test=this->dualImageSpace(), _vector=ie ) +=
176 _expr=expr*( -grad( this->dualImageSpace()->element() )*vf::N() +
177 M_gamma / vf::hFace() *
id( this->dualImageSpace()->element() ) ),
178 _quad=quad, _quad1=quad1, _geomap=geomap );
182 if ( M_proj_type == DIFF )
184 form1( _test=this->dualImageSpace(), _vector=ie ) +=
186 _expr=expr*M_epsilon*( -grad( this->domainSpace()->element() )*vf::N() +
187 M_gamma / vf::hFace() *
id( this->dualImageSpace()->element() ) ),
188 _quad=quad, _quad1=quad1, _geomap=geomap );
193 M_matrixFull = M_backend->newMatrix( _trial=this->domainSpace(), _test=this->dualImageSpace(), _pattern=Pattern::EXTENDED );
194 auto bilinearForm = form2( _trial=this->domainSpace(), _test=this->dualImageSpace(), _matrix=M_matrixFull );
196 if ( ( M_proj_type == LIFT ) && ( M_dir == WEAK ) )
200 ( -trans(
id( this->dualImageSpace()->element() ) )*gradt( this->domainSpace()->element() )*vf::N()
201 -trans( idt( this->domainSpace()->element() ) )* grad( this->dualImageSpace()->element() )*vf::N()
202 + M_gamma * trans( idt( this->domainSpace()->element() ) )
203 *
id( this->dualImageSpace()->element() ) / vf::hFace()
204 ), _quad=quad, _quad1=quad1, _geomap=geomap );
207 M_matrixFull->close();
208 M_matrixFull->addMatrix( 1., M_matrix );
210 if ( ( M_proj_type == LIFT ) && ( M_dir == STRONG ) )
211 this->applyOn(range, expr);
213 M_backend->solve( M_matrixFull, sol, ie );
218 template<
typename RhsExpr>
220 operator()( RhsExpr
const& rhs_expr )
222 return this->
project( rhs_expr );
225 template<
typename RhsExpr>
227 operator()( domain_element_type& de, RhsExpr
const& rhs_expr )
229 de = this->
project( rhs_expr );
233 operator()( image_element_type
const& ie )
235 domain_element_type de = this->domainSpace()->element();
236 M_backend->solve( M_matrix, de, ie );
241 operator()( domain_element_type &de, image_element_type
const& ie )
243 M_backend->solve( M_matrix, de, ie );
246 template<
typename Range,
typename Expr>
248 operator()( Range
const& range ,Expr
const& expr )
250 return this->
project( expr, range );
254 apply( domain_element_type& de,
255 image_element_type
const& ie )
257 M_backend->solve( M_matrix, de, ie );
260 template<
typename RhsExpr>
262 apply( domain_element_type& de,
263 RhsExpr
const& rhs_expr )
269 template<
typename Expr>
270 domain_element_type derivate( Expr expr )
272 de = this->domainSpace()->element();
273 ie = M_backend->newVector( this->dualImageSpace() );
275 form1( _test=this->dualImageSpace(), _vector=ie, _init=
true );
277 form1( _test=this->dualImageSpace(), _vector=ie ) +=
279 _expr = - trace( expr * trans( grad( this->dualImageSpace()->element() ) ) ) );
281 form1(_test=this->dualImageSpace(), _vector=ie) +=
283 _expr = trans(
id( this->dualImageSpace()->element() ) ) * expr * vf::N() );
287 M_backend->solve( M_matrix, de, ie );
300 auto a = form2 ( _trial=this->domainSpace(),
301 _test=this->dualImageSpace(),
305 switch ( M_proj_type )
310 trans( idt( this->domainSpace()->element() ) )
311 *
id( this->dualImageSpace()->element() )
319 trans( idt( this->domainSpace()->element() ) )
320 *
id( this->dualImageSpace()->element() )
322 trace( gradt( this->domainSpace()->element() )
323 * trans( grad( this->dualImageSpace()->element() ) ) )
331 trans( idt( this->domainSpace()->element() ) )
332 *
id( this->dualImageSpace()->element() )
335 trace( gradt( this->domainSpace()->element() )
336 * trans( grad( this->dualImageSpace()->element() ) ) )
340 M_epsilon*( -trans(
id( this->dualImageSpace()->element() ) )*gradt( this->domainSpace()->element() )*vf::N() ) );
342 M_epsilon*( -trans( idt( this->domainSpace()->element() ) )* grad( this->dualImageSpace()->element() )*vf::N() ) );
344 M_epsilon*( M_gamma * trans( idt( this->domainSpace()->element() ) )
345 *
id( this->dualImageSpace()->element() ) / vf::hFace()
353 trans( idt( this->domainSpace()->element() ) )
354 *
id( this->dualImageSpace()->element() )
356 ( divt( this->domainSpace()->element() ) *
357 div( this->dualImageSpace()->element() ) )
365 trans( idt( this->domainSpace()->element() ) )
366 *
id( this->dualImageSpace()->element() )
369 curlzt( this->domainSpace()->element() )
370 * curlz( this->dualImageSpace()->element() )
378 trace( gradt( this->domainSpace()->element() )
379 * trans( grad( this->dualImageSpace()->element() ) ) )
388 trans( idt( this->domainSpace()->element() ) )
389 *
id( this->dualImageSpace()->element() )
393 M_gamma * hFace() * hFace()
394 * trans(jumpt( gradt(this->domainSpace()->element()) ))
395 * jump( grad(this->dualImageSpace()->element()) )
407 template<
typename Range,
typename Expr>
408 void applyOn( Range range, Expr expr, mpl::int_<MESH_ELEMENTS> ){}
410 template<
typename Range,
typename Expr>
411 void applyOn( Range range, Expr expr, mpl::int_<MESH_FACES> )
413 form2 ( _trial=this->domainSpace(),
414 _test=this->dualImageSpace(),
415 _matrix=M_matrixFull ) += on( _range=range , _element=de, _rhs=ie, _expr=expr );
418 backend_ptrtype M_backend;
419 const double M_epsilon;
420 const double M_gamma;
421 const ProjectorType M_proj_type;
423 matrix_ptrtype M_matrix;
424 matrix_ptrtype M_matrixFull;
425 domain_element_type de;
439 template<
typename TDomainSpace,
typename TDualImageSpace>
440 boost::shared_ptr< Projector<TDomainSpace, TDualImageSpace> >
441 projector( boost::shared_ptr<TDomainSpace>
const& domainspace,
442 boost::shared_ptr<TDualImageSpace>
const& imagespace,
443 typename Projector<TDomainSpace, TDualImageSpace>::backend_ptrtype
const& backend =
Backend<double>::build( BACKEND_PETSC ),
444 ProjectorType proj_type=L2,
double epsilon=0.01,
double gamma = 20, DirichletType dirichlet_type = WEAK)
447 boost::shared_ptr<Proj_type> proj(
new Proj_type( domainspace, imagespace, backend, proj_type, epsilon, gamma, dirichlet_type ) );
451 BOOST_PARAMETER_FUNCTION( (
typename Feel::detail::projector_args<Args>::return_type ),
455 ( domainSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
456 ( imageSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
459 ( type, (ProjectorType), L2 )
460 ( penaldir, *( boost::is_arithmetic<mpl::_> ), 20. )
461 ( backend, *, Backend<double>::build( BACKEND_PETSC ) )
464 return projector( domainSpace,imageSpace, backend, type, 0.01, penaldir );
467 BOOST_PARAMETER_FUNCTION( (
typename Feel::detail::lift_args<Args>::lift_return_type ),
471 ( domainSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
474 ( type, (DirichletType), WEAK )
475 ( penaldir, *( boost::is_arithmetic<mpl::_> ), 20. )
476 ( backend, *, Backend<double>::build( BACKEND_PETSC ) )
479 return projector( domainSpace, domainSpace, backend, ProjectorType::LIFT, 0.01 , penaldir, type );
Poly::value_type integrate(Polynomial< Poly, PsetType > const &p)
integrate a polynomial over a convex
Definition: operations.hpp:51
elements_type const & elements() const
Definition: elements.hpp:355
Polynomial< Pset, Scalar > project(Pset const &pset, Func const &f, IM const &im)
Definition: operations.hpp:436
base class for all linear algebra backends
Definition: backend.hpp:133
Projection made easy.
Definition: projector.hpp:38
boost::shared_ptr< Projector< TDomainSpace, TDualImageSpace > > projector(boost::shared_ptr< TDomainSpace > const &domainspace, boost::shared_ptr< TDualImageSpace > const &imagespace, typename Projector< TDomainSpace, TDualImageSpace >::backend_ptrtype const &backend=Backend< double >::build(BACKEND_PETSC), ProjectorType proj_type=L2, double epsilon=0.01, double gamma=20, DirichletType dirichlet_type=WEAK)
Definition: projector.hpp:441
boost::tuple< mpl::size_t< MESH_FACES >, typename MeshTraits< MeshType >::location_face_const_iterator, typename MeshTraits< MeshType >::location_face_const_iterator > boundaryfaces(MeshType const &mesh)
Definition: filters.hpp:1158
boost::tuple< mpl::size_t< MESH_FACES >, typename MeshTraits< MeshType >::location_face_const_iterator, typename MeshTraits< MeshType >::location_face_const_iterator > internalfaces(MeshType const &mesh)
Definition: filters.hpp:1175