fe_ex3.cpp

Go to the documentation of this file.

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

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

int main() {

    initSyFi(2); 

    Triangle t(lst(0,0), lst(1,0), lst(0,1));  
    int order = 2; 
    ex polynom; 
    lst variables; 
    lst basis_functions; 

    ex polynom_space = pol(order, 2, "a"); 
    // the polynomial spaces on the form: 
    // first item:     a0 + a1*x + a2*y + a3*x^2 + a4*x*y ...     the polynom
    // second item:    a0, a1, a2, ...                            the coefficents 
    // third  item     1, x, y, x^2                               the basis  
    // Could also do:
    // GiNaC::ex polynom_space = bernstein(order, t, "a"); 

    polynom = polynom_space.op(0); 
    variables = ex_to<lst>(polynom_space.op(1));
    // the variables a0,a1,a2 ..

    ex Nj; 
    // The bezier ordinates (in which the basis function should be either 0 or 1)
    lst points = bezier_ordinates(t,order); 

    // Loop over all basis functions Nj and all points. 
    // Each basis function Nj is determined by a set of linear equations: 
    //   Nj(xi) = dirac(i,j) 
    // This system of equations is then solved by lsolve
    for (unsigned int j=1; j <= points.nops(); j++) {
        lst equations; 
        for (unsigned int i=1; i<= points.nops() ; i++ ) { 
            // The point xi 
            ex point = points.op(i-1); 
            // The equation Nj(x) = dirac(i,j)  
            ex eq = polynom == dirac(i,j); 
            // Substitute x = xi and y = yi and appended the equation 
            // to the list of equations
            equations.append(eq.subs(lst(x == point.op(0) , y == point.op(1))));  
        }


        // We solve the linear system 
        //      cout <<"equations "<<equations<<endl; 
        //      cout <<"variables "<<variables<<endl; 
        ex subs = lsolve(equations, variables); 
        // Substitute to get the N_j 
        Nj = polynom.subs(subs);   
        basis_functions.append(Nj); 
        cout <<"Nj "<<Nj<<endl; 
    }



    // regression test

    archive ar; 
    for (unsigned int i=0; i< basis_functions.nops(); i++) {
            ar.archive_ex(basis_functions[i], istr("N",i).c_str()); 
    }

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

    return 0; 

}