nljacobian_ex.cpp

Go to the documentation of this file.

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

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

void compute_poisson_element_matrix(
        FE& fe, 
        Dof& dof, 
        std::map<std::pair<unsigned int,unsigned int>, ex>& A) 
{
    std::pair<unsigned int,unsigned int> index;
    Polygon& domain = fe.get_polygon(); 

    ex ujs = symbolic_matrix(1,fe.nbf(), "u");  
    ex u;  
    for (unsigned int k=0; k< fe.nbf(); k++)  {
        u += ujs.op(k)*fe.N(k);
    }


    for (unsigned int i=0; i< fe.nbf() ; i++) {
        index.first = i; 
        ex Fi = inner(grad(u), grad(fe.N(i)));   
        for (unsigned int j=0; j< fe.nbf() ; j++) {
            index.second = j; 
            symbol uj = ex_to<symbol>(ujs.op(j)); 
            ex nabla = Fi.diff(uj,1);  
            ex Aij = domain.integrate(nabla);   
            A[index] += Aij;  
        }
    }
}

void compute_nlconvdiff_element_matrix(
        FE& fe, 
        Dof& dof, 
        std::map<std::pair<unsigned int,unsigned int>, ex>& A) 
{
    std::pair<unsigned int,unsigned int> index;
    Polygon& domain = fe.get_polygon(); 

    // insert the local dofs into the global Dof object
    for (unsigned int i=0; i< fe.nbf() ; i++) { 
        dof.insert_dof(1,i,fe.dof(i)); 
    }

    // create the local U field: U = sum_k u_k N_k 
    ex UU = matrix(2,1,lst(0,0)); 
    ex ujs = symbolic_matrix(1,fe.nbf(), "u");  
    for (unsigned int k=0; k< fe.nbf(); k++)  {
        UU +=ujs.op(k)*fe.N(k);   // U += u_k N_k  
    }

    //Get U represented as a matrix 
    matrix U = ex_to<matrix>(UU.evalm());

    for (unsigned int i=0; i< fe.nbf() ; i++) {
        index.first = dof.glob_dof(fe.dof(i));           // fetch global dof associated with i 

        // First: the diffusion term in Fi
        ex gradU = grad(U);                              // compute the gradient of U  
        ex Fi_diffusion  = inner(gradU, grad(fe.N(i)));  // inner product of grad(U) and grad(Ni) 

        // Second: the convection term in Fi 
        ex Ut = U.transpose();                           // get the transposed of U 
        ex UgradU = (Ut*gradU).evalm();                  // compute U*grad(U)     
        ex Fi_convection = inner(UgradU, fe.N(i), true); // compute U*grad(U)*Ni    

        // add together terms for convection and diffusion 
        ex Fi = Fi_convection + Fi_diffusion;            


        // Loop over all uj and differentiate Fi with respect
        // to uj to get the Jacobian Jij 
        for (unsigned int j=0; j< fe.nbf() ; j++) {
            index.second = dof.glob_dof(fe.dof(j));    // fetch global dof associated with j 
            symbol uj = ex_to<symbol>(ujs.op(j));    // cast uj to a symbol  
            ex Jij = Fi.diff(uj,1);                    // differentiate Fi with respect to uj  
            ex Aij = domain.integrate(Jij);            // intergrate the Jacobian Jij  
            A[index] += Aij;                           // update the global matrix  
        }
    }
}




int main() {

    initSyFi(2); 

    Triangle T(lst(0,0), lst(1,0), lst(0,1), "t"); 
    int order = 2; 

    // First we compute a standard Poisson problem, i.e., 
    // The differentiation of F(u_i) with respect to u_j 
    // should give the standard Poisson problem. 
    Lagrange fe; 
    fe.set_order(order); 
    fe.set_polygon(T); 
    fe.compute_basis_functions(); 

    Dof dof1; 
    std::map<std::pair<unsigned int,unsigned int>, ex> A1;
    compute_poisson_element_matrix(fe,dof1,A1); 
    print(A1); 

    // Second we compute a nonlinear convection
    // diffusion problem. 
    VectorLagrange vfe; 
    vfe.set_order(order); 
    vfe.set_size(2); 
    vfe.set_polygon(T); 
    vfe.compute_basis_functions(); 
    usage(vfe); 

    Dof dof2; 
    std::map<std::pair<unsigned int,unsigned int>, ex> A2;
    compute_nlconvdiff_element_matrix(vfe,dof2,A2); 
    cout <<"standard output"<<endl; 
    print(A2); 
    cout <<"LaTeX output"<<endl; 
    cout <<latex; 
    print(A2); 
    cout <<"Python output"<<endl; 
    cout <<python; 
    print(A2); 
    cout <<"C output"<<endl; 
    cout <<csrc; 
    print(A2); 



    // regression test 
    
    archive ar; 
    map<std::pair<unsigned int,unsigned int>,ex>::iterator iter; 
    for (iter = A1.begin(); iter != A1.end() ; iter++) {
      ar.archive_ex((*iter).second, istr("A1_", 
                    (*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("nljacobian_ex.gar.v"); 
    vfile << ar; vfile.close(); 
    if(!compare_archives("nljacobian_ex.gar.v", "nljacobian_ex.gar.r")) { 
            cerr << "Failure!" << endl;
            return -1;
    }

    return 0; 
}