ArnoldFalkWintherWeakSym.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 "ArnoldFalkWintherWeakSym.h"
#include "Nedelec2Hdiv.h"
#include "DiscontinuousLagrange.h"
#include "P0.h"
#include "utilities.h"

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

namespace SyFi
{

        //------------ sigma element

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

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

        void ArnoldFalkWintherWeakSymSigma:: compute_basis_functions()
        {

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

                if ( order < 1 )
                {
                        throw(std::logic_error("Arnold-Falk-Winther elements must be of order 1 or higher."));
                }

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

                description = istr("ArnoldFalkWintherWeakSymSigma_", order) + "_3D";

                Nedelec2Hdiv fe;
                fe.set_order(order);
                fe.set_polygon(*p);
                fe.compute_basis_functions();

                for (int d=0; d<3; d++)
                {
                        for (unsigned int i=0; i<fe.nbf(); i++)
                        {
                                GiNaC::matrix Nmat = GiNaC::matrix(3,3);
                                Nmat(d,0) = fe.N(i).op(0);
                                Nmat(d,1) = fe.N(i).op(1);
                                Nmat(d,2) = fe.N(i).op(2);
                                Ns.insert(Ns.end(), Nmat);
                                dofs.insert(dofs.end(),GiNaC::lst(fe.dof(i).op(0), fe.dof(i).op(1), d));
                        }
                }
        }

        //------------ u element

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

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

        void ArnoldFalkWintherWeakSymU:: compute_basis_functions()
        {

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

                if ( order < 1 )
                {
                        throw(std::logic_error("Arnold-Falk-Winther elements must be of order 1 or higher."));
                }

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

                description = istr("ArnoldFalkWintherWeakSymU_", order) + "_3D";

                if ( order  > 1 )
                {
                        VectorDiscontinuousLagrange fe;
                        fe.set_order(order-1);
                        fe.set_size(3);
                        fe.set_polygon(*p);
                        fe.compute_basis_functions();

                        for (unsigned int i=0; i<fe.nbf(); i++)
                        {
                                GiNaC::lst Ni = GiNaC::lst(fe.N(i).op(0), fe.N(i).op(1), fe.N(i).op(2));
                                GiNaC::ex Nmat = GiNaC::matrix(3,1,Ni);
                                Ns.insert(Ns.end(), Nmat);
                                dofs.insert(dofs.end(),GiNaC::lst(fe.dof(i)));
                        }
                }
                else if ( order == 1 )
                {
                        VectorP0 fe;
                        fe.set_order(order-1);
                        fe.set_size(3);
                        fe.set_polygon(*p);
                        fe.compute_basis_functions();

                        for (unsigned int i=0; i<fe.nbf(); i++)
                        {
                                GiNaC::lst Ni = GiNaC::lst(fe.N(i).op(0), fe.N(i).op(1), fe.N(i).op(2));
                                GiNaC::ex Nmat = GiNaC::matrix(3,1,Ni);
                                Ns.insert(Ns.end(), Nmat);
                                GiNaC::ex d = GiNaC::lst(fe.dof(i), 0);
                                dofs.insert(dofs.end(),GiNaC::lst(d));
                        }
                }

        }

        //------------ p element

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

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

        void ArnoldFalkWintherWeakSymP:: compute_basis_functions()
        {

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

                if ( order < 1 )
                {
                        cout <<"Arnold-Falk-Winther elements must be of order 1 or higher."<<endl;
                        return;
                }

                if ( p == NULL )
                {
                        cout <<"You need to set a polygon before the basisfunctions can be computed"<<endl;
                        return;
                }

                description = istr("ArnoldFalkWintherWeakSymP_", order) + "_3D";

                if ( order > 1 )
                {

                        VectorDiscontinuousLagrange fe;
                        fe.set_order(order);
                        fe.set_size(3);
                        fe.set_polygon(*p);
                        fe.compute_basis_functions();

                        for (unsigned int i=0; i<fe.nbf(); i++)
                        {
                                GiNaC::matrix Nmat = GiNaC::matrix(3,3);
                                Nmat(1,2) = -fe.N(i).op(0);
                                Nmat(2,1) =  fe.N(i).op(0);

                                Nmat(0,2) =  fe.N(i).op(1);
                                Nmat(2,0) = -fe.N(i).op(1);

                                Nmat(0,1) = -fe.N(i).op(2);
                                Nmat(1,0) =  fe.N(i).op(2);

                                Ns.insert(Ns.end(), Nmat);
                                dofs.insert(dofs.end(),GiNaC::lst(fe.dof(i)));
                        }
                }
                else if  ( order == 1 )
                {

                        VectorP0 fe;
                        fe.set_order(order);
                        fe.set_size(3);
                        fe.set_polygon(*p);
                        fe.compute_basis_functions();

                        for (unsigned int i=0; i<fe.nbf(); i++)
                        {
                                GiNaC::matrix Nmat = GiNaC::matrix(3,3);
                                Nmat(1,2) = -fe.N(i).op(0);
                                Nmat(2,1) =  fe.N(i).op(0);

                                Nmat(0,2) =  fe.N(i).op(1);
                                Nmat(2,0) = -fe.N(i).op(1);

                                Nmat(0,1) = -fe.N(i).op(2);
                                Nmat(1,0) =  fe.N(i).op(2);

                                Ns.insert(Ns.end(), Nmat);
                                GiNaC::ex d = GiNaC::lst(fe.dof(i), 0);
                                dofs.insert(dofs.end(),GiNaC::lst(d));
                        }
                }

        }

}                                                                // namespace SyFi