Rheolef  7.1
an efficient C++ finite element environment
basis_symbolic_cxx.cc
Go to the documentation of this file.
1 #include "basis_symbolic.h"
22 #include <sstream>
23 using namespace rheolef;
24 using namespace std;
25 using namespace GiNaC;
26 
27 static
28 void
29 put_gpl (ostream& out)
30 {
31  out << "///" << endl
32  << "/// This file is part of Rheolef." << endl
33  << "///" << endl
34  << "/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>" << endl
35  << "///" << endl
36  << "/// Rheolef is free software; you can redistribute it and/or modify" << endl
37  << "/// it under the terms of the GNU General Public License as published by" << endl
38  << "/// the Free Software Foundation; either version 2 of the License, or" << endl
39  << "/// (at your option) any later version." << endl
40  << "///" << endl
41  << "/// Rheolef is distributed in the hope that it will be useful," << endl
42  << "/// but WITHOUT ANY WARRANTY; without even the implied warranty of" << endl
43  << "/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the" << endl
44  << "/// GNU General Public License for more details." << endl
45  << "///" << endl
46  << "/// You should have received a copy of the GNU General Public License" << endl
47  << "/// along with Rheolef; if not, write to the Free Software" << endl
48  << "/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA" << endl
49  << "///" << endl
50  << "/// =========================================================================" << endl
51  ;
52 }
53 ex
55 {
56  size_type d = dimension();
57  ex expr = expr0;
58  // first optimize;
59  if (d > 0) expr = collect(expr,x);
60  if (d > 1) expr = collect(expr,y);
61  if (d > 2) expr = collect(expr,z);
62  // then subst symbols:
63  if (d > 0) expr = expr.subs(x == symbol("hat_x[0]"));
64  if (d > 1) expr = expr.subs(y == symbol("hat_x[1]"));
65  if (d > 2) expr = expr.subs(z == symbol("hat_x[2]"));
66  return expr;
67 }
68 void
70 {
71  if (size() == 0) return;
72  stringstream class_name;
73  class_name << "basis_" << name() << "_" << _hat_K.name();
74  out << "template<class T>" << endl
75  << "class " << class_name.str() << " {" << endl
76  << "public:" << endl
77  << " typedef basis_rep<T> base;" << endl
78  << " typedef typename base::size_type size_type;" << endl
79  << " static void evaluate (const point_basic<T>& hat_x, Eigen::Matrix<T,Eigen::Dynamic,1>& values);" << endl
80  << " static void grad_evaluate (const point_basic<T>& hat_x, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& values);" << endl
81  << " static void hat_node (Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&);" << endl
82  << "};" << endl;
83 }
84 void
86 {
87  if (size() == 0) return;
88  size_type d = dimension();
89  stringstream class_name;
90  class_name << "basis_" << name() << "_" << _hat_K.name();
91  // --------------------------------------------------
92  // evaluate
93  // --------------------------------------------------
94  out << "template<class T>" << endl
95  << "void" << endl
96  << class_name.str() << "<T>::evaluate (" << endl
97  << " const point_basic<T>& hat_x," << endl
98  << " Eigen::Matrix<T,Eigen::Dynamic,1>& values)" << endl
99  << "{" << endl
100  << " values.resize(" << _basis.size() << ");" << endl;
101  for (size_type i = 0; i < _basis.size(); i++) {
102  ex expr = indexed_symbol (_basis[i]);
103  stringstream idx_stream;
104  out << " values[" << i << "] = ";
105  expr.print(print_csrc_double(cout));
106  out << ";" << endl;
107  }
108  out << "}" << endl;
109  // --------------------------------------------------
110  // grad_evaluate
111  // --------------------------------------------------
112  out << "template<class T>" << endl
113  << "void" << endl
114  << class_name.str() << "<T>::grad_evaluate (" << endl
115  << " const point_basic<T>& hat_x," << endl
116  << " Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& values)" << endl
117  << "{" << endl
118  << " values.resize(" << _basis.size() << ");" << endl;
119  for (size_type i = 0; i < _basis.size(); i++) {
120  for (size_type j = 0; j < d; j++) {
121  ex expr = indexed_symbol (_grad_basis[i][j]);
122  out << " values[" << i << "][" << j << "] = ";
123  expr.print(print_csrc_double(cout));
124  out << ";" << endl;
125  }
126  }
127  out << "}" << endl;
128  // --------------------------------------------------
129  // hat_node
130  // --------------------------------------------------
131  out << "template<class T>" << endl
132  << "void" << endl
133  << class_name.str() << "<T>::hat_node (Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x)" << endl
134  << "{" << endl
135  << " x.resize(" << _basis.size() << ");" << endl;
136  for (size_type i = 0; i < _basis.size(); i++) {
137  out << " x[" << i << "] = point_basic<T>(";
138  for (size_type j = 0; j < d; j++) {
139  ex expr = indexed_symbol (_node[i][j]);
140  expr.print(print_csrc_double(cout));
141  if (j != d-1) out << ", ";
142  }
143  out << ");" << endl;
144  }
145  out << "}" << endl;
146 }
147 void
149 {
150  stringstream class_name;
151  class_name << "basis_" << name();
152 
153  out << "// file automatically generated by \"" << __FILE__ << "\"" << endl;
154  put_gpl(out);
155  out << "#ifndef _RHEOLEF_" << name() << "_H" << endl
156  << "#define _RHEOLEF_" << name() << "_H" << endl
157  << "#include \"rheolef/basis.h\"" << endl
158  << "#include \"rheolef/tensor.h\"" << endl;
159  if (have_index_parameter() || name() == "bubble") { // P0, P1, bubble
160  out << "#include \"Pk_get_local_idof_on_side.icc\"" << endl
161  << "#include \"basis_fem_Pk_lagrange.h\"" << endl;
162  }
163  out << "namespace rheolef {" << endl
164  << endl
165  << "template<class T>" << endl
166  << "class " << class_name.str() << ": public basis_rep<T> {" << endl
167  << "public:" << endl
168  << " typedef basis_rep<T> base;" << endl
169  << " typedef typename base::size_type size_type;" << endl;
170  if (!have_continuous_feature()) { // P0, bubble, P1qd, ...
171  out << " " << class_name.str() << "();" << endl;
172  } else {
173  out << " " << class_name.str() << " (const basis_option& sopt);" << endl;
174  }
175  out << " ~" << class_name.str() << "();" << endl;
176  if (!have_index_parameter()) { // bubble, P1qd, ...
177  out << " std::string name() const { return family_name(); }" << endl
178  << " bool have_index_parameter() const { return false; }" << endl
179  << " std::string family_name() const { return \"" << name() << "\"; }" << endl;
180  } else {
181  out << " bool have_index_parameter() const { return true; }" << endl
182  << " std::string family_name() const { return \"" << family_name() << "\"; }" << endl
183  << " void local_idof_on_side (" << endl
184  << " reference_element hat_K," << endl
185  << " const side_information_type& sid," << endl
186  << " Eigen::Matrix<size_type,Eigen::Dynamic,1>& loc_idof) const" << endl
187  << " {" << endl
188  << " details::Pk_get_local_idof_on_side (hat_K, sid, degree(), loc_idof);" << endl
189  << " }" << endl;
190  }
191  if (name() == "bubble") { // P0, P1qd, ...
192  out << " bool have_compact_support_inside_element() const { return true; }" << endl;
193  }
194  out << " size_type degree() const;" << endl
195  << " bool is_nodal() const { return true; }" << endl
196 #ifdef TO_CLEAN
197  << " size_type ndof (reference_element hat_K) const;" << endl
198 #endif // TO_CLEAN
199  << " void evaluate (" << endl
200  << " reference_element hat_K," << endl
201  << " const point_basic<T>& hat_x," << endl
202  << " Eigen::Matrix<T,Eigen::Dynamic,1>& values) const;" << endl
203  << " void grad_evaluate (" << endl
204  << " reference_element hat_K," << endl
205  << " const point_basic<T>& hat_x," << endl
206  << " Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& values) const;" << endl
207  << " void _compute_dofs (" << endl
208  << " reference_element hat_K," << endl
209  << " const Eigen::Matrix<T,Eigen::Dynamic,1>& f_xnod," << endl
210  << " Eigen::Matrix<T,Eigen::Dynamic,1>& dof) const;" << endl
211  << " void _initialize_cstor_sizes() const;" << endl
212  << " void _initialize_data (reference_element hat_K) const;" << endl
213  << " const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& hat_node(reference_element hat_K) const {" << endl
214  << " base::_initialize_data_guard(hat_K);" << endl
215  << " return _hat_node [hat_K.variant()]; }" << endl
216  << " mutable std::array<" << endl
217  << " Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>," << endl
218  << " reference_element::max_variant> _hat_node;" << endl
219  << "};" << endl
220  << "} // namespace rheolef" << endl
221  << "#endif // _RHEOLEF_" << name() << "_H" << endl
222  << endl;
223 }
224 void
226 {
227  out << "// file automatically generated by \"" << __FILE__ << "\"" << endl;
228  put_gpl(out);
229  out << "#include \"" << name() << ".h\"" << endl
230  << "#include \"piola_fem_lagrange.h\"" << endl
231  << "namespace rheolef {" << endl
232  << "using namespace std;" << endl;
233  for (size_type i = 0; i < reference_element::max_variant; i++) {
234  operator[](i).put_cxx_header(out);
235  }
236  for (size_type i = 0; i < reference_element::max_variant; i++) {
237  operator[](i).put_cxx_body(out);
238  }
239  stringstream class_name;
240  class_name << "basis_" << name();
241 
242  // --------------------------------------------------
243  // allocator
244  // --------------------------------------------------
245  out << "template<class T>" << endl;
246  if (!have_continuous_feature()) { // P0, bubble, P1qd, ...
247  out << class_name.str() << "<T>::" << class_name.str() << "()" << endl
248  << " : base(basis_option()), _hat_node()" << endl
249  << "{" << endl;
250  if (name() != "bubble") { // P0, P1qd, ...
251  out << " base::_sopt.set_continuous(false);" << endl;
252  }
253  } else {
254  out << class_name.str() << "<T>::" << class_name.str() << " (const basis_option& sopt)" << endl
255  << " : base(sopt), _hat_node()" << endl
256  << "{" << endl;
257  }
258  out << " _initialize_cstor_sizes();" << endl;
259  if (name() == "P1") {
260  out << " base::_name = base::standard_naming (family_name(), 1, base::_sopt);" << endl;
261  } else {
262  out << " base::_name = \"" << name() << "\";" << endl;
263  }
264  out << " base::_piola_fem.piola_fem<T>::base::operator= (new_macro(piola_fem_lagrange<T>));" << endl
265  << "}" << endl;
266  // --------------------------------------------------
267  // destructor
268  // --------------------------------------------------
269  out << "template<class T>" << endl
270  << class_name.str() << "<T>::~" << class_name.str() << "()" << endl
271  << "{" << endl
272  << "}" << endl;
273  // --------------------------------------------------
274  // degree
275  // --------------------------------------------------
276  out << "template<class T>" << endl
277  << "typename " << class_name.str() << "<T>::size_type" << endl
278  << class_name.str() << "<T>::degree () const" << endl
279  << "{" << endl
280  << " return " << degree() << ";" << endl
281  << "}" << endl;
282 #ifdef TO_CLEAN
283  // --------------------------------------------------
284  // ndof
285  // --------------------------------------------------
286  out << "template<class T>" << endl
287  << "typename " << class_name.str() << "<T>::size_type" << endl
288  << class_name.str() << "<T>::ndof (" << endl
289  << " reference_element hat_K) const" << endl
290  << "{" << endl
291  << " switch (hat_K.variant()) {" << endl;
292  for (size_type i = 0; i < reference_element::max_variant; i++) {
293  const basis_symbolic_nodal_on_geo& b = operator[](i);
294  if (b.size() == 0) continue;
295  out << " case reference_element::" << b.hat_K().name() << ": {" << endl
296  << " return " << b.size() << ";" << endl
297  << " }" << endl;
298  }
299  out << " default : {" << endl
300  << " error_macro (\"size: unsupported `\" << hat_K.name() << \"' element type\");" << endl
301  << " return 0;" << endl
302  << " }" << endl
303  << " }" << endl
304  << "}" << endl;
305 #endif // TO_CLEAN
306  // --------------------------------------------------
307  // evaluate
308  // --------------------------------------------------
309  out << "template<class T>" << endl
310  << "void" << endl
311  << class_name.str() << "<T>::evaluate (" << endl
312  << " reference_element hat_K," << endl
313  << " const point_basic<T>& hat_x," << endl
314  << " Eigen::Matrix<T,Eigen::Dynamic,1>& values) const" << endl
315  << "{" << endl
316  << " switch (hat_K.variant()) {" << endl;
317  for (size_type i = 0; i < reference_element::max_variant; i++) {
318  const basis_symbolic_nodal_on_geo& b = operator[](i);
319  if (b.size() == 0) continue;
320  out << " case reference_element::" << b.hat_K().name() << ": {" << endl
321  << " return " << class_name.str() << "_" << b.hat_K().name() << "<T>::evaluate (hat_x, values);" << endl
322  << " }" << endl;
323  }
324  out << " default : {" << endl
325  << " error_macro (\"evaluate: unsupported `\" << hat_K.name() << \"' element type\");" << endl
326  << " }" << endl
327  << " }" << endl
328  << "}" << endl;
329  // --------------------------------------------------
330  // grad_evaluate
331  // --------------------------------------------------
332  out << "template<class T>" << endl
333  << "void" << endl
334  << class_name.str() << "<T>::grad_evaluate (" << endl
335  << " reference_element hat_K," << endl
336  << " const point_basic<T>& hat_x," << endl
337  << " Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& values) const" << endl
338  << "{" << endl
339  << " switch (hat_K.variant()) {" << endl;
340  for (size_type i = 0; i < reference_element::max_variant; i++) {
341  const basis_symbolic_nodal_on_geo& b = operator[](i);
342  if (b.size() == 0) continue;
343  out << " case reference_element::" << b.hat_K().name() << ": {" << endl
344  << " return " << class_name.str() << "_" << b.hat_K().name() << "<T>::grad_evaluate (hat_x, values);" << endl
345  << " }" << endl;
346  }
347  out << " default : {" << endl
348  << " error_macro (\"grad_evaluate: unsupported `\" << hat_K.name() << \"' element type\");" << endl
349  << " }" << endl
350  << " }" << endl
351  << "}" << endl;
352  // --------------------------------------------------
353  // interpolate
354  // --------------------------------------------------
355  out << "template<class T>" << endl
356  << "void" << endl
357  << class_name.str() << "<T>::_compute_dofs (" << endl
358  << " reference_element hat_K," << endl
359  << " const Eigen::Matrix<T,Eigen::Dynamic,1>& f_xnod," << endl
360  << " Eigen::Matrix<T,Eigen::Dynamic,1>& dof) const" << endl
361  << "{" << endl
362  << " dof = f_xnod;" << endl
363  << "}" << endl;
364  // --------------------------------------------------
365  // init sizes
366  // --------------------------------------------------
367  out << "template<class T>" << endl
368  << "void" << endl
369  << class_name.str() << "<T>::_initialize_cstor_sizes() const" << endl
370  << "{" << endl;
371  if (have_index_parameter()) { // P0, P1
372  out << " basis_fem_Pk_lagrange<T>::initialize_local_first (" << endl
373  << " " << degree() << "," << endl
374  << " base::is_continuous()," << endl
375  << " base::_ndof_on_subgeo," << endl
376  << " base::_nnod_on_subgeo," << endl
377  << " base::_first_idof_by_dimension," << endl
378  << " base::_first_inod_by_dimension);" << endl;
379  } else if (name() == "bubble") {
380  out << " basis_fem_Pk_lagrange<T>::initialize_local_first (" << endl
381  << " 0," << endl
382  << " false," << endl
383  << " base::_ndof_on_subgeo," << endl
384  << " base::_nnod_on_subgeo," << endl
385  << " base::_first_idof_by_dimension," << endl
386  << " base::_first_inod_by_dimension);" << endl;
387  } else {
388  out << " fatal_macro(\"::_initialize_cstor_sizes: not yet\"); /* TODO */" << endl;
389  }
390  out << "}" << endl;
391  // --------------------------------------------------
392  // hat_node (vectorial)
393  // --------------------------------------------------
394  out << "template<class T>" << endl
395  << "void" << endl
396  << class_name.str() << "<T>::_initialize_data(" << endl
397  << " reference_element hat_K) const" << endl
398  << "{" << endl
399  << " switch (hat_K.variant()) {" << endl;
400  for (size_type i = 0; i < reference_element::max_variant; i++) {
401  const basis_symbolic_nodal_on_geo& b = operator[](i);
402  if (b.size() == 0) continue;
403  out << " case reference_element::" << b.hat_K().name() << ": {" << endl
404  << " return " << class_name.str() << "_" << b.hat_K().name()
405  << "<T>::hat_node (_hat_node[hat_K.variant()]);" << endl
406  << " }" << endl;
407  }
408  out << " default : {" << endl
409  << " error_macro (\"hat_node: unsupported `\" << hat_K.name() << \"' element type\");" << endl
410  << " }" << endl
411  << " }" << endl
412  << "}" << endl;
413  // --------------------------------------------------
414  // call to a constructor
415  // --------------------------------------------------
416  out << "// instantiation in library:" << endl
417  << "template class basis_" << name() << "<Float>;" << endl
418  << "} // namespace rheolef" << endl;
419 }
420 void
421 basis_symbolic_nodal::put_cxx_main (int argc, char**argv) const
422 {
423  if (argc <= 1 || string(argv[1]) == "-h") {
424  put_cxx_header (cout);
425  } else {
426  put_cxx_body (cout);
427  }
428 }
rheolef::io::out
@ out
Definition: rheostream.h:167
mkgeo_ball.expr
expr
Definition: mkgeo_ball.sh:361
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
rheolef::basis_symbolic_nodal::size_type
basis_symbolic_nodal_on_geo::size_type size_type
Definition: basis_symbolic.h:123
rheolef::basis_symbolic_nodal_on_geo
Definition: basis_symbolic.h:39
rheolef::basis_symbolic_nodal::put_cxx_body
void put_cxx_body(std::ostream &out) const
Definition: basis_symbolic_cxx.cc:225
rheolef::basis_symbolic_nodal_on_geo::put_cxx_header
void put_cxx_header(std::ostream &out) const
Definition: basis_symbolic_cxx.cc:69
dimension
const size_t dimension
Definition: edge.icc:64
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::reference_element::max_variant
static const variant_type max_variant
Definition: reference_element.h:82
rheolef::basis_symbolic_nodal_on_geo::indexed_symbol
polynom_type indexed_symbol(const polynom_type &expr0) const
Definition: basis_symbolic_cxx.cc:54
basis_symbolic.h
rheolef::basis_symbolic_nodal::put_cxx_main
void put_cxx_main(int argc, char **argv) const
Definition: basis_symbolic_cxx.cc:421
rheolef::basis_symbolic_nodal_on_geo::put_cxx_body
void put_cxx_body(std::ostream &out) const
Definition: basis_symbolic_cxx.cc:85
rheolef::basis_symbolic_nodal::put_cxx_header
void put_cxx_header(std::ostream &out) const
Definition: basis_symbolic_cxx.cc:148
rheolef::basis_symbolic_nodal_on_geo::size_type
std::vector< int >::size_type size_type
Definition: basis_symbolic.h:44
rheolef::std
Definition: vec_expr_v2.h:402
mkgeo_contraction.name
string name
Definition: mkgeo_contraction.sh:133