CrouzeixRaviart.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 "CrouzeixRaviart.h"

#include "tools.h"

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

namespace SyFi
{

        CrouzeixRaviart:: CrouzeixRaviart() : StandardFE()
        {
                order = 1;
        }

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

        void CrouzeixRaviart:: 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 (order != 1)
                {
                        throw(std::logic_error("Only Crouziex-Raviart elements of order 1 is possible"));
                }

                // see e.g. Brezzi and Fortin book page 116 for the definition

                if ( p->str().find("ReferenceLine") != string::npos )
                {
                        cout <<"Can not define the Raviart-Thomas element on a line"<<endl;
                }
                else if ( p->str().find("Triangle") != string::npos )
                {

                        description = "CrouzeixRaviart_2D";

                        Triangle& triangle = (Triangle&)(*p);

                        // create the polynomial space
                        GiNaC::ex polynom_space = bernstein(1, triangle, "a");
                        GiNaC::ex polynom = polynom_space.op(0);
                        GiNaC::lst variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
                        GiNaC::ex basis = polynom_space.op(2);

                        // create the dofs
                        GiNaC::symbol t("t");
                        for (int j=0; j< 3; j++)
                        {

                                // solve the linear system to compute
                                // each of the basis functions
                                GiNaC::lst equations;
                                for (int i=0; i< 3; i++)
                                {

                                        Line line = triangle.line(i);
                                        //                          GiNaC::ex dofi = line.integrate(polynom);
                                        GiNaC::lst midpoint = GiNaC::lst(
                                                (line.vertex(0).op(0) +  line.vertex(1).op(0))/2,
                                                (line.vertex(0).op(1) +  line.vertex(1).op(1))/2);
                                        dofs.insert(dofs.end(), midpoint);

                                        //                          GiNaC::ex dofi = polynom.subs( x == midpoint.op(0)).subs( y == midpoint.op(1));
                                        GiNaC::ex dofi = line.integrate(polynom);
                                        equations.append( dofi == dirac(i,j));

                                        if (j == 1)
                                        {
                                                //                                  GiNaC::lst d = GiNaC::lst(line.vertex(0) ,  line.vertex(1));
                                                dofs.insert(dofs.end(), midpoint);
                                        }

                                }
                                GiNaC::ex sub = lsolve(equations, variables);
                                GiNaC::ex Ni = polynom.subs(sub);
                                Ns.insert(Ns.end(),Ni);
                        }

                }
                else if ( p->str().find("Tetrahedron") != string::npos )
                {

                        description = "CrouzeixRaviart_3D";

                        Tetrahedron& tetrahedron = (Tetrahedron&)(*p);
                        GiNaC::ex polynom_space = bernstein(1, tetrahedron, "a");
                        GiNaC::ex polynom = polynom_space.op(0);
                        GiNaC::lst variables = GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1));
                        GiNaC::ex basis = polynom_space.op(2);

                        GiNaC::ex bernstein_pol;

                        GiNaC::symbol t("t");
                        // dofs related to edges
                        for (int j=0; j< 4; j++)
                        {

                                GiNaC::lst equations;
                                for (int i=0; i< 4; i++)
                                {
                                        Triangle triangle = tetrahedron.triangle(i);
                                        GiNaC::lst midpoint = GiNaC::lst(
                                                (triangle.vertex(0).op(0) + triangle.vertex(1).op(0) + triangle.vertex(2).op(0))/3,
                                                (triangle.vertex(0).op(1) + triangle.vertex(1).op(1) + triangle.vertex(2).op(1))/3,
                                                (triangle.vertex(0).op(2) + triangle.vertex(1).op(2) + triangle.vertex(2).op(2))/3
                                                );

                                        GiNaC::ex dofi = polynom.subs(x == midpoint.op(0)).subs(y == midpoint.op(1)).subs(z == midpoint.op(2));
                                        equations.append( dofi == dirac(i,j));

                                        if ( j == 1 )
                                        {
                                                //                                  GiNaC::lst d = GiNaC::lst(triangle.vertex(0), triangle.vertex(1), triangle.vertex(2));
                                                dofs.insert(dofs.end(), midpoint);
                                        }
                                }
                                GiNaC::ex sub = lsolve(equations, variables);
                                GiNaC::ex Ni = polynom.subs(sub);
                                Ns.insert(Ns.end(),Ni);

                        }
                }
        }

        // ------------VectorCrouzeixRaviart ---

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

        VectorCrouzeixRaviart:: VectorCrouzeixRaviart() : StandardFE()
        {
                description = "VectorCrouzeixRaviart";
                order = 1;
        }

        void VectorCrouzeixRaviart:: compute_basis_functions()
        {

                if (order != 1)
                {
                        throw(std::logic_error("Only Crouziex-Raviart elements of order 1 is possible"));
                }

                CrouzeixRaviart fe;
                fe.set_polygon(*p);
                fe.compute_basis_functions();

                description = "Vector" + fe.str();

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

                for (unsigned int s=0; s< size ; s++)
                {
                        for (unsigned int i=0; i< fe.nbf() ; i++)
                        {
                                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);
                        }
                }
        }

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

}                                                                // namespace SyFi