ElementComputations.cpp

Go to the documentation of this file.

// Copyright (C) 2006-2009 Kent-Andre Mardal and Simula Research Laboratory
//
// This file is part of SyFi.
//
// SyFi is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 2 of the License, or
// (at your option) any later version.
//
// SyFi is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with SyFi. If not, see <http://www.gnu.org/licenses/>.

#include "ElementComputations.h"
#include "tools.h"

using std::cout;
using std::endl;

namespace SyFi
{

        void usage(FE& fe)
        {
                for (unsigned int i=0; i< fe.nbf(); i++)
                {
                        cout <<"fe.N("<<i<<")         =   "<<fe.N(i)<<endl;
                        cout <<"grad(fe.N("<<i<<"))   =   "<<grad(fe.N(i))<<endl;
                        cout <<"fe.dof("<<i<<")       =   "<<fe.dof(i)<<endl;
                }
        }

        void usage(FE& v_fe, FE& p_fe)
        {
                for (unsigned int i=0; i< v_fe.nbf(); i++)
                {
                        cout <<"v_fe.N("<<i<<")         =   "<<v_fe.N(i)<<endl;
                        cout <<"grad(v_fe.N("<<i<<"))   =   "<<grad(v_fe.N(i))<<endl;
                        cout <<"v_fe.dof("<<i<<")       =   "<<v_fe.dof(i)<<endl;
                }
                for (unsigned int i=0; i< p_fe.nbf(); i++)
                {
                        cout <<"p_fe.N("<<i<<")=   "<<p_fe.N(i)<<endl;
                        cout <<"p_fe.dof("<<i<<")= "<<p_fe.dof(i)<<endl;
                }
        }

        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++)
                {
                                                                 // fetch the global dof for Ni
                        index.first = dof.glob_dof(fe.dof(i));
                        for (unsigned int j=0; j< fe.nbf(); j++)
                        {
                                                                 // fetch the global dof for Nj
                                index.second = dof.glob_dof(fe.dof(j));
                                                                 // compute the integrand
                                GiNaC::ex nabla = inner(grad(fe.N(i)),
                                        grad(fe.N(j)));
                                                                 // compute the integral
                                GiNaC::ex Aij = domain.integrate(nabla);
                                A[index] += Aij; // add to global matrix
                        }
                }
        }

        void compute_Stokes_element_matrix(
                FE& v_fe,
                FE& p_fe,
                Dof& dof,
                std::map<std::pair<unsigned int,unsigned int>, GiNaC::ex>& A)
        {
                std::pair<unsigned int,unsigned int> index;
                std::pair<unsigned int,unsigned int> index2;

                // FIXME: need to check that p_fe
                // contains the same domain
                Polygon& domain = v_fe.get_polygon();

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

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

                // The term -(div u, q)
                for (unsigned int i=0; i< p_fe.nbf(); i++)
                {
                                                                 // fetch the global dof for p_i
                        index.first = dof.glob_dof(p_fe.dof(i));
                        for (unsigned int j=0; j< v_fe.nbf(); j++)
                        {
                                                                 // fetch the global dof for v_j
                                index.second=dof.glob_dof(v_fe.dof(j));
                                                                 // compute the integrand
                                GiNaC::ex divV= -p_fe.N(i)*div(v_fe.N(j));
                                                                 // compute the integral
                                GiNaC::ex Aij = domain.integrate(divV);
                                A[index] += Aij; // add to global matrix

                                // Do not need to compute the term (grad(p),v), since the system is
                                // symmetric. We simply set Aji = Aij
                                index2.first = index.second;
                                index2.second = index.first;
                                A[index2] += Aij;
                        }
                }
        }

        void compute_mixed_Poisson_element_matrix(
                FE& v_fe,
                FE& p_fe,
                Dof& dof,
                std::map<std::pair<unsigned int,unsigned int>, GiNaC::ex>& A)
        {
                std::pair<unsigned int,unsigned int> index;
                std::pair<unsigned int,unsigned int> index2;

                // FIXME: need to check that p_fe
                // contains the same domain
                Polygon& domain = v_fe.get_polygon();

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

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

                // The term -(div u, q)
                for (unsigned int i=0; i< p_fe.nbf(); i++)
                {
                                                                 // fetch the global dof for p_i
                        index.first = dof.glob_dof(p_fe.dof(i));
                        for (unsigned int j=0; j< v_fe.nbf(); j++)
                        {
                                                                 // fetch the global dof for v_j
                                index.second=dof.glob_dof(v_fe.dof(j));
                                                                 // compute the integrand
                                GiNaC::ex divV= -p_fe.N(i)*div(v_fe.N(j));
                                                                 // compute the integral
                                GiNaC::ex Aij = domain.integrate(divV);
                                A[index] += Aij; // add to global matrix

                                // Do not need to compute the term (grad(p),v), since the system is
                                // symmetric we simply set Aji = Aij
                                index2.first = index.second;
                                index2.second = index.first;
                                A[index2] += Aij;
                        }
                }
        }

}