DGFunction.hpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2008 Stéphane Del Pino

//  This program 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, or (at your option)
//  any later version.

//  This program 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 this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: DGFunction.hpp,v 1.1 2008/01/31 21:38:15 delpinux Exp $

#ifndef DG_FUNCTION_HPP
#define DG_FUNCTION_HPP

#include <FiniteElementTraits.hpp>
#include <DGFunctionBase.hpp>

template <typename MeshType,
          typename FiniteElementTraits>
class DGFunction
  : public DGFunctionBase
{
private:
  typedef typename MeshType::CellType 
  CellType;                     
  typedef
  typename FiniteElementTraits::Transformation
  Transformation;               
  typedef
  typename FiniteElementTraits::Type
  FiniteElementType;            
  ConstReferenceCounting<MeshType>
  __mesh;                       
  DGFunction(const DGFunction<MeshType,FiniteElementTraits>& f);
public:
  real_t operator()(const TinyVector<3,real_t>& x) const
  {
    typename MeshType::const_iterator icell = __mesh->find(x);
    if (icell.end()) {
      return __outsideValue;
    }

    const CellType& K = *icell;
    Transformation T(K);

    TinyVector<3, real_t> xhat;
    T.invertT(x, xhat);

    real_t value = 0;
    for (size_t l=0; l<FiniteElementType::numberOfDegreesOfFreedom; ++l) {
      value += __values[(*__dofPositionsSet)(icell.number(),l)]*FiniteElementType::instance().W(l,xhat);
    }

    return value;
  }

  void operator=(const ScalarFunctionBase& f)
  {
    // during an affectation, outside values are set to 0 (the dg
    // function is "pure")
    __outsideValue = 0;

    if (f.type() == this->type()) {
      const DGFunctionBase& dgBase = dynamic_cast<const DGFunctionBase&>(f);
      if ((dgBase.discretizationType() == this->discretizationType()) and
          (dgBase.baseMesh() == this->baseMesh())) {
        // if the function is of the same kind: just copy values
        const DGFunction<MeshType,FiniteElementTraits>& dg
          = dynamic_cast<const DGFunction<MeshType,FiniteElementTraits>&>(dgBase);
        __values = dg.__values;
        return;
      }
    }
    for (size_t i=0; i<__values.size(); i++) {
      const TinyVector<3>& x = __dofPositionsSet->vertex(i);
      __values[i] = f(x);
    }
  }

  TinyVector<3,real_t>
  gradient(const TinyVector<3,real_t>& x) const
  {
    typename MeshType::const_iterator icell = __mesh->find(x);
    if (icell.end()) {
      return 0;
    }

    const CellType& K = *icell;
    Transformation T(K);

    TinyVector<3, real_t> xhat;
    T.invertT(x, xhat);

    TinyVector<3, real_t> referenceGradient = 0;
    for (size_t l=0; l<FiniteElementType::numberOfDegreesOfFreedom; ++l) {
      const real_t value = __values[(*__dofPositionsSet)(icell.number(),l)];
      referenceGradient[0] += value*FiniteElementType::instance().dxW(l,xhat);
      referenceGradient[1] += value*FiniteElementType::instance().dyW(l,xhat);
      referenceGradient[2] += value*FiniteElementType::instance().dzW(l,xhat);
    }

    TinyMatrix<3,3, real_t> J;
    {
      TinyVector<3, real_t> temp;

      T.dx(x,temp);
      for(size_t i=0; i<3; ++i) {
        J(0,i) = temp[i];
      }

      T.dy(x,temp);
      for(size_t i=0; i<3; ++i) {
        J(1,i) = temp[i];
      }

      T.dz(x,temp);
      for(size_t i=0; i<3; ++i) {
        J(2,i) = temp[i];
      }
    }
    // now we use 
    TinyVector<3, real_t> result;

    gaussPivot(J, referenceGradient, result);
    return result;
  }

  real_t dx(const TinyVector<3>& x) const
  {
    return gradient(x)[0];
  }

  real_t dy(const TinyVector<3>& x) const
  {
    return gradient(x)[1];
  }

  real_t dz(const TinyVector<3>& x) const
  {
    return gradient(x)[2];
  }

  DGFunction(ConstReferenceCounting<MeshType> mesh)
    : DGFunctionBase(mesh,
                      ScalarDiscretizationTypeBase::Type(FiniteElementTraits::ScalarDiscretizationTypeBase)),
      __mesh(mesh)
  {
    ;
  }

  DGFunction(ConstReferenceCounting<MeshType> mesh,
              const ScalarFunctionBase& f)
    : DGFunctionBase(mesh,
                      ScalarDiscretizationTypeBase::Type(FiniteElementTraits::ScalarDiscretizationTypeBase)),
      __mesh(mesh)
  {
    (*this) = f;
  }

  DGFunction(ConstReferenceCounting<MeshType> mesh,
              const real_t& d)
    : DGFunctionBase(mesh, ScalarDiscretizationTypeBase::Type(FiniteElementTraits::ScalarDiscretizationTypeBase)),
      __mesh(mesh)
  {
    __values = d;
  }

  DGFunction(ConstReferenceCounting<MeshType> mesh,
              const Vector<real_t>& values)
    : DGFunctionBase(mesh, ScalarDiscretizationTypeBase::Type(FiniteElementTraits::ScalarDiscretizationTypeBase)),
      __mesh(mesh)
  {
    ASSERT(__values.size() == values.size());
    __values = values;
  }

  ~DGFunction()
  {
    ;
  }
};

#endif // DG_FUNCTION_HPP

---