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
modelcrbbase.hpp
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): Stephane Veys <stephane.veys@imag.fr>
6  Date: 2013-02-22
7 
8  Copyright (C) 2008-2012 Universite Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef ModelCrbBase_H
30 #define ModelCrbBase_H 1
31 
32 #include <feel/feel.hpp>
33 #include <feel/feelcrb/eim.hpp>
34 
35 namespace Feel
36 {
37 
38 
39 
40 class ParameterDefinitionBase
41 {
42 public :
43  typedef ParameterSpace<1> parameterspace_type ;
44 };
45 
46 class FunctionSpaceDefinitionBase
47 {
48 public :
49  /*mesh*/
50  typedef Simplex<1> entity_type ;
51  typedef Mesh<entity_type > mesh_type ;
52 
53  /*basis*/
54  typedef Lagrange<1,Scalar> basis_type ;
55 
56  /*space*/
57  typedef FunctionSpace<mesh_type , basis_type > space_type ;
58 };
59 
60 template <typename ParameterDefinition, typename FunctionSpaceDefinition>
61 class EimDefinitionBase
62 {
63 
64 public :
65  typedef typename ParameterDefinition::parameterspace_type parameterspace_type;
66  typedef typename FunctionSpaceDefinition::space_type space_type;
67 
68  /* EIM */
69  typedef EIMFunctionBase<space_type , space_type , parameterspace_type > fun_type ;
70  typedef EIMFunctionBase<space_type , space_type , parameterspace_type > fund_type ;
71 
72 };
73 
74 
75 
76 template <typename ParameterDefinition, typename FunctionSpaceDefinition, typename EimDefinition = EimDefinitionBase<ParameterDefinition,FunctionSpaceDefinition> >
77 class ModelCrbBase : public ModelCrbBaseBase
78 {
79 
80 public :
81  typedef typename EimDefinition::fun_type fun_type;
82  typedef typename EimDefinition::fund_type fund_type;
83 
84  typedef typename ParameterDefinition::parameterspace_type parameterspace_type;
85  typedef typename parameterspace_type::element_type parameter_type;
86 
87  typedef typename FunctionSpaceDefinition::space_type space_type;
88 
89  typedef boost::shared_ptr<fun_type> fun_ptrtype;
90  typedef std::vector<fun_ptrtype> funs_type;
91 
92  typedef boost::shared_ptr<fund_type> fund_ptrtype;
93  typedef std::vector<fund_ptrtype> funsd_type;
94 
95  typedef Eigen::VectorXd vectorN_type;
96 
97  typedef std::vector< std::vector< double > > beta_vector_type;
98 
99  typedef OperatorLinear< space_type , space_type > operator_type;
100  typedef boost::shared_ptr<operator_type> operator_ptrtype;
101 
102  typedef OperatorLinearComposite< space_type , space_type > operatorcomposite_type;
103  typedef boost::shared_ptr<operatorcomposite_type> operatorcomposite_ptrtype;
104 
105  typedef FsFunctionalLinearComposite< space_type > functionalcomposite_type;
106  typedef boost::shared_ptr<functionalcomposite_type> functionalcomposite_ptrtype;
107 
108  typedef FsFunctionalLinear< space_type > functional_type;
109  typedef boost::shared_ptr<functional_type> functional_ptrtype;
110 
111  typedef backend_type::vector_type vector_type;
112  typedef backend_type::vector_ptrtype vector_ptrtype;
113 
114  typedef backend_type::sparse_matrix_type sparse_matrix_type;
115  typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
116 
117  ModelCrbBase()
118  :
119  M_is_initialized( false )
120  {
121  }
122 
123  void setInitialized( const bool & b )
124  {
125  M_is_initialized = b ;
126  }
127 
128  bool isInitialized()
129  {
130  return M_is_initialized;
131  }
132 
133  virtual funs_type scalarContinuousEim () const
134  {
135  return M_funs;
136  }
137 
138  virtual funsd_type scalarDiscontinuousEim () const
139  {
140  return M_funs_d;
141  }
142 
143  virtual void initModel() = 0;
144 
145  virtual operatorcomposite_ptrtype operatorCompositeA()
146  {
147  return M_compositeA;
148  }
149  virtual operatorcomposite_ptrtype operatorCompositeM()
150  {
151  return M_compositeM;
152  }
153  virtual std::vector< functionalcomposite_ptrtype > functionalCompositeF()
154  {
155  return M_compositeF;
156  }
157 
158 
159  //for linear steady models, mass matrix does not exist
160  //non-linear steady models need mass matrix for the initial guess
161  //this class can provide the function operatorCompositeM() to ensure compilation
162  //but we need to know if the model can provide an operator composite for mass matrix
163  //because non-linear problems have to do these operators
164  //because CRBModel is not able to construct such operators in the general case.
165  //Elements of function space are members of CRBModel
166  //then for composite spaces, we need a view of these elements
167  //BUT we can't have a view of an element of a non-composite space
168  //this function returns true if the model provides an operator composite for mass matrix
169  //false by default
170  virtual bool constructOperatorCompositeM()
171  {
172  return false;
173  }
174 
175 
176 
182  virtual int QInitialGuess() const
183  {
184  return 1;
185  }
186 
187  virtual int mMaxInitialGuess( int q ) const
188  {
189  if( q == 0 )
190  return 1;
191  else
192  throw std::logic_error( "[ModelCrbBase::mMaxInitialGuess(q)] ERROR wrong index q, should be 0" );
193  }
194 
195  virtual beta_vector_type computeBetaInitialGuess( parameter_type const& mu )
196  {
197  beta_vector_type beta;
198  beta.resize(1); //q=1
199  beta[0].resize(1); //m=1
200  beta[0][0]=0;
201  return beta;
202  }
203 
204 
205 
211  virtual double scalarProductForMassMatrix( vector_ptrtype const& X, vector_ptrtype const& Y )
212  {
213  throw std::logic_error("Your model is time-dependant so you MUST implement scalarProductForMassMatrix function");
214  return 0;
215  }
216  virtual double scalarProductForMassMatrix( vector_type const& x, vector_type const& y )
217  {
218  throw std::logic_error("Your model is time-dependant so you MUST implement scalarProductForMassMatrix function");
219  return 0;
220  }
221 
226  virtual sparse_matrix_ptrtype const& innerProductForMassMatrix () const
227  {
228  throw std::logic_error("Your model is time-dependant so you MUST implement innerProductForMassMatrix function");
229  return M;
230  }
231  virtual sparse_matrix_ptrtype innerProductForMassMatrix ()
232  {
233  throw std::logic_error("Your model is time-dependant so you MUST implement innerProductForMassMatrix function");
234  return M;
235  }
236 
242  vectorN_type computeStatistics( Eigen::VectorXd vector , std::string name )
243  {
244  double min=0,max=0,mean=0,mean1=0,mean2=0,standard_deviation=0,variance=0;
245  Eigen::MatrixXf::Index index;
246  Eigen::VectorXd square;
247 
248  if( vector.size() > 0 )
249  {
250  bool force = option("eim.use-dimension-max-functions").template as<bool>();
251  int Neim=0;
252  if( force )
253  Neim = option("eim.dimension-max").template as<int>();
254 
255  int N = vector.size();
256 
257  if( force )
258  LOG( INFO ) <<" statistics for "<<name<<" ( was called "<< N << " times with "<<Neim<<" basis functions )";
259  else
260  LOG( INFO ) <<" statistics for "<<name<<" ( was called "<< N << " times )";
261 
262  min = vector.minCoeff(&index);
263  max = vector.maxCoeff(&index);
264  mean = vector.mean();
265  mean1 = mean * mean;
266  square = vector.array().pow(2);
267  mean2 = square.mean();
268  standard_deviation = math::sqrt( mean2 - mean1 );
269  LOG(INFO)<<"min : "<<min<<" - max : "<<max<<" mean : "<<mean<<" standard deviation : "<<standard_deviation;
270  }
271  vectorN_type result(4);
272  result(0)=min;
273  result(1)=max;
274  result(2)=mean;
275  result(3)=standard_deviation;
276  return result;
277  }
278 
279 
280 protected :
281 
282  funs_type M_funs;
283  funsd_type M_funs_d;
284  bool M_is_initialized;
285 
286  sparse_matrix_ptrtype M;
287 
288  operatorcomposite_ptrtype M_compositeA;
289  operatorcomposite_ptrtype M_compositeM;
290  std::vector< functionalcomposite_ptrtype > M_compositeF;
291 
292 };
293 
294 }//Feel
295 #endif /* __Model_H */

Generated on Sun Oct 20 2013 08:25:02 for Feel++ by doxygen 1.8.4