SpaceTimeElement.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 "symbol_factory.h"
#include "diff_tools.h"
#include "Lagrange.h"
#include "P0.h"
#include "SpaceTimeElement.h"

namespace SyFi
{

        SpaceTimeDomain::SpaceTimeDomain(Line& time_line_, Polygon& polygon_)
        {
                time_line = time_line_.copy();
                polygon = polygon_.copy();
        }

        SpaceTimeDomain* SpaceTimeDomain::copy() const
        {
                return new SpaceTimeDomain(*this);
        }

        SpaceTimeDomain::SpaceTimeDomain(const SpaceTimeDomain& domain)
        {

                if (time_line)
                {
                        delete time_line;
                }
                if (polygon)
                {
                        delete polygon;
                }

                time_line = domain.get_time_domain().copy();
                polygon   = domain.get_space_domain().copy();
        }

        const std::string SpaceTimeDomain::str() const
        {
                return "Time" + polygon->str();
        }

        unsigned int SpaceTimeDomain:: no_space_dim() const
        {
                return polygon->no_space_dim() +1;
        }

        Line SpaceTimeDomain ::line(unsigned int i) const
        {
                //FIXME
                // Could use the convention that the time line is the first line, the
                // next lines are the lines in the polygon
                return Line();
        }

        GiNaC::ex SpaceTimeDomain:: repr(Repr_format format)  const
        {
                return GiNaC::lst(time_line->repr(t, format), polygon->repr(format));
        }

        GiNaC::ex SpaceTimeDomain::integrate(GiNaC::ex f,  Repr_format format)
        {
                GiNaC::ex intf;

                // integrate in space
                intf = polygon->integrate(f, format);

                // integrate in time ((x,y,z) are now integrated away)
                intf = intf.subs( t == x );
                intf = time_line->integrate(intf , format);

                return intf;

        }

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

        SpaceTimeElement:: SpaceTimeElement(Line* time_line_, unsigned int order_, StandardFE* fe_)
        {
                time_line = time_line_;
                order = order_;
                fe = fe_;
                compute_basis_functions();
        }

        void SpaceTimeElement:: set_time_domain(Line* line_)
        {
                time_line = line_;
        }

        void SpaceTimeElement:: set_order_in_time(unsigned int order_)
        {
                order = order_;
        }

        void SpaceTimeElement:: set_spatial_element(StandardFE* fe_)
        {
                fe = fe_;
        }

        void SpaceTimeElement:: compute_basis_functions()
        {

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

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

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

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

                StandardFE* time_element;
                if ( order == 0)
                {
                        time_element = new P0(*time_line);
                }
                else
                {
                        time_element = new Lagrange(*time_line, order);
                }

                for (unsigned int j = 0; j < fe->nbf(); j++)
                {
                        GiNaC::ex Nj = fe->N(j);
                        for (unsigned int i = 0; i < (*time_element).nbf(); i++)
                        {
                                GiNaC::ex Ni = (*time_element).N(i);
                                Ni = Ni.subs(x == t);
                                GiNaC::ex N = Nj*Ni;
                                Ns.insert(Ns.end(), N);
                                dofs.insert(dofs.end(), GiNaC::lst((*time_element).dof(i), fe->dof(j)));
                        }
                }

                description = time_element->str() + "_" +  fe->str();
                delete time_element;
        }

}                                                                // namespace SyFi