PDEOperatorExpression.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: PDEOperatorExpression.cpp,v 1.6 2007/02/10 17:20:53 delpinux Exp $

#include <PDEOperatorExpression.hpp>
#include <PDEOperator.hpp>

#include <FunctionExpressionUnaryMinus.hpp>
#include <ScalarFunctionBase.hpp>

#include <DivMuGrad.hpp>
#include <SecondOrderOperator.hpp>

#include <FirstOrderOperator.hpp>

#include <MassOperator.hpp>

PDEOperatorExpression::
PDEOperatorExpression(const PDEOperatorExpression& e)
  : Expression(e),
    __pdeOperatorType(e.__pdeOperatorType)
{
  ;
}

PDEOperatorExpression::
PDEOperatorExpression(PDEOperatorType t)
  : Expression(Expression::pdeOperator),
    __pdeOperatorType(t)
{
  ;
}

PDEOperatorExpression::~PDEOperatorExpression()
{
  ;
}

void PDEOperatorExpression::execute()
{
  this->__execute();
  if (this->__checkBoundaryExpression()) {
    throw ErrorHandler(__FILE__, __LINE__,
                       "cannot evaluate the expression \""
                       +stringify(*this)+
                       "\" in the volume",
                       ErrorHandler::normal);
  }
}

ReferenceCounting<PDEOperator>
PDEScalarOperatorExpression::pdeOperator()
{
  return __pdeOperator;
}

PDEScalarOperatorExpression::PDEScalarOperatorExpression(const PDEScalarOperatorExpression& e)
  : PDEOperatorExpression(e),
    __unknownName(e.__unknownName),
    __pdeOperator(e.__pdeOperator),
    __pdeScalarOperatorType(e.__pdeScalarOperatorType)
{
  ;
}

PDEScalarOperatorExpression::
PDEScalarOperatorExpression(const std::string& unknownName,
                            PDEScalarOperatorType t)
  : PDEOperatorExpression(PDEOperatorExpression::scalar),
    __unknownName(unknownName),
    __pdeScalarOperatorType(t)
{
  ;
}

PDEScalarOperatorExpression::~PDEScalarOperatorExpression()
{
  ;
}


bool PDEScalarOperatorExpressionDivMuGrad::__checkBoundaryExpression() const
{
  return __mu->hasBoundaryExpression();
}

void PDEScalarOperatorExpressionDivMuGrad::__execute()
{
  __mu->execute();
  ScalarFunctionBuilder functionBuilder;
  functionBuilder.setFunction(__mu->function());
  functionBuilder.setUnaryMinus();

  __pdeOperator = new DivMuGrad(functionBuilder.getBuiltFunction());
}

bool PDEScalarOperatorExpressionOrderZero::__checkBoundaryExpression() const
{
  return __Alpha->hasBoundaryExpression();
}

void PDEScalarOperatorExpressionOrderZero::__execute()
{
  __Alpha->execute();
  __pdeOperator = new MassOperator(__Alpha->function());
}

bool PDEScalarOperatorExpressionOrderOne::__checkBoundaryExpression() const
{
  for (size_t i=0; i<3; ++i) {
    if(__nu[i] != 0) {
      if (__nu[i]->hasBoundaryExpression()) {
        return true;
      }
    }
  }
  return false;
}

void PDEScalarOperatorExpressionOrderOne::__execute()
{
  this->preexec();

  ReferenceCounting<TinyVector<3, ConstReferenceCounting<ScalarFunctionBase> > > nu
    = new TinyVector<3, ConstReferenceCounting<ScalarFunctionBase> >;

  for (size_t i=0; i<3; ++i) {
    if(__nu[i] != 0) {
      (*nu)[i] = __nu[i]->function();
    } else {
      (*nu)[i] = 0;
    }
  }

  __pdeOperator = new FirstOrderOperator(nu);
}


PDEVectorialOperatorExpression::
PDEVectorialOperatorExpression(const PDEVectorialOperatorExpression::
                               PDEScalarOperatorType& t)
  : PDEOperatorExpression(PDEOperatorExpression::vectorial),
    __pdeOperatorType(t)
{
  ;
}

PDEVectorialOperatorExpression::
PDEVectorialOperatorExpression(const PDEVectorialOperatorExpression& e)
  : PDEOperatorExpression(e),
    __pdeOperatorList(e.__pdeOperatorList),
    __pdeOperatorType(e.__pdeOperatorType)
{
  ;
}

PDEVectorialOperatorExpression::
~PDEVectorialOperatorExpression()
{
  ;
}

bool PDEVectorialOperatorExpressionOrderOne::
__checkBoundaryExpression() const
{
  for(FirstOrderSum::const_iterator i = __sum.begin();
      i != __sum.end(); ++i) {
    if ((*(*i).second).__checkBoundaryExpression()) {
      return true;
    }
  }
  for(FirstOrderSum::const_iterator i = __diff.begin();
      i != __diff.end(); ++i) {
    if ((*(*i).second).__checkBoundaryExpression()) {
      return true;
    }
  }
  return false;
}

bool PDEVectorialOperatorExpressionOrderTwo::
__checkBoundaryExpression() const
{
  typedef
    TinyMatrix<3,3,ReferenceCounting<FunctionExpression> > 
    FunctionMatrix;

  typedef
    std::map<FunctionVariable*,FunctionMatrix > 
    FunctionMatrixSet;
  FunctionMatrixSet matrixSet;
  
  typedef
    PDEVectorialOperatorExpressionOrderOne::FirstOrderSum::const_iterator
    PDEVectorialOperatorIterator;

  for (FirstOperatorList::const_iterator i = __firstOrderOp.begin();
       i != __firstOrderOp.end(); ++i) {
    const PDEVectorialOperatorExpressionOrderOne& pdeOpOrderOne
      = (*(*i).second);
    
    if (pdeOpOrderOne.__checkBoundaryExpression()) {
      return true;
    }
  }

  return false;
}

void PDEVectorialOperatorExpressionOrderTwo::__execute()
{
  typedef
    TinyMatrix<3,3,ReferenceCounting<FunctionExpression> > 
    FunctionMatrix;

  typedef
    std::map<std::string,FunctionMatrix > 
    FunctionMatrixSet;
  FunctionMatrixSet matrixSet;
  
  typedef
    PDEVectorialOperatorExpressionOrderOne::FirstOrderSum::iterator
    PDEVectorialOperatorIterator;

  for (FirstOperatorList::iterator i = __firstOrderOp.begin();
       i != __firstOrderOp.end(); ++i) {
    int I = (*i).first;
    PDEVectorialOperatorExpressionOrderOne& pdeOpOrderOne
      = (*(*i).second);
    pdeOpOrderOne.preexec();
    
    for (PDEVectorialOperatorIterator j = pdeOpOrderOne.__sum.begin();
         j != pdeOpOrderOne.__sum.end(); ++j) {
      PDEScalarOperatorExpressionOrderOne& operatorOrderOne = *j->second;

      FunctionMatrix& A = matrixSet[(*j).first];

      for (size_t J=0; J<3; ++J) {
        if (operatorOrderOne.nu(J) != 0) {
          if (A(I,J) == 0) {
            A(I,J) = operatorOrderOne.nu(J);
          } else {
            A(I,J)
              = new FunctionExpressionBinaryOperation(BinaryOperation::sum,A(I,J),operatorOrderOne.nu(J));
          }
        }
      }
    }

    for (PDEVectorialOperatorIterator j = pdeOpOrderOne.__diff.begin();
         j != pdeOpOrderOne.__diff.end(); ++j) {
      PDEScalarOperatorExpressionOrderOne& operatorOrderOne = (*(*j).second);

      FunctionMatrix& A = matrixSet[(*j).first];

      for (size_t J=0; J<3; ++J) {
        if (operatorOrderOne.nu(J) != 0) {
          if (A(I,J) == 0) {
            A(I,J) = new FunctionExpressionUnaryMinus(operatorOrderOne.nu(J));
          } else {
            A(I,J)
              = new FunctionExpressionBinaryOperation(BinaryOperation::difference, A(I,J),operatorOrderOne.nu(J));
          }
        }
      }
    }
  }
  for (FunctionMatrixSet::iterator i = matrixSet.begin();
       i != matrixSet.end(); ++i) {
    FunctionMatrix& a = (*i).second;
    ReferenceCounting<SecondOrderOperator::Matrix> A
      = new SecondOrderOperator::Matrix;
    for (size_t m=0; m<3; ++m)
      for (size_t n=0; n<3; ++n) {
        if (a(m,n) != 0) {
          ScalarFunctionBuilder functionBuilder;
          functionBuilder.setFunction(a(m,n)->function());
          functionBuilder.setUnaryMinus();
          (*A)(m,n) = functionBuilder.getBuiltFunction();
        }
      }
    SecondOrderOperator* secondOrderOp = new SecondOrderOperator(A);
    __pdeOperatorList.insert(std::make_pair(i->first,secondOrderOp));
  }
}

PDEVectorialOperatorExpressionOrderTwo::
PDEVectorialOperatorExpressionOrderTwo()
  : PDEVectorialOperatorExpression(PDEVectorialOperatorExpression::
                                   secondOrder)
{
  ;
}

PDEVectorialOperatorExpressionOrderTwo::
~PDEVectorialOperatorExpressionOrderTwo()
{
  ;
}

---