code_gen.cpp

Go to the documentation of this file.

#include <SyFi.h>
#include <fstream>

using namespace GiNaC; 
using namespace std; 
using namespace SyFi; 

void compute_Poisson_element_matrix(
        FE& fe, 
        Dof& dof, 
        std::map<std::pair<unsigned int,unsigned int>, GiNaC::ex>& A)
{
    std::pair<unsigned int,unsigned int> index; 

    // Insert the local degrees of freedom into the global Dof
    for (unsigned int i=0; i< fe.nbf(); i++) {
        dof.insert_dof(1,i,fe.dof(i)); 
    }

    Polygon& domain = fe.get_polygon(); 

    // The term (grad u, grad v)  
    for (unsigned int i=0; i< fe.nbf(); i++) {
        index.first = dof.glob_dof(fe.dof(i));               // fetch the global dof for Ni 
        for (unsigned int j=0; j< fe.nbf(); j++) {
            index.second = dof.glob_dof(fe.dof(j));            // fetch the global dof for Nj 
            GiNaC::ex nabla = inner(grad(fe.N(i)),grad(fe.N(j))); 
            GiNaC::ex Aij = domain.integrate(nabla);   // compute the integral  
            A[index] += Aij;                           // add to global matrix 
        }
    }
}

void code_gen2D(FE& fe){ 
    cout <<csrc; 
    for (unsigned int i=0; i< fe.nbf(); i++) {
        cout <<"double N"<<i<<"(double x, double y){"<<endl; 
        cout <<"  return "<<fe.N(i)<<";"<<endl; 
        cout <<"}"<<endl; 
    }

    for (unsigned int i=0; i< fe.nbf(); i++) {
        cout <<"double dN"<<i<<"dx(double x, double y){"<<endl; 
        cout <<"  return "<<diff(fe.N(i),x)<<";"<<endl; 
        cout <<"}"<<endl; 
    }

    for (unsigned int i=0; i< fe.nbf(); i++) {
        cout <<"double dN"<<i<<"dy(double x, double y){"<<endl; 
        cout <<"  return "<<diff(fe.N(i),y)<<";"<<endl; 
        cout <<"}"<<endl; 
    }

    cout <<"double N(int i, double x, double y){"<<endl; 
    cout <<"  switch(i) {"<<endl; 
    for (unsigned int i=0; i< fe.nbf(); i++) {
        cout <<"    case "<<i<<" :  return N"<<i<<"(x,y);"<<endl;  
    }
    cout <<"  }"<<endl; 
    cout <<"}"<<endl; 

    cout <<"double dNdx(int i, double x, double y){"<<endl; 
    cout <<"  switch(i) {"<<endl; 
    for (unsigned int i=0; i< fe.nbf(); i++) {
        cout <<"    case "<<i<<" :  return dN"<<i<<"dx(x,y);"<<endl;  
    }
    cout <<"  }"<<endl; 
    cout <<"}"<<endl; 

    cout <<"double dNdy(int i, double x, double y){"<<endl; 
    cout <<"  switch(i) {"<<endl; 
    for (unsigned int i=0; i< fe.nbf(); i++) {
        cout <<"    case "<<i<<" :  return dN"<<i<<"dy(x,y);"<<endl;  
    }
    cout <<"  }"<<endl; 
    cout <<"}"<<endl; 
}

int main() {

    initSyFi(2); 

    Triangle triangle(lst(0,0), lst(1,0), lst(0,1));   
    Lagrange fe; 
    fe.set_order(2); 
    fe.set_polygon(triangle); 
    fe.compute_basis_functions(); 


    code_gen2D(fe); 


    Dof dof; 
    std::map<std::pair<unsigned int,unsigned int>, ex> A; 
    ::compute_Poisson_element_matrix(fe, dof, A); 
    cout <<"C code format on output "<<endl; 
    cout <<csrc; 
    print(A); 


    archive ar; 
    for (unsigned int i=0; i<fe.nbf(); i++) {
      ar.archive_ex(fe.N(i), istr("N",i).c_str()); 

    }
    pair<unsigned int,unsigned int> index; 
    for (unsigned int i=0; i<fe.nbf(); i++) {
      index.first = i; 
      for (unsigned int j=0; j<fe.nbf(); j++) {
        index.second = j; 
        ar.archive_ex(A[index], istr("A",i,j).c_str()); 
      }

    }
    ofstream vfile("code_gen.gar.v"); 
    vfile << ar; vfile.close(); 
    if(!compare_archives("code_gen.gar.v", "code_gen.gar.r")) { 
      cerr << "Failure!" << endl;
      return -1;
    }

    return 0; 
}