Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
basis_symbolic_cxx.cc
Go to the documentation of this file.
1
21#include "basis_symbolic.h"
22#include <sstream>
23using namespace rheolef;
24using namespace std;
25using namespace GiNaC;
26
27static
28void
29put_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}
53ex
55{
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}
68void
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}
84void
86{
87 if (size() == 0) return;
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}
147void
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}
224void
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_internal," << endl
376 << " base::_ndof_on_subgeo," << endl
377 << " base::_nnod_on_subgeo_internal," << endl
378 << " base::_nnod_on_subgeo," << endl
379 << " base::_first_idof_by_dimension_internal," << endl
380 << " base::_first_idof_by_dimension," << endl
381 << " base::_first_inod_by_dimension_internal," << endl
382 << " base::_first_inod_by_dimension);" << endl;
383 } else if (name() == "bubble") {
384 out << " basis_fem_Pk_lagrange<T>::initialize_local_first (" << endl
385 << " 0," << endl
386 << " false," << endl
387 << " base::_ndof_on_subgeo_internal," << endl
388 << " base::_ndof_on_subgeo," << endl
389 << " base::_nnod_on_subgeo_internal," << endl
390 << " base::_nnod_on_subgeo," << endl
391 << " base::_first_idof_by_dimension_internal," << endl
392 << " base::_first_idof_by_dimension," << endl
393 << " base::_first_inod_by_dimension_internal," << endl
394 << " base::_first_inod_by_dimension);" << endl;
395 } else {
396 out << " fatal_macro(\"::_initialize_cstor_sizes: not yet\"); /* TODO */" << endl;
397 }
398 out << "}" << endl;
399 // --------------------------------------------------
400 // hat_node (vectorial)
401 // --------------------------------------------------
402 out << "template<class T>" << endl
403 << "void" << endl
404 << class_name.str() << "<T>::_initialize_data(" << endl
405 << " reference_element hat_K) const" << endl
406 << "{" << endl
407 << " switch (hat_K.variant()) {" << endl;
408 for (size_type i = 0; i < reference_element::max_variant; i++) {
409 const basis_symbolic_nodal_on_geo& b = operator[](i);
410 if (b.size() == 0) continue;
411 out << " case reference_element::" << b.hat_K().name() << ": {" << endl
412 << " return " << class_name.str() << "_" << b.hat_K().name()
413 << "<T>::hat_node (_hat_node[hat_K.variant()]);" << endl
414 << " }" << endl;
415 }
416 out << " default : {" << endl
417 << " error_macro (\"hat_node: unsupported `\" << hat_K.name() << \"' element type\");" << endl
418 << " }" << endl
419 << " }" << endl
420 << "}" << endl;
421 // --------------------------------------------------
422 // call to a constructor
423 // --------------------------------------------------
424 out << "// instantiation in library:" << endl
425 << "template class basis_" << name() << "<Float>;" << endl
426 << "} // namespace rheolef" << endl;
427}
428void
429basis_symbolic_nodal::put_cxx_main (int argc, char**argv) const
430{
431 if (argc <= 1 || string(argv[1]) == "-h") {
432 put_cxx_header (cout);
433 } else {
434 put_cxx_body (cout);
435 }
436}
polynom_type indexed_symbol(const polynom_type &expr0) const
std::vector< point_basic< polynom_type > > _grad_basis
std::vector< point_basic< GiNaC::ex > > _node
std::vector< polynom_type > _basis
void put_cxx_body(std::ostream &out) const
std::vector< int >::size_type size_type
void put_cxx_header(std::ostream &out) const
void put_cxx_main(int argc, char **argv) const
basis_symbolic_nodal_on_geo::size_type size_type
void put_cxx_body(std::ostream &out) const
std::string family_name() const
void put_cxx_header(std::ostream &out) const
static const variant_type max_variant
This file is part of Rheolef.
STL namespace.