basis_symbolic_cxx.cc
Go to the documentation of this file.
#include "basis_symbolic.h"
#include <sstream>
using namespace rheolef;
using namespace std;
using namespace GiNaC;
static
void
put_gpl (ostream& out)
{
out << "///" << endl
<< "/// This file is part of Rheolef." << endl
<< "///" << endl
<< "/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>" << endl
<< "///" << endl
<< "/// Rheolef is free software; you can redistribute it and/or modify" << endl
<< "/// it under the terms of the GNU General Public License as published by" << endl
<< "/// the Free Software Foundation; either version 2 of the License, or" << endl
<< "/// (at your option) any later version." << endl
<< "///" << endl
<< "/// Rheolef is distributed in the hope that it will be useful," << endl
<< "/// but WITHOUT ANY WARRANTY; without even the implied warranty of" << endl
<< "/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the" << endl
<< "/// GNU General Public License for more details." << endl
<< "///" << endl
<< "/// You should have received a copy of the GNU General Public License" << endl
<< "/// along with Rheolef; if not, write to the Free Software" << endl
<< "/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA" << endl
<< "///" << endl
<< "/// =========================================================================" << endl
;
}
ex
basis_symbolic_nodal_on_geo::indexed_symbol (const ex& expr0) const
{
size_type d = dimension();
ex expr = expr0;
// first optimize;
if (d > 0) expr = collect(expr,x);
if (d > 1) expr = collect(expr,y);
if (d > 2) expr = collect(expr,z);
// then subst symbols:
if (d > 0) expr = expr.subs(x == symbol("hat_x[0]"));
if (d > 1) expr = expr.subs(y == symbol("hat_x[1]"));
if (d > 2) expr = expr.subs(z == symbol("hat_x[2]"));
return expr;
}
void
basis_symbolic_nodal_on_geo::put_cxx_header(ostream& out) const
{
if (size() == 0) return;
stringstream class_name;
class_name << "basis_" << name() << "_" << _hat_K.name();
out << "template<class T>" << endl
<< "class " << class_name.str() << " {" << endl
<< "public:" << endl
<< " typedef basis_rep<T> base;" << endl
<< " typedef typename base::size_type size_type;" << endl
<< " static void evaluate (const point_basic<T>& hat_x, Eigen::Matrix<T,Eigen::Dynamic,1>& values);" << endl
<< " static void grad_evaluate (const point_basic<T>& hat_x, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& values);" << endl
<< " static void hat_node (Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&);" << endl
<< "};" << endl;
}
void
basis_symbolic_nodal_on_geo::put_cxx_body(ostream& out) const
{
if (size() == 0) return;
size_type d = dimension();
stringstream class_name;
class_name << "basis_" << name() << "_" << _hat_K.name();
// --------------------------------------------------
// evaluate
// --------------------------------------------------
out << "template<class T>" << endl
<< "void" << endl
<< class_name.str() << "<T>::evaluate (" << endl
<< " const point_basic<T>& hat_x," << endl
<< " Eigen::Matrix<T,Eigen::Dynamic,1>& values)" << endl
<< "{" << endl
<< " values.resize(" << _basis.size() << ");" << endl;
for (size_type i = 0; i < _basis.size(); i++) {
ex expr = indexed_symbol (_basis[i]);
stringstream idx_stream;
out << " values[" << i << "] = ";
expr.print(print_csrc_double(cout));
out << ";" << endl;
}
out << "}" << endl;
// --------------------------------------------------
// grad_evaluate
// --------------------------------------------------
out << "template<class T>" << endl
<< "void" << endl
<< class_name.str() << "<T>::grad_evaluate (" << endl
<< " const point_basic<T>& hat_x," << endl
<< " Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& values)" << endl
<< "{" << endl
<< " values.resize(" << _basis.size() << ");" << endl;
for (size_type i = 0; i < _basis.size(); i++) {
for (size_type j = 0; j < d; j++) {
ex expr = indexed_symbol (_grad_basis[i][j]);
out << " values[" << i << "][" << j << "] = ";
expr.print(print_csrc_double(cout));
out << ";" << endl;
}
}
out << "}" << endl;
// --------------------------------------------------
// hat_node
// --------------------------------------------------
out << "template<class T>" << endl
<< "void" << endl
<< class_name.str() << "<T>::hat_node (Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x)" << endl
<< "{" << endl
<< " x.resize(" << _basis.size() << ");" << endl;
for (size_type i = 0; i < _basis.size(); i++) {
out << " x[" << i << "] = point_basic<T>(";
for (size_type j = 0; j < d; j++) {
ex expr = indexed_symbol (_node[i][j]);
expr.print(print_csrc_double(cout));
if (j != d-1) out << ", ";
}
out << ");" << endl;
}
out << "}" << endl;
}
void
basis_symbolic_nodal::put_cxx_header(ostream& out) const
{
stringstream class_name;
class_name << "basis_" << name();
out << "// file automatically generated by \"" << __FILE__ << "\"" << endl;
put_gpl(out);
out << "#ifndef _RHEOLEF_" << name() << "_H" << endl
<< "#define _RHEOLEF_" << name() << "_H" << endl
<< "#include \"rheolef/basis.h\"" << endl
<< "#include \"rheolef/tensor.h\"" << endl;
if (have_index_parameter() || name() == "bubble") { // P0, P1, bubble
out << "#include \"Pk_get_local_idof_on_side.icc\"" << endl
<< "#include \"basis_fem_Pk_lagrange.h\"" << endl;
}
out << "namespace rheolef {" << endl
<< endl
<< "template<class T>" << endl
<< "class " << class_name.str() << ": public basis_rep<T> {" << endl
<< "public:" << endl
<< " typedef basis_rep<T> base;" << endl
<< " typedef typename base::size_type size_type;" << endl;
if (!have_continuous_feature()) { // P0, bubble, P1qd, ...
out << " " << class_name.str() << "();" << endl;
} else {
out << " " << class_name.str() << " (const basis_option& sopt);" << endl;
}
out << " ~" << class_name.str() << "();" << endl;
if (!have_index_parameter()) { // bubble, P1qd, ...
out << " std::string name() const { return family_name(); }" << endl
<< " bool have_index_parameter() const { return false; }" << endl
<< " std::string family_name() const { return \"" << name() << "\"; }" << endl;
} else {
out << " bool have_index_parameter() const { return true; }" << endl
<< " std::string family_name() const { return \"" << family_name() << "\"; }" << endl
<< " void local_idof_on_side (" << endl
<< " reference_element hat_K," << endl
<< " const side_information_type& sid," << endl
<< " Eigen::Matrix<size_type,Eigen::Dynamic,1>& loc_idof) const" << endl
<< " {" << endl
<< " details::Pk_get_local_idof_on_side (hat_K, sid, degree(), loc_idof);" << endl
<< " }" << endl;
}
if (name() == "bubble") { // P0, P1qd, ...
out << " bool have_compact_support_inside_element() const { return true; }" << endl;
}
out << " size_type degree() const;" << endl
<< " bool is_nodal() const { return true; }" << endl
#ifdef TO_CLEAN
<< " size_type ndof (reference_element hat_K) const;" << endl
#endif // TO_CLEAN
<< " void evaluate (" << endl
<< " reference_element hat_K," << endl
<< " const point_basic<T>& hat_x," << endl
<< " Eigen::Matrix<T,Eigen::Dynamic,1>& values) const;" << endl
<< " void grad_evaluate (" << endl
<< " reference_element hat_K," << endl
<< " const point_basic<T>& hat_x," << endl
<< " Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& values) const;" << endl
<< " void _compute_dofs (" << endl
<< " reference_element hat_K," << endl
<< " const Eigen::Matrix<T,Eigen::Dynamic,1>& f_xnod," << endl
<< " Eigen::Matrix<T,Eigen::Dynamic,1>& dof) const;" << endl
<< " void _initialize_cstor_sizes() const;" << endl
<< " void _initialize_data (reference_element hat_K) const;" << endl
<< " const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& hat_node(reference_element hat_K) const {" << endl
<< " base::_initialize_data_guard(hat_K);" << endl
<< " return _hat_node [hat_K.variant()]; }" << endl
<< " mutable std::array<" << endl
<< " Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>," << endl
<< " reference_element::max_variant> _hat_node;" << endl
<< "};" << endl
<< "} // namespace rheolef" << endl
<< "#endif // _RHEOLEF_" << name() << "_H" << endl
<< endl;
}
void
basis_symbolic_nodal::put_cxx_body(ostream& out) const
{
out << "// file automatically generated by \"" << __FILE__ << "\"" << endl;
put_gpl(out);
out << "#include \"" << name() << ".h\"" << endl
<< "#include \"piola_fem_lagrange.h\"" << endl
<< "namespace rheolef {" << endl
<< "using namespace std;" << endl;
for (size_type i = 0; i < reference_element::max_variant; i++) {
operator[](i).put_cxx_header(out);
}
for (size_type i = 0; i < reference_element::max_variant; i++) {
operator[](i).put_cxx_body(out);
}
stringstream class_name;
class_name << "basis_" << name();
// --------------------------------------------------
// allocator
// --------------------------------------------------
out << "template<class T>" << endl;
if (!have_continuous_feature()) { // P0, bubble, P1qd, ...
out << class_name.str() << "<T>::" << class_name.str() << "()" << endl
<< " : base(basis_option()), _hat_node()" << endl
<< "{" << endl;
if (name() != "bubble") { // P0, P1qd, ...
out << " base::_sopt.set_continuous(false);" << endl;
}
} else {
out << class_name.str() << "<T>::" << class_name.str() << " (const basis_option& sopt)" << endl
<< " : base(sopt), _hat_node()" << endl
<< "{" << endl;
}
out << " _initialize_cstor_sizes();" << endl;
if (name() == "P1") {
out << " base::_name = base::standard_naming (family_name(), 1, base::_sopt);" << endl;
} else {
out << " base::_name = \"" << name() << "\";" << endl;
}
out << " base::_piola_fem.piola_fem<T>::base::operator= (new_macro(piola_fem_lagrange<T>));" << endl
<< "}" << endl;
// --------------------------------------------------
// destructor
// --------------------------------------------------
out << "template<class T>" << endl
<< class_name.str() << "<T>::~" << class_name.str() << "()" << endl
<< "{" << endl
<< "}" << endl;
// --------------------------------------------------
// degree
// --------------------------------------------------
out << "template<class T>" << endl
<< "typename " << class_name.str() << "<T>::size_type" << endl
<< class_name.str() << "<T>::degree () const" << endl
<< "{" << endl
<< " return " << degree() << ";" << endl
<< "}" << endl;
#ifdef TO_CLEAN
// --------------------------------------------------
// ndof
// --------------------------------------------------
out << "template<class T>" << endl
<< "typename " << class_name.str() << "<T>::size_type" << endl
<< class_name.str() << "<T>::ndof (" << endl
<< " reference_element hat_K) const" << endl
<< "{" << endl
<< " switch (hat_K.variant()) {" << endl;
for (size_type i = 0; i < reference_element::max_variant; i++) {
const basis_symbolic_nodal_on_geo& b = operator[](i);
if (b.size() == 0) continue;
out << " case reference_element::" << b.hat_K().name() << ": {" << endl
<< " return " << b.size() << ";" << endl
<< " }" << endl;
}
out << " default : {" << endl
<< " error_macro (\"size: unsupported `\" << hat_K.name() << \"' element type\");" << endl
<< " return 0;" << endl
<< " }" << endl
<< " }" << endl
<< "}" << endl;
#endif // TO_CLEAN
// --------------------------------------------------
// evaluate
// --------------------------------------------------
out << "template<class T>" << endl
<< "void" << endl
<< class_name.str() << "<T>::evaluate (" << endl
<< " reference_element hat_K," << endl
<< " const point_basic<T>& hat_x," << endl
<< " Eigen::Matrix<T,Eigen::Dynamic,1>& values) const" << endl
<< "{" << endl
<< " switch (hat_K.variant()) {" << endl;
for (size_type i = 0; i < reference_element::max_variant; i++) {
const basis_symbolic_nodal_on_geo& b = operator[](i);
if (b.size() == 0) continue;
out << " case reference_element::" << b.hat_K().name() << ": {" << endl
<< " return " << class_name.str() << "_" << b.hat_K().name() << "<T>::evaluate (hat_x, values);" << endl
<< " }" << endl;
}
out << " default : {" << endl
<< " error_macro (\"evaluate: unsupported `\" << hat_K.name() << \"' element type\");" << endl
<< " }" << endl
<< " }" << endl
<< "}" << endl;
// --------------------------------------------------
// grad_evaluate
// --------------------------------------------------
out << "template<class T>" << endl
<< "void" << endl
<< class_name.str() << "<T>::grad_evaluate (" << endl
<< " reference_element hat_K," << endl
<< " const point_basic<T>& hat_x," << endl
<< " Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& values) const" << endl
<< "{" << endl
<< " switch (hat_K.variant()) {" << endl;
for (size_type i = 0; i < reference_element::max_variant; i++) {
const basis_symbolic_nodal_on_geo& b = operator[](i);
if (b.size() == 0) continue;
out << " case reference_element::" << b.hat_K().name() << ": {" << endl
<< " return " << class_name.str() << "_" << b.hat_K().name() << "<T>::grad_evaluate (hat_x, values);" << endl
<< " }" << endl;
}
out << " default : {" << endl
<< " error_macro (\"grad_evaluate: unsupported `\" << hat_K.name() << \"' element type\");" << endl
<< " }" << endl
<< " }" << endl
<< "}" << endl;
// --------------------------------------------------
// interpolate
// --------------------------------------------------
out << "template<class T>" << endl
<< "void" << endl
<< class_name.str() << "<T>::_compute_dofs (" << endl
<< " reference_element hat_K," << endl
<< " const Eigen::Matrix<T,Eigen::Dynamic,1>& f_xnod," << endl
<< " Eigen::Matrix<T,Eigen::Dynamic,1>& dof) const" << endl
<< "{" << endl
<< " dof = f_xnod;" << endl
<< "}" << endl;
// --------------------------------------------------
// init sizes
// --------------------------------------------------
out << "template<class T>" << endl
<< "void" << endl
<< class_name.str() << "<T>::_initialize_cstor_sizes() const" << endl
<< "{" << endl;
if (have_index_parameter()) { // P0, P1
out << " basis_fem_Pk_lagrange<T>::initialize_local_first (" << endl
<< " " << degree() << "," << endl
<< " base::is_continuous()," << endl
<< " base::_ndof_on_subgeo," << endl
<< " base::_nnod_on_subgeo," << endl
<< " base::_first_idof_by_dimension," << endl
<< " base::_first_inod_by_dimension);" << endl;
} else if (name() == "bubble") {
out << " basis_fem_Pk_lagrange<T>::initialize_local_first (" << endl
<< " 0," << endl
<< " false," << endl
<< " base::_ndof_on_subgeo," << endl
<< " base::_nnod_on_subgeo," << endl
<< " base::_first_idof_by_dimension," << endl
<< " base::_first_inod_by_dimension);" << endl;
} else {
out << " fatal_macro(\"::_initialize_cstor_sizes: not yet\"); /* TODO */" << endl;
}
out << "}" << endl;
// --------------------------------------------------
// hat_node (vectorial)
// --------------------------------------------------
out << "template<class T>" << endl
<< "void" << endl
<< class_name.str() << "<T>::_initialize_data(" << endl
<< " reference_element hat_K) const" << endl
<< "{" << endl
<< " switch (hat_K.variant()) {" << endl;
for (size_type i = 0; i < reference_element::max_variant; i++) {
const basis_symbolic_nodal_on_geo& b = operator[](i);
if (b.size() == 0) continue;
out << " case reference_element::" << b.hat_K().name() << ": {" << endl
<< " return " << class_name.str() << "_" << b.hat_K().name()
<< "<T>::hat_node (_hat_node[hat_K.variant()]);" << endl
<< " }" << endl;
}
out << " default : {" << endl
<< " error_macro (\"hat_node: unsupported `\" << hat_K.name() << \"' element type\");" << endl
<< " }" << endl
<< " }" << endl
<< "}" << endl;
// --------------------------------------------------
// call to a constructor
// --------------------------------------------------
out << "// instantiation in library:" << endl
<< "template class basis_" << name() << "<Float>;" << endl
<< "} // namespace rheolef" << endl;
}
void
basis_symbolic_nodal::put_cxx_main (int argc, char**argv) const
{
if (argc <= 1 || string(argv[1]) == "-h") {
put_cxx_header (cout);
} else {
put_cxx_body (cout);
}
}
rheolef::io::out
@ out
Definition: rheostream.h:167
mkgeo_ball.expr
expr
Definition: mkgeo_ball.sh:361
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
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
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:391
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133