P0.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 "P0.h"
#include "tools.h"

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

namespace SyFi
{

        P0::P0() : StandardFE()
        {
                description = "P0";
        }

        P0:: P0(Polygon& p, unsigned int order) : StandardFE (p,order)
        {
                compute_basis_functions();
        }

        void P0:: compute_basis_functions()
        {

                // remove previously computed basis functions and dofs
                Ns.clear();
                dofs.clear();

                if ( p == NULL )
                {
                        throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
                }

                // insert basis function
                Ns.insert(Ns.end(), GiNaC::numeric(1));

                GiNaC::lst midpoint = GiNaC::lst();
                // create and insert dof
                                                                 // p is a lst
                if (GiNaC::is_a<GiNaC::lst>(p->vertex(0)))
                {
                        for (unsigned int d=0; d< p->vertex(1).nops(); d++)
                        {
                                midpoint.append(GiNaC::numeric(0));
                        }
                        for (unsigned int i=0; i< p->no_vertices(); i++)
                        {
                                int nops;
                                nops = p->vertex(i).nops();
                                for (int d=0; d< nops; d++)
                                {
                                        midpoint.let_op(d)  += p->vertex(i).op(d);
                                }
                        }
                        for (unsigned int d=0; d< p->vertex(1).nops(); d++)
                        {
                                midpoint.let_op(d)  = midpoint.op(d)/p->no_vertices();
                        }
                }
                else
                {
                        midpoint.append(GiNaC::numeric(0));
                        for (unsigned int i=0; i< p->no_vertices(); i++)
                        {
                                midpoint.let_op(0)  += p->vertex(i);
                        }
                        midpoint.let_op(0)  = midpoint.op(0)/p->no_vertices();
                }

                dofs.insert(dofs.end(), midpoint);

                description = istr("P0_"  , midpoint.nops()) + "D";
        }

        // ------------VectorP0 ---

        VectorP0::VectorP0() : StandardFE()
        {
                description = "VectorP0";
        }

        VectorP0::VectorP0(Polygon& p, unsigned int order, unsigned int size_) : StandardFE(p, order)
        {
                size = size_ < 0 ? nsd: size_;
                compute_basis_functions();
        }

        void VectorP0:: compute_basis_functions()
        {

                // remove previously computed basis functions and dofs
                Ns.clear();
                dofs.clear();

                if ( p == NULL )
                {
                        throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
                }

                if ( size == 0)
                {
                        throw(std::logic_error("You need to set the size of the vector before the basisfunctions can be computed"));
                }

                P0 fe(*p);
                GiNaC::lst zero_list;
                for (unsigned int s=1; s<= size ; s++)
                {
                        zero_list.append(0);
                }

                for (unsigned int i=0; i< fe.nbf() ; i++)
                {
                        for (unsigned int s=0; s< size ; s++)
                        {
                                GiNaC::lst Nis = zero_list;
                                Nis.let_op(s) = fe.N(i);
                                GiNaC::ex Nmat = GiNaC::matrix(size,1,Nis);
                                Ns.insert(Ns.end(), Nmat);

                                GiNaC::lst dof = GiNaC::lst(fe.dof(i), s) ;
                                dofs.insert(dofs.end(), dof);
                        }
                }
                description = "Vector" + fe.str();
        }

        void VectorP0:: set_size(unsigned int size_)
        {
                size = size_;
        }

        // ------------TensorP0 ---

        TensorP0::TensorP0() : StandardFE()
        {
                description = "TensorP0";
        }

        TensorP0::TensorP0(Polygon& p, unsigned int order, unsigned int size_) : StandardFE(p, order)
        {
                size = size_ < 0 ? nsd: size_;
                compute_basis_functions();
        }

        void TensorP0:: compute_basis_functions()
        {

                // remove previously computed basis functions and dofs
                Ns.clear();
                dofs.clear();

                if ( p == NULL )
                {
                        throw(std::logic_error("You need to set a polygon before the basisfunctions can be computed"));
                }

                if ( size == 0)
                {
                        throw(std::logic_error("You need to set the size of the vector before the basisfunctions can be computed"));
                }

                P0 fe(*p);
                GiNaC::lst zero_list;
                for (unsigned int s=1; s<= size*size ; s++)
                {
                        zero_list.append(0);
                }

                for (unsigned int i=0; i< fe.nbf() ; i++)
                {
                        for (unsigned int r=0; r< size ; r++)
                        {
                                for (unsigned int s=0; s< size ; s++)
                                {
                                        GiNaC::lst Nis = zero_list;
                                        Nis.let_op((size)*r + s) = fe.N(i);
                                        GiNaC::ex Nmat = GiNaC::matrix(size,size,Nis);
                                        Ns.insert(Ns.end(), Nmat);

                                        GiNaC::lst dof = GiNaC::lst(fe.dof(i), r, s) ;
                                        dofs.insert(dofs.end(), dof);
                                }
                        }
                }
                description = "Tensor" + fe.str();
        }

        void TensorP0:: set_size(unsigned int size_)
        {
                size = size_;
        }

}