code_gen.cpp
Go to the documentation of this file.00001 #include <SyFi.h>
00002 #include <fstream>
00003
00004 using namespace GiNaC;
00005 using namespace std;
00006 using namespace SyFi;
00007
00008 void compute_Poisson_element_matrix(
00009 FE& fe,
00010 Dof& dof,
00011 std::map<std::pair<unsigned int,unsigned int>, GiNaC::ex>& A)
00012 {
00013 std::pair<unsigned int,unsigned int> index;
00014
00015
00016 for (unsigned int i=0; i< fe.nbf(); i++) {
00017 dof.insert_dof(1,i,fe.dof(i));
00018 }
00019
00020 Polygon& domain = fe.get_polygon();
00021
00022
00023 for (unsigned int i=0; i< fe.nbf(); i++) {
00024 index.first = dof.glob_dof(fe.dof(i));
00025 for (unsigned int j=0; j< fe.nbf(); j++) {
00026 index.second = dof.glob_dof(fe.dof(j));
00027 GiNaC::ex nabla = inner(grad(fe.N(i)),grad(fe.N(j)));
00028 GiNaC::ex Aij = domain.integrate(nabla);
00029 A[index] += Aij;
00030 }
00031 }
00032 }
00033
00034 void code_gen2D(FE& fe){
00035 cout <<csrc;
00036 for (unsigned int i=0; i< fe.nbf(); i++) {
00037 cout <<"double N"<<i<<"(double x, double y){"<<endl;
00038 cout <<" return "<<fe.N(i)<<";"<<endl;
00039 cout <<"}"<<endl;
00040 }
00041
00042 for (unsigned int i=0; i< fe.nbf(); i++) {
00043 cout <<"double dN"<<i<<"dx(double x, double y){"<<endl;
00044 cout <<" return "<<diff(fe.N(i),x)<<";"<<endl;
00045 cout <<"}"<<endl;
00046 }
00047
00048 for (unsigned int i=0; i< fe.nbf(); i++) {
00049 cout <<"double dN"<<i<<"dy(double x, double y){"<<endl;
00050 cout <<" return "<<diff(fe.N(i),y)<<";"<<endl;
00051 cout <<"}"<<endl;
00052 }
00053
00054 cout <<"double N(int i, double x, double y){"<<endl;
00055 cout <<" switch(i) {"<<endl;
00056 for (unsigned int i=0; i< fe.nbf(); i++) {
00057 cout <<" case "<<i<<" : return N"<<i<<"(x,y);"<<endl;
00058 }
00059 cout <<" }"<<endl;
00060 cout <<"}"<<endl;
00061
00062 cout <<"double dNdx(int i, double x, double y){"<<endl;
00063 cout <<" switch(i) {"<<endl;
00064 for (unsigned int i=0; i< fe.nbf(); i++) {
00065 cout <<" case "<<i<<" : return dN"<<i<<"dx(x,y);"<<endl;
00066 }
00067 cout <<" }"<<endl;
00068 cout <<"}"<<endl;
00069
00070 cout <<"double dNdy(int i, double x, double y){"<<endl;
00071 cout <<" switch(i) {"<<endl;
00072 for (unsigned int i=0; i< fe.nbf(); i++) {
00073 cout <<" case "<<i<<" : return dN"<<i<<"dy(x,y);"<<endl;
00074 }
00075 cout <<" }"<<endl;
00076 cout <<"}"<<endl;
00077 }
00078
00079 int main() {
00080
00081 initSyFi(2);
00082
00083 int order =1;
00084 Triangle triangle(lst(0,0), lst(1,0), lst(0,1));
00085 Lagrange fe;
00086 fe.set_order(2);
00087 fe.set_polygon(triangle);
00088 fe.compute_basis_functions();
00089
00090
00091 code_gen2D(fe);
00092
00093
00094 Dof dof;
00095 std::map<std::pair<unsigned int,unsigned int>, ex> A;
00096 ::compute_Poisson_element_matrix(fe, dof, A);
00097 cout <<"C code format on output "<<endl;
00098 cout <<csrc;
00099 print(A);
00100
00101
00102 archive ar;
00103 for (unsigned int i=0; i<fe.nbf(); i++) {
00104 ar.archive_ex(fe.N(i), istr("N",i).c_str());
00105
00106 }
00107 pair<unsigned int,unsigned int> index;
00108 for (unsigned int i=0; i<fe.nbf(); i++) {
00109 index.first = i;
00110 for (unsigned int j=0; j<fe.nbf(); j++) {
00111 index.second = j;
00112 ar.archive_ex(A[index], istr("A",i,j).c_str());
00113 }
00114
00115 }
00116 ofstream vfile("code_gen.gar.v");
00117 vfile << ar; vfile.close();
00118 if(!compare_archives("code_gen.gar.v", "code_gen.gar.r")) {
00119 cerr << "Failure!" << endl;
00120 return -1;
00121 }
00122
00123 return 0;
00124 }