BrezziDouglasMarini.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 "BrezziDouglasMarini.h"
#include <fstream>

#include "tools.h"

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

namespace SyFi
{

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

        BrezziDouglasMarini:: BrezziDouglasMarini(Polygon& p, int order, bool pointwise_) : StandardFE(p, order)
        {
                pointwise = pointwise_;
                compute_basis_functions();
        }

        void BrezziDouglasMarini:: compute_basis_functions()
        {

                if ( order < 1 )
                {
                        throw(std::logic_error("Brezzi-Douglas-Marini 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"));
                }

                GiNaC::ex nsymb = GiNaC::symbol("n");
                if (pointwise)
                {

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

                                cout <<"Can not define the Brezzi-Douglas-Marini element on a line"<<endl;

                        }

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

                                cout <<"Can not define the Brezzi-Douglas-Marini element on a tetrahedron, see Nedelec2Hdiv"<<endl;

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

                                description = istr("BrezziDouglasMarini_", order) + "_2D";

                                Triangle& triangle = (Triangle&)(*p);
                                GiNaC::lst equations;
                                GiNaC::lst variables;
                                GiNaC::lst polynom_space = bernsteinv(2,order, triangle, "b");
                                GiNaC::ex pspace = polynom_space.op(0);

                                variables = collapse(GiNaC::ex_to<GiNaC::lst>(polynom_space.op(1)));

                                GiNaC::symbol t("t");
                                GiNaC::ex dofi;
                                // dofs related to edges
                                for (int i=0; i< 3; i++)
                                {
                                        Line line = triangle.line(i);
                                        GiNaC::lst normal_vec = normal(triangle, i);
                                        GiNaC::ex Vn = inner(pspace, normal_vec);
                                        GiNaC::lst points = interior_coordinates(line, order);

                                        GiNaC::ex point;
                                        for (unsigned int j=0; j< points.nops(); j++)
                                        {
                                                point = points.op(j);
                                                dofi = Vn.subs(x == point.op(0)).subs(y == point.op(1));
                                                GiNaC::ex eq = dofi == GiNaC::numeric(0);
                                                equations.append(eq);
                                                dofs.insert(dofs.end(), GiNaC::lst(point,nsymb));
                                        }
                                }

//                              std::cout <<"no variables "<<variables.nops()<<std::endl;
//                              std::cout <<"no equations "<<equations.nops()<<std::endl;

                                // invert the matrix:
                                // GiNaC has a bit strange way to invert a matrix.
                                // It solves the system AA^{-1} = Id.
                                // It seems that this way is the only way to do
                                // properly with the solve_algo::gauss flag.
                                //
                                GiNaC::matrix b; GiNaC::matrix A;
                                matrix_from_equations(equations, variables, A, b);

                                unsigned int ncols = A.cols();
                                GiNaC::matrix vars_sq(ncols, ncols);

                                // matrix of symbols
                                for (unsigned r=0; r<ncols; ++r)
                                        for (unsigned c=0; c<ncols; ++c)
                                                vars_sq(r, c) = GiNaC::symbol();

                                GiNaC::matrix id(ncols, ncols);

                                // identity
                                const GiNaC::ex _ex1(1);
                                for (unsigned i=0; i<ncols; ++i)
                                        id(i, i) = _ex1;

                                // invert the matrix
                                GiNaC::matrix m_inv(ncols, ncols);
                                m_inv = A.solve(vars_sq, id, GiNaC::solve_algo::gauss);

                                for (unsigned int i=0; i<dofs.size(); i++)
                                {
                                        b.let_op(i) = GiNaC::numeric(1);
                                        GiNaC::ex xx = m_inv.mul(GiNaC::ex_to<GiNaC::matrix>(b));

                                        GiNaC::lst subs;
                                        for (unsigned int ii=0; ii<xx.nops(); ii++)
                                        {
                                                subs.append(variables.op(ii) == xx.op(ii));
                                        }
                                        GiNaC::ex Nj1 = pspace.op(0).subs(subs);
                                        GiNaC::ex Nj2 = pspace.op(1).subs(subs);
                                        Ns.insert(Ns.end(), GiNaC::matrix(2,1,GiNaC::lst(Nj1,Nj2)));
                                        b.let_op(i) = GiNaC::numeric(0);
                                }
                        }
                }

        }

}                                                                // namespace SyFi