arnoldfalkwinther.cpp

Go to the documentation of this file.

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



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



int main() {


    archive ar; 

    initSyFi(3); 
    ReferenceTetrahedron tetrahedron; 
    int order = 2; 

    ArnoldFalkWintherWeakSymSigma sigma_fe; 
    sigma_fe.set_order(order); 
    sigma_fe.set_polygon(tetrahedron); 
    sigma_fe.compute_basis_functions(); 
    for (unsigned int i=0; i<sigma_fe.nbf(); i++) {
        cout <<"sigma_fe.N("<<i<<")="<<sigma_fe.N(i)<<endl; ; 
        cout <<"div(sigma_fe.N("<<i<<"))="<<div(sigma_fe.N(i))<<endl; ; 
        cout <<"sigma_fe.dof("<<i<<"))="<<sigma_fe.dof(i)<<endl; ; 

        ar.archive_ex(sigma_fe.N(i), istr("sN", i).c_str()); 

    }

    map<pair<int,int>,ex> A; 
    pair<int,int> index; 
    for (unsigned int i=0; i<sigma_fe.nbf(); i++) {
        index.first = i; 
        for (unsigned int j=0; j<sigma_fe.nbf(); j++) {
            index.second = j; 
            ex integrand = inner(sigma_fe.N(i), sigma_fe.N(j)); 
            A[index] = tetrahedron.integrate(integrand); 
            cout <<"A["<<i<<","<<j<<"]="<<A[index]<<endl; 

            ar.archive_ex(A[index], istr("A", i,j).c_str()); 
        }
    }




    ArnoldFalkWintherWeakSymU u_fe; 
    u_fe.set_order(order); 
    u_fe.set_polygon(tetrahedron); 
    u_fe.compute_basis_functions(); 
    for (unsigned int i=0; i<u_fe.nbf(); i++) {
        cout <<"u_fe.N("<<i<<")="<<u_fe.N(i)<<endl; ; 
        ar.archive_ex(u_fe.N(i), istr("uN", i).c_str()); 
    }

    map<pair<int,int>,ex> B; 
    for (unsigned int i=0; i<sigma_fe.nbf(); i++) {
        index.first = i; 
        for (unsigned int j=0; j<u_fe.nbf(); j++) {
            index.second = j; 
            ex integrand = inner(div(sigma_fe.N(i)), u_fe.N(j)); 
            B[index] = tetrahedron.integrate(integrand); 
            cout <<"B["<<i<<","<<j<<"]="<<B[index]<<endl; 
        }
    }




    ArnoldFalkWintherWeakSymP p_fe; 
    p_fe.set_order(order); 
    p_fe.set_polygon(tetrahedron); 
    p_fe.compute_basis_functions(); 
    for (unsigned int i=0; i<p_fe.nbf(); i++) {
        cout <<"p_fe.N("<<i<<")="<<p_fe.N(i)<<endl; ; 
        ar.archive_ex(p_fe.N(i), istr("pN", i).c_str()); 
    }


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