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
nirb.hpp
1 //
2 // nirb.hpp
3 //
4 
5 
6 #ifndef _nirb_hpp
7 #define _nirb_hpp
8 
9 
10 
11 #include <feel/options.hpp>
12 #include <feel/feelalg/backend.hpp>
13 #include <feel/feeldiscr/functionspace.hpp>
16 #include <feel/feeldiscr/mesh.hpp>
17 #include <feel/feelfilters/exporter.hpp>
18 #include <feel/feelvf/vf.hpp>
19 
20 //#include <feel/feelalg/solvereigen.hpp>
21 #include <cstdlib>
22 #include <string>
23 #include <sstream>
24 #include <fstream>
25 #include <iostream>
26 
27 #include <Eigen/Core>
28 #include <Eigen/LU>
29 #include <Eigen/Dense>
30 #include <Eigen/Eigenvalues>
31 
32 #include <vector>
33 #include <algorithm>
34 
35 #include<boost/range/algorithm/max_element.hpp>
36 #include<boost/filesystem.hpp>
37 
39 using namespace Feel;
40 using namespace Feel::vf;
41 
42 using namespace std;
43 
44 template<int PolynomialOrder> class NIRBTEST
45  :
46 public Simget
47 {
48  typedef Simget super;
49 public:
50 
51  //static const uint16_type Order = 3;
52 
54  typedef double value_type;
55 
59  typedef boost::shared_ptr<backend_type> backend_ptrtype;
61  typedef typename backend_type::vector_type vector_type;
63  typedef typename backend_type::vector_ptrtype vector_ptrtype;
65  typedef typename backend_type::sparse_matrix_type sparse_matrix_type;
67  typedef typename backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
68 
69 
71  typedef Simplex<2> convex_type;
73  typedef Mesh<convex_type> mesh_type;
75  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
76 
77 
79  typedef bases<Lagrange<PolynomialOrder,Scalar> > basis_type;
80 
82  typedef FunctionSpace<mesh_type, basis_type> space_type;
84  typedef boost::shared_ptr<space_type> space_ptrtype;
85  typedef typename space_type::element_type element_type;
86 
88  typedef Exporter<mesh_type> export_type;
90  typedef boost::shared_ptr<export_type> export_ptrtype;
91  typedef std::vector<element_type> vector_of_element_type;
92 
96  NIRBTEST()
97  :
98  super(),
99  M_backend( backend_type::build( this->vm() ) ),
100  CoarseMeshSize( this->vm()["hcoarsesize"].template as<double>() ),
101  FineMeshSize( this->vm()["hfinsize"].template as<double>() ),
102  ReadingMeshes( this->vm()["ReadingMeshes"].template as<int>() ),
103  NbSnapshot( this->vm()["NbSnapshot"].template as<int>() ),
104  sizeRB( this->vm()["sizeRB"].template as<int>() ),
105  muMin( this->vm()["muMin"].template as<double>() ),
106  muMax( this->vm()["muMax"].template as<double>() ),
107  mu( this->vm()["mu"].template as<double>() ),
108  Sampling( this->vm()["Sampling"].template as<int>() ),
109  SamplingCoarse( this->vm()["SamplingCoarse"].template as<int>() ),
110  Offline( this->vm()["Offline"].template as<int>() ),
111  ComputeError( this->vm()["ComputeError"].template as<int>() )
112  {
113  }
114 
115  void run();
116 
117  void run( const double* X, unsigned long P, double* Y, unsigned long N );
118 
119 
120  // computation of the "blackbox" solution
121  element_type blackbox( space_ptrtype Xh, double param );
122 
123 
124  //Construction of NIRB basis
125 
126  void ComputeSnapshot( space_ptrtype Xh, std::string filename );
127 
128  void ChooseRBFunction( space_ptrtype Xh,
129  vector_of_element_type &VuBasis,
130  sparse_matrix_ptrtype const &MassMatrix,
131  sparse_matrix_ptrtype const &StiffMatrix );
132 
133  void OrthogonalisationRBFunction( space_ptrtype Xh,
134  vector_of_element_type &Vu,
135  vector_of_element_type &M_Vu,
136  sparse_matrix_ptrtype const &StiffMatrix,
137  sparse_matrix_ptrtype const &MassMatrix );
138 
139  void OrthogonalisationRBFunctionL2GrammSchmidt( space_ptrtype Xh,
140  vector_of_element_type &Vu,
141  sparse_matrix_ptrtype const &MassMatrix );
142 
143  double OrthogonalisationRBFunctionL2GrammSchmidt( space_ptrtype Xh,
144  vector_of_element_type &Vu,
145  int n,
146  sparse_matrix_ptrtype const &MassMatrix );
147 
148 
149  Eigen::MatrixXd ConstructStiffMatrixSnapshot( space_ptrtype Xh
150  ,sparse_matrix_ptrtype const &StiffMatrix );
151 
152  Eigen::MatrixXd ConstructStiffMatrixRB( space_ptrtype Xh,
153  vector_of_element_type &Vu,
154  sparse_matrix_ptrtype const &StiffMatrix );
155 
156  Eigen::MatrixXd ConstructMassMatrixRB( space_ptrtype Xh,
157  vector_of_element_type &Vu,
158  sparse_matrix_ptrtype const &MassMatrix );
159 
160  void ConstructNIRB ( space_ptrtype Xh,
161  vector_of_element_type &M_VNirbBasis,
162  vector_of_element_type &VNirbBasis );
163 
164 
165  //Construction of NIRB solution
166  element_type BuildNirbSolutionWithoutPostProcess( space_ptrtype XhFine,
167  element_type uCoarseInterpolate,
168  vector_of_element_type const &M_VNirbBasis,
169  vector_of_element_type const &VNirbBasis,
170  Eigen::VectorXd & BetaiH );
171  element_type BuildCoarseInterpolation( space_ptrtype XhFine,
172  space_ptrtype XhCoarse,
173  double param );
174 
175  element_type BuildNirbSolution( space_ptrtype XhFine,
176  space_ptrtype XhCoarse,
177  vector_of_element_type const &M_VNirbBasis,
178  vector_of_element_type const &VNirbBasis,
179  double param,
180  Eigen::VectorXd BetaiH );
181 
182  //Post-processing of NIRB solution
183  Eigen::MatrixXd BuildBetaH( space_ptrtype XhFine,
184  space_ptrtype XhCoarse,
185  vector_of_element_type &M_VNirbBasis );
186 
187  Eigen::MatrixXd BuildBetah( space_ptrtype XhFine,
188  vector_of_element_type &M_VNirbBasis );
189 
190  Eigen::MatrixXd BuildPostProcessMatrix( space_ptrtype XhFine,
191  space_ptrtype XhCoarse,
192  vector_of_element_type &M_VNirbBasis );
193 
194  element_type BuildNirbSolutionWithPostProcess( space_ptrtype XhFine,
195  space_ptrtype XhCoarse,
196  vector_of_element_type &M_VNirbBasis,
197  vector_of_element_type &VNirbBasis,
198  Eigen::VectorXd &BetaiuH,
199  element_type &uCoarseInterpolate );
200 
201 
202  //Create geometry
203  gmsh_ptrtype createGeo( double hsize,std::string MeshFileName );
204 
205 private:
206 
208  backend_ptrtype M_backend;
209 
211  double CoarseMeshSize;
212  double FineMeshSize;
213  int ReadingMeshes; // By Default = 0
214 
216  int NbSnapshot,sizeRB;
217  double muMin,muMax,mu;
218 
219  // Paramater to set OFF or ON the "offline" procedure
220 
221  int Sampling;
222  int SamplingCoarse;
223  // if sampling = 1 ==> compute sampling (by default)
224  // if sampling = 0 ==> read sampling from a file
225 
226  int Offline;
227  // if Offline = 1 then ConstructNIRB is called
228  // if Offline = 0 then ConstructNIRB is not called, the Nirb Basis function are read from a file
229 
230  // Paramater to set OFF or ON the computation of the error between the F.E and the NIRB Methods
231  int ComputeError;
232  //if ComputeError = 1 == ON
233  //if ComputeError = 0 == OFF
234 
235  int polynomialOrder;
236 }; // NIRBTEST
237 
238 
239 //-----------------------------------------
240 //-----------------------------------------
241 //const uint16_type NIRBTEST::Order;
242 template< int PolynomialOrder>
243 void NIRBTEST<PolynomialOrder>::run()
244 {
245 
246  std::cout << "------------------------------------------------------------\n";
247  std::cout << "Execute NIRBTEST<" << PolynomialOrder << ">\n";
248  std::vector<double> X( 12 );
249  X[0] = FineMeshSize;
250  X[1] = CoarseMeshSize;
251  X[2] = ReadingMeshes;
252  X[3] = NbSnapshot;
253  X[4] = sizeRB;
254  X[5] = muMin;
255  X[6] = muMax;
256  X[7] = mu;
257  X[8] = Sampling;
258  X[9] = SamplingCoarse;
259  X[10] = Offline;
260  X[11] = ComputeError;
261  std::vector<double> Y( 1 );
262  run( X.data(), X.size(), Y.data(), Y.size() );
263 }
264 //-----------------------------------------
265 //-----------------------------------------
266 template< int PolynomialOrder>
267 void NIRBTEST<PolynomialOrder>::run( const double* X, unsigned long P, double* Y, unsigned long N )
268 {
269 
270  std::cout<< "Size reduced basis : " << sizeRB <<std::endl;
271  std::cout<< "Polynomial order : P" << PolynomialOrder <<std::endl;
272 
273  if ( !this->vm().count( "nochdir" ) )
274  Environment::changeRepository( boost::format( "%1%/P%2%/" )
275  % this->about().appName()
276  % PolynomialOrder );
277 
278 
279  mesh_ptrtype meshExtraCoarse, meshCoarse, meshFine ;
280 
281  if ( ReadingMeshes )
282  {
283  meshCoarse = loadGMSHMesh( _mesh=new mesh_type,
284  _filename="/home/chakir/Mesh/Mesh_H1.msh",
285  _update=MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
286  meshFine = loadGMSHMesh( _mesh=new mesh_type,
287  _filename="/home/chakir/Mesh/Mesh_H2.msh",
288  _update=MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
289  std::cout << "Meshes read in file : done " <<std::endl;
290  }
291 
292  else
293  {
294  std::cout << "Meshes build using 'createGMSHMesh' " <<std::endl;
295  std::cout << "Coarse mesh size : " << X[1] << std::endl;
296  meshCoarse = createGMSHMesh( _mesh= new mesh_type,
297  _desc = createGeo( X[1],"Mesh_Coarse" ),
298  _refine=1 );
299 
300  if ( X[1]> X[0] )
301  {
302 
303  std::cout << "Fine mesh size : " << X[0] << std::endl;
304  meshFine = createGMSHMesh( _mesh= new mesh_type,
305  _desc = createGeo( X[0],"Mesh_Fine" ),
306  _refine=1 );
307  }
308 
309  else
310  {
311 
312  std::cout << "WARNING : 'hfinesize > 'hcoarsesize' --> Fine mesh size set to hcoarsesize/2 : " << X[1]/2 << std::endl;
313  meshFine = createGMSHMesh( _mesh= new mesh_type,
314  _desc = createGeo( X[1]/2,"Mesh_Fine" ),
315  _refine=1 );
316  }
317 
318  }
319 
320 
321 
322  space_ptrtype XhFine = space_type::New( meshFine );
323  //Fine F.E space used to build the NIRB basis
324 
325  space_ptrtype XhCoarse = space_type::New( meshCoarse );
326 
327 
328 
329  //STEP ONE : Construction of the "non intruisive reduced basis (nirb) functions"
330  //Computation of the X[3] snapshots solution on Xhfine
331  //Selection of X[2] fonctions to build the "reduced basis" using a POD technique
332  //Orthogonalisation in L2 and H1 norm of "reduced basis function"
333  //Save this final functions
334 
335 
336  //STEP ONE : Construction of the "non intruisive reduced basis (nirb) functions"
337  //Computation of the X[3] snapshots solution on Xhfine
338  //Selection of X[2] fonctions to build the "reduced basis" using a POD technique
339  //Orthogonalisation in L2 and H1 norm of "reduced basis function"
340  //Save this final functions
341 
342 
343  auto ui = XhFine->element();
344  ui.zero();
345  vector_of_element_type VNirbBasis( sizeRB,ui );
346  vector_of_element_type MassMat_x_VuNirb( sizeRB,ui );
347 
348  if ( !Offline )
349  {
350  std :: string path;
351 
352  //Checking if the NIRB basis files exist
353  //WARNING : We also must check if the basis have been computed on the right FE space (Not done yet)
354  for ( int i = 0 ; i < sizeRB; i++ )
355  {
356  path = ( boost::format( "./RB%1%NIRB_BasisFile_%2%" )%sizeRB%i ).str();
357  ui.load( _path=path );
358 
359  if ( ui.l2Norm()== 0. )
360  {
361  std::cout <<"WARNING Error : RB " << sizeRB << "NIRB_BasisFile_" << i << " is equal to zero => OFFline parameter set to 1 " <<std::endl;
362  Offline = 1;
363  break;
364  }
365 
366  VNirbBasis[i] = ui;
367  path =( boost::format( "./RB%1%D_Nirb_Basis%2%File" ) %sizeRB%i ).str() ;
368  ui.load( _path=path );
369 
370  if ( ui.l2Norm()== 0. )
371  {
372  std::cout <<"WARNING Error : RB" << sizeRB << "D_NIRB_Basis" << i << "File is equal to zero => Re-calculating M*Ui" <<std::endl;
373  auto Mui = M_backend->newVector( XhFine );
374  auto Ui = M_backend->newVector( XhFine );
375 
376  auto uj = XhFine->element();
377 
378  sparse_matrix_ptrtype MassMatrix = M_backend->newMatrix( _test=XhFine, _trial=XhFine ); //Sparse FE mass matrix
379  form2( _test=XhFine , _trial=XhFine, _matrix= MassMatrix ) =
380  integrate( _range=elements( XhFine->mesh() ), _expr=id( uj )*idt( ui ) );
381 
382  *Ui = VNirbBasis[i];
383  Ui->close();
384  MassMatrix->multVector( Ui,Mui );
385  MassMat_x_VuNirb[i] = *Mui;
386  }
387 
388  else
389  {
390  MassMat_x_VuNirb[i] = ui;
391  }
392  }
393  }
394 
395  if ( Offline )
396  {
397  ConstructNIRB( XhFine,MassMat_x_VuNirb,VNirbBasis );
398 
399  if ( SamplingCoarse )
400  {
401  std::cout << "Computation of Coarse snapshot " <<std::endl;
402  boost::timer ti;
403  ComputeSnapshot( XhCoarse,"_Coarse_" );
404 
405  double Time_snapshot_Coarse = ti.elapsed();
406  std::cout << "Computation of the " << NbSnapshot << " snapshots : done -- ";
407  std::cout << "Time per snapshot: " << Time_snapshot_Coarse/NbSnapshot << " sec " <<std::endl;
408  SamplingCoarse = 0;
409  }
410  }
411 
412  //STEP TWO : Approximation of the solution using the "nirb" functions for a choosen mu
413 
414 
415 
416 
417  double p = mu;
418 
419 
420  std :: cout << "OFFLINE Procedure has been skipped" <<std::endl;
421 
422 
423  auto uCoarseInterpolate = XhFine->element();
424  auto uNirbCoarse = XhFine->element(); // Fine/Coarse Grid NIRB solution
425  Eigen::VectorXd BetaiH( sizeRB );
426  double TimeCoarse,TimeCoarsePostProcess,TimeFine;
427  boost::timer ti;
428  ti.restart();
429  uCoarseInterpolate = BuildCoarseInterpolation( XhFine,XhCoarse,p );
430  TimeCoarse = ti.elapsed();
431  std :: cout << "Calculation of uCoarse and it's interpolation on XhFine :" << TimeCoarse << " sec" <<std::endl;
432  uNirbCoarse = BuildNirbSolutionWithoutPostProcess( XhFine,uCoarseInterpolate,MassMat_x_VuNirb,VNirbBasis,BetaiH );
433  //uNirbCoarse = BuildNirbSolution(XhFine,XhCoarse,MassMat_x_VuNirb,VNirbBasis,p,BetaiH);
434  TimeCoarse = ti.elapsed();
435  std :: cout << "Construction of NIRB solution (uNirbCoarse) - Fine/Coarse Grid (saved in nirb2GridCoarse):" <<std::endl;
436  std::cout << "Time to build solution " << TimeCoarse << " sec" <<std::endl;
437 
438  export_ptrtype exporter2GridCoarse( export_type::New( this->vm(), "nirb2GridCoarse" ) );
439  exporter2GridCoarse->step( 0 )->setMesh( meshFine );
440  exporter2GridCoarse->step( 0 )->add( "uNirbCoarse",uNirbCoarse );
441  exporter2GridCoarse->save();
442 
443  auto uNirbCoarsePostProcess = XhFine->element(); // Fine/Coarse Grid NIRB solution
444  ti.restart();
445  if (SamplingCoarse)
446  {
447  std::cout << "Computation of Coarse snapshot " <<std::endl;
448  boost::timer ti;
449  ComputeSnapshot( XhCoarse,"_Coarse_" );
450 
451  double Time_snapshot_Coarse = ti.elapsed();
452  std::cout << "Computation of the " << NbSnapshot << " snapshots : done -- ";
453  std::cout << "Time per snapshot: " << Time_snapshot_Coarse/NbSnapshot << " sec " <<std::endl;
454  }
455  uNirbCoarsePostProcess = BuildNirbSolutionWithPostProcess( XhFine, XhCoarse, MassMat_x_VuNirb,
456  VNirbBasis,BetaiH,uCoarseInterpolate );
457 
458  TimeCoarsePostProcess = ti.elapsed();
459  std :: cout << "Post-processing of NIRB solution (uNirbCoarsePostProcess) - Fine/Coarse Grid (saved in nirb2GridCoarse):" <<std::endl;
460  std::cout << "Time " << TimeCoarsePostProcess << " sec" <<std::endl;
461 
462  export_ptrtype exporter2GridCoarsePostProcess( export_type::New( this->vm(), "nirb2GridCoarsePostProcess" ) );
463  exporter2GridCoarsePostProcess->step( 0 )->setMesh( meshFine );
464  exporter2GridCoarsePostProcess->step( 0 )->add( "uNirbCoarsePostProcess",uNirbCoarsePostProcess );
465  exporter2GridCoarsePostProcess->save();
466 
467  if ( ComputeError )
468  {
469  std :: cout << "Error calculation " <<std::endl;
470  Eigen::VectorXd Betaih( sizeRB );
471 
472  auto u1Grid = XhFine->element();
473  ti.restart();
474  u1Grid = BuildNirbSolution( XhFine,XhFine,MassMat_x_VuNirb,VNirbBasis,p,Betaih );
475  TimeFine = ti.elapsed();
476  std::cout << "Construction of uNirbFine - Fine/Fine Grid (saved in nirb1Grid): done "<<std::endl;
477  std::cout << "Time to build " << TimeFine << " sec" <<std::endl;
478 
479 
480 
481  mesh_ptrtype meshRef;
482 
483  if ( ReadingMeshes )
484  {
485  meshRef = loadGMSHMesh( _mesh=new mesh_type,
486  _filename="/home/chakir/Mesh/Mesh_H3.msh",
487  _update=MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
488  }
489 
490  else
491  {
492  if ( X[1] > X[0] )
493  {
494  meshRef = createGMSHMesh( _mesh= new mesh_type,
495  _desc = createGeo( X[0]/2,"Mesh_Fine" ),
496  _refine=1 );
497  }
498 
499  else
500  {
501  meshRef = createGMSHMesh( _mesh= new mesh_type,
502  _desc = createGeo( X[1]/4,"Mesh_Fine" ),
503  _refine=1 );
504  }
505  }
506 
507 
508 
509  space_ptrtype XhRef = space_type::New( meshRef );
510 
511 
512 
513  //Reference F.E space used to compute the "reference" solution
514  //which supposed to be accurate enough
515 
516 
517  auto uRef = XhRef->element();
518 
519  uRef = blackbox( XhRef, p );
520  std :: cout << "Computation of FE solution (uRef) - Ref Grid (not saved): done" <<std::endl;
521 
522 
523  //Computation of the H1 norm of uRef
524  double L2NormUref = std :: sqrt( integrate( _range=elements( XhRef->mesh() ),
525  _expr=( idv( uRef )*idv( uRef ) ) ).evaluate()( 0,0 ) );
526  double SemiH1NormUref = std :: sqrt( integrate( _range=elements( XhRef->mesh() ),
527  _expr=( gradv( uRef )*trans( gradv( uRef ) ) ) ).evaluate()( 0,0 ) );
528  double H1NormUref = std :: sqrt( SemiH1NormUref *SemiH1NormUref + L2NormUref*L2NormUref );
529 
530  //std::cout << "L2NormURef = " << L2NormUref << " -- H1NormUref = " << H1NormUref <<std::endl;
531  //std::cout <<std::endl;
532 
533 
534 
535  // //Computation of relative error measured in H1 norm
536 
537 
538  // // auto DFine = M_backend->newMatrix( _test=XhRef, _trial=XhFine );//Sparse F.E mass matrix D
539  // // form2( _test=XhRef, _trial=XhFine, _matrix=DFine ) =
540  // // integrate( _range=elements(XhFine->mesh()), _expr=idt(uRef)*id(uFine));
541  // // auto DCoarse = M_backend->newMatrix( _test=XhRef, _trial=XhCoarse );//Sparse F.E mass matrix D
542  // // form2( _test=XhRef, _trial=XhCoarse, _matrix=DCoarse ) =
543  // // integrate( _range=elements(XhFine->mesh()), _expr=idt(uRef)*id(uCoarse));
544 
545 
546 
547 
548  double ErrH1uNirbCoarse,ErrH1uNirbCoarsePostProcess,ErrH1u1Grid;
549 
550 
551 
552  // uRef - uNirbCoarse
553  ErrH1uNirbCoarsePostProcess = integrate( _range=elements( XhRef->mesh() ),
554  _expr=( ( gradv( uRef )-gradv( uNirbCoarsePostProcess ) )*trans( gradv( uRef )-gradv( uNirbCoarsePostProcess ) )
555  + ( idv( uRef )-idv( uNirbCoarsePostProcess ) )*( idv( uRef )-idv( uNirbCoarsePostProcess ) ) ) ).evaluate()( 0,0 );
556  ErrH1uNirbCoarsePostProcess = sqrt( ErrH1uNirbCoarsePostProcess )/H1NormUref;
557 
558 
559  // uRef - uNirbCoarse
560  ErrH1uNirbCoarse = integrate( _range=elements( XhRef->mesh() ),
561  _expr=( ( gradv( uRef )-gradv( uNirbCoarse ) )*trans( gradv( uRef )-gradv( uNirbCoarse ) )
562  + ( idv( uRef )-idv( uNirbCoarse ) )*( idv( uRef )-idv( uNirbCoarse ) ) ) ).evaluate()( 0,0 );
563  ErrH1uNirbCoarse = sqrt( ErrH1uNirbCoarse )/H1NormUref;
564 
565  //uRef - uNirbFine (= uRef - u1Grid)
566  ErrH1u1Grid = integrate( _range=elements( XhRef->mesh() ),
567  _expr=( ( gradv( uRef )-gradv( u1Grid ) )*trans( gradv( uRef )-gradv( u1Grid ) )
568  + ( idv( uRef )-idv( u1Grid ) )*( idv( uRef )-idv( u1Grid ) ) ) ).evaluate()( 0,0 );
569  ErrH1u1Grid = sqrt( ErrH1u1Grid )/H1NormUref;
570 
571  std :: cout << "H1-norm NIRB Error " <<std::endl;
572  std :: cout << "||u_ref - u_NirbCoarse ||_{H1} = " << ErrH1uNirbCoarse <<std::endl;
573  std :: cout << "||u_ref - u_NirbCoarsePostProcess ||_{H1} = " << ErrH1uNirbCoarsePostProcess <<std::endl;
574  std :: cout << "||u_ref - u_NirbFine ||_{H1} = " << ErrH1u1Grid<<std::endl;
575  std :: cout << std::endl;
576 
577 
578 
579 
580 
581  export_ptrtype exporter1Grid( export_type::New( this->vm(), "nirbFine" ) );
582  exporter1Grid->step( 0 )->setMesh( meshFine );
583  exporter1Grid->step( 0 )->add( "uNirbFine",u1Grid );
584  exporter1Grid->save();
585 
586 
587 
588  auto uFine = XhFine->element();
589  auto uCoarse = XhCoarse->element();
590 
591  uCoarse = blackbox( XhCoarse, p );
592  std :: cout << "Computation of FE solution (uCoarse) - Coarse Grid (not saved):"<<std::endl;
593 
594  uFine = blackbox( XhFine, p );
595  std :: cout << "Computation of FE solution (uFine) - Fine Grid (saved in EFFine): done" <<std::endl;
596 
597  /*
598 
599  // uref - uFine
600  double ErrSemiH1uFine,ErrL2uFine,ErrH1uFine;
601 
602 
603  ErrL2uFine = std::sqrt(integrate(_range=elements(XhRef->mesh()),
604  _expr=( (idv(uRef)-idv(uFine))*(idv(uRef)-idv(uFine)) )).evaluate()(0,0));
605 
606 
607  ErrSemiH1uFine = std::sqrt(integrate(_range=elements(XhRef->mesh()),
608  _expr=( (gradv(uRef)-gradv(uFine))*trans(gradv(uRef)-gradv(uFine)))).evaluate()(0,0));
609 
610 
611  ErrH1uFine = std::sqrt(ErrSemiH1uFine *ErrSemiH1uFine + ErrL2uFine*ErrL2uFine);
612  ErrH1uFine = ErrH1uFine/H1NormUref;
613  ErrL2uFine = ErrL2uFine/L2NormUref;
614 
615 
616 
617  // uref - uCoarse
618  double ErrL2uCoarse ,ErrSemiH1uCoarse ,ErrH1uCoarse;
619  ErrL2uCoarse = std::sqrt(integrate(_range=elements(XhRef->mesh()),
620  _expr=( (idv(uRef)-idv(uCoarse))*(idv(uRef)-idv(uCoarse)) )).evaluate()(0,0));
621 
622 
623  ErrSemiH1uCoarse = std::sqrt(integrate(_range=elements(XhRef->mesh()),
624  _expr=( (gradv(uRef)-gradv(uCoarse))*trans(gradv(uRef)-gradv(uCoarse)))).evaluate()(0,0));
625 
626  ErrH1uCoarse = std::sqrt(ErrSemiH1uCoarse * ErrSemiH1uCoarse + ErrL2uCoarse*ErrL2uCoarse);
627  ErrH1uCoarse = ErrH1uCoarse/H1NormUref;
628  ErrL2uCoarse = ErrL2uCoarse/L2NormUref;
629 
630 
631  std :: cout << "H1-Norm F.E Error " <<std::endl;
632  std :: cout << "||u_ref - u_coarse ||_{H1} = " << ErrH1uCoarse <<std::endl;
633  std :: cout << "||u_ref - u_fine ||_{H1} = " << ErrH1uFine <<std::endl;
634  std :: cout <<std::endl;
635 
636  */
637  export_ptrtype exporterFine( export_type::New( this->vm(), "EF_Fine" ) );
638  export_ptrtype exporterCoarse( export_type::New( this->vm(), "EF_Coarse" ) );
639  exporterFine->step( 0 )->setMesh( meshFine );
640  exporterCoarse->step( 0 )->setMesh( meshCoarse );
641 
642  exporterFine->step( 0 )->add( "uFine", uFine );
643  exporterCoarse->step( 0 )->add( "uCoarse", uCoarse );
644 
645  exporterFine->save();
646  exporterCoarse->save();
647 
648  }
649 
650 
651 
652 
653 
654 
655 
656 }
657 
658 
659 //-----------------------------------------
660 //-----------------------------------------
661 
662 template< int PolynomialOrder>
663 gmsh_ptrtype NIRBTEST<PolynomialOrder>::createGeo( double hsize, std::string MeshFileName )
664 {
665  double hCornersize1 = hsize/15.;
666  double hCornersize2 = hsize/30.;
667  std::ostringstream ostr;
668  ostr << "Mesh.MshFileVersion = 2.2;" << "\n";
669  ostr << "Mesh.CharacteristicLengthExtendFromBoundary=1;" <<std::endl;
670  ostr << "Mesh.CharacteristicLengthFromPoints=1;" <<std::endl;
671  ostr << "Mesh.ElementOrder=1;" <<std::endl;
672  ostr << "Mesh.SecondOrderIncomplete = 0;" <<std::endl;
673  ostr << "Mesh.Algorithm = 6;" <<std::endl;
674  ostr << "Mesh.OptimizeNetgen=1;" <<std::endl;
675  ostr << "// partitioning data" <<std::endl;
676  ostr << "Mesh.Partitioner=1;" <<std::endl;
677  ostr << "Mesh.NbPartitions=1;" <<std::endl;
678  ostr << "Mesh.MshFilePartitioned=0;" <<std::endl;
679  ostr << "Point(1) = {0,0,0.0," << hsize << "};" <<std::endl;
680  ostr << "Point(2) = {1,0,0.0," << hCornersize2 << "};" <<std::endl;
681  ostr << "Point(3) = {1,1,0.0," << hCornersize1 << "};"<<std::endl;
682  ostr << "Point(4) = {0,1,0.0," << hCornersize2 << "};"<<std::endl;
683  ostr << "Line(1) = {4,1};" <<std::endl;
684  ostr << "Line(2) = {1,2};" <<std::endl;
685  ostr << "Line(3) = {2,3};" <<std::endl;
686  ostr << "Line(4) = {3,4};" <<std::endl;
687  ostr << "Line Loop(5) = {1,2,3,4};" <<std::endl;
688  ostr << "Plane Surface(6) = {5};" <<std::endl;
689  ostr << "Physical Line(\"Dirichlet1\") = {1,2};" <<std::endl;
690  ostr << "Physical Line(\"Dirichlet2\") = {3};" <<std::endl;
691  ostr << "Physical Line(\"Dirichlet3\") = {4};" <<std::endl;
692  ostr << "Physical Surface(\"Mat1\") = {6};" <<std::endl;
693  std::ostringstream nameStr;
694  nameStr.precision( 3 );
695  nameStr << MeshFileName;
696  //nameStr << MeshFileName.c_str();
697  gmsh_ptrtype gmshp( new Gmsh );
698  gmshp->setPrefix( nameStr.str() );
699  gmshp->setDescription( ostr.str() );
700  return gmshp;
701 
702 }
703 
704 //-----------------------------------------
705 //-----------------------------------------
706 
707 template< int PolynomialOrder>
708 void NIRBTEST<PolynomialOrder>::ComputeSnapshot( space_ptrtype Xh,std::string filename )
709 {
710 
711  int XhNdof = Xh->nLocalDof();
712  double pi = 3.14159265;
713 
714  auto ui = Xh->element();
715 
716  for ( int i =0; i<NbSnapshot; i++ )
717  {
718  double theta = ( i+1 )*( muMax-muMin )/( 113. );
719  ui = blackbox( Xh,theta );
720  if (ui.l2Norm()==0.)
721  {
722  std::cerr << "In ComputeSnapshot : ERROR IN USING BLACKBOX -> ui = 0" << endl;
723  exit(0);
724  }
725  std::string path = "./Sol" + filename + ( boost::format( "%1%" ) %i ).str() ;
726  ui.save( _path=path );
727  //boost::filesystem::path full_path_Sol( boost::filesystem::current_path() );
728  //cout << "Path where the sol are save " << full_path_Sol << endl;
729 
730  }
731 
732 }
733 
734 //-----------------------------------------
735 //-----------------------------------------
736 template< int PolynomialOrder>
737 Eigen::MatrixXd NIRBTEST<PolynomialOrder> :: ConstructStiffMatrixSnapshot( space_ptrtype Xh,sparse_matrix_ptrtype const & StiffMatrix )
738 {
739 
740 
741  auto ui = Xh->element();
742  auto uj = Xh->element();
743  /*
744  auto StiffMatrix = M_backend->newMatrix( _test=Xh, _trial=Xh );//Sparse F.E stiffness matrix D
745  form2( _test=Xh, _trial=Xh, _matrix=D ) =
746  integrate( _range=elements(Xh->mesh()), _expr=gradt(ui)*trans(grad(uj)));
747  */
748  Eigen::MatrixXd S ( NbSnapshot,NbSnapshot ); //Dense RB stiffness matrix S
749  cout << " Dans ConstructStiffMatrixSnapshot " << endl;
750  //boost::filesystem::path full_path_Sol( boost::filesystem::current_path() );
751  //cout << "Path where the sol are load " << full_path_Sol << endl;
752 
753  for ( int i = 0; i< NbSnapshot; i++ )
754  {
755  ui.zero();
756  std::string path = ( boost::format( "./Sol_%1%" ) %i ).str() ;
757  //std::string suffixe = "-";
758  //ui.load( _path=path,_suffix=suffixe );
759 
760  ui.load( _path=path);
761 
762  if ( ui.l2Norm()==0. )
763  {
764  std::cerr << "In 'ConstructStiffMatrixSnapshot': ERROR IN LOADING FILE " << path <<std::endl;
765  auto t = Xh->element();
766  t.load(_path=path);
767  cout << "Reloading of ui done -> L2norm(ui)= " << t.l2Norm() << endl;
768  exit( 0 );
769  }
770 
771  for ( int j=0; j< i; j++ )
772  {
773  uj.zero();
774  path = ( boost::format( "./Sol_%1%" ) %j ).str() ;
775  uj.load( _path=path );
776 
777  if ( uj.l2Norm()==0. )
778  {
779  std::cerr << "In 'ConstructStiffMatrixSnapshot': ERROR IN LOADING FILE " << path <<std::endl;
780  exit( 0 );
781  }
782 
783  double Sij = StiffMatrix->energy( ui,uj );
784  S( i,j ) = Sij;
785  S( j,i ) = Sij;
786  }
787 
788  double Sii = StiffMatrix->energy( ui,ui );
789  S( i,i ) = Sii;
790  }
791 
792  std::string filename = ( boost::format( "StiffMatrixP%1%" )%PolynomialOrder ).str();
793  std::ofstream fileMat( filename.c_str() );
794  fileMat << NbSnapshot <<std::endl;
795 
796  for ( int i = 0; i<NbSnapshot; i++ )
797  {
798  for ( int j =0; j<NbSnapshot; j++ )
799  {
800  fileMat << S( i,j ) <<std::endl;
801  }
802  }
803 
804  fileMat.close();
805 
806  return S;
807 }
808 
809 
810 //--------------------------------------------------------------
811 //--------------------------------------------------------------
812 template< int PolynomialOrder>
813 Eigen::MatrixXd NIRBTEST<PolynomialOrder> :: ConstructStiffMatrixRB( space_ptrtype Xh,vector_of_element_type &Vu,sparse_matrix_ptrtype const & StiffMatrix )
814 {
815 
816  Eigen::MatrixXd S ( sizeRB,sizeRB ); //Dense RB stiffness matrix
817 
818  for ( int i = 0; i< sizeRB; i++ )
819  {
820  for ( int j=0; j< i; j++ )
821  {
822  double Sij = StiffMatrix->energy( Vu[i],Vu[j] );
823  S( i,j ) = Sij;
824  S( j,i ) = Sij;
825  }
826 
827  double Sii = StiffMatrix->energy( Vu[i],Vu[i] );
828  S( i,i ) = Sii;
829  }
830 
831  return S;
832 }
833 //--------------------------------------------------------------
834 //--------------------------------------------------------------
835 template< int PolynomialOrder>
836 Eigen::MatrixXd NIRBTEST<PolynomialOrder> :: ConstructMassMatrixRB( space_ptrtype Xh,vector_of_element_type &Vu,sparse_matrix_ptrtype const &MassMatrix )
837 {
838 
839  Eigen::MatrixXd S ( sizeRB,sizeRB ); //Dense RB mass matrix
840 
841  for ( int i = 0; i< sizeRB; i++ )
842  {
843  for ( int j=0; j< i; j++ )
844  {
845  double Sij = MassMatrix->energy( Vu[i],Vu[j] );
846  S( i,j ) = Sij;
847  S( j,i ) = Sij;
848  }
849 
850  double Sii = MassMatrix->energy( Vu[i],Vu[i] );
851  S( i,i ) = Sii;
852  }
853 
854  return S;
855 
856 }
857 
858 //-----------------------------------------
859 //-----------------------------------------
860 template< int PolynomialOrder>
861 void NIRBTEST<PolynomialOrder> ::ChooseRBFunction( space_ptrtype Xh,vector_of_element_type &VuBasis,sparse_matrix_ptrtype const &MassMatrix,sparse_matrix_ptrtype const &StiffMatrix )
862 {
863 
864  //std :: cout << "Inside ChooseRBFunction" <<std::endl;
865  Eigen::MatrixXd S ( NbSnapshot,NbSnapshot ); //Dense Stiffness Matrix
866 
867  if ( !Sampling )
868  {
869  std::string str = ( boost::format( "StiffMatrixP%1%" )%PolynomialOrder ).str();
870  std::ifstream fileMat( str.c_str() );
871 
872  if ( !fileMat )
873  {
874  std::cerr << "Error in reading file " << ( boost::format( "StiffMatrixP%1%" )%PolynomialOrder ).str() << "does not exist " << std::endl;
875  std::cout << "Need to re-build the Stiffness Matrix -> SAMPLING SET TO 1 " << std::endl;
876  Sampling = 1;
877  }
878 
879  else
880  {
881  int Itemp;
882  fileMat >> Itemp;
883 
884  if ( Itemp != NbSnapshot )
885  {
886  std :: cerr << "Error in reading file" << ( boost::format( "StiffMatrixP%1%" )%PolynomialOrder ).str() << " -- " << Itemp << " != NbSnapshot(" << NbSnapshot << ")" << std::endl;
887  std::cout << "Need to re-build the Stiffness Matrix -> SAMPLING SET TO 1 " << std::endl;
888  Sampling = 1;
889  }
890 
891  else
892  {
893  for ( int i=0; i<NbSnapshot; i++ )
894  {
895  for ( int j =0; j<NbSnapshot; j++ )
896  {
897  double Dtemp;
898  fileMat >> Dtemp;
899  S( i,j ) = Dtemp;
900  }
901  }
902  }
903  }
904 
905  }
906 
907  if ( Sampling )
908  {
909  S = ConstructStiffMatrixSnapshot( Xh,StiffMatrix );
910  }
911 
912 
913  Eigen::EigenSolver <Eigen::MatrixXd> eigen_solver( S );
914 
915  //int nb_EigenValue = eigen_solver.eigenvalues().size();
916  //eigen_solver.eigenvectors().col(i) = eigenvector #i
917 
918  std :: cout << "Computation of the eigenvalues of the stiffness matrix S (NbSnapshot,NbSnapshot) : done" <<std::endl;
919 
920  //Sorting N = "SizeBR" F.E solutions uh(mu_k) the more represented in the eigenvectors
921  //associated to the N largest eigenvalues
922 
923 
924 
925 
926  // saving the index's number to identify the NIRB basis functions
927  std::string path = ( boost::format( "./IndBR%1%" ) %sizeRB ).str() ;
928  std::ofstream find( path.c_str() );
929 
930  if ( !find )
931  {
932  std :: cerr <<" 'ChooseRBFunction routine' - ERROR IN OPENING FILE:" << path <<std::endl;
933  exit( 0 );
934  }
935 
936  std::vector<int> Uind( NbSnapshot );
937 
938  for ( int i=0; i<NbSnapshot; i++ )
939  {
940  Uind[i] = 0; //Set to zero if u(Uind[i]) is not NIRB basis function.
941  }
942 
943 
944  auto ui = Xh->element();
945 
946 
947  double Tol = 1e-30;
948  //std::ofstream fout("fout");
949  //Sorting eigenvector #i
950  std::vector<double>Vi( NbSnapshot );
951 
952  for ( int i=0; i<sizeRB; i++ )
953  {
954  for ( int ri = 0; ri < NbSnapshot; ri++ )
955  {
956  //If u(Uind[ri]) is a NIRB basis function, the "Uind[ri]" component of Vi is set to zero
957  //to avoid it.
958  if ( Uind[ri] == 0 )
959  {
960  Vi[ri] = std::abs( real( eigen_solver.eigenvectors().col( i )[ri] ) );
961  }
962 
963  else
964  {
965  Vi[ri] = 0;
966  }
967  }
968 
969  //Using max_element to find the position of the maximum component of Vi
970 
971  int IndMax = std::distance( Vi.begin(),boost::max_element( Vi ) );
972 
973  //copy u[IndMax] in Mu
974  ui.zero();
975  path = ( boost::format( "./Sol_%1%" ) %IndMax ).str();
976  ui.load( _path=path );
977 
978  //cout << "LOAD DONE " <<std::endl;
979  if ( ui.l2Norm()==0. )
980  {
981  std::cerr <<"'ChooseRBFunction routine': ERROR IN LOADING FILE " << path <<std::endl;
982  exit( 0 );
983  }
984 
985  VuBasis[i]= ui;
986 
987  if ( i>0 )
988  {
989  double normL2 = OrthogonalisationRBFunctionL2GrammSchmidt( Xh,VuBasis,i,MassMatrix );
990 
991  //std :: cout << "i = " << i << " --- NormU(" << IndMax << ") = " << normL2 <<std::endl;
992  //fout << "i = " << i << " --- NormU(" << IndMax << ") = " << normL2 <<std::endl;
993  if ( normL2 > Tol )
994  {
995  Uind[IndMax] = i;
996  find << IndMax <<std::endl;
997  }
998 
999  else
1000  {
1001  Uind[IndMax] = -1;
1002  i--;
1003  }
1004  }
1005 
1006  else
1007  {
1008  Uind[IndMax] = i;
1009  find << IndMax <<std::endl;
1010  }
1011 
1012  }
1013 
1014  find.close();
1015  //fout.close();
1016 
1017 
1018  double L2normVui;
1019 
1020  //save in a file the matrix A
1021  for ( int i=0; i<sizeRB; i++ )
1022  {
1023  L2normVui = std::sqrt( MassMatrix->energy( VuBasis[i],VuBasis[i] ) );
1024  VuBasis[i].scale( 1./L2normVui );
1025  VuBasis[i].save( _path=( boost::format( "./Sol_OR%1%" ) %i ).str() );
1026  }
1027 
1028 
1029 }
1030 
1031 
1032 //-----------------------------------------
1033 //-----------------------------------------
1034 template< int PolynomialOrder>
1035 double NIRBTEST <PolynomialOrder> :: OrthogonalisationRBFunctionL2GrammSchmidt( space_ptrtype Xh, vector_of_element_type &Vu,int n,sparse_matrix_ptrtype const &MassMatrix )
1036 {
1037 
1038  if ( MassMatrix->linftyNorm()==0. )
1039  {
1040  std::cout << "ERROR in OrthogonalisationRBFunctionL2GrammSchmidt : MassMatrix is null " <<std::endl;
1041  exit( 0. );
1042  }
1043 
1044  double Dtemp1,Dtemp2,Dtemp3,Dtemp4;
1045  //We suppose that the n-1 first function has been already orthogonalized
1046 
1047  for ( int k=0; k<n; k++ )
1048  {
1049 
1050  Dtemp1 = MassMatrix->energy( Vu[n],Vu[k] );
1051  Dtemp2 = MassMatrix->energy( Vu[k],Vu[k] );
1052  Dtemp3 = -Dtemp1/Dtemp2;
1053  Vu[n].add( Dtemp3,Vu[k] );
1054 
1055  }
1056 
1057  Dtemp4 = MassMatrix->energy( Vu[n],Vu[n] );
1058  return Dtemp4;
1059 }
1060 
1061 //-----------------------------------------
1062 //-----------------------------------------
1063 
1064 //-----------------------------------------
1065 template< int PolynomialOrder>
1066 void NIRBTEST<PolynomialOrder> :: OrthogonalisationRBFunctionL2GrammSchmidt
1067 ( space_ptrtype Xh,vector_of_element_type &Vu, sparse_matrix_ptrtype const &MassMatrix )
1068 {
1069 
1070  double Dtemp,Dtemp1,Dtemp2,Dtemp3, L2norm;
1071 
1072  for ( int i=1; i<sizeRB; i++ )
1073  {
1074  for ( int k=0; k<i; k++ )
1075  {
1076  Dtemp1 = MassMatrix->energy( Vu[i],Vu[k] );
1077  Dtemp2 = MassMatrix->energy( Vu[k],Vu[k] );
1078  Dtemp3 = - Dtemp1/Dtemp2;
1079  Vu[i].add( Dtemp3,Vu[k] );
1080  }
1081 
1082  }
1083 
1084  for ( int i=0; i<sizeRB; i++ )
1085  {
1086  L2norm = std::sqrt( MassMatrix->energy( Vu[i],Vu[i] ) );
1087  Vu[i].scale( 1./L2norm ) ;
1088  }
1089 
1090 }
1091 
1092 //-----------------------------------------
1093 //-----------------------------------------
1094 template< int PolynomialOrder>
1095 void NIRBTEST<PolynomialOrder> ::OrthogonalisationRBFunction( space_ptrtype Xh,vector_of_element_type &Vu,vector_of_element_type &M_Vu,sparse_matrix_ptrtype const &StiffMatrix,sparse_matrix_ptrtype const &MassMatrix )
1096 {
1097  //First step : L2 -Pre-orthogonalization and normalization using a Gramm-schmidt method
1098 
1099  std::string path;
1100  auto ui = Xh->element();
1101  ui.zero();
1102 
1103  //Second step : Solving a generalized eigenvalue problem
1104  //A*Epsilon_i = \lambda_i*B*Epsilon_i
1105  //A : Dense stiffness matrix in the pre-orthogonalized nirb basis
1106  //B : Dense mass matrix in the pre-orthogonalized nirb basis = Idendity matrix
1107 
1108  //Construction of the matrix A and B
1109 
1110  // M.resize(sizeRB,sizeRB);
1111  Eigen :: MatrixXd A ( sizeRB,sizeRB );
1112  Eigen :: MatrixXd B ( sizeRB,sizeRB );
1113  Eigen :: MatrixXd C ( sizeRB,sizeRB );
1114 
1115  A = ConstructStiffMatrixRB( Xh,Vu,StiffMatrix );
1116  B = ConstructMassMatrixRB( Xh,Vu,MassMatrix );
1117  // Eigen:: GeneralizedSelfAdjointEigenSolver<Eigen::MatrixXd> es(A,B);
1118  C = B.inverse()*A;
1119 
1120  Eigen:: EigenSolver<Eigen::MatrixXd> es;
1121 
1122  es.compute( C );
1123  vector_of_element_type VuNirb( sizeRB,ui );
1124 
1125 
1126  for ( int i =0; i<sizeRB; i++ )
1127  {
1128  for ( int k=0; k<sizeRB; k++ )
1129  {
1130  VuNirb[i].add( real( es.eigenvectors().col( i )[k] ),Vu[k] );
1131  }
1132  }
1133 
1134 
1135  OrthogonalisationRBFunctionL2GrammSchmidt( Xh,VuNirb,MassMatrix );
1136 
1137  auto Mui = M_backend->newVector( Xh );
1138  auto Ui = M_backend->newVector( Xh );
1139 
1140  //save in a file the final NIRB basis functions
1141  for ( int i=0; i<sizeRB; i++ )
1142  {
1143  path = ( boost::format( "./RB%1%NIRB_BasisFile_%2%" ) %sizeRB%i ).str() ;
1144  VuNirb[i].save( _path=path );
1145  Vu[i] = VuNirb[i];
1146  *Ui = VuNirb[i];
1147  Ui->close();
1148  MassMatrix->multVector( Ui,Mui );
1149  M_Vu[i] = *Mui;
1150 
1151  if ( M_Vu[i].l2Norm()>1e-20 )
1152  {
1153  path = ( boost::format( "./RB%1%D_Nirb_Basis%2%File" )%sizeRB%i ).str();
1154  M_Vu[i].save( _path=path );
1155  }
1156 
1157  else
1158  {
1159  std::cerr << "ERROR: in computation of D*NIRB Basis function (" << i << ") : Result equal to zero " <<std::endl;
1160  exit( 0 );
1161  }
1162 
1163  }
1164 }
1165 
1166 //---------------------------------------------------
1167 //---------------------------------------------------
1168 
1169 template< int PolynomialOrder>
1170 typename NIRBTEST<PolynomialOrder>::element_type NIRBTEST<PolynomialOrder> :: BuildNirbSolutionWithoutPostProcess ( space_ptrtype XhFine, element_type uCoarseInterpolate,vector_of_element_type const &M_VNirbBasis, vector_of_element_type const &VNirbBasis, Eigen::VectorXd & BetaiH )
1171 {
1172 
1173 
1174 
1175  auto uNirb = XhFine->element();
1176 
1177  auto UcoarseInterpol = M_backend->newVector( XhFine );
1178  *UcoarseInterpol = uCoarseInterpolate;
1179 
1180  //Computation of the coefficiant \BetaiH = \int uCoarse*\Epsilon_i
1181  //with \Epsilon_i being the final Nirb basis functions
1182 
1183  uNirb.zero();
1184  auto Mui = M_backend->newVector( XhFine );
1185 
1186  for ( int i =0; i<sizeRB; i++ )
1187  {
1188  *Mui=M_VNirbBasis[i];
1189  BetaiH[i] = inner_product( UcoarseInterpol,Mui );
1190  uNirb.add( BetaiH[i], VNirbBasis[i] );
1191  }
1192 
1193 
1194  return uNirb;
1195 
1196 
1197 }
1198 //---------------------------------------------------
1199 //---------------------------------------------------
1200 
1201 template< int PolynomialOrder>
1202 typename NIRBTEST<PolynomialOrder>::element_type NIRBTEST<PolynomialOrder> ::BuildCoarseInterpolation( space_ptrtype XhFine,space_ptrtype XhCoarse, double param )
1203 {
1204  //boost::timer ti;
1205  //ti.restart();
1206  auto uCoarse = XhCoarse->element();
1207  auto uCoarseInterpolate = XhFine->element();
1208  auto ui = XhFine->element();
1209  uCoarse = blackbox( XhCoarse,param );
1210  //std :: cout << "Calculation of uCoarse " << ti.elapsed() << " sec" <<std::endl;
1211  //ti.restart();
1213 
1214  //auto opI = opInterpolation(_domainSpace=XhCoarse,_imageSpace=XhFine);
1215  //std :: cout << "Construction of the operator" << ti.elapsed() << " sec" <<std::endl;
1216  //ti.restart();
1217  //opI->apply(uCoarse,uCoarseInterpolate);
1218  //std :: cout << "Application of the operator" << ti.elapsed() << " sec" <<std::endl;
1219  interpolate ( XhFine,uCoarse,uCoarseInterpolate );
1220  return uCoarseInterpolate;
1221 
1222 }
1223 
1224 
1225 // //---------------------------------------------------
1227 
1228 template< int PolynomialOrder>
1229 typename NIRBTEST<PolynomialOrder>::element_type NIRBTEST<PolynomialOrder> :: BuildNirbSolution( space_ptrtype XhFine,space_ptrtype XhCoarse, vector_of_element_type const &M_VNirbBasis, vector_of_element_type const &VNirbBasis, double param,Eigen::VectorXd BetaiH )
1230 {
1231 
1232  //Eigen::VectorXd BetaiH(sizeRB);
1233  auto uCoarseInterpolate = XhFine->element();
1234  auto uNirb = XhFine->element();
1235  uCoarseInterpolate = BuildCoarseInterpolation( XhFine,XhCoarse,param );
1236  uNirb = BuildNirbSolutionWithoutPostProcess( XhFine,uCoarseInterpolate,M_VNirbBasis,VNirbBasis,BetaiH );
1237  return uNirb;
1238 
1239 }
1240 
1241 //-----------------------------------------
1242 template< int PolynomialOrder>
1243 void NIRBTEST<PolynomialOrder> :: ConstructNIRB( space_ptrtype Xh, vector_of_element_type &M_VNirbBasis, vector_of_element_type &VNirbBasis )
1244 {
1245 
1246  std :: cout << "OFFLINE PROCEDURE : Construction of the 'non intruisive' reduced basis (nirb) " <<std::endl;
1247  boost::timer ti;
1248  double Time_snapshot = 0.;
1249 
1250  if ( Sampling )
1251  {
1252  std::cout << "Sampling Procedure :" <<std::endl;
1253  ComputeSnapshot( Xh,"_" );
1254  Time_snapshot = ti.elapsed();
1255  std::cout << "Computation of the " << NbSnapshot << " snapshots : done -- ";
1256  std::cout << "Time per snapshot: " << Time_snapshot/NbSnapshot << " sec " <<std::endl;
1257  }
1258 
1259  ti.restart();
1260 
1261  //construction of elementary mass and stiffness matrix
1262  auto ui = Xh->element();
1263  auto uj = Xh->element();
1264 
1265  sparse_matrix_ptrtype MassMatrix = M_backend->newMatrix( _test=Xh, _trial=Xh ); //Sparse FE mass matrix
1266  form2( _test=Xh , _trial=Xh, _matrix= MassMatrix ) =
1267  integrate( _range=elements( Xh->mesh() ), _expr=id( uj )*idt( ui ) );
1268 
1269  sparse_matrix_ptrtype StiffMatrix = M_backend->newMatrix( _test=Xh, _trial=Xh ); //Sparse FE mass matrix
1270  form2( _test=Xh , _trial=Xh, _matrix= StiffMatrix ) =
1271  integrate( _range=elements( Xh->mesh() ), _expr=gradt( ui )*trans( grad( uj ) ) );
1272  double Time_Construct_Elementary_Matrix = ti.elapsed();
1273  std::cout << "Time to build elementary F.E matrix = " << Time_Construct_Elementary_Matrix << " sec " <<std::endl;
1274 
1275  ti.restart();
1276  ChooseRBFunction( Xh,VNirbBasis,MassMatrix,StiffMatrix );
1277  double TimeChooseRB = ti.elapsed();
1278  std::cout << "Choice of " << sizeRB << " reduced basis functions : done -- ";
1279  std::cout << "Time: " << TimeChooseRB << " sec "<<std::endl;
1280  ti.restart();
1281 
1282  //Orthogonalisation de Gram-schmidt
1283 
1284  OrthogonalisationRBFunction( Xh,VNirbBasis,M_VNirbBasis,StiffMatrix,MassMatrix );
1285  //OrthogonalisationRBFunction(Xh);
1286  double Time_buildNIRBbasis = ti.elapsed();
1287  std::cout << "H1 and L2 orthogonalisation of the reduced basis functions : done -- ";
1288  std::cout << "Time: " << Time_buildNIRBbasis << " sec " <<std::endl;
1289  double TotalTime = Time_buildNIRBbasis +TimeChooseRB + Time_Construct_Elementary_Matrix + Time_snapshot;
1290  std::cout <<"Total time for the OFFLINE procedure: " << TotalTime << " sec "<<std::endl;
1291 
1292 
1293 
1294 }
1295 
1296 
1297 
1298 //---------------------------------------------------
1299 //---------------------------------------------------
1300 template< int PolynomialOrder>
1301 typename NIRBTEST<PolynomialOrder>::element_type NIRBTEST<PolynomialOrder>::blackbox( space_ptrtype Xh, double param )
1302 {
1303 
1304  auto u = Xh->element();
1305  auto v = Xh->element();
1306  // std::cout << "Xh ndof=" << Xh->nLocalDof() << "\n";
1307  value_type pi = M_PI;
1308  double penaltyTerm = 30.;
1309 
1311  auto velocity = vec( cst( cos( param ) ),cst( sin( param ) ) );
1312  /*
1313  auto g = ( chi(abs(Px()-1) < 1e-10 )*Px()*Px()+
1314  chi(abs(Py()-1) < 1e-10 )*Py()*Py() );
1315  */
1316 
1317  auto g = ( Px()*Px()*Py()*Py() );
1318  auto F = M_backend->newVector( Xh );
1319  /*
1320  form1( _test=Xh, _vector=F ) =
1321  integrate( _range=boundaryfaces(Xh->mesh()),
1322  _expr=g*(-0.01*dn(v)+penaltyTerm*id(v)/hFace()) );
1323  */
1324  auto D = M_backend->newMatrix( _test=Xh, _trial=Xh );
1325  form2( _test=Xh, _trial=Xh, _matrix=D ) =
1326  integrate( _range=elements( Xh->mesh() ), _expr=0.01*gradt( u )*trans( grad( v ) )+ ( gradt( u )*velocity )*id( v ) );
1327 
1328  /*
1329  form2( _test=Xh, _trial=Xh, _matrix=D ) +=
1330  integrate( boundaryfaces(Xh->mesh()),
1331  -0.01*dnt(u)*id(v)-0.01*dn(v)*idt(u)+penaltyTerm*id(v)*idt(u)/hFace());
1332  */
1333 
1334  form2( _test=Xh, _trial=Xh, _matrix=D ) += on( boundaryfaces( Xh->mesh() ), u, F,g );
1335  backend_type::build()->solve( _matrix=D, _solution=u, _rhs=F );
1336 
1337  return u;
1338 }
1339 
1340 
1341 //---------------------------------------------------
1342 //---------------------------------------------------
1343 //Post-processing of NIRB solution
1344 template< int PolynomialOrder>
1345 Eigen::MatrixXd NIRBTEST<PolynomialOrder> :: BuildBetaH( space_ptrtype XhFine,
1346  space_ptrtype XhCoarse,vector_of_element_type &M_VNirbBasis )
1347 {
1348 
1349  //boost::timer ti;
1350  //ti.restart();
1351  //auto opI = opInterpolation(_domainSpace=XhCoarse,_imageSpace=XhFine);
1352  //std :: cout << "Construction of the operator" << ti.elapsed() << " sec" <<std::endl;
1353  //ti.restart();
1354  auto uCoarse = XhCoarse->element();
1355  auto uCoarseInterpolation = XhFine->element();
1356 
1357  std::vector<int> IndTab( sizeRB );
1358  Eigen::MatrixXd BetaH( sizeRB,sizeRB );
1359  std::string path = ( boost::format( "./IndBR%1%" ) %sizeRB ).str() ;
1360  std::ifstream find( path.c_str() );
1361 
1362  if ( !find )
1363  {
1364  std::cerr << "ERROR : Problem in opening 'IndBR'" << sizeRB << " file " <<std::endl;
1365  exit( 0 );
1366 
1367 
1368  }
1369 
1370  for ( int i =0 ; i < sizeRB; i++ )
1371  {
1372  int Itemp;
1373  find >> Itemp;
1374  IndTab[i] = Itemp;
1375  }
1376 
1377  path = ( boost::format( "./AlphaH_%1%" ) %sizeRB ).str();
1378  std::ofstream fB( path.c_str() );
1379  fB << sizeRB << " " << sizeRB <<std::endl;
1380 
1381  for ( int i =0; i<sizeRB; i++ )
1382  {
1383  uCoarse.zero();
1384  uCoarseInterpolation.zero();
1385  path = ( boost::format( "./Sol_Coarse_%1%" ) %IndTab[i] ).str();
1386  uCoarse.load( _path=path );
1387 
1388  if ( uCoarse.l2Norm()==0. )
1389  {
1390  std::cerr <<"'BuildBetaH routine': ERROR IN LOADING FILE " << path <<std::endl;
1391  std::cerr <<"Coarse sampling has to be done " << std::endl;
1392  ComputeSnapshot( XhCoarse,"_Coarse_" );
1393  i--;
1394 
1395  }
1396 
1397  //opI->apply(uCoarse,uCoarseInterpolation);
1398  interpolate( XhFine,uCoarse,uCoarseInterpolation );
1399  auto UcoarseInterpol_i = M_backend->newVector( XhFine );
1400  *UcoarseInterpol_i = uCoarseInterpolation;
1401  auto Muj = M_backend->newVector( XhFine );
1402 
1403  for ( int j =0; j<sizeRB; j++ )
1404  {
1405  *Muj=M_VNirbBasis[j];
1406  BetaH( j,i ) = inner_product( UcoarseInterpol_i,Muj ); //beta_j(\mu_i)
1407  fB << BetaH( j,i ) <<std::endl;
1408  }
1409  }
1410 
1411  fB.close();
1412  return BetaH;
1413 
1414 }
1415 //---------------------------------------------------
1416 //---------------------------------------------------
1417 
1418 template< int PolynomialOrder>
1419 Eigen::MatrixXd NIRBTEST<PolynomialOrder> :: BuildBetah( space_ptrtype XhFine,
1420  vector_of_element_type &M_VNirbBasis )
1421 {
1422 
1423  auto uFine = XhFine->element();
1424 
1425  std::vector<int> IndTab( sizeRB );
1426  Eigen::MatrixXd BetaH( sizeRB,sizeRB );
1427  std::string path = ( boost::format( "./IndBR%1%" ) %sizeRB ).str() ;
1428  std::ifstream find( path.c_str() );
1429 
1430  if ( !find )
1431  {
1432  std::cerr << "ERROR : Problem in opening 'IndBR'" << sizeRB << " file " <<std::endl;
1433  exit( 0 );
1434  }
1435 
1436  for ( int i =0 ; i < sizeRB; i++ )
1437  {
1438  int Itemp;
1439  find >> Itemp;
1440  IndTab[i] = Itemp;
1441  }
1442 
1443  path = ( boost::format( "./Betah_%1%" ) %sizeRB ).str();
1444  std::ofstream fB( path.c_str() );
1445  fB << sizeRB << " " << sizeRB <<std::endl;
1446 
1447  for ( int i =0; i<sizeRB; i++ )
1448  {
1449  uFine.zero();
1450  path = ( boost::format( "./Sol_%1%" ) %IndTab[i] ).str();
1451  uFine.load( _path=path );
1452 
1453  if ( uFine.l2Norm()==0. )
1454  {
1455  std::cerr <<"'ChooseRBFunction routine': ERROR IN LOADING FILE " << path <<std::endl;
1456  exit( 0 );
1457  }
1458 
1459  auto Ufine_i = M_backend->newVector( XhFine );
1460  *Ufine_i = uFine;
1461  auto Muj = M_backend->newVector( XhFine );
1462 
1463  for ( int j =0; j<sizeRB; j++ )
1464  {
1465  *Muj=M_VNirbBasis[j];
1466  BetaH( j,i ) = inner_product( Ufine_i,Muj ); // beta_j(\mu_i)
1467  fB << BetaH( j,i ) <<std::endl;
1468  }
1469  }
1470 
1471  fB.close();
1472 
1473  return BetaH;
1474 
1475 }
1476 //---------------------------------------------------
1477 //---------------------------------------------------
1478 
1479 template< int PolynomialOrder>
1480 Eigen::MatrixXd NIRBTEST<PolynomialOrder> ::BuildPostProcessMatrix( space_ptrtype XhFine, space_ptrtype XhCoarse,vector_of_element_type &M_VNirbBasis )
1481 {
1482 
1483  Eigen::MatrixXd BetaH( sizeRB,sizeRB );
1484  Eigen::MatrixXd Betah( sizeRB,sizeRB );
1485  Eigen::MatrixXd T( sizeRB,sizeRB );
1486 
1487  int FileOkA = 1,FileOkB = 1 ;
1488 
1489  if ( !Offline )
1490  {
1491  std::string path = ( boost::format( "./AlphaH_%1%" ) %sizeRB ).str();
1492  std::ifstream fA ( path.c_str() );
1493 
1494  if ( !fA )
1495  {
1496  FileOkA = 0;
1497  BetaH = BuildBetaH( XhFine,XhCoarse,M_VNirbBasis );
1498  }
1499 
1500  path = ( boost::format( "./betah_%1%" ) %sizeRB ).str();
1501  std::ifstream fB ( path.c_str() );
1502 
1503  if ( !fB )
1504  {
1505  FileOkB = 0;
1506  Betah = BuildBetah( XhFine,M_VNirbBasis );
1507  }
1508 
1509  int Itemp1,Itemp2;
1510  double dTempA,dTempB;
1511 
1512  if ( FileOkA && FileOkB )
1513  {
1514  fA >> Itemp1 >> Itemp2;
1515  fB >> Itemp1 >> Itemp2;
1516 
1517  for ( int i = 0; i< sizeRB; i++ )
1518  {
1519  for ( int j=0; j<sizeRB; j++ )
1520  {
1521  fA >> dTempA;
1522  fB >> dTempB;
1523  BetaH( j,i ) = dTempA;
1524  Betah( j,i ) = dTempB;
1525  }
1526  }
1527  }
1528 
1529  else
1530  {
1531  if ( FileOkA )
1532  {
1533  fA >> Itemp1 >> Itemp2;
1534 
1535  for ( int i = 0; i< sizeRB; i++ )
1536  {
1537  for ( int j=0; j<sizeRB; j++ )
1538  {
1539  fA >> dTempA;
1540  BetaH( j,i ) = dTempA;
1541  }
1542  }
1543  }
1544 
1545  if ( FileOkB )
1546  {
1547  fB >> Itemp1 >> Itemp2;
1548 
1549  for ( int i = 0; i< sizeRB; i++ )
1550  {
1551  for ( int j=0; j<sizeRB; j++ )
1552  {
1553  fB >> dTempB;
1554  Betah( j,i ) = dTempB;
1555  }
1556  }
1557  }
1558  }
1559 
1560  fA.close();
1561  fB.close();
1562 
1563  }
1564 
1565  if ( Offline )
1566  {
1567  BetaH = BuildBetaH( XhFine,XhCoarse,M_VNirbBasis );
1568  Betah = BuildBetah( XhFine,M_VNirbBasis );
1569  }
1570 
1571  T = Betah* BetaH.inverse() ;
1572 
1573  return T;
1574 
1575 
1576 }
1577 
1578 
1579 // //---------------------------------------------------
1581 
1582 template< int PolynomialOrder>
1583 typename NIRBTEST<PolynomialOrder>::element_type NIRBTEST<PolynomialOrder> ::
1584 BuildNirbSolutionWithPostProcess( space_ptrtype XhFine,
1585  space_ptrtype XhCoarse,
1586  vector_of_element_type &M_VNirbBasis,
1587  vector_of_element_type &VNirbBasis,
1588  Eigen::VectorXd &BetaiuH,
1589  element_type &uCoarseInterpolate )
1590 {
1591 
1592  Eigen::MatrixXd T ( sizeRB,sizeRB );
1593  T = BuildPostProcessMatrix( XhFine,XhCoarse,M_VNirbBasis );
1594 
1595  //Eigen::VectorXd AlphaiH (sizeRB);
1596  //AlphaiH = T*BetaiuH;
1597  double AlphaiH;
1598  auto uNirb = XhFine->element();
1599  uNirb.zero();
1600 
1601  for ( int i =0; i<sizeRB; i++ )
1602  {
1603  AlphaiH =0.;
1604 
1605  for ( int k=0; k<sizeRB; k++ )
1606  {
1607  AlphaiH += T( i,k )*BetaiuH[k];
1608  }
1609 
1610  uNirb.add( AlphaiH, VNirbBasis[i] );
1611  }
1612 
1613 
1614  return uNirb;
1615 }
1616 
1617 
1618 
1619 #endif

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