DiscretizedOperators.hpp

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: DiscretizedOperators.hpp,v 1.7 2007/11/02 20:02:21 delpinux Exp $

#ifndef DISCRETIZED_OPERATORS_HPP
#define DISCRETIZED_OPERATORS_HPP

#include <TinyMatrix.hpp>

#include <ReferenceCounting.hpp>

#include <Problem.hpp>
#include <ScalarFunctionBase.hpp>

#include <ElementaryMatrixSet.hpp>

#include <PDESystem.hpp>

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


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

#include <ErrorHandler.hpp>

#include <map>

template <typename ElementaryMatrixType>
class DiscretizedOperators
{
public:
  class FunctionAndPosition
  {
  private:
    const size_t __i;           
    const size_t __j;           
    ConstReferenceCounting<ScalarFunctionBase>
    __function;                 
  public:
    inline const size_t& i() const
    {
      return __i;
    }

    inline const size_t& j() const
    {
      return __j;
    }

    inline real_t operator()(const TinyVector<3, real_t>& x) const
    {
      return (*__function)(x);
    }

    FunctionAndPosition(const FunctionAndPosition& fap)
      : __i(fap.__i),
        __j(fap.__j),
        __function(fap.__function)
    {
      ;
    }

    FunctionAndPosition(const size_t& line,
                        const size_t& column,
                        ConstReferenceCounting<ScalarFunctionBase> function)
      : __i(line),
        __j(column),
        __function(function)
    {
      ;
    }

    ~FunctionAndPosition()
    {
      ;
    }
  };

  // ================================
  // Now we describe the class itself 
  // ================================

  typedef std::pair<const ElementaryMatrixType*,
                    FunctionAndPosition > ListPair;
  typedef std::multimap<const ElementaryMatrixType*,
                        FunctionAndPosition > List;
  typedef typename List::iterator iterator;
  typedef typename List::const_iterator const_iterator;

  typename DiscretizedOperators::iterator
  begin()
  {
    return __list.begin();
  }

  const typename DiscretizedOperators::const_iterator
  begin() const
  {
    return __list.begin();
  }

  typename DiscretizedOperators::iterator
  end()
  {
    return __list.end();
  }

  const typename DiscretizedOperators::const_iterator
  end() const
  {
    return __list.end();
  }

private:
  const ElementaryMatrixSet<ElementaryMatrixType>&
  __elementaryMatrixSet;        
  List __list;                  
  DiscretizedOperators(const DiscretizedOperators& c);

public:
  DiscretizedOperators(const ElementaryMatrixSet<ElementaryMatrixType>& e,
                       const Problem& problem)
    : __elementaryMatrixSet(e)
  {
    switch (problem.type()) {
    case Problem::pde: {
      const PDESystem& pdeSystem
        = dynamic_cast<const PDESystem&>(problem);

      for (size_t i=0; i<problem.numberOfUnknown(); ++i) {
        const PDE& pde = pdeSystem[i].pde();
        for (size_t j=0; j<pdeSystem.numberOfEquations(); ++j) {
          const PDEOperatorSum& pdeOpSum = *pde[j];
          for (size_t k=0; k<pdeOpSum.numberOfOperators(); ++k) {
            const PDEOperator& pdeOperator = *pdeOpSum[k];
            switch (pdeOperator.type()) {
            case PDEOperator::firstorderop: {
              const FirstOrderOperator& firstOrderOperator
                = dynamic_cast<const FirstOrderOperator&>(pdeOperator);
              for (size_t m=0; m<3; ++m) {
                if (firstOrderOperator.isSet(m)) {
                  __list.insert(ListPair(&__elementaryMatrixSet.firstOrderOperatorDxUV(m),
                                         FunctionAndPosition(i,j,firstOrderOperator.nu(m))));
                }
              }
              break;
            }
            case PDEOperator::divmugrad: {
              const DivMuGrad& divMuGrad
                = dynamic_cast<const DivMuGrad&>(pdeOperator);
              __list.insert(ListPair(&__elementaryMatrixSet.divMuGrad(),
                                     FunctionAndPosition(i,j,divMuGrad.mu())));
              break;
            }
            case PDEOperator::secondorderop: {
              const SecondOrderOperator& secondOrderOperator
                = dynamic_cast<const SecondOrderOperator&>(pdeOperator);
              for (size_t m=0; m<3; ++m)
                for (size_t n=0; n<3; ++n) {
                  if (secondOrderOperator.isSet(m,n)) {
                    __list.insert(ListPair(&__elementaryMatrixSet.secondOrderOperator(m,n),
                                           FunctionAndPosition(i,j,secondOrderOperator.A(m,n))));
                  }
                }
              break;
            }
            case PDEOperator::massop: {
              const MassOperator& massOperator
                = dynamic_cast<const MassOperator&>(pdeOperator);
              __list.insert(ListPair(&__elementaryMatrixSet.massOperator(),
                                     FunctionAndPosition(i,j,massOperator.alpha())));
              break;
            }
            default: {
              throw ErrorHandler(__FILE__,__LINE__,
                                 "unknown operator",
                                 ErrorHandler::unexpected);
            }
            }
          }
        }
      }
      break;
    }
    case Problem::variational: {
      const VariationalProblem& P
        = dynamic_cast<const VariationalProblem&>(problem);
      for (VariationalProblem::bilinearOperatorConst_iterator
             i = P.beginBilinearOperator();
           i != P.endBilinearOperator(); ++i) {
        switch ((*(*i)).type()) {
        case VariationalBilinearOperator::muGradUGradV: {
          const VariationalMuGradUGradVOperator& O
            = dynamic_cast<const VariationalMuGradUGradVOperator&>(*(*i));

          __list.insert(ListPair(&__elementaryMatrixSet.divMuGrad(),
                                 FunctionAndPosition(O.unknownNumber(),
                                                     O.testFunctionNumber(),
                                                     O.mu())));

          break;
        }
        case VariationalBilinearOperator::alphaDxUDxV: {
          const VariationalAlphaDxUDxVOperator& O
            = dynamic_cast<const VariationalAlphaDxUDxVOperator&>(*(*i));

          const size_t i = O.i();
          const size_t j = O.j();
        
          __list.insert(ListPair(&__elementaryMatrixSet.secondOrderOperator(i,j),
                                 FunctionAndPosition(O.unknownNumber(),
                                                     O.testFunctionNumber(),
                                                     O.alpha())));
          break;
        }
        case VariationalBilinearOperator::nuUdxV: {
          const VariationalNuUdxVOperator& O
            = dynamic_cast<const VariationalNuUdxVOperator&>(*(*i));

          const size_t n = O.i();
          __list.insert(ListPair(&__elementaryMatrixSet.firstOrderOperatorUdxV(n),
                                 FunctionAndPosition(O.unknownNumber(),
                                                     O.testFunctionNumber(),
                                                     O.nu())));
          break;
        }
        case VariationalBilinearOperator::nuDxUV: {
          const VariationalNuDxUVOperator& O
            = dynamic_cast<const VariationalNuDxUVOperator&>(*(*i));

          const size_t n = O.i();
          __list.insert(ListPair(&__elementaryMatrixSet.firstOrderOperatorDxUV(n),
                                 FunctionAndPosition(O.unknownNumber(),
                                                     O.testFunctionNumber(),
                                                     O.nu())));
          break;
        }
        case VariationalBilinearOperator::alphaUV: {
          const VariationalAlphaUVOperator& O
            = dynamic_cast<const VariationalAlphaUVOperator&>(*(*i));

          __list.insert(ListPair(&__elementaryMatrixSet.massOperator(),
                                 FunctionAndPosition(O.unknownNumber(),
                                                     O.testFunctionNumber(),
                                                     O.alpha())));

          break;
        }
        default: {
          throw ErrorHandler(__FILE__,__LINE__,
                             "unknown variational form",
                             ErrorHandler::unexpected);
        }
        }
      }
      break;
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
                         "Unknown problem type",
                         ErrorHandler::unexpected);
    }
    }
  }
};

#endif // DISCRETIZED_OPERATORS_HPP

---