RealExpression.cpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 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: RealExpression.cpp,v 1.29 2008/09/29 23:28:46 delpinux Exp $

#include <RealExpression.hpp>
#include <BooleanExpression.hpp>

#include <FunctionExpression.hpp>
#include <Vector3Expression.hpp>

#include <MeshExpression.hpp>

#include <Variable.hpp>

#include <Structured3DMesh.hpp>
#include <MeshOfHexahedra.hpp>
#include <MeshOfTetrahedra.hpp>
#include <SurfaceMeshOfTriangles.hpp>
#include <SurfaceMeshOfQuadrangles.hpp>

#include <SpectralMesh.hpp>
#include <SpectralConformTransformation.hpp>

#include <GaussLobattoManager.hpp>

#include <Interval.hpp>
#include <Information.hpp>

#include <Hexahedron.hpp>
#include <Tetrahedron.hpp>
#include <Triangle.hpp>
#include <Quadrangle.hpp>

#include <ConformTransformation.hpp>
#include <FiniteElementTraits.hpp>

#include <ScalarFunctionBase.hpp>
#include <ScalarFunctionMaxComputer.hpp>
#include <ScalarFunctionMinComputer.hpp>

#include <VariableRepository.hpp>

#include <limits>

ReferenceCounting<RealExpression>
RealExpressionVariable::
value()
{
  return __expression->value();
}

void RealExpressionVariable::execute()
{
  RealVariable* realVariable = VariableRepository::instance().findVariable<RealVariable>(__variableName);

  __expression = realVariable->expression();
  __realValue = __expression->realValue();
}

std::istream& RealExpressionVariable::
_get(std::istream& is)
{
  if (not (is >> __realValue)) {
    throw ErrorHandler(__FILE__,__LINE__,
                       "Cannot read real value",
                       ErrorHandler::normal);
  }
  RealVariable* realVariable = VariableRepository::instance().findVariable<RealVariable>(__variableName);

  (*realVariable) = new RealExpressionValue(__realValue);

  return is;
}


RealExpressionVariable::
RealExpressionVariable(const std::string& variableName)
  : __variableName(variableName),
    __expression(0)
{
  ;
}

RealExpressionVariable::RealExpressionVariable(const RealExpressionVariable& e)
  : __realVariable(e.__realVariable),
    __expression(e.__expression)
{
  ;
}

RealExpressionVariable::~RealExpressionVariable()
{
  ;
}

/* Function Evaluate */

ReferenceCounting<RealExpression>
RealExpressionFunctionEvaluate::
value()
{
  return new RealExpressionValue(__realValue);
}

void RealExpressionFunctionEvaluate::execute()
{
  (*__realFunction).execute();
  if ((*__realFunction).hasBoundaryExpression()) {
    throw ErrorHandler(__FILE__,__LINE__,
                       "Cannot evaluate \""
                       +stringify(*__realFunction)+
                       "\": missing boundary",
                       ErrorHandler::normal);
  }
  TinyVector<3,real_t> X;
  if (__v == 0) {
    (*__x).execute();
    (*__y).execute();
    (*__z).execute();

    X[0] = (*__x).realValue();
    X[1] = (*__y).realValue();
    X[2] = (*__z).realValue();
  } else {
    (*__v).execute();
    for (size_t i=0; i<3; ++i) {
      X[i] = (*__v).value(i);
    }
  }
  const ScalarFunctionBase& f = *__realFunction->function();
  __realValue = f(X);
}

RealExpressionFunctionEvaluate::
RealExpressionFunctionEvaluate (ReferenceCounting<FunctionExpression> f,
                                ReferenceCounting<RealExpression> x,
                                ReferenceCounting<RealExpression> y,
                                ReferenceCounting<RealExpression> z)
  : __realFunction(f),
    __v(0),
    __x(x),
    __y(y),
    __z(z)
{
  ;
}


RealExpressionFunctionEvaluate::
RealExpressionFunctionEvaluate(ReferenceCounting<FunctionExpression> f,
                               ReferenceCounting<Vector3Expression> v)
  : __realFunction(f),
    __v(v),
    __x(0),
    __y(0),
    __z(0)
{
  ;
}

RealExpressionFunctionEvaluate::
RealExpressionFunctionEvaluate(const RealExpressionFunctionEvaluate& e)
  : __realFunction(e.__realFunction),
    __v(e.__v),
    __x(e.__x),
    __y(e.__y),
    __z(e.__z)
{
  ;
}

RealExpressionFunctionEvaluate::~RealExpressionFunctionEvaluate()
{
  ;
}

/* Integrate */

ReferenceCounting<RealExpression>
RealExpressionIntegrate::
value()
{
  return new RealExpressionValue(__realValue);
}

template <typename MeshType, typename QuadratureType>
real_t RealExpressionIntegrate::__integrate(const MeshType& M,
                                            const QuadratureType& Q,
                                            FunctionExpression& fe)
{
  const ScalarFunctionBase& f = *fe.function();
  real_t integral=0;
  typedef typename MeshType::CellType CellType;
  for (typename MeshType::const_iterator iCell(M);
       not(iCell.end()); ++iCell) {
    const CellType& C = *iCell;
    typename FiniteElementTraits<CellType,
                                 ScalarDiscretizationTypeBase::lagrangianFEM1>::Transformation T(C);
    typename FiniteElementTraits<CellType,
                                 ScalarDiscretizationTypeBase::lagrangianFEM1>::JacobianTransformation J(T);

    TinyVector<3, real_t> X;
    for (size_t i=0; i<QuadratureType::numberOfQuadraturePoints; ++i) {
      T.value(Q[i],X);
      integral += J.jacobianDet()*Q.weight(i)*f(X);
    }
  }
  return integral;
}


template <typename MeshType>
real_t RealExpressionIntegrate::__integrate(const MeshType& M,
                                            FunctionExpression& f)
{
  switch(this->__discretizationType) {
  case ScalarDiscretizationTypeBase::lagrangianFEM0: {
    typedef typename FiniteElementTraits<typename MeshType::CellType,
                                         ScalarDiscretizationTypeBase::lagrangianFEM0>::Type
      FiniteElementType;
    typedef typename FiniteElementType::QuadratureType QuadratureType;
    const QuadratureType& Q  = QuadratureType::instance();
    return __integrate(M,Q,f);
  }
  case ScalarDiscretizationTypeBase::lagrangianFEM1: {
    typedef typename FiniteElementTraits<typename MeshType::CellType,
                                         ScalarDiscretizationTypeBase::lagrangianFEM1>::Type
      FiniteElementType;
    typedef typename FiniteElementType::QuadratureType QuadratureType;
    const QuadratureType& Q  = QuadratureType::instance();
    return __integrate(M,Q,f);
  }
  case ScalarDiscretizationTypeBase::lagrangianFEM2: {
    typedef typename FiniteElementTraits<typename MeshType::CellType,
                                         ScalarDiscretizationTypeBase::lagrangianFEM2>::Type
      FiniteElementType;
    typedef typename FiniteElementType::QuadratureType QuadratureType;
    const QuadratureType& Q  = QuadratureType::instance();
    return __integrate(M,Q,f);
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "unexpected quadrature type",
                       ErrorHandler::unexpected);    
    return 0;
  }
  }
}

void RealExpressionIntegrate::execute()
{
  __mesh->execute();
  Information::instance().setMesh(__mesh->mesh());

  const Mesh& M = *__mesh->mesh();

  __realFunction->execute();
  FunctionExpression& f = (*__realFunction);
  
  if (M.family() ==Mesh::volume) {
    if (f.hasBoundaryExpression()) {
      throw ErrorHandler(__FILE__,__LINE__,
                         "Cannot evaluate \""
                         +stringify(f)+
                         "\": missing boundary",
                         ErrorHandler::normal);
    }
  }

  real_t& integral = __realValue;
  integral = 0;

  switch (M.type()) {
  case Mesh::cartesianHexahedraMesh: {
    const Structured3DMesh& m = static_cast<const Structured3DMesh&>(M);
    integral = __integrate(m, f);
    break;
  }
  case Mesh::hexahedraMesh: {
    const MeshOfHexahedra& m = static_cast<const MeshOfHexahedra&>(M);
    integral = __integrate(m, f);
    break;
  }
  case Mesh::tetrahedraMesh: {
    const MeshOfTetrahedra& m = static_cast<const MeshOfTetrahedra&>(M);
    integral = __integrate(m, f);
    break;
  }
  case Mesh::surfaceMeshTriangles: {
    const SurfaceMeshOfTriangles& m
      = static_cast<const SurfaceMeshOfTriangles&>(M);
    integral = __integrate(m, f);
    break;
  }
  case Mesh::surfaceMeshQuadrangles: {
    const SurfaceMeshOfQuadrangles& m
      = static_cast<const SurfaceMeshOfQuadrangles&>(M);
    integral = __integrate(m, f);
    break;
  }
  case Mesh::spectralMesh: {
    const SpectralMesh& m 
      = static_cast<const SpectralMesh&>(M);

    const ScalarFunctionBase& integratedFunction = *f.function();
    TinyVector<3, real_t> X;
    Interval intervalX(m.shape().a()[0],m.shape().b()[0]);
    Interval intervalY(m.shape().a()[1],m.shape().b()[1]);
    Interval intervalZ(m.shape().a()[2],m.shape().b()[2]);
    
    SpectralConformTransformation  transformX(intervalX);
    SpectralConformTransformation  transformY(intervalY);
    SpectralConformTransformation  transformZ(intervalZ);
    
    const GaussLobatto& gaussLobattoX = GaussLobattoManager::instance().get(m.degree(0)+1);
    const GaussLobatto& gaussLobattoY = GaussLobattoManager::instance().get(m.degree(1)+1);
    const GaussLobatto& gaussLobattoZ = GaussLobattoManager::instance().get(m.degree(2)+1);

    Vector<real_t>  nodesX(gaussLobattoX.numberOfPoints());
    for (size_t i=0; i<gaussLobattoX.numberOfPoints(); ++i) {
      nodesX[i] = transformX(gaussLobattoX(i));
    }
    Vector<real_t> nodesY(gaussLobattoY.numberOfPoints());
    for (size_t i=0; i<gaussLobattoY.numberOfPoints(); ++i) {
      nodesY[i] = transformY(gaussLobattoY(i));
    }
    Vector<real_t> nodesZ(gaussLobattoZ.numberOfPoints());
    for (size_t i=0; i<gaussLobattoZ.numberOfPoints(); ++i) {
      nodesZ[i] = transformZ(gaussLobattoZ(i));
    }
    
    for (size_t i=0; i < gaussLobattoX.numberOfPoints(); ++i) {
      X[0]= nodesX[i];
      const real_t& wi = gaussLobattoX.weight(i);
      for (size_t j=0; j < gaussLobattoY.numberOfPoints(); ++j) {
        X[1] = nodesY[j];
        const real_t& wj = gaussLobattoY.weight(j);
        for (size_t k=0; k < gaussLobattoZ.numberOfPoints(); ++k) {       
          X[2]= nodesZ[k];
          const real_t& wk = gaussLobattoZ.weight(k);
          integral += wi * wj * wk *integratedFunction(X);
        }
      }
    }
    integral =  integral *1/(transformX.inverseDeterminant()
                             * transformY.inverseDeterminant()
                             *transformZ.inverseDeterminant()
                             );   
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
                       "unexpected mesh type",
                       ErrorHandler::unexpected);
  }
  }
  
  Information::instance().unsetMesh();
}

RealExpressionIntegrate::RealExpressionIntegrate(ReferenceCounting<FunctionExpression> f,
                                                 ReferenceCounting<MeshExpression> m,
                                                 const ScalarDiscretizationTypeBase::Type& discretizationType)
  : __realFunction(f),
    __mesh(m),
    __discretizationType(discretizationType)
{
  ;
}

RealExpressionIntegrate::RealExpressionIntegrate(const RealExpressionIntegrate& e)
  : __realFunction(e.__realFunction),
    __mesh(e.__mesh),
    __discretizationType(e.__discretizationType)
{
}

RealExpressionIntegrate::~RealExpressionIntegrate()
{
  ;
}


/* Max */

ReferenceCounting<RealExpression>
RealExpressionMinMax::
value()
{
  return new RealExpressionValue(__realValue);
}

void RealExpressionMinMax::execute()
{
  __mesh->execute();
  __realFunction->execute();
  Information::instance().setMesh(__mesh->mesh());

  if (__operatorName == "max") {

    ScalarFunctionMaxComputer computer(__mesh->mesh(), 
                                       __realFunction->function());
    __realValue = computer.getValue();
  } else if (__operatorName == "min") {
    ScalarFunctionMinComputer computer(__mesh->mesh(), 
                                       __realFunction->function());
    __realValue = computer.getValue();
  } else {
    throw ErrorHandler(__FILE__,__LINE__,
                       "unknown operator: "+__operatorName,
                       ErrorHandler::unexpected);
  }

  Information::instance().unsetMesh();
}

RealExpressionMinMax::RealExpressionMinMax(const std::string& operatorName,
                                           ReferenceCounting<FunctionExpression> f,
                                           ReferenceCounting<MeshExpression> m)
  : __operatorName(operatorName),
    __realFunction(f),
    __mesh(m)
{
  ;
}

RealExpressionMinMax::RealExpressionMinMax(const RealExpressionMinMax& e)
  : __operatorName(e.__operatorName),
    __realFunction(e.__realFunction),
    __mesh(e.__mesh)
{
}

RealExpressionMinMax::~RealExpressionMinMax()
{
  ;
}


/* Boolean */

ReferenceCounting<RealExpression>
RealExpressionBoolean::
value()
{
  if((*__booleanExpression).boolValue()) {
    return new RealExpressionValue(1.);
  } else {
    return new RealExpressionValue(0.);
  }
}

void RealExpressionBoolean::execute()
{
  (*__booleanExpression).execute();
}

RealExpressionBoolean::
RealExpressionBoolean(ReferenceCounting<BooleanExpression> be)
  : __booleanExpression(be)
{
  ;
}

RealExpressionBoolean::RealExpressionBoolean(const RealExpressionBoolean& re)
  : RealExpression(re),
    __booleanExpression(re.__booleanExpression)
{
  ;
}

RealExpressionBoolean::~RealExpressionBoolean()
{
  ;
}

---