dune-fem  2.6-git
codegen.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_BASEFUNCTIONSETS_CODEGEN_HH
2 #define DUNE_FEM_BASEFUNCTIONSETS_CODEGEN_HH
3 
4 #include <cassert>
5 #include <cstdlib>
6 #include <iostream>
7 #include <fstream>
8 #include <sstream>
9 #include <vector>
10 #include <set>
11 #include <map>
12 
13 #include <dune/common/exceptions.hh>
14 #include <dune/fem/io/io.hh>
15 #include <dune/fem/io/parameter.hh>
16 
17 namespace Dune
18 {
19 
20  namespace Fem
21  {
22 
23  class CodegenInfoFinished : public Dune :: Exception {};
24 
26  {
27  std::string path_;
28  std::string outFileName_;
29 
30  int nopMax_;
31  int nopMin_;
32 
33  int baseMax_;
34  int baseMin_;
35 
36  typedef std::set< int > DimRangeSetType;
37  typedef std::set< void* > BaseSetPointerType;
38  DimRangeSetType savedimRanges_;
39  mutable DimRangeSetType dimRanges_;
40  BaseSetPointerType savedBaseSets_;
41 
42  typedef void codegenfunc_t (std::ostream& out,
43  const int dim,
44  const int dimRange,
45  const size_t numRows,
46  const size_t numCols );
47 
48  typedef std::pair< std::string, int > EntryType;
49  std::vector< EntryType > filenames_;
50 
51  typedef std::vector< int > KeyType;
52  typedef std::set< KeyType > KeyStorageType;
53 
54  typedef std::pair< int , int > EvalPairType ;
55  typedef std::set< EvalPairType > EvalSetType ;
56 
57  std::map< codegenfunc_t* , KeyStorageType > entryMap_;
58  EvalSetType evalSet_;
59 
60  CodegenInfo()
61  : path_( Dune::Fem::Parameter::commonInputPath() + "/autogeneratedcode"),
62  outFileName_( "/autogeneratedcode.hh" ),
63  nopMax_(0), nopMin_(0), baseMax_(0), baseMin_(0),
64  dimRanges_(), savedBaseSets_(), filenames_(),
65  entryMap_(),
66  evalSet_()
67  {
68  }
69 
70  public:
72  {
73  static CodegenInfo info;
74  return info;
75  }
76 
78  void clear() {
79  savedBaseSets_.clear();
80  filenames_.clear();
81  dimRanges_.clear();
82  entryMap_.clear();
83  evalSet_.clear();
84  nopMax_ = 0;
85  nopMin_ = 0;
86  baseMax_ = 0;
87  baseMin_ = 0;
88  }
89 
91  void setPath(const std::string& path ) { path_ = path ; }
92 
94  void setFileName(const std::string& filename ) { outFileName_ = filename ; }
95 
96  template <class BaseSet>
97  void addDimRange(const BaseSet* baseSet,
98  const int dimRange)
99  {
100  void* ptr = (void *) baseSet;
101  if( savedBaseSets_.find( ptr ) == savedBaseSets_.end() )
102  {
103  savedBaseSets_.insert( ptr );
104  std::cout << "Add dimRange " << dimRange << std::endl;
105  dimRanges_.insert( dimRange ) ;
106  }
107  }
108 
109  void addEntry(const std::string& fileprefix,
110  codegenfunc_t* codegenfunc,
111  const int dim, const int dimRange, const int quadNop, const int numBase )
112  {
113  KeyStorageType& keyMap = entryMap_[ codegenfunc ];
114 
115  typedef KeyStorageType :: iterator iterator;
116  KeyType key( 4 );
117  key[ 0 ] = dim;
118  key[ 1 ] = dimRange;
119  key[ 2 ] = quadNop;
120  key[ 3 ] = numBase;
121 
122  // search for key, if already exists, then do nothing
123  iterator it = keyMap.find( key );
124  if( it != keyMap.end() ) return ;
125 
126  // make sure dimRange was registered
127  assert( dimRanges_.find( dimRange ) != dimRanges_.end() );
128 
129  // create directory to store files
130  createDirectory( path_ );
131 
132  std::stringstream filename;
133  filename << fileprefix << dimRange << "_" << quadNop << "_" << numBase << ".hh";
134 
135  // second check, if file exists, do nothing
136  int pos = exists( filename.str() );
137  if( pos >= 0 ) return ;
138 
139  std::string filenameToWrite( path_ + "/" + filename.str() );
140  std::ofstream file( filenameToWrite.c_str(), std::ios::out );
141 
142  // call code generation function
143  codegenfunc( file, dim, dimRange, quadNop, numBase );
144  std::cout << "Generate code " << fileprefix << " for (" << dimRange << ","
145  << quadNop << "," << numBase << ")" << std::endl;
146 
147  // insert evaluate pair
148  EvalPairType evalPair ( quadNop, numBase );
149  evalSet_.insert( evalPair );
150 
151  if( baseMin_ == 0 ) baseMin_ = numBase;
152  if( nopMin_ == 0 ) nopMin_ = quadNop;
153 
154  EntryType entry ( filename.str() , dimRange );
155  filenames_.push_back( entry );
156  nopMax_ = std::max( quadNop, nopMax_ );
157  nopMin_ = std::min( quadNop, nopMin_ );
158  baseMax_ = std::max( numBase, baseMax_ );
159  baseMin_ = std::min( numBase, baseMin_ );
160  }
161 
162  void finalize () const
163  {
164  if( checkAbort() )
165  {
166  dumpInfo();
167  std::cerr << "All automated code generated, bye, bye !! " << std::endl;
168  DUNE_THROW(CodegenInfoFinished,"All automated code generated, bye, bye !!");
169  }
170  }
171 
172 
173  void dumpInfo() const
174  {
175  {
176  std::ofstream file( outFileName_.c_str() );
177  file << "#include \"" <<path_<< "/" << outFileName_ << "\"";
178  }
179 
180  {
181  std::string filename( path_ + "/" + outFileName_ );
182  std::ofstream file( filename.c_str() );
183  file << "#ifdef CODEGEN_INCLUDEMAXNUMS" << std::endl;
184  file << "#ifndef CODEGEN_INCLUDEMAXNUMS_INCLUDED" << std::endl;
185  file << "#define CODEGEN_INCLUDEMAXNUMS_INCLUDED" << std::endl << std::endl;
186  file << "#define MAX_NUMBER_OF_QUAD_POINTS " << nopMax_ << std::endl;
187  file << "#define MIN_NUMBER_OF_QUAD_POINTS " << nopMin_ << std::endl;
188  file << "#define MAX_NUMBER_OF_BASE_FCT " << baseMax_ << std::endl;
189  file << "#define MIN_NUMBER_OF_BASE_FCT " << baseMin_ << std::endl << std::endl;
190  file << "/* include all headers with inner loop extern declarations */" << std::endl;
191  file << "#define CODEGEN_COMPILE_INNERLOOPS 1" << std::endl;
192  for( size_t i = 0; i < filenames_.size(); ++i )
193  {
194  file << "#include \""<< filenames_[ i ].first << "\"" << std::endl;
195  }
196  file << "#undef CODEGEN_COMPILE_INNERLOOPS" << std::endl << std::endl;
197  file << "#include \"" << filename << ".c\"" <<std::endl;
198 
199  file << "#endif // CODEGEN_INCLUDEMAXNUMS_INCLUDED" << std::endl << std::endl;
200  file << "#elif defined CODEGEN_INCLUDEEVALCALLERS" << std::endl;
201  file << "#ifndef CODEGEN_EVALCALLERS_INCLUDED" << std::endl;
202  file << "#define CODEGEN_EVALCALLERS_INCLUDED" << std::endl << std::endl;
203  typedef EvalSetType :: iterator iterator ;
204  const iterator endit = evalSet_.end();
205  for( iterator it = evalSet_.begin(); it != endit; ++it )
206  {
207  file << " template <class Traits>" << std::endl;
208  file << " struct EvaluateImplementation< Traits, " << it->first << " , " << it->second << " >" << std::endl;
209  file << " : public EvaluateRealImplementation< Traits, " << it->first << " , " << it->second << " >" << std::endl;
210  file << " {" << std::endl;
211  file << " typedef EvaluateRealImplementation< Traits, " << it->first << " , " << it->second << " > BaseType;" << std::endl;
212  file << " typedef typename BaseType :: RangeVectorType RangeVectorType;" << std::endl;
213  file << " EvaluateImplementation( const RangeVectorType& rv ) : BaseType ( rv ) {}" << std::endl;
214  file << " };"<< std::endl;
215  file << std::endl;
216  }
217  file << "#endif // CODEGEN_EVALCALLERS_INCLUDED" << std::endl << std::endl;
218  file << "#else" << std::endl << std::endl ;
219  file << "#ifndef CODEGEN_INCLUDE_IMPLEMENTATION" << std::endl;
220  file << "#define CODEGEN_INCLUDE_IMPLEMENTATION" << std::endl;
221  file << "#undef CODEGEN_COMPILE_INNERLOOPS" << std::endl;
222  for( size_t i = 0; i < filenames_.size(); ++i )
223  {
224  file << "#include \""<< filenames_[ i ].first << "\"" << std::endl;
225  }
226  file << "#endif // CODEGEN_INCLUDE_IMPLEMENTATION" << std::endl << std::endl;
227  file << "#endif // CODEGEN_INCLUDEMAXNUMS" << std::endl;
228 
229  // write C file with implementation of inner loop functions
230  filename += ".c";
231  std::ofstream Cfile( filename.c_str() );
232 
233  Cfile << "#include <stdlib.h>" << std::endl;
234  Cfile << "/* include all headers with inner loop implementation */" << std::endl;
235  Cfile << "#define CODEGEN_COMPILE_INNERLOOPS 2" << std::endl;
236  for( size_t i = 0; i < filenames_.size(); ++i )
237  {
238  Cfile << "#include \""<< filenames_[ i ].first << "\"" << std::endl;
239  }
240  }
241  }
242  protected:
243  int exists( const std::string& filename ) const
244  {
245  for( size_t i = 0; i < filenames_.size(); ++i )
246  {
247  if( filename == filenames_[ i ].first )
248  {
249  return i;
250  }
251  }
252  return -1;
253  }
254 
255  bool checkAbort() const
256  {
257  DimRangeSetType found ;
258  bool canAbort = true ;
259  for( size_t i = 0; i < filenames_.size(); ++i )
260  {
261  found.insert( filenames_[ i ].second );
262  if ( filenames_[ i ].second > 0 )
263  {
264  canAbort = false ;
265  }
266  }
267  typedef DimRangeSetType :: iterator iterator ;
268  for( iterator it = found.begin(); it != found.end(); ++it )
269  {
270  dimRanges_.erase( *it );
271  }
272 
273  if( canAbort && dimRanges_.size() == 0 )
274  return true ;
275  else
276  return false ;
277  }
278  };
279 
282  {
283  static void evaluateCodegen(std::ostream& out,
284  const int dim,
285  const int dimRange,
286  const size_t numRows,
287  const size_t numCols )
288  {
289  out << "#if ! CODEGEN_COMPILE_INNERLOOPS" << std::endl;
290  out << "template <class BaseFunctionSet>" << std::endl;
291  out << "struct EvaluateRanges<BaseFunctionSet, Fem :: EmptyGeometry, " << dimRange << ", " << numRows << ", " << numCols << ">" << std::endl;
292  out << "{" << std::endl;
293  out << " template< class QuadratureType,"<< std::endl;
294  out << " class RangeVectorType," << std::endl;
295  out << " class RangeFactorType," << std::endl;
296  out << " class LocalDofVectorType>" << std::endl;
297  out << " static void eval( const QuadratureType& quad," << std::endl;
298  out << " const RangeVectorType& rangeStorage," << std::endl;
299  out << " const LocalDofVectorType& dofStorage," << std::endl;
300  out << " RangeFactorType &rangeVector)" << std::endl;
301  out << " {" << std::endl;
302  out << " typedef typename ScalarRangeType :: field_type field_type;" << std::endl;
303  out << " typedef typename RangeVectorType :: value_type value_type;" << std::endl;
304  // make length sse conform
305  //const int sseDimRange = dimRange + (dimRange % 2);
306  const int sseDimRange = dimRange; // + (dimRange % 2);
307  out << " const field_type dofs[ " << numCols * dimRange << " ] = { " << std::endl;
308  for( size_t col = 0, colR = 0; col < numCols; ++ col )
309  {
310  out << " ";
311  for( int r = 0; r < dimRange; ++ r , ++colR )
312  {
313  out << "dofStorage[ " << colR << " ]";
314 
315  if( colR < (numCols * dimRange) - 1 )
316  {
317  out << ",";
318  if( r == dimRange - 1) out << std::endl;
319  }
320  else out << " };" << std::endl;
321  }
322  }
323  out << " for( size_t row=0; row<"<<numRows<<" ; ++row )"<<std::endl;
324  out << " {" << std::endl;
325  out << " const value_type& rangeStorageRow = rangeStorage[ quad.cachingPoint( row ) ];" << std::endl;
326  out << " field_type result [ " << sseDimRange << " ] = { ";
327  for( int r = 0; r < sseDimRange-1; ++r ) out << " 0 ,";
328  out << " 0 };" << std::endl;
329  for( size_t col = 0, colR = 0; col < numCols; ++col )
330  {
331  out << " const field_type phi"<< col << " = rangeStorageRow[ " << col << " ][ 0 ];" << std::endl;
332  for( int r = 0; r < dimRange; ++r , ++colR )
333  {
334  out << " result[ " << r << " ] += dofs[ " << colR << " ] * phi" << col << ";" << std::endl;
335  if( sseDimRange != dimRange && r == dimRange - 1 )
336  {
337  if( (colR + 1) < numCols * dimRange )
338  out << " result[ " << r+1 << " ] += dofs[ " << colR + 1 << " ] * phi" << col << ";" << std::endl;
339  else
340  out << " result[ " << r+1 << " ] += 0.5 * phi" << col << ";" << std::endl;
341  }
342  }
343  }
344  out << " // store result in vector" << std::endl;
345  out << " RangeType& realResult = rangeVector[ row ];" << std::endl;
346  for( int r = 0; r < dimRange; ++r)
347  {
348  out << " realResult[ " << r << " ] = result[ " << r << " ];" << std::endl;
349  }
350  out << " }" << std::endl;
351  out << " }" << std::endl << std::endl;
352  out << "};" << std::endl;
353  out << "#endif" << std::endl;
354  }
355 
356  static void axpyCodegen(std::ostream& out,
357  const int dim, const int dimRange, const size_t numRows, const size_t numCols )
358  {
359  out << "#if ! CODEGEN_COMPILE_INNERLOOPS" << std::endl;
360  out << "template <class BaseFunctionSet>" << std::endl;
361  out << "struct AxpyRanges<BaseFunctionSet, Fem :: EmptyGeometry, " << dimRange << ", " << numRows << ", " << numCols << ">" << std::endl;
362  out << "{" << std::endl;
363  out << " template< class QuadratureType,"<< std::endl;
364  out << " class RangeVectorType," << std::endl;
365  out << " class RangeFactorType," << std::endl;
366  out << " class LocalDofVectorType>" << std::endl;
367  out << " static void axpy( const QuadratureType& quad," << std::endl;
368  out << " const RangeVectorType& rangeStorage," << std::endl;
369  out << " const RangeFactorType &rangeFactors," << std::endl;
370  out << " LocalDofVectorType& dofs)" << std::endl;
371  out << " {" << std::endl;
372  for( size_t row=0; row< numRows; ++row )
373  {
374  out << " const RangeType& factor" << row << " = rangeFactors[ " << row << " ];" << std::endl;
375  }
376  out << std::endl ;
377 
378  out << " typedef typename RangeVectorType :: value_type value_type;" << std::endl;
379  for( size_t row = 0; row<numRows ; ++ row )
380  {
381  out << " const value_type& rangeStorage"<<row<<" = rangeStorage[ quad.cachingPoint( " << row << " ) ];" << std::endl;
382  }
383  out << std::endl;
384  out << " typedef typename ScalarRangeType :: field_type field_type;" << std::endl;
385  const int sseDimRange = dimRange + (dimRange % 2);
386  out << " for( size_t col = 0, dof = 0; col <"<<numCols<< " ; ++col )" << std::endl;
387  {
388  out << " {" << std::endl;
389  out << " field_type result [ " << sseDimRange << " ] = { ";
390  for( int r = 0; r < sseDimRange-1; ++r ) out << " 0 ,";
391  out << " 0 };" << std::endl << std::endl;
392 
393  out << " const field_type phi[ " << numRows << " ] = {" << std::endl;
394  for( size_t row=0; row< numRows; ++row )
395  {
396  out << " rangeStorage" << row << "[ col ][ 0 ]";
397  if( row < numRows - 1)
398  out << " ," << std::endl;
399  else out << " };" << std::endl;
400  }
401  for( size_t row=0; row< numRows; ++row )
402  {
403  for( int r = 0; r < dimRange; ++r )
404  {
405  out << " result[ " << r << " ] += factor" << row << "[ " << r << " ] * phi[ " << row << " ];" << std::endl;
406  if( sseDimRange != dimRange && r == dimRange - 1 )
407  out << " result[ " << r+1 << " ] += 0.5 * phi[ " << row << " ];" << std::endl;
408  }
409  }
410  for( int r = 0; r < dimRange; ++r ) //, ++colR )
411  out << " dofs[ dof++ ] += result[ " << r << " ];" << std::endl;
412 
413  out << " }" << std::endl;
414  }
415  out << " }" << std::endl << std::endl;
416  out << "};" << std::endl;
417  out << "#endif" << std::endl;
418  }
419 
420  static void evaluateJacobiansCodegen(std::ostream& out,
421  const int dim, const int dimRange, const size_t numRows, const size_t numCols )
422  {
423  out << "#if ! CODEGEN_COMPILE_INNERLOOPS" << std::endl;
424  out << "template <class BaseFunctionSet>" << std::endl;
425  out << "struct EvaluateJacobians<BaseFunctionSet, Fem :: EmptyGeometry, " << dimRange << ", " << numRows << ", " << numCols << ">" << std::endl;
426  out << "{" << std::endl;
427  out << " template< class QuadratureType,"<< std::endl;
428  out << " class JacobianRangeVectorType," << std::endl;
429  out << " class LocalDofVectorType," << std::endl;
430  out << " class JacobianRangeFactorType>" << std::endl;
431  out << " static void eval( const QuadratureType&," << std::endl;
432  out << " const Fem :: EmptyGeometry&," << std::endl;
433  out << " const JacobianRangeVectorType&," << std::endl;
434  out << " const LocalDofVectorType&," << std::endl;
435  out << " JacobianRangeFactorType &)" << std::endl;
436  out << " {" << std::endl;
437  out << " std::cerr << \"ERROR: wrong code generated for VectorialBaseFunctionSet::axpyJacobians\" << std::endl;" << std::endl;
438  out << " abort();" << std::endl;
439  out << " }" << std::endl;
440  out << "};" << std::endl << std::endl;
441  out << "template <class BaseFunctionSet, class Geometry>" << std::endl;
442  out << "struct EvaluateJacobians<BaseFunctionSet, Geometry, " << dimRange << ", " << numRows << ", " << numCols << ">" << std::endl;
443  out << "{" << std::endl;
444  out << " template< class QuadratureType,"<< std::endl;
445  out << " class JacobianRangeVectorType," << std::endl;
446  out << " class LocalDofVectorType," << std::endl;
447  out << " class JacobianRangeFactorType>" << std::endl;
448  out << " static void eval( const QuadratureType& quad," << std::endl;
449  out << " const Geometry& geometry," << std::endl;
450  out << " const JacobianRangeVectorType& jacStorage," << std::endl;
451  out << " const LocalDofVectorType& dofs," << std::endl;
452  out << " JacobianRangeFactorType& jacFactors)" << std::endl;
453  out << " {" << std::endl;
454  out << " evalJac( quad, geometry, jacStorage, dofs, jacFactors, jacFactors[ 0 ] );" << std::endl;
455  out << " }" << std::endl;
456  out << "private:" << std::endl;
457  out << " template< class QuadratureType,"<< std::endl;
458  out << " class JacobianRangeVectorType," << std::endl;
459  out << " class LocalDofVectorType," << std::endl;
460  out << " class JacobianRangeFactorType," << std::endl;
461  out << " class GlobalJacobianRangeType>" << std::endl;
462  out << " static void evalJac( const QuadratureType& quad," << std::endl;
463  out << " const Geometry& geometry," << std::endl;
464  out << " const JacobianRangeVectorType& jacStorage," << std::endl;
465  out << " const LocalDofVectorType& dofs," << std::endl;
466  out << " JacobianRangeFactorType& jacVector," << std::endl;
467  out << " const GlobalJacobianRangeType& )" << std::endl;
468  out << " {" << std::endl;
469  out << " typedef typename JacobianRangeVectorType :: value_type value_type;" << std::endl;
470  out << " for( size_t row = 0; row < " << numRows << " ; ++ row )" << std::endl;
471  out << " {" << std::endl;
472  out << " const value_type& jacStorageRow = jacStorage[ quad.cachingPoint( row ) ];" << std::endl;
473  out << " typedef typename Geometry::JacobianInverseTransposed GeometryJacobianType;" << std::endl;
474  out << " // use reference to GeometryJacobianType to make code compile with SPGrid Geometry" << std::endl;
475  out << " const GeometryJacobianType& gjit = geometry.jacobianInverseTransposed( quad.point( row ) );" << std::endl << std::endl;
476  out << " GlobalJacobianRangeType& result = jacVector[ row ];" << std::endl;
477  out << " result = 0;" << std::endl;
478  out << " typedef typename GlobalJacobianRangeType :: row_type JacobianRangeType;" << std::endl;
479  out << " JacobianRangeType gradPhi;" << std::endl;
480  for( size_t col = 0, colR = 0; col < numCols; ++col )
481  {
482  out << " gjit.mv( jacStorageRow[ "<< col << " ][ 0 ], gradPhi );" << std::endl;
483  for( int r = 0; r < dimRange; ++r, ++colR )
484  {
485  out << " result[ " << r << " ].axpy( dofs[ "<< colR << " ], gradPhi );" << std::endl;
486  }
487  }
488  out << " }" << std::endl;
489  out << " }" << std::endl;
490  out << "};" << std::endl;
491  out << "#endif" << std::endl;
492  }
493 
494  static void axpyJacobianCodegen(std::ostream& out,
495  const int dim, const int dimRange, const size_t numRows, const size_t numCols )
496  {
497  out << "#if ! CODEGEN_COMPILE_INNERLOOPS" << std::endl;
498  out << "template <class BaseFunctionSet>" << std::endl;
499  out << "struct AxpyJacobians<BaseFunctionSet, Fem :: EmptyGeometry, " << dimRange << ", " << numRows << ", " << numCols << ">" << std::endl;
500  out << "{" << std::endl;
501  out << " template< class QuadratureType,"<< std::endl;
502  out << " class JacobianRangeVectorType," << std::endl;
503  out << " class JacobianRangeFactorType," << std::endl;
504  out << " class LocalDofVectorType>" << std::endl;
505  out << " static void axpy( const QuadratureType&," << std::endl;
506  out << " const Fem :: EmptyGeometry&," << std::endl;
507  out << " const JacobianRangeVectorType&," << std::endl;
508  out << " const JacobianRangeFactorType&," << std::endl;
509  out << " LocalDofVectorType&)" << std::endl;
510  out << " {" << std::endl;
511  out << " std::cerr << \"ERROR: wrong code generated for VectorialBaseFunctionSet::axpyJacobians\" << std::endl;" << std::endl;
512  out << " abort();" << std::endl;
513  out << " }" << std::endl;
514  out << "};" << std::endl << std::endl;
515  out << "template <class BaseFunctionSet, class Geometry>" << std::endl;
516  out << "struct AxpyJacobians<BaseFunctionSet, Geometry, " << dimRange << ", " << numRows << ", " << numCols << ">" << std::endl;
517  out << "{" << std::endl;
518  out << " template< class QuadratureType,"<< std::endl;
519  out << " class JacobianRangeVectorType," << std::endl;
520  out << " class JacobianRangeFactorType," << std::endl;
521  out << " class LocalDofVectorType>" << std::endl;
522  out << " static void axpy( const QuadratureType& quad," << std::endl;
523  out << " const Geometry& geometry," << std::endl;
524  out << " const JacobianRangeVectorType& jacStorage," << std::endl;
525  out << " const JacobianRangeFactorType& jacFactors," << std::endl;
526  out << " LocalDofVectorType& dofs)" << std::endl;
527  out << " {" << std::endl;
528  const int sseDimRange = dimRange + (dimRange % 2);
529  //const int sseDimRange = dimRange;// + (dimRange % 2);
530  out << " typedef typename JacobianRangeVectorType :: value_type value_type;" << std::endl;
531  out << " typedef typename JacobianRangeType :: field_type field_type;" << std::endl;
532  const size_t dofs = sseDimRange * numCols ;
533  out << " field_type result [ " << dofs << " ] = {";
534  for( size_t dof = 0 ; dof < dofs-1 ; ++ dof ) out << " 0,";
535  out << " 0 };" << std::endl << std::endl;
536  out << " for( size_t row = 0; row < " << numRows << " ; ++ row )" << std::endl;
537  out << " {" << std::endl;
538  out << " const value_type& jacStorageRow = jacStorage[ quad.cachingPoint( row ) ];" << std::endl;
539  out << " typedef typename Geometry::JacobianInverseTransposed GeometryJacobianType;" << std::endl;
540  out << " // use reference to GeometryJacobianType to make code compile with SPGrid Geometry" << std::endl;
541  out << " const GeometryJacobianType& gjit = geometry.jacobianInverseTransposed( quad.point( row ) );" << std::endl << std::endl;
542  out << " JacobianRangeType jacFactorTmp;" << std::endl;
543  out << " for( int r = 0; r < " << dimRange << " ; ++r )" << std::endl;
544  out << " {"<<std::endl;
545  out << " gjit.mtv( jacFactors[ row ][ r ], jacFactorTmp[ r ] );" << std::endl;
546  out << " }" << std::endl << std::endl;
547  out << " // calculate updates" << std::endl;
548  out << " // rearrange values to have linear memory walk through" << std::endl;
549  out << " const field_type jacFactorInv[ " << dim * sseDimRange << " ] = {" << std::endl;
550  std::string delimiter(" ");
551  for( int i =0; i < dim; ++i )
552  {
553  for( int r = 0; r < dimRange; ++ r )
554  {
555  out << delimiter;
556  delimiter = ", ";
557  out << "jacFactorTmp[ " << r << " ][ " << i << " ]";
558  }
559  if( dimRange < sseDimRange )
560  {
561  out << ", 0.5";
562  }
563  delimiter = ",\n ";
564  }
565  out << " };" << std::endl;
566  for( size_t col = 0, colR = 0; col < numCols; ++col )
567  {
568  out << " {" << std::endl;
569  out << " const field_type phi[ " << dim * sseDimRange << " ] = {" << std::endl;
570  delimiter = std::string(" ");
571  for( int d =0; d < dim; ++d )
572  {
573  for( int r = 0; r < dimRange; ++ r )
574  {
575  out << delimiter;
576  delimiter = ", ";
577  out << " jacStorageRow[ " << col << " ][ 0 ][ " << d << " ]";
578  //FactorTmp[ " << r << " ][ " << i << " ]";
579  }
580  if( dimRange < sseDimRange )
581  {
582  out << ", 0.5";
583  }
584  delimiter = ",\n ";
585  }
586  out << std::endl;
587  out << " };" << std::endl;
588  for( int d = 0, dr = 0 ; d < dim ; ++ d )
589  {
590  for( int r = 0; r < sseDimRange; ++r, ++dr )
591  {
592  out << " result[ " << colR+r << " ] += phi[ " << dr << " ] * jacFactorInv[ " << dr << " ];" << std::endl;
593  }
594  }
595  colR += sseDimRange;
596  out << " }" << std::endl;
597  }
598  /*
599  for( size_t col = 0, colR = 0; col < numCols; ++col )
600  {
601  out << " {" << std::endl;
602  for( int d = 0 ; d < dim ; ++ d )
603  {
604  out << " const field_type& phi"<< d << " = jacStorageRow[ " << col << " ][ 0 ][ " << d << " ];" << std::endl;
605  }
606  for( int r = 0; r < sseDimRange; ++r )
607  {
608  for( int d = 0 ; d < dim ; ++ d )
609  {
610  out << " result[ " << colR+r << " ] += phi" << d << " * jacFactorTmp[ " << r << " ][ " << d << " ];" << std::endl;
611  }
612  }
613  colR += sseDimRange;
614  out << " }" << std::endl;
615  }
616  */
617  out << " }" << std::endl << std::endl;
618 
619  for( size_t col = 0, colD = 0, colR = 0; col < numCols; ++ col )
620  {
621  for(int r = 0; r < dimRange; ++ r, ++colD, ++colR )
622  {
623  out << " dofs[ " << colD << " ] += result[ " << colR << " ];" << std::endl;
624  }
625  if( sseDimRange != dimRange ) ++ colR ;
626  }
627  out << " }" << std::endl << std::endl;
628  out << "};" << std::endl;
629  out << "#endif" << std::endl;
630  }
631  };
632 
633  // if this pre processor variable is defined then
634  // we assume that CODEGENERATOR_REPLACEMENT is CodeGenerator of choice
635 #ifndef FEM_CODEGENERATOR_REPLACEMENT
637 #else
638  typedef FEM_CODEGENERATOR_REPLACEMENT CodeGeneratorType;
639 #endif
640 
641  } // namespace Fem
642 
643 } // namespace Dune
644 
645 #endif // #ifndef DUNE_FEM_BASEFUNCTIONSETS_CODEGEN_HH
std::string path
Definition: readioparams.cc:156
Definition: bindguard.hh:11
bool createDirectory(const std::string &inName)
create a directory
Definition: io.cc:19
DefaultCodeGenerator CodeGeneratorType
Definition: codegen.hh:636
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76
static constexpr T max(T a)
Definition: utility.hh:77
static constexpr T min(T a)
Definition: utility.hh:93
static std::string commonInputPath()
obtain common input path
Definition: io/parameter.hh:431
Definition: grcommon.hh:31
Definition: codegen.hh:23
Definition: codegen.hh:26
void finalize() const
Definition: codegen.hh:162
int exists(const std::string &filename) const
Definition: codegen.hh:243
void clear()
clear all registered entries
Definition: codegen.hh:78
void addEntry(const std::string &fileprefix, codegenfunc_t *codegenfunc, const int dim, const int dimRange, const int quadNop, const int numBase)
Definition: codegen.hh:109
void setFileName(const std::string &filename)
overwrite filename
Definition: codegen.hh:94
void setPath(const std::string &path)
overwrite path
Definition: codegen.hh:91
bool checkAbort() const
Definition: codegen.hh:255
void dumpInfo() const
Definition: codegen.hh:173
static CodegenInfo & instance()
Definition: codegen.hh:71
void addDimRange(const BaseSet *baseSet, const int dimRange)
Definition: codegen.hh:97
default code generator methods
Definition: codegen.hh:282
static void axpyJacobianCodegen(std::ostream &out, const int dim, const int dimRange, const size_t numRows, const size_t numCols)
Definition: codegen.hh:494
static void axpyCodegen(std::ostream &out, const int dim, const int dimRange, const size_t numRows, const size_t numCols)
Definition: codegen.hh:356
static void evaluateJacobiansCodegen(std::ostream &out, const int dim, const int dimRange, const size_t numRows, const size_t numCols)
Definition: codegen.hh:420
static void evaluateCodegen(std::ostream &out, const int dim, const int dimRange, const size_t numRows, const size_t numCols)
Definition: codegen.hh:283