30 #ifndef __OperatorLinearComposite_H
31 #define __OperatorLinearComposite_H 1
33 #include <feel/feel.hpp>
34 #include <feel/feelcore/pslogger.hpp>
39 template<
class DomainSpace,
class DualImageSpace>
40 class OperatorLinearComposite :
public OperatorLinear<DomainSpace, DualImageSpace>
43 typedef OperatorLinear<DomainSpace,DualImageSpace> super;
44 typedef boost::shared_ptr<super> super_ptrtype;
48 typedef typename super::domain_space_type domain_space_type;
49 typedef typename super::dual_image_space_type dual_image_space_type;
50 typedef typename super::domain_space_ptrtype domain_space_ptrtype;
51 typedef typename super::dual_image_space_ptrtype dual_image_space_ptrtype;
53 typedef typename super::backend_type backend_type;
54 typedef typename super::backend_ptrtype backend_ptrtype;
56 typedef typename super::matrix_type matrix_type;
57 typedef boost::shared_ptr<matrix_type> matrix_ptrtype;
59 typedef FsFunctionalLinear<DualImageSpace> image_element_type;
60 typedef typename super::domain_element_type domain_element_type;
61 typedef typename super::dual_image_element_type dual_image_element_type;
63 typedef typename super::domain_element_range_type domain_element_range_type;
64 typedef typename super::domain_element_slice_type domain_element_slice_type;
65 typedef typename super::dual_image_element_range_type dual_image_element_range_type;
66 typedef typename super::dual_image_element_slice_type dual_image_element_slice_type;
68 typedef OperatorLinearComposite<DomainSpace, DualImageSpace> this_type;
69 typedef boost::shared_ptr<this_type> this_ptrtype;
71 OperatorLinearComposite (domain_space_ptrtype domainSpace,
72 dual_image_space_ptrtype dualImageSpace,
73 backend_ptrtype backend,
76 super( domainSpace,dualImageSpace,backend , false ),
86 int size = M_operators1.size();
87 M_operators1.insert( std::make_pair( size , op ) ) ;
91 void addElement(
int q, super_ptrtype
const& op )
93 M_operators1.insert( std::make_pair( q , op ) ) ;
96 void addList( std::vector< super_ptrtype >
const& vec )
99 for(
int q=0; q<q_max; q++)
103 void displayOperatorsNames()
105 int size1 = M_operators1.size();
106 int size2 = M_operators2.size();
108 LOG( INFO ) <<
" the composite operator linear "<<this->name()<<
" has following operators : ";
112 auto end = M_operators1.end();
113 for(
auto it=M_operators1.begin(); it!=end; it++)
114 LOG(INFO)<<it->second->name();
118 auto end = M_operators2.end();
119 for(
auto it=M_operators2.begin(); it!=end; it++)
120 LOG(INFO)<<it->second->name();
126 int size1 = M_operators1.size();
127 int size2 = M_operators2.size();
133 void addElement( boost::tuple<int,int>
const& qm , super_ptrtype
const& op )
135 M_operators2.insert( std::make_pair( qm , op ) );
138 void addList( std::vector< std::vector< super_ptrtype > >
const& vec )
140 int q_max = vec.size();
141 for(
int q=0; q<q_max; q++)
143 int m_max = vec[q].size();
144 for(
int m=0; m<m_max; m++)
146 auto tuple = boost::make_tuple( q , m );
152 void setScalars( std::vector<double> scalars )
154 M_scalars1 = scalars;
157 void setScalars( std::vector< std::vector< double > > scalars )
159 M_scalars2 = scalars;
162 void sumAllMatrices(matrix_ptrtype & matrix,
bool use_scalar_one=
false )
const
164 int size1 = M_operators1.size();
165 int size2 = M_operators2.size();
166 bool size_error=
false;
167 if( size1 > 0 && size2 > 0 )
169 if( (size1 + size2) == 0 )
172 FEELPP_ASSERT( !size_error )( size1 )( size2 ).error(
"OperatorLinearComposite has no elements, or both maps have elements" );
175 sumAllMatrices1( matrix, use_scalar_one );
177 sumAllMatrices2( matrix, use_scalar_one );
183 void sumAllMatrices1( matrix_ptrtype & matrix,
bool use_scalar_one=
false )
const
185 int size1 = M_operators1.size();
186 int size2 = M_operators2.size();
188 FEELPP_ASSERT( size1 > 0 )( size1 )( size2 ).error(
"OperatorLinearComposite has no elements" );
191 auto temp_matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
193 auto end = M_operators1.end();
194 for(
auto it=M_operators1.begin(); it!=end; ++it)
196 int position = it->first;
198 if( ! use_scalar_one )
199 scalar = M_scalars1[position];
200 it->second->matPtr(temp_matrix);
201 matrix->addMatrix( scalar , temp_matrix );
209 void sumAllMatrices2( matrix_ptrtype & matrix,
bool use_scalar_one=
false )
const
211 int size1 = M_operators1.size();
212 int size2 = M_operators2.size();
214 FEELPP_ASSERT( size2 > 0 )( size1 )( size2 ).error(
"OperatorLinearComposite has no elements" );
217 auto temp_matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
219 auto end = M_operators2.end();
220 for(
auto it=M_operators2.begin(); it!=end; ++it)
222 auto position = it->first;
223 int q = position.template get<0>();
224 int m = position.template get<1>();
226 if( ! use_scalar_one )
227 scalar = M_scalars2[q][m];
228 it->second->matPtr(temp_matrix);
229 matrix->addMatrix( scalar , temp_matrix );
237 void matrixPtr(
int q , matrix_ptrtype& matrix)
239 int q_max = M_operators1.size();
240 FEELPP_ASSERT( q < q_max )( q_max )( q ).error(
"OperatorLinearComposite has not enough elements" );
242 M_operators1.template at(q)->matPtr(matrix);
245 void matrixPtr(
int q , matrix_ptrtype& matrix)
const
247 int q_max = M_operators1.size();
248 FEELPP_ASSERT( q < q_max )( q_max )( q ).error(
"OperatorLinearComposite has not enough elements" );
250 M_operators1.template at(q)->matPtr(matrix);
254 void matrixPtr(
int q,
int m , matrix_ptrtype& matrix)
257 auto tuple = boost::make_tuple(q,m);
258 M_operators2.template at(tuple)->matPtr(matrix);
261 void matrixPtr(
int q,
int m , matrix_ptrtype& matrix)
const
263 auto tuple = boost::make_tuple(q,m);
264 M_operators2.template at(tuple)->matPtr(matrix);
269 super_ptrtype & operatorlinear(
int q)
271 int q_max = M_operators1.size();
272 FEELPP_ASSERT( q < q_max )( q_max )( q ).error(
"OperatorLinearComposite has not enough elements" );
273 return M_operators1.template at(q);
276 super_ptrtype & operatorlinear(
int q,
int m)
278 auto tuple = boost::make_tuple(q,m);
279 return M_operators2.template at(tuple);
284 apply(
const domain_element_type& de,
285 image_element_type& ie )
const
288 auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
289 sumAllMatrices( matrix,
true );
292 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
293 *_v1 = de;_v1->close();
294 vector_ptrtype _v2( M_backend->newVector( _test=ie.space() ) );
295 M_backend->prod( matrix, _v1, _v2 );
296 ie.container() = *_v2;
301 energy(
const typename domain_space_type::element_type & de,
302 const typename dual_image_space_type::element_type & ie )
const
305 auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
306 sumAllMatrices( matrix,
true );
309 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
310 *_v1 = de;_v1->close();
311 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
313 vector_ptrtype _v3( M_backend->newVector( _test=ie.functionSpace() ) );
314 M_backend->prod( matrix, _v1, _v3 );
320 apply(
const typename domain_space_type::element_type & de,
321 typename dual_image_space_type::element_type & ie )
323 auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
324 sumAllMatrices( matrix,
true );
327 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
328 *_v1 = de;_v1->close();
329 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
330 M_backend->prod( matrix, _v1, _v2 );
331 ie.container() = *_v2;
336 apply(
const domain_element_range_type & de,
337 typename dual_image_space_type::element_type & ie )
339 auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
340 sumAllMatrices(matrix,
true );
343 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
344 *_v1 = de;_v1->close();
345 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
346 M_backend->prod( matrix, _v1, _v2 );
347 ie.container() = *_v2;
351 apply(
const typename domain_space_type::element_type & de,
352 dual_image_element_range_type & ie )
354 auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
355 sumAllMatrices( matrix,
true );
358 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
359 *_v1 = de;_v1->close();
360 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
361 M_backend->prod( matrix, _v1, _v2 );
362 ie.container() = *_v2;
366 apply(
const domain_element_range_type & de,
367 dual_image_element_range_type & ie )
369 auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
370 sumAllMatrices( matrix,
true );
373 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
374 *_v1 = de;_v1->close();
375 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
376 M_backend->prod( matrix, _v1, _v2 );
377 ie.container() = *_v2;
381 apply(
const domain_element_slice_type & de,
382 typename dual_image_space_type::element_type & ie )
385 auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
386 sumAllMatrices( matrix,
true );
389 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
390 *_v1 = de;_v1->close();
391 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
392 M_backend->prod( matrix, _v1, _v2 );
393 ie.container() = *_v2;
398 apply(
const typename domain_space_type::element_type & de,
399 dual_image_element_slice_type & ie )
401 auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
402 sumAllMatrices( matrix,
true );
405 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
406 *_v1 = de;_v1->close();
407 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
408 M_backend->prod( matrix, _v1, _v2 );
409 ie.container() = *_v2;
413 apply( domain_element_slice_type de,
414 dual_image_element_slice_type ie )
416 auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
417 sumAllMatrices( matrix,
true );
420 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
421 *_v1 = de;_v1->close();
422 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
423 M_backend->prod( matrix, _v1, _v2 );
424 ie.container() = *_v2;
430 apply(
const domain_element_range_type & de,
431 dual_image_element_slice_type & ie )
433 auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
434 sumAllMatrices( matrix,
true );
437 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
438 *_v1 = de;_v1->close();
439 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
440 M_backend->prod( matrix, _v1, _v2 );
441 ie.container() = *_v2;
445 apply(
const domain_element_slice_type & de,
446 dual_image_element_range_type & ie )
448 auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
449 sumAllMatrices( matrix,
true );
452 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
453 *_v1 = de;_v1->close();
454 vector_ptrtype _v2( M_backend->newVector( _test=ie.functionSpace() ) );
455 M_backend->prod( matrix, _v1, _v2 );
456 ie.container() = *_v2;
461 applyInverse( domain_element_type& de,
462 const image_element_type& ie )
464 auto matrix = M_backend->newMatrix( _test=this->dualImageSpace(), _trial=this->domainSpace(), _pattern=M_pattern );
465 sumAllMatrices( matrix,
true );
468 vector_ptrtype _v1( M_backend->newVector( _test=de.functionSpace() ) );
469 vector_ptrtype _v2( M_backend->newVector( _test=ie.space() ) );
470 *_v2 = ie.container();
471 M_backend->solve( matrix, matrix, _v1, _v2 );
475 this_type&
operator=( this_type
const& m )
477 M_backend = m.M_backend;
478 M_pattern = m.M_pattern;
479 M_operators1 = m.M_operators1;
480 M_operators2 = m.M_operators2;
481 M_scalars1 = m.M_scalars1;
482 M_scalars2 = m.M_scalars2;
490 std::map< int, super_ptrtype > M_operators1;
491 std::map< boost::tuple<int,int> , super_ptrtype > M_operators2;
492 backend_ptrtype M_backend;
494 std::vector< double > M_scalars1;
495 std::vector< std::vector<double> > M_scalars2;
502 template<
typename Args>
503 struct compute_opLinearComposite_return
505 typedef typename boost::remove_reference<typename parameter::binding<Args, tag::domainSpace>::type>::type::element_type domain_space_type;
506 typedef typename boost::remove_reference<typename parameter::binding<Args, tag::imageSpace>::type>::type::element_type image_space_type;
508 typedef OperatorLinearComposite<domain_space_type, image_space_type> type;
509 typedef boost::shared_ptr<OperatorLinearComposite<domain_space_type, image_space_type> > ptrtype;
513 BOOST_PARAMETER_FUNCTION(
514 (
typename Feel::detail::compute_opLinearComposite_return<Args>::ptrtype ),
518 ( domainSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
519 ( imageSpace, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ) )
522 ( backend, *, Backend<
typename Feel::detail::compute_opLinearComposite_return<Args>::domain_space_type::value_type>::build() )
523 ( pattern, *, (
size_type)Pattern::COUPLED )
528 Feel::detail::ignore_unused_variable_warning( args );
529 typedef typename Feel::detail::compute_opLinearComposite_return<Args>::type oplincomposite_type;
530 typedef typename Feel::detail::compute_opLinearComposite_return<Args>::ptrtype oplincomposite_ptrtype;
531 return oplincomposite_ptrtype (
new oplincomposite_type( domainSpace,imageSpace,backend,pattern) );