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

#include <VariationalProblemExpression.hpp>

#include <VariationalProblem.hpp>

#include <ScalarFunctionBase.hpp>

#include <VariationalOperatorMuGradUGradV.hpp>
#include <VariationalOperatorAlphaDxUDxV.hpp>
#include <VariationalOperatorAlphaUV.hpp>
#include <VariationalOperatorNuUdxV.hpp>
#include <VariationalOperatorNuDxUV.hpp>

#include <VariationalBorderOperatorAlphaUV.hpp>
#include <VariationalBorderOperatorFV.hpp>

#include <VariationalOperatorFV.hpp>
#include <VariationalOperatorFdxGV.hpp>
#include <VariationalOperatorFdxV.hpp>
#include <VariationalOperatorFgradGgradV.hpp>

#include <Information.hpp>

#include <BoundaryConditionExpressionDirichlet.hpp>

#include <BoundaryList.hpp>

class VariationalProblemExpression::AlgebraicPlus
{
public:
  static
  ConstReferenceCounting<ScalarFunctionBase>
  convert(ConstReferenceCounting<FunctionExpression> f)
  {
    return f->function();
  }
};

class VariationalProblemExpression::AlgebraicMinus
{
public:
  static
  ConstReferenceCounting<ScalarFunctionBase>
  convert(ReferenceCounting<FunctionExpression> f)
  {
    ScalarFunctionBuilder functionBuilder;
    functionBuilder.setFunction(f->function());
    functionBuilder.setUnaryMinus();
    return functionBuilder.getBuiltFunction();
  }
};

void
VariationalProblemExpression::
__splitBilinearBoundaryList(const size_t& unknownNumber,
                            const size_t& testNumber,
                            const ScalarFunctionBase* alpha,
                            const Boundary* b)
{
  const BoundaryList& boundaryList = dynamic_cast<const BoundaryList&>(*b);

  for (BoundaryList::List::const_iterator i = boundaryList.list().begin();
       i != boundaryList.list().end(); ++i) {
    if ((**i).type() == Boundary::list) {
      this->__splitBilinearBoundaryList(unknownNumber, testNumber, alpha, *i);
    } else {
      __variationalProblem->add(new VariationalBorderOperatorAlphaUV(unknownNumber,
                                                                     testNumber,
                                                                     alpha,
                                                                     *i));
    }
  }
}

void
VariationalProblemExpression::
__splitLinearBoundaryList(const size_t& testNumber,
                          const ScalarFunctionBase* f,
                          const Boundary* b)
{
  const BoundaryList& boundaryList = dynamic_cast<const BoundaryList&>(*b);

  for (BoundaryList::List::const_iterator i = boundaryList.list().begin();
       i != boundaryList.list().end(); ++i) {
    if ((*i)->type() == Boundary::list) {
      this->__splitLinearBoundaryList(testNumber, f, *i);
    } else {
      __variationalProblem->add(new VariationalBorderOperatorFV(testNumber, f, *i));
    }
  }
}

template <typename AlgebraicOperator>
void VariationalProblemExpression::
__internalSetBilinearBC(VariationalFormulaExpression::BilinearOperatorList& bilinearList)
{
  for (VariationalFormulaExpression::BilinearOperatorList::iterator
         i = bilinearList.begin(); i != bilinearList.end(); ++i) {
    switch ((*i)->operatorType()) {
    case VariationalBilinearOperatorExpression::alphaUV: {
      VariationalAlphaUVExpression& I
        = dynamic_cast<VariationalAlphaUVExpression&>(**i);

      ConstReferenceCounting<ScalarFunctionBase> alpha
        = AlgebraicOperator::convert(I.alpha());

      const size_t unknownNumber = __unknownList->number(I.unknownName());
      const size_t testNumber = __testFunctionList->number(I.testFunctionName());

      ReferenceCounting<BoundaryExpression> boundary = I.border();

      if (boundary->boundary()->type() == Boundary::list) {
        this->__splitBilinearBoundaryList(unknownNumber, testNumber, alpha, boundary->boundary());
      } else {
        __variationalProblem->add(new VariationalBorderOperatorAlphaUV(unknownNumber,
                                                                       testNumber,
                                                                       alpha,
                                                                       boundary->boundary()));
      }
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unexpected boundary operator type",
                         ErrorHandler::unexpected);

    }
    }
  }
}



template <typename AlgebraicOperator>
void VariationalProblemExpression::
__internalSetBilinear(VariationalFormulaExpression::BilinearOperatorList& bilinearList)
{
  for (VariationalFormulaExpression::BilinearOperatorList::iterator
         i = bilinearList.begin(); i != bilinearList.end(); ++i) {
    switch ((*i)->operatorType()) {
    case VariationalBilinearOperatorExpression::mugradUgradV: {
      VariationalMuGradUGradVExpression& I
        = dynamic_cast<VariationalMuGradUGradVExpression&>(**i);

      ConstReferenceCounting<ScalarFunctionBase> mu = AlgebraicOperator::convert(I.mu());
      const size_t unknownNumber = __unknownList->number(I.unknownName());
      const size_t testNumber = __testFunctionList->number(I.testFunctionName());

      __variationalProblem->add(new VariationalMuGradUGradVOperator(unknownNumber,
                                                                    I.unknownProperty(),
                                                                    testNumber,
                                                                    I.testFunctionProperty(),
                                                                    mu));
      break;
    }
    case VariationalBilinearOperatorExpression::alphaDxUDxV: {
      VariationalAlphaDxUDxVExpression& I
        = dynamic_cast<VariationalAlphaDxUDxVExpression&>(*(*i));

      ConstReferenceCounting<ScalarFunctionBase> alpha
        = AlgebraicOperator::convert(I.alpha());

      const size_t unknownNumber = __unknownList->number(I.unknownName());
      const size_t testNumber = __testFunctionList->number(I.testFunctionName());

      __variationalProblem->add(new VariationalAlphaDxUDxVOperator(unknownNumber,
                                                                   I.unknownProperty(),
                                                                   testNumber,
                                                                   I.testFunctionProperty(),
                                                                   alpha,
                                                                   I.i(),
                                                                   I.j()));
      break;
    }
    case VariationalBilinearOperatorExpression::nuUdxV: {
      VariationalNuUdxVExpression& I
        = dynamic_cast<VariationalNuUdxVExpression&>(*(*i));

      ConstReferenceCounting<ScalarFunctionBase> nu
        = AlgebraicOperator::convert(I.nu());

      const size_t unknownNumber = __unknownList->number(I.unknownName());
      const size_t testNumber = __testFunctionList->number(I.testFunctionName());

      __variationalProblem->add(new VariationalNuUdxVOperator(unknownNumber,
                                                              I.unknownProperty(),
                                                              testNumber,
                                                              I.testFunctionProperty(),
                                                              nu,
                                                              I.number()));
      break;
    }
    case VariationalBilinearOperatorExpression::nuDxUV: {
      VariationalNuDxUVExpression& I
        = dynamic_cast<VariationalNuDxUVExpression&>(*(*i));

      ConstReferenceCounting<ScalarFunctionBase> nu
        = AlgebraicOperator::convert(I.nu());

      const size_t unknownNumber = __unknownList->number(I.unknownName());
      const size_t testNumber = __testFunctionList->number(I.testFunctionName());

      __variationalProblem->add(new VariationalNuDxUVOperator(unknownNumber,
                                                              I.unknownProperty(),
                                                              testNumber,
                                                              I.testFunctionProperty(),
                                                              nu,
                                                              I.number()));
      break;
    }
    case VariationalBilinearOperatorExpression::alphaUV: {
      VariationalAlphaUVExpression& I
        = dynamic_cast<VariationalAlphaUVExpression&>(*(*i));

      ConstReferenceCounting<ScalarFunctionBase> alpha
        = AlgebraicOperator::convert(I.alpha());

      const size_t& unknownNumber = __unknownList->number(I.unknownName());
      const size_t& testNumber = __testFunctionList->number(I.testFunctionName());

      __variationalProblem->add(new VariationalAlphaUVOperator(unknownNumber,
                                                               I.unknownProperty(),
                                                               testNumber,
                                                               I.testFunctionProperty(),
                                                               alpha));
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unexpected operator type",
                         ErrorHandler::unexpected);
    }
    }
  }
}


template <typename AlgebraicOperator>
void VariationalProblemExpression::
__internalSetLinearBC(VariationalFormulaExpression::LinearOperatorList& linearList)
{
  for (VariationalFormulaExpression::LinearOperatorList::iterator
         i = linearList.begin(); i != linearList.end(); ++i) {
    switch ((*i)->operatorType()) {
    case VariationalLinearOperatorExpression::FV: {
      VariationalFVExpression& I
        = dynamic_cast<VariationalFVExpression&>(**i);

      const size_t testNumber = __testFunctionList->number(I.testFunctionName());
      ConstReferenceCounting<ScalarFunctionBase> f = AlgebraicOperator::convert(I.f());

      ReferenceCounting<BoundaryExpression> boundary = I.border();
      ASSERT(boundary != 0);

      if (boundary->boundary()->type() == Boundary::list) {
        this->__splitLinearBoundaryList(testNumber, f, boundary->boundary());
      } else {
        __variationalProblem->add(new VariationalBorderOperatorFV(testNumber,
                                                                  f,
                                                                  boundary->boundary()));
      }
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unexpected operator type",
                         ErrorHandler::unexpected);
    }
    }
  }
}


template <typename AlgebraicOperator>
void VariationalProblemExpression::
__internalSetLinear(VariationalFormulaExpression::LinearOperatorList& linearList)
{
  for (VariationalFormulaExpression::LinearOperatorList::iterator
         i = linearList.begin(); i != linearList.end(); ++i) {

    VariationalOperator::Property testFunctionProperty = VariationalOperator::normal;

    switch ((*i)->operatorType()) {
    case VariationalLinearOperatorExpression::FV: {
      VariationalFVExpression& I
        = dynamic_cast<VariationalFVExpression&>(**i);
      I.execute();

      const size_t testFunctionNumber = __testFunctionList->number(I.testFunctionName());
      ConstReferenceCounting<ScalarFunctionBase> f = AlgebraicOperator::convert(I.f());

      __variationalProblem->add(new VariationalOperatorFV(testFunctionNumber,
                                                          testFunctionProperty,
                                                          f));
      break;
    }
    case VariationalLinearOperatorExpression::FdxGV: {
      VariationalFdxGVExpression& I
        = dynamic_cast<VariationalFdxGVExpression&>(**i);

      I.execute();

      const size_t testFunctionNumber = __testFunctionList->number(I.testFunctionName());
      ConstReferenceCounting<ScalarFunctionBase> f = AlgebraicOperator::convert(I.f());
      ConstReferenceCounting<ScalarFunctionBase> g = I.g()->function();

      __variationalProblem->add(new VariationalOperatorFdxGV(testFunctionNumber,
                                                             testFunctionProperty,
                                                             f,g,
                                                             I.number()));

      break;
    }
    case VariationalLinearOperatorExpression::FdxV: {
      VariationalFdxVExpression& I
        = dynamic_cast<VariationalFdxVExpression&>(**i);

      I.execute();

      const size_t testFunctionNumber = __testFunctionList->number(I.testFunctionName());
      ConstReferenceCounting<ScalarFunctionBase> f = AlgebraicOperator::convert(I.f());

      __variationalProblem->add(new VariationalOperatorFdxV(testFunctionNumber,
                                                            testFunctionProperty,
                                                            f,I.number()));
      break;
    }
    case VariationalLinearOperatorExpression::FgradGgradV: {
      VariationalFgradGgradVExpression& I
        = dynamic_cast<VariationalFgradGgradVExpression&>(*(*i));

      I.execute();

      const size_t testFunctionNumber = __testFunctionList->number(I.testFunctionName());
      ConstReferenceCounting<ScalarFunctionBase> f = AlgebraicOperator::convert(I.f());
      ConstReferenceCounting<ScalarFunctionBase> g = I.g()->function();

      __variationalProblem->add(new VariationalOperatorFgradGgradV(testFunctionNumber,
                                                                   testFunctionProperty,
                                                                   f,g));

      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "unexpected operator type",
                         ErrorHandler::unexpected);
    }
    }
  }
}

void
VariationalProblemExpression::
__splitBoundaryList(BoundaryConditionSet& bcSet, const PDECondition* pde, const Boundary* b)
{
  const BoundaryList& boundaryList = dynamic_cast<const BoundaryList&>(*b);

  for (BoundaryList::List::const_iterator i = boundaryList.list().begin();
       i != boundaryList.list().end(); ++i) {
    if ((*i)->type() == Boundary::list) {
      this->__splitBoundaryList(bcSet, pde, *i);
    } else {
      bcSet.addBoundaryCondition(new BoundaryCondition(pde, *i));
    }
  }
}

bool VariationalProblemExpression::
hasPOVBoundary() const
{
  if (__dirichletList->hasPOVBoundary()) return true;
  if (__variationalFormula->hasPOVBoundary()) return true;
  return false;
}

bool VariationalProblemExpression::
hasPredefinedBoundary() const
{
  if (__dirichletList->hasPredefinedBoundary()) return true;
  if (__variationalFormula->hasPredefinedBoundary()) return true;
  return false;
}

void VariationalProblemExpression::
execute()
{
  __variationalFormula->execute();
  __dirichletList->execute();
  __testFunctionList->execute();

  __unknownList = Information::instance().getUnknownList();

  std::vector<ReferenceCounting<BoundaryConditionSet> > listOfBCSet;
  listOfBCSet.resize(__testFunctionList->size());

  for (size_t i=0; i<listOfBCSet.size(); ++i) {
    listOfBCSet[i] = new BoundaryConditionSet();
  }

  for (VariationalDirichletListExpression::iterator
         i = __dirichletList->begin();
       i != __dirichletList->end(); ++i) {
    size_t number = __unknownList->number((*i)->unknownName());
    BoundaryConditionSet& bcSet = (*listOfBCSet[number]);
    const BoundaryCondition& b = * (*i)->boundaryCondition();
    if (b.boundary()->type() == Boundary::list) {
      const Boundary* bc = b.boundary();
      const PDECondition* pde = b.condition();

      this->__splitBoundaryList(bcSet, pde, bc);
    } else {
      bcSet.addBoundaryCondition((*i)->boundaryCondition());
    }
  }

  // Now Boundary conditions will not change:
  std::vector<ConstReferenceCounting<BoundaryConditionSet> > listOfConstBCSet;
  listOfConstBCSet.resize(listOfBCSet.size());
  for (size_t i=0; i<listOfBCSet.size(); ++i) {
    listOfConstBCSet[i] = listOfBCSet[i];
  }

  __variationalProblem =  new VariationalProblem(listOfConstBCSet);

  __internalSetBilinear<AlgebraicPlus>(__variationalFormula->__bilinearPlus);
  __internalSetBilinear<AlgebraicMinus>(__variationalFormula->__bilinearMinus);

  __internalSetLinear<AlgebraicPlus>(__variationalFormula->__linearPlus);
  __internalSetLinear<AlgebraicMinus>(__variationalFormula->__linearMinus);

  __internalSetBilinearBC<AlgebraicPlus>(__variationalFormula->__bilinearBorderPlus);
  __internalSetBilinearBC<AlgebraicMinus>(__variationalFormula->__bilinearBorderMinus);

  __internalSetLinearBC<AlgebraicPlus>(__variationalFormula->__linearBorderPlus);
  __internalSetLinearBC<AlgebraicMinus>(__variationalFormula->__linearBorderMinus);
}


VariationalProblemExpression::
VariationalProblemExpression(ReferenceCounting<VariationalFormulaExpression> v,
                             ReferenceCounting<VariationalDirichletListExpression> d,
                             ReferenceCounting<TestFunctionExpressionList> t)
  : ProblemExpression(ProblemExpression::variationalProblem),
    __variationalFormula(v),
    __dirichletList(d),
    __testFunctionList(t)
{
  ;
}


VariationalProblemExpression::
VariationalProblemExpression(const VariationalProblemExpression& vp)
  : ProblemExpression(vp),
    __variationalFormula(vp.__variationalFormula),
    __dirichletList(vp.__dirichletList),
    __unknownList(vp.__unknownList),
    __testFunctionList(vp.__testFunctionList),
    __variationalProblem(vp.__variationalProblem)
{
  ;
}


VariationalProblemExpression::
~VariationalProblemExpression()
{
  ;
}


ReferenceCounting<VariationalProblem>
VariationalProblemExpression::variationalProblem()
{
  return __variationalProblem;
}

---