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
oseen.hpp
Go to the documentation of this file.
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): Christoph Winkelmann <christoph.winkelmann@epfl.ch>
6  Date: 2006-10-03
7 
8  Copyright (C) 2006 EPFL
9  Copyright (C) 2011 Université Joseph Fourier
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
31 /* Oseen solver with interior penalty stabilization, and
32  weak Dirichlet and homogeneous Neumann boundary conditions [1],
33  building u-p-coupled matrix solved by an adaptive gmm backend
34 
35  equation solved:
36  sigma u + (grad u)*beta - div(2 nu sgrad u) + grad p = f on domain
37  div u + epsilon p = c on domain
38  u = g on D-bdry
39  2 nu sgrad u n - p n = 0 on N-bdry
40 
41  where sgrad is the symmetric gradient tensor.
42 
43  References:
44  [1] E. Burman, M. Fernandez and P. Hansbo:
45  Continuous interior penalty finite element method for Oseen's equations
46  SIAM J. Numer. Anal. Vol. 44, No. 3, pp. 1248-1274
47 */
48 #include <set>
49 
50 #include <boost/program_options/parsers.hpp>
51 #include <boost/program_options/variables_map.hpp>
52 
53 #include <feel/feeldiscr/functionspace.hpp>
54 #include <feel/feelpoly/im.hpp>
55 
56 #include <feel/feelalg/backend.hpp>
57 
58 #include <feel/feelvf/vf.hpp>
60 
61 namespace Feel
62 {
63 
64 
65 template<class Space,
66  uint16_type imOrder = 3*boost::fusion::result_of::value_at_c<typename Space::basis_type, 0>::basis_type::nOrder-1,
67  template<uint16_type,uint16_type,uint16_type> class Entity = Simplex>
68 class Oseen
69 {
70 public:
71 
72  // -- TYPEDEFS --
73  typedef Space space_type;
74 
75  static const uint16_type Dim = space_type::nDim;
76  typedef typename space_type::value_type value_type;
77 
78  /* backend */
79  typedef Backend<value_type> backend_type;
80 
81 
82  typedef boost::shared_ptr<backend_type> backend_ptrtype;
83  typedef typename backend_type::sparse_matrix_type sparse_matrix_type;
84  typedef typename backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
85  typedef typename backend_type::vector_type vector_type;
86  typedef typename backend_type::vector_ptrtype vector_ptrtype;
87 
88  /* mesh */
89  typedef typename space_type::mesh_type mesh_type;
90  typedef typename space_type::mesh_ptrtype mesh_ptrtype;
91 
92  /* P0 basis */
93  typedef bases<Lagrange<0,Scalar,Discontinuous> > basis_i_type;
94 
95  /* spaces */
96  typedef FunctionSpace<mesh_type, basis_i_type, Discontinuous,value_type > space_i_type;
97  typedef boost::shared_ptr<space_type> space_ptrtype;
98  typedef boost::shared_ptr<space_i_type> space_i_ptrtype;
99  typedef typename space_type::element_type element_type;
100  typedef typename element_type::template sub_element<0>::type element_u_type;
101  typedef typename element_type::template sub_element<1>::type element_p_type;
102  typedef typename space_i_type::element_type element_i_type;
103 
104  /* quadrature */
105  static const uint16_type uOrder = boost::fusion::result_of::value_at_c<typename Space::basis_type, 0>::type::nOrder;
106  static const uint16_type pOrder = boost::fusion::result_of::value_at_c<typename Space::basis_type, 1>::type::nOrder;
107  template<int Order>
108  struct im
109  {
110  typedef IM<Dim, Order, value_type, Entity> type;
111  };
112 
113  Oseen( const space_ptrtype& space,
114  const backend_ptrtype& backend,
115  const std::set<flag_type>& dirichletFlags,
116  const std::set<flag_type>& neumannFlags,
117  bool weak_dirichlet = true,
118  bool export_matlab = false
119  );
120 
121  Oseen( const space_ptrtype& space,
122  const backend_ptrtype& backend,
123  const std::set<flag_type>& dirichletFlags,
124  const std::set<flag_type>& neumannFlags,
125  po::variables_map const& vm,
126  std::string const& prefix = ""
127  );
128 
129  // setting of options
130  void setBCCoeffDiff( value_type bccoeffdiff )
131  {
132  M_bcCoeffDiff = bccoeffdiff;
133  }
134 
135  void setBCCoeffConv( value_type bccoeffconv )
136  {
137  M_bcCoeffConv = bccoeffconv;
138  }
139 
140  void setStabCoeffDiv( value_type stabcoeffdiv )
141  {
142  M_stabCoeffDiv = stabcoeffdiv;
143  }
144 
145  void setStabCoeffP( value_type stabcoeffp )
146  {
147  M_stabCoeffP = stabcoeffp;
148  }
149 
150  void decouplePstab( element_p_type& pStab, value_type theta )
151  {
152  M_pStab = &pStab;
153  M_thetaPstab = theta;
154  }
155 
156  void couplePstab()
157  {
158  M_pStab = 0;
159  M_thetaPstab = -1.0;
160  }
161 
162  void setEpsCompress( value_type epscompress )
163  {
164  M_epsCompress = epscompress;
165  }
166 
167  void setDivDivCoeff( value_type divdivcoeff );
168 
169  // update operator and rhs with given expressions
170  template<typename ItRange, typename EsigmaInc,
171  typename EnuInc, typename EnuAbs,
172  typename Ebeta, typename Ef, typename Ec, typename Eg,
173  typename EnoSlip>
174  void update( const ItRange& itRange,
175  const EsigmaInc& sigmaInc,
176  const EnuInc& nuInc,
177  const EnuAbs& nuAbs,
178  const Ebeta& beta,
179  const Ef& f,
180  const Ec& c,
181  const Eg& g,
182  const EnoSlip& noSlip,
183  bool updateStabilization = true );
184 
185  // solve system, call not needed, but possible
186  void solve();
187 
188  // results
189  const element_u_type& velocity()
190  {
191  solve();
192  return u;
193  }
194  const element_p_type& pressure()
195  {
196  solve();
197  return p;
198  }
199 
200  element_type& solution()
201  {
202  return M_sol;
203  }
204  const element_type& solution() const
205  {
206  return M_sol;
207  }
208 
209  value_type stabilizationEnergy() const;
210 
211 private:
212 
216  void solveNonSym( sparse_matrix_ptrtype const& D, element_type& u, vector_ptrtype const& F );
217 
218 private:
219 
220  mesh_ptrtype M_mesh;
221 
222  space_ptrtype M_space;
223 
224  element_type M_sol;
225  element_u_type u;
226  element_p_type p;
227 
228  sparse_matrix_ptrtype M_matrixFull;
229  sparse_matrix_ptrtype M_matrixAu;
230  sparse_matrix_ptrtype M_matrixStab;
231  vector_ptrtype M_vectorRhsFull;
232  vector_ptrtype M_vectorSolFull;
233 
234  value_type M_bcCoeffDiff;
235  value_type M_bcCoeffConv;
236 
237  value_type M_stabCoeffDiv;
238  value_type M_stabCoeffP;
239  element_p_type* M_pStab;
240  value_type M_thetaPstab;
241 
242  value_type M_epsCompress;
243  value_type M_divDivCoeff;
244 
245  bool M_updated;
246 
247  backend_ptrtype M_backend;
248 
249  std::set<flag_type> M_dirichletFlags;
250  std::set<flag_type> M_neumannFlags;
251 
252  space_i_ptrtype M_space_i;
253  element_i_type M_stabWeightP;
254 
255  bool M_weak_dirichlet;
256  bool M_export_matlab;
257 
258 }; // class Oseen
259 
260 template<class Space, uint16_type imOrder,
261  template<uint16_type,uint16_type,uint16_type> class Entity>
262 Oseen<Space, imOrder, Entity>::Oseen( const space_ptrtype& space,
263  const backend_ptrtype& backend,
264  const std::set<flag_type>& dirichletFlags,
265  const std::set<flag_type>& neumannFlags,
266  bool weak_dirichlet,
267  bool export_matlab )
268  :
269  M_mesh( space->mesh() ),
270  M_space( space ),
271  M_sol( space, "sol" ),
272  u( M_sol.template element<0>() ),
273  p( M_sol.template element<1>() ),
274  M_matrixFull( backend->newMatrix( M_space, M_space ) ),
275  M_matrixAu( backend->newMatrix( M_space, M_space ) ),
276  M_matrixStab( backend->newMatrix( M_space, M_space ) ),
277  M_vectorRhsFull( backend->newVector( M_space ) ),
278  M_vectorSolFull( backend->newVector( M_space ) ),
279  M_bcCoeffDiff( 100 ),
280  M_bcCoeffConv( 100 ),
281  M_stabCoeffDiv( 0.0 ),
282  M_stabCoeffP( 0.0 ),
283  M_pStab( 0 ),
284  M_thetaPstab( -1.0 ),
285  M_epsCompress( 0.0 ),
286  M_divDivCoeff( 0.0 ),
287  M_updated( false ),
288  M_backend( backend ),
289  M_dirichletFlags( dirichletFlags ),
290  M_neumannFlags( neumannFlags ),
291  M_space_i( new space_i_type( space->mesh() ) ),
292  M_stabWeightP( M_space_i, "cs" ),
293  M_weak_dirichlet( weak_dirichlet ),
294  M_export_matlab( export_matlab )
295 {
296  using namespace Feel::vf;
297 
298  element_u_type& v = u;
299  element_p_type& q = p;
300 
301  DVLOG(2) << "[Oseen::Oseen] -(p,div v) + (q,div u)\n";
302  form2( _test=M_space, _trial=M_space, _matrix=M_matrixAu, _init=true ) =
303  integrate( elements( M_mesh ),
304  - div( v ) * idt( p )
305  )+
306  integrate( elements( M_mesh ),
307  + id( q ) * divt( u )
308  );
309 
310 
311  //M_matrixStab.close();
312 
313  // --- Set rhs to zero (pressure component is not set otherwise)
314  M_vectorRhsFull->zero();
315 
316 } // Oseen constructor
317 
318 template<class Space, uint16_type imOrder,
319  template<uint16_type,uint16_type,uint16_type> class Entity>
320 Oseen<Space, imOrder, Entity>::Oseen( const space_ptrtype& space,
321  const backend_ptrtype& backend,
322  const std::set<flag_type>& dirichletFlags,
323  const std::set<flag_type>& neumannFlags,
324  po::variables_map const& vm,
325  std::string const& prefix )
326  :
327  M_mesh( space->mesh() ),
328  M_space( space ),
329  M_sol( space, "sol" ),
330  u( M_sol.template element<0>() ),
331  p( M_sol.template element<1>() ),
332  M_matrixFull( backend->newMatrix( M_space, M_space ) ),
333  M_matrixAu( backend->newMatrix( M_space, M_space ) ),
334  M_matrixStab( backend->newMatrix( M_space, M_space ) ),
335  M_vectorRhsFull( backend->newVector( M_space ) ),
336  M_vectorSolFull( backend->newVector( M_space ) ),
337  M_bcCoeffDiff( 0.0 ),
338  M_bcCoeffConv( 0.0 ),
339  M_stabCoeffDiv( 0.0 ),
340  M_stabCoeffP( 0.0 ),
341  M_epsCompress( 0.0 ),
342  M_divDivCoeff( 0.0 ),
343  M_updated( false ),
344  M_backend( backend ),
345  M_dirichletFlags( dirichletFlags ),
346  M_neumannFlags( neumannFlags ),
347  M_space_i( new space_i_type( space->mesh() ) ),
348  M_stabWeightP( M_space_i, "cs" ),
349  M_weak_dirichlet( true ),
350  M_export_matlab( true )
351 {
352  std::string _prefix = prefix;
353 
354  if ( !_prefix.empty() )
355  _prefix += "-";
356 
357  M_bcCoeffDiff = vm[_prefix+"oseen-bc-coeff-diff"].template as<double>();
358  M_bcCoeffConv = vm[_prefix+"oseen-bc-coeff-conv"].template as<double>();
359  M_stabCoeffDiv = vm[_prefix+"oseen-stab-coeff-div"].template as<double>();
360  M_stabCoeffP = vm[_prefix+"oseen-stab-coeff-p"].template as<double>();
361  M_epsCompress = vm[_prefix+"oseen-eps-compress"].template as<double>();
362  M_divDivCoeff = vm[_prefix+"oseen-divdiv-coeff"].template as<double>();
363  M_weak_dirichlet = vm[_prefix+"oseen-weak-dirichlet"].template as<bool>();
364  M_export_matlab = vm[_prefix+"oseen-export-matlab"].template as<bool>();
365 
366  using namespace Feel::vf;
367 
368  element_u_type& v = u;
369  element_p_type& q = p;
370 
371  DVLOG(2) << "[Oseen::Oseen] -(p,div v) + (q,div u)\n";
372  form2( _test=M_space, _trial=M_space, _matrix=M_matrixAu, _init=true ) =
373  integrate( elements( M_mesh ),
374  - div( v ) * idt( p )
375  )+
376  integrate( elements( M_mesh ),
377  + id( q ) * divt( u )
378  );
379 
380  //M_matrixStab.close();
381 
382  // --- Set rhs to zero (pressure component is not set otherwise)
383  M_vectorRhsFull->zero();
384 
385 } // Oseen constructor
386 
387 template<class Space, uint16_type imOrder,
388  template<uint16_type,uint16_type,uint16_type> class Entity>
389 void
390 Oseen<Space, imOrder, Entity>::setDivDivCoeff( value_type divdivcoeff )
391 {
392  using namespace Feel::vf;
393 
394  element_u_type& v = u;
395 
396  value_type delta = divdivcoeff - M_divDivCoeff;
397 
398  if ( delta != 0.0 )
399  {
400  DVLOG(2) << "[Oseen::set_divdivcoeff] (h gamma div u,div v)\n";
401  form2( _test=M_space, _trial=M_space, _matrix=M_matrixAu ) +=
402  integrate( elements( M_mesh ),
403  h()*delta*div( v )*divt( u )
404  );
405  M_divDivCoeff = divdivcoeff;
406  }
407 }
408 template<class Space, uint16_type imOrder,
409  template<uint16_type,uint16_type,uint16_type> class Entity>
410 typename Oseen<Space, imOrder, Entity>::value_type
411 Oseen<Space, imOrder, Entity>::stabilizationEnergy() const
412 {
413  vector_ptrtype temp( M_backend->newVector( M_space ) );
414  //backend_type::applyMatrix( M_matrixStab, M_sol.container(), temp );
415  M_backend->prod( *M_matrixStab, M_sol.container(), *temp );
416  return M_backend->dot( *temp, M_sol.container() );
417 }
418 
419 
420 template<class Space, uint16_type imOrder,
421  template<uint16_type,uint16_type,uint16_type> class Entity>
422 template<typename ItRange, typename EsigmaInc,
423  typename EnuInc, typename EnuAbs,
424  typename Ebeta, typename Ef, typename Ec, typename Eg,
425  typename EnoSlip>
426 void Oseen<Space, imOrder, Entity>::update( const ItRange& itRange,
427  const EsigmaInc& sigmaInc,
428  const EnuInc& nuInc,
429  const EnuAbs& nuAbs,
430  const Ebeta& beta,
431  const Ef& f,
432  const Ec& c,
433  const Eg& g,
434  const EnoSlip& noSlip,
435  bool updateStabilization )
436 {
437  element_u_type& v = u;
438  element_p_type& q = p;
439 
440  M_updated = true;
441 
442  using namespace Feel::vf;
443 
444  auto bcCoeffVeloFull = ( nuAbs )*( noSlip )*M_bcCoeffDiff/hFace()-chi( ( trans( beta )*N() ) < 0 ) * trans( beta )*N();
445  auto bcCoeffVeloNorm = M_bcCoeffConv * vf::max( vf::sqrt( trans( beta )*( beta ) ), ( nuAbs )/hFace() );
446 
447  // --- right hand side rhs
448 
449  // rhs volume terms
450  DVLOG(2) << "[Oseen::update] (f,v)+(c,q)\n";
451  form1( M_space, M_vectorRhsFull, _init=true ) =
452  integrate( elements( M_mesh ), val( trans( f ) )*id( v ) + val( c )*id( q ) );
453  DVLOG(2) << "[Oseen::update] (f,v)+(c,q) done\n";
454 
455  if ( M_weak_dirichlet )
456  {
457  // rhs boundary terms
458  for ( auto diriIter = M_dirichletFlags.begin(),
459  diriEnd = M_dirichletFlags.end();
460  diriIter != diriEnd; ++diriIter )
461  {
462  DVLOG(2) << "[Oseen::update] <bc(g),v>_Gamma_"
463  << *diriIter << "\n";
464  form1( M_space, M_vectorRhsFull ) +=
465  integrate( markedfaces( M_mesh, *diriIter ),
466  ( -val( noSlip )*( val( nuAbs )*trans( N() )
467  *( grad( v )+trans( grad( v ) ) )
468  + id( q )*trans( N() ) )
469  + val( bcCoeffVeloFull )*trans( id( v ) )
470  + val( bcCoeffVeloNorm )*( trans( id( v ) )*N() )*trans( N() ) )
471  * val( g )
472  );
473  DVLOG(2) << "[Oseen::update] <bc(g),v>_Gamma_"
474  << *diriIter << " done\n";
475  }
476  } // M_weak_dirichlet
477 
478  // --- full matrix
479 
480  // incremental volume terms
481  DVLOG(2) << "[Oseen::update] adding (nu D(u),D(v)) + (sigma u,v)\n";
482  form2( M_space, M_space, _matrix=M_matrixAu ) +=
483  integrate( itRange,
484  val( nuInc )*trace( grad( v )*( gradt( u )+trans( gradt( u ) ) ) )
485  + val( sigmaInc )*trans( id( v ) )*idt( u )
486  );
487  DVLOG(2) << "[Oseen::update] adding (nu D(u),D(v)) + (sigma u,v) done\n";
488 
489  //M_matrixFull = M_matrixAu;
490  //M_matrixFull->addMatrix( 1.0, M_matrixAu );
491 
492  // convective terms
493  DVLOG(2) << "[Oseen::update] ((beta*grad)u,v)\n";
494  form2( M_space, M_space, M_matrixFull, _init=true ) =
495  integrate( elements( M_mesh ),
496  trans( id( v ) ) * ( gradt( u )*val( beta ) )
497  );
498 #if 0
499 
500  // --- stabilization matrices
501  if ( ( M_stabCoeffDiv != 0.0 ) ||
502  ( M_stabCoeffP != 0.0 ) ||
503  ( M_epsCompress != 0.0 ) )
504  {
505  if ( updateStabilization )
506  {
507  form2( M_space, M_space, M_matrixStab, _init=true );
508 #if 1
509 
510  if ( M_stabCoeffDiv != 0.0 )
511  {
512  DVLOG(2) << "[Oseen::update] <gamma_div [div u],[div v]>_Gamma_int\n";
513  form( M_space, M_space, M_matrixStab ) +=
514  integrate( internalfaces( M_mesh ), typename im<imOrder-1>::type(),
515  val( M_stabCoeffDiv * vf::pow( hFace(),2.0 ) *
516  vf::sqrt( trans( beta )*( beta ) ) )
517  * trans( jump( div( v ) ) ) * jumpt( divt( u ) )
518  );
519  }
520 
521  if ( M_stabCoeffP != 0.0 )
522  {
523  DVLOG(2) << "[Oseen::update] <gamma_p [grad p],[grad q]>_Gamma_int\n";
524  M_stabWeightP =
525  vf::project( M_space_i, elements( M_mesh ),
526  M_stabCoeffP * vf::pow( h(),3.0 ) /
527  max( h() * vf::sqrt( trans( beta ) * ( beta ) ),
528  ( nuAbs ) )
529  );
530 
531  if ( M_pStab )
532  {
533  form( M_space, M_space, M_matrixStab ) +=
534  integrate( internalfaces( M_mesh ),
535  typename im<2*pOrder-2>::type(),
536  val( M_thetaPstab
537  * maxface( idv( M_stabWeightP ) ) )
538  * ( leftface ( grad ( q )*N() ) *
539  leftfacet ( gradt( p )*N() ) +
540  rightface ( grad ( q )*N() ) *
541  rightfacet( gradt( p )*N() ) )
542  );
543  }
544 
545  else
546  {
547  form( M_space, M_space, M_matrixStab ) +=
548  integrate( internalfaces( M_mesh ),
549  typename im<2*pOrder-2>::type(),
550  maxface( idv( M_stabWeightP ) )
551  * jump( grad( q ) ) * jumpt( gradt( p ) )
552  );
553  }
554  }
555 
556 #endif
557 
558  if ( M_epsCompress > 0 )
559  {
560  VLOG(1) << "[Oseen] adding pseudo compressibility term with coeff= " << M_epsCompress << "\n";
561  form( M_space, M_space, M_matrixStab ) +=
562  integrate( elements( M_mesh ), typename im<2*pOrder>::type(),
563  M_epsCompress * id( q ) * idt( p )
564  );
565  VLOG(1) << "[Oseen] adding pseudo compressibility term with coeff= " << M_epsCompress << " done\n";
566  }
567 
568  M_matrixStab->close();
569  }
570 
571  }
572 
573 #endif
574 
575 #if 0
576 
577  if ( M_pStab && ( M_stabCoeffP != 0.0 ) )
578  {
579  DVLOG(2) << "[Oseen::update] added rhs stab\n";
580  const int staborder = 2*pOrder-2;
581  form( M_space, *M_vectorRhsFull ) +=
582  integrate( internalfaces( M_mesh ), typename im<staborder>::type(),
583  val( maxface( idv( M_stabWeightP ) ) )
584  * ( - ( leftface ( grad ( q )*N() ) *
585  rightfacev( gradv( *M_pStab )*N() ) +
586  rightface ( grad ( q )*N() ) *
587  leftfacev ( gradv( *M_pStab )*N() ) )
588  + ( M_thetaPstab - 1.0 )
589  * ( ( leftface ( grad ( q )*N() ) *
590  leftfacev ( gradv( *M_pStab )*N() ) +
591  rightface ( grad ( q )*N() ) *
592  rightfacev ( gradv( *M_pStab )*N() ) )
593  )
594  )
595  );
596  DVLOG(2) << "[Oseen::update] added rhs stab done\n";
597  }
598 
599 #endif
600 
601  // boundary terms
602  for ( auto diriIter = M_dirichletFlags.begin(),
603  diriEnd = M_dirichletFlags.end();
604  diriIter != diriEnd; ++diriIter )
605  {
606  if ( M_weak_dirichlet )
607  {
608  DVLOG(2) << "[Oseen::update] <bc(u),v>_Gamma_"
609  << *diriIter << "\n";
610  form2( M_space, M_space, M_matrixFull ) +=
611  integrate( markedfaces( M_mesh, *diriIter ),
612  -val( noSlip )*( val( ( nuAbs )*trans( N() ) )
613  *( ( grad( v )+trans( grad( v ) ) ) *idt( u ) +
614  ( gradt( u )+trans( gradt( u ) ) )*id( v ) )
615  + id( q )*trans( N() )*idt( u )
616  - idt( p )*trans( N() )*id( v ) )
617  + val( bcCoeffVeloFull )*( trans( id( v ) )*idt( u ) )
618  + val( bcCoeffVeloNorm )
619  * ( trans( id( v ) )*N() ) * ( trans( N() )*idt( u ) )
620  );
621  DVLOG(2) << "[Oseen::update] <bc(u),v>_Gamma_"
622  << *diriIter << " done\n";
623  }
624 
625  else
626  {
627  M_matrixFull->close();
628  M_vectorRhsFull->close();
629  DVLOG(2) << "[Oseen::update] on(u)_Gamma_"
630  << *diriIter << "\n";
631  form2( M_space, M_space, M_matrixFull ) +=
632  on( markedfaces( *M_mesh, *diriIter ),
633  u, *M_vectorRhsFull, g );
634  DVLOG(2) << "[Oseen::update] on(u)_Gamma_"
635  << *diriIter << "done \n";
636  }
637  }
638 
639  //
640  // close global matrix and vector to make them usable by the
641  // solvers
642  //
643  M_matrixFull->close();
644  M_vectorRhsFull->close();
645 
646  // added pressure and viscous terms. \attention addMatrix can be
647  // called only on closed matrices
648  //M_matrixFull->addMatrix( 1.0, M_matrixStab );
649  VLOG(1) << "[Oseen] added stabilisation matrix\n";
650  M_matrixFull->addMatrix( 1.0, M_matrixAu );
651  VLOG(1) << "[Oseen] added standard oseen matrix\n";
652 
653  if ( M_export_matlab )
654  {
655  M_matrixFull->printMatlab( "M_oseen.m" );
656  M_vectorRhsFull->printMatlab( "F_oseen.m" );
657  }
658 
659 } // update
660 
661 template<class Space, uint16_type imOrder,
662  template<uint16_type,uint16_type,uint16_type> class Entity>
663 void Oseen<Space, imOrder, Entity>::solve()
664 {
665  // -- make sure solve is needed
666  if ( !M_updated )
667  {
668  return;
669  }
670 
671  M_updated = false;
672 
673 #if 0
674  // -- solve
675  solveNonSym( M_matrixFull, M_vectorSolFull, M_vectorRhsFull );
676 
677  // -- extract solution
678  std::copy( M_vectorSolFull.begin(),
679  M_vectorSolFull.end(),
680  M_sol.begin() );
681 #else
682 
683  // -- solve
684  solveNonSym( M_matrixFull, M_sol, M_vectorRhsFull );
685 
686 #endif
687 
688 
689 } // Oseen::solve
690 
691 
692 template<class Space, uint16_type imOrder,
693  template<uint16_type,uint16_type,uint16_type> class Entity>
694 void
695 Oseen<Space, imOrder, Entity>::solveNonSym( sparse_matrix_ptrtype const& D,
696  element_type & u,
697  vector_ptrtype const& F )
698 {
699  M_backend->solve( _matrix=D, _solution=u, _rhs=F );
700 } // Oseen::solveNonSym
701 
702 
703 
704 } // Feel

Generated on Fri Oct 25 2013 14:24:22 for Feel++ by doxygen 1.8.4