VariationalFormulaExpression.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: VariationalFormulaExpression.cpp,v 1.13 2007/12/19 00:45:45 delpinux Exp $

#include <VariationalFormulaExpression.hpp>

#include <FunctionExpressionConstant.hpp>

#include <PDESystem.hpp>

bool VariationalFormulaExpression::
hasPOVBoundary() const
{
  return __givenMultiLinearExpression->hasPOVBoundary();
}

bool VariationalFormulaExpression::
hasPredefinedBoundary() const
{
  return __givenMultiLinearExpression->hasPOVBoundary();
}

void VariationalFormulaExpression::
__getBandLForms2D(MultiLinearExpressionSum::iterator begin,
                  MultiLinearExpressionSum::iterator end,
                  BilinearOperatorList& bilinearPlus,
                  BilinearOperatorList& bilinearMinus,
                  LinearOperatorList& linearPlus,
                  LinearOperatorList& linearMinus,
                  ReferenceCounting<BoundaryExpression> borderExpression)
{
  for (MultiLinearExpressionSum::iterator j = begin;
       j != end; ++j) {
    ReferenceCounting<FunctionExpression> f
      = (*j)->getFunction();

    switch ((*j)->operatorType()) {
    case MultiLinearExpression::UV: {
      const std::string& unknownName = (*j)->getUnknownName();
      const std::string& testFunctionName = (*j)->getTestFunctionName();

      const VariationalOperator::Property unknownProperty = (*j)->getUnknownProperty();
      const VariationalOperator::Property testFunctionProperty = (*j)->getTestFunctionProperty();

      bilinearPlus.push_back(new VariationalAlphaUVExpression(f,
                                                              unknownName,
                                                              unknownProperty,
                                                              testFunctionName,
                                                              testFunctionProperty,
                                                              borderExpression));
      break;
    }
    case MultiLinearExpression::V: {
      const std::string& testFunctionName = (*j)->getTestFunctionName();
      const VariationalOperator::Property testFunctionProperty = (*j)->getTestFunctionProperty();

      linearMinus.push_back(new VariationalFVExpression(f,
                                                        testFunctionName,
                                                        testFunctionProperty,
                                                        borderExpression));
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unexpected operator type",
                         ErrorHandler::unexpected);
    }
    }
  }
}


template <typename ListType>
std::ostream& 
VariationalFormulaExpression::__putList(const ListType& listType,
                                        char delimiter,
                                        bool& first,
                                        std::ostream& os) const
{
  for (typename ListType::const_iterator i = listType.begin();
       i != listType.end(); ++i) {
    os << "\t  ";
    if (not(first) or (delimiter != '+')) {
      os << delimiter;
    } else {
      os << ' ';
    }
    first = false;
    os << "int";
    if ((*i)->border() != 0) {
      os << '[' << *(*i)->border() << ']';
    }
    os << '(' <<  **i << ")\n";
  }
  return os;
}


std::ostream& VariationalFormulaExpression::
put(std::ostream& os) const
{
  bool first = true;
  os << '\n';
  __putList(__bilinearPlus, '+' , first, os);
  __putList(__bilinearMinus, '-', first , os);

  __putList(__bilinearBorderPlus, '+', first , os);
  __putList(__bilinearBorderMinus, '-', first , os);

  os << "\t   =";

  if ( __linearPlus.size()
      +__linearMinus.size()
      +__linearBorderPlus.size()
       +__linearBorderMinus.size()==0) {
    os << " 0";
  } else {
    os << '\n';
    first = true;

    __putList(__linearPlus, '+', first , os);
    __putList(__linearMinus, '-', first , os);

    __putList(__linearBorderPlus, '+', first , os);
    __putList(__linearBorderMinus,  '-', first , os);

    os << ";\n";
  }
  return os;
}


void VariationalFormulaExpression::
__getBandLForms3D(MultiLinearExpressionSum::iterator begin,
                  MultiLinearExpressionSum::iterator end,
                  BilinearOperatorList& bilinearPlus,
                  BilinearOperatorList& bilinearMinus,
                  LinearOperatorList& linearPlus,
                  LinearOperatorList& linearMinus)
{
  for (MultiLinearExpressionSum::iterator j = begin;
       j != end; ++j) {

    ReferenceCounting<FunctionExpression> f
      = (*j)->getFunction();

    switch ((*j)->operatorType()) {
    case MultiLinearExpression::gradUgradV: {
      const std::string& unknownName = (*j)->getUnknownName();
      const std::string& testFunctionName = (*j)->getTestFunctionName();

      const VariationalOperator::Property unknownProperty = (*j)->getUnknownProperty();
      const VariationalOperator::Property testFunctionProperty = (*j)->getTestFunctionProperty();

      bilinearPlus.push_back(new VariationalMuGradUGradVExpression(f,
                                                                   unknownName,
                                                                   unknownProperty,
                                                                   testFunctionName,
                                                                   testFunctionProperty));
      break;
    }
    case MultiLinearExpression::dxUdxV: {
      const std::string& unknownName = (*j)->getUnknownName();
      const std::string& testFunctionName = (*j)->getTestFunctionName();

      IntegratedOperatorExpression::OperatorType t
        = (*j)->getTestFunctionOperator();

      size_t m = 0;
      switch (t) {
      case (IntegratedOperatorExpression::dx): {
        m=0;
        break;
      }
      case (IntegratedOperatorExpression::dy): {
        m=1;
        break;
      }
      case (IntegratedOperatorExpression::dz): {
        m=2;
        break;
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "unexpected derivative number",
                           ErrorHandler::unexpected);
      }
      }

      t = (*j)->getUnknownOperator();

      size_t n = 0;
      switch (t) {
      case (IntegratedOperatorExpression::dx): {
        n=0;
        break;
      }
      case (IntegratedOperatorExpression::dy): {
        n=1;
        break;
      }
      case (IntegratedOperatorExpression::dz): {
        n=2;
        break;
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "unexpected derivative number",
                           ErrorHandler::unexpected);
      }
      }

      const VariationalOperator::Property unknownProperty = (*j)->getUnknownProperty();
      const VariationalOperator::Property testFunctionProperty = (*j)->getTestFunctionProperty();

      bilinearPlus.push_back(new VariationalAlphaDxUDxVExpression(f,m,n,
                                                                  unknownName,
                                                                  unknownProperty,
                                                                  testFunctionName,
                                                                  testFunctionProperty));
      break;
    }
    case MultiLinearExpression::UdxV: {
      const std::string& unknownName = (*j)->getUnknownName();
      const std::string& testFunctionName = (*j)->getTestFunctionName();

      IntegratedOperatorExpression::OperatorType t
        = (*j)->getTestFunctionOperator();

      size_t n = 0;
      switch (t) {
      case (IntegratedOperatorExpression::dx): {
        n=0;
        break;
      }
      case (IntegratedOperatorExpression::dy): {
        n=1;
        break;
      }
      case (IntegratedOperatorExpression::dz): {
        n=2;
        break;
      }
      case (IntegratedOperatorExpression::gradient): {
        throw ErrorHandler(__FILE__,__LINE__,
                           "\""+stringify(**j)+"\" is not an allowed construction",
                           ErrorHandler::normal);       
        break;
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "unexpected derivative number",
                           ErrorHandler::unexpected);
      }
      }

      const VariationalOperator::Property unknownProperty = (*j)->getUnknownProperty();
      const VariationalOperator::Property testFunctionProperty = (*j)->getTestFunctionProperty();

      bilinearPlus.push_back(new VariationalNuUdxVExpression(f,n,
                                                             unknownName,
                                                             unknownProperty,
                                                             testFunctionName,
                                                             testFunctionProperty));
      break;
    }
    case MultiLinearExpression::dxUV: {
      const std::string& unknownName = (*j)->getUnknownName();
      const std::string& testFunctionName = (*j)->getTestFunctionName();

      IntegratedOperatorExpression::OperatorType t
        = (*j)->getUnknownOperator();

      size_t n = 0;
      switch (t) {
      case (IntegratedOperatorExpression::dx): {
        n=0;
        break;
      }
      case (IntegratedOperatorExpression::dy): {
        n=1;
        break;
      }
      case (IntegratedOperatorExpression::dz): {
        n=2;
        break;
      }
      case (IntegratedOperatorExpression::gradient): {
        throw ErrorHandler(__FILE__,__LINE__,
                           "\""+stringify(**j)+"\" is not an allowed construction",
                           ErrorHandler::normal);       
        break;
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "unexpected derivative number",
                           ErrorHandler::unexpected);
      }
      }
      const VariationalOperator::Property unknownProperty = (*j)->getUnknownProperty();
      const VariationalOperator::Property testFunctionProperty = (*j)->getTestFunctionProperty();

      bilinearPlus.push_back(new VariationalNuDxUVExpression(f,n,
                                                             unknownName,
                                                             unknownProperty,
                                                             testFunctionName,
                                                             testFunctionProperty));
      break;
    }
    case MultiLinearExpression::UV: {
      const std::string& unknownName = (*j)->getUnknownName();
      const std::string& testFunctionName = (*j)->getTestFunctionName();

      const VariationalOperator::Property unknownProperty = (*j)->getUnknownProperty();
      const VariationalOperator::Property testFunctionProperty = (*j)->getTestFunctionProperty();

      bilinearPlus.push_back(new VariationalAlphaUVExpression(f,
                                                              unknownName,
                                                              unknownProperty,
                                                              testFunctionName,
                                                              testFunctionProperty));
      break;
    }
    case MultiLinearExpression::gradV: {
      const std::string& testFunctionName = (*j)->getTestFunctionName();
      throw ErrorHandler(__FILE__,__LINE__,
                         "cannot discretize 'int("+stringify(*f)+"*grad("+testFunctionName+"))'",
                         ErrorHandler::normal);
      break;
    }
    case MultiLinearExpression::dxV: {
      const std::string& testFunctionName = (*j)->getTestFunctionName();

      IntegratedOperatorExpression::OperatorType t
        = (*j)->getTestFunctionOperator();

      size_t n = 0;
      switch (t) {
      case (IntegratedOperatorExpression::dx): {
        n=0;
        break;
      }
      case (IntegratedOperatorExpression::dy): {
        n=1;
        break;
      }
      case (IntegratedOperatorExpression::dz): {
        n=2;
        break;
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "unexpected derivative number",
                           ErrorHandler::unexpected);
      }
      }

      const VariationalOperator::Property testFunctionProperty = (*j)->getTestFunctionProperty();

      linearMinus.push_back(new VariationalFdxVExpression(f,
                                                          testFunctionName,
                                                          testFunctionProperty,
                                                          n));
      break;
    }
    case MultiLinearExpression::V: {
      const std::string& testFunctionName = (*j)->getTestFunctionName();

      const VariationalOperator::Property testFunctionProperty = (*j)->getTestFunctionProperty();

      linearMinus.push_back(new VariationalFVExpression(f,
                                                        testFunctionName,
                                                        testFunctionProperty));
      break;
    }
    case MultiLinearExpression::dxFV: {
      ReferenceCounting<FunctionExpression> g
        = (*j)->getFunctionExpression();

      const std::string& testFunctionName = (*j)->getTestFunctionName();

      IntegratedOperatorExpression::OperatorType t
        = (*j)->getFunctionOperator();

      size_t n = 0;
      switch (t) {
      case (IntegratedOperatorExpression::dx): {
        n=0;
        break;
      }
      case (IntegratedOperatorExpression::dy): {
        n=1;
        break;
      }
      case (IntegratedOperatorExpression::dz): {
        n=2;
        break;
      }
      default: {
        throw ErrorHandler(__FILE__,__LINE__,
                           "unexpected derivative number",
                           ErrorHandler::unexpected);
      }
      }

      const VariationalOperator::Property testFunctionProperty = (*j)->getTestFunctionProperty();

      linearMinus.push_back(new VariationalFdxGVExpression(f,g,
                                                           testFunctionName,
                                                           testFunctionProperty,
                                                           n));
      break;
    }
    case MultiLinearExpression::gradFgradV: {
      ReferenceCounting<FunctionExpression> g
        = (*j)->getFunctionExpression();
      const std::string& testFunctionName = (*j)->getTestFunctionName();

      const VariationalOperator::Property testFunctionProperty = (*j)->getTestFunctionProperty();

      linearMinus.push_back(new VariationalFgradGgradVExpression(f,g,
                                                                 testFunctionName,
                                                                 testFunctionProperty));
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unexpected operator type",
                         ErrorHandler::unexpected);
    }
    }
  }
}

void VariationalFormulaExpression::
__getBilinearAndLinearForms(MultiLinearFormSumExpression::iterator begin,
                            MultiLinearFormSumExpression::iterator end,
                            BilinearOperatorList& bilinearPlus,
                            BilinearOperatorList& bilinearMinus,
                            BilinearOperatorList& bilinearBorderPlus,
                            BilinearOperatorList& bilinearBorderMinus,
                            LinearOperatorList& linearPlus,
                            LinearOperatorList& linearMinus,
                            LinearOperatorList& linearBorderPlus,
                            LinearOperatorList& linearBorderMinus)
{
  ReferenceCounting<RealExpression> r
    = new RealExpressionValue(1);

  ReferenceCounting<FunctionExpression> one
    = new FunctionExpressionConstant(r);

  for(MultiLinearFormSumExpression::iterator i = begin;
      i != end; ++i) {

    MultiLinearExpressionSum& m = *(*i)->multiLinearExpression();

    if ((*i)->formType() == MultiLinearFormExpression::threeD) {
      __getBandLForms3D(m.beginPlus(), m.endPlus(),
                        bilinearPlus, bilinearMinus,
                        linearPlus, linearMinus);

      __getBandLForms3D(m.beginMinus(), m.endMinus(),
                        bilinearMinus, bilinearPlus,
                        linearMinus, linearPlus);
    } else {
      __getBandLForms2D(m.beginPlus(), m.endPlus(),
                        bilinearBorderPlus, bilinearBorderMinus,
                        linearBorderPlus, linearBorderMinus,
                        (*i)->boundaryExpression());

      __getBandLForms2D(m.beginMinus(), m.endMinus(),
                        bilinearBorderMinus, bilinearBorderPlus,
                        linearBorderMinus, linearBorderPlus,
                        (*i)->boundaryExpression());
    }
  }
}


void VariationalFormulaExpression::
execute()
{
  for (BilinearOperatorList::iterator i = __bilinearPlus.begin();
       i != __bilinearPlus.end(); ++i) {
    (*i)->execute();
  }

  for (BilinearOperatorList::iterator i = __bilinearMinus.begin();
       i != __bilinearMinus.end(); ++i) {
    (*i)->execute();
  }

  for (LinearOperatorList::iterator i = __linearPlus.begin();
       i != __linearPlus.end(); ++i) {
    (*i)->execute();
  }

  for (LinearOperatorList::iterator i = __linearMinus.begin();
       i != __linearMinus.end(); ++i) {
    (*i)->execute();
  }

  for (BilinearOperatorList::iterator i = __bilinearBorderPlus.begin();
       i != __bilinearBorderPlus.end(); ++i) {
    (*i)->execute();
  }

  for (BilinearOperatorList::iterator i = __bilinearBorderMinus.begin();
       i != __bilinearBorderMinus.end(); ++i) {
    (*i)->execute();
  }

  for (LinearOperatorList::iterator i = __linearBorderPlus.begin();
       i != __linearBorderPlus.end(); ++i) {
    (*i)->execute();
  }

  for (LinearOperatorList::iterator i = __linearBorderMinus.begin();
       i != __linearBorderMinus.end(); ++i) {
    (*i)->execute();
  }
}

VariationalFormulaExpression::
VariationalFormulaExpression(ReferenceCounting<MultiLinearFormSumExpression> m)
  : Expression(Expression::variationalFormula),
    __givenMultiLinearExpression(m)
{
  this->__getBilinearAndLinearForms(__givenMultiLinearExpression->plusBegin(),
                                    __givenMultiLinearExpression->plusEnd(),
                                    __bilinearPlus,
                                    __bilinearMinus,
                                    __bilinearBorderPlus,
                                    __bilinearBorderMinus,
                                    __linearPlus,
                                    __linearMinus,
                                    __linearBorderPlus,
                                    __linearBorderMinus);

  this->__getBilinearAndLinearForms(__givenMultiLinearExpression->minusBegin(),
                                    __givenMultiLinearExpression->minusEnd(),
                                    __bilinearMinus,
                                    __bilinearPlus,
                                    __bilinearBorderMinus,
                                    __bilinearBorderPlus,
                                    __linearMinus,
                                    __linearPlus,
                                    __linearBorderMinus,
                                    __linearBorderPlus);
}

VariationalFormulaExpression::
VariationalFormulaExpression(const VariationalFormulaExpression& v)
  : Expression(v),
    __givenMultiLinearExpression(v.__givenMultiLinearExpression),
    __bilinearPlus(v.__bilinearPlus),
    __bilinearMinus(v.__bilinearMinus),
    __linearPlus(v.__linearPlus),
    __linearMinus(v.__linearMinus)
{
  ;
}

VariationalFormulaExpression::
~VariationalFormulaExpression()
{
  ;
}

---