mxpoisson_ex.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); 

    ReferenceTriangle domain; 

    // First the lowest order. 
    // The lowest order is special, because we must use P0 
    RaviartThomas v_fe;
    v_fe.set_order(1);
    v_fe.set_polygon(domain);
    v_fe.compute_basis_functions();

    P0 p_fe; 
    p_fe.set_polygon(domain);
    p_fe.compute_basis_functions();

    usage(v_fe, p_fe); 

    Dof dof; 
    std::map<std::pair<unsigned int,unsigned int>, ex> A; 
    compute_mixed_Poisson_element_matrix(v_fe, p_fe, dof, A);  
    print(A); 


    // Then for instance order 3
    int order = 3; 
    RaviartThomas v_fe2;
    v_fe2.set_order(order);
    v_fe2.set_polygon(domain);
    v_fe2.compute_basis_functions();

    DiscontinuousLagrange p_fe2; 
    p_fe2.set_order(order); 
    p_fe2.set_polygon(domain);
    p_fe2.compute_basis_functions();

    usage(v_fe2, p_fe2); 

    Dof dof2; 
    std::map<std::pair<unsigned int,unsigned int>, ex> A2; 
    compute_mixed_Poisson_element_matrix(v_fe2, p_fe2, dof2, A2);  
    print(A2); 


    /* comment out test for now

    // regression test 
    
    archive ar; 
    map<std::pair<unsigned int,unsigned int>,ex>::iterator iter; 
    for (iter = A.begin(); iter != A.end() ; iter++) {
        ar.archive_ex((*iter).second, istr("A_", (*iter).first.first, (*iter).first.second).c_str()); 
    }

    for (iter = A2.begin(); iter != A2.end() ; iter++) {
        ar.archive_ex((*iter).second, istr("A2_", (*iter).first.first, (*iter).first.second).c_str()); 
    }
    ofstream vfile("mxpoisson_ex.gar.v"); 
    vfile << ar; vfile.close(); 
    if(!compare_archives("mxpoisson_ex.gar.v", "mxpoisson_ex.gar.r")) { 
            cerr << "Failure!" << endl;
            return -1;
    }

    */

    return 0; 


}