ExtendedProblemStat.h

/******************************************************************************
 *
 * Extension of AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: Simon Praetorius et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/


#ifndef EXTENSIONS_EXTENDED_PROBLEM_STAT_H
#define EXTENSIONS_EXTENDED_PROBLEM_STAT_H

#include "AMDiS.h"
#include "SingularDirichletBC2.h"
#if defined NONLIN_PROBLEM
  #include "nonlin/ProblemNonLin.h"
#endif

using namespace AMDiS;

const Flag INIT_EXACT_SOLUTION = 0X10000L;
const Flag UPDATE_PERIODIC_BC  = 0X20000L;
const Flag UPDATE_DIRICHLET_BC = 0X40000L;

#if defined NONLIN_PROBLEM
typedef ProblemNonLin ProblemStat_;
#else
typedef ProblemStat ProblemStat_;
#endif
class ExtendedProblemStat : public ProblemStat_
{
public:
  
  ExtendedProblemStat(std::string nameStr, ProblemIterationInterface *problemIteration = NULL)
  :
  #if defined NONLIN_PROBLEM
  ProblemStat_(nameStr)
  #else
  ProblemStat_(nameStr, problemIteration)
  #endif
  , oldMeshChangeIdx(0)
  {
    exactSolutions.resize(nComponents);
    for (int i = 0; i < nComponents; ++i)
      exactSolutions[i] = NULL;
  }
  
  template<typename SubProblemType>
  ExtendedProblemStat(std::string nameStr, SubProblemType *subProblem)
  :
  ProblemStat_(nameStr, subProblem)
  , oldMeshChangeIdx(0)
  {
    exactSolutions.resize(nComponents);
    for (int i = 0; i < nComponents; ++i)
      exactSolutions[i] = NULL;
  }
  
  
  //////////////////////////////////////////////////////////////////////////////
  ~ExtendedProblemStat()
  {
    for (int i = 0; i < nComponents; ++i) {
      if (exactSolutions[i]) {
	delete exactSolutions[i];
	exactSolutions[i] = NULL;
      }
    }
    
    for (size_t i = 0; i < DirichletBcDataList.size(); i++)
      delete DirichletBcDataList[i];
    
  }

  
  //////////////////////////////////////////////////////////////////////////////
  void initialize(Flag initFlag,
		  ProblemStatSeq *adoptProblem = NULL,
		  Flag adoptFlag = INIT_NOTHING)
  {    
    ProblemStat_::initialize(initFlag, adoptProblem, adoptFlag);
    for (int i = 0; i < nComponents; ++i)
      exactSolutions[i] = new DOFVector< double >(getFeSpace(i), "exact_solution");
  }
  

  inline DOFVector<double>* getRhsVector(int i = 0)
  {
    return rhs->getDOFVector(i);
  }

  
  //////////////////////////////////////////////////////////////////////////////
  void setExactSolution(DOFVector<double> *dof, int component) 
  {
    TEST_EXIT(exactSolutions[component] != NULL)
      ("You have to initialize the exactSolutions vector! Set the initFlag INIT_EXACT_SOLUTION for the problem.\n");
    exactSolutions[component]->copy(*dof);
  }

  
  //////////////////////////////////////////////////////////////////////////////
  inline DOFVector<double> *getExactSolution(int i = 0) 
  {
    TEST_EXIT(exactSolutions[i] != NULL)
      ("You have to initialize the exactSolutions vector! Set the initFlag INIT_EXACT_SOLUTION for the problem.\n");
    return exactSolutions[i];
  }

  
  //////////////////////////////////////////////////////////////////////////////
  void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
			  bool asmMatrix, bool asmVector)
  {  
    ProblemStat_::buildAfterCoarsen(adaptInfo, flag, asmMatrix, asmVector);

    // update periodic data
    if (oldMeshChangeIdx != getMesh()->getChangeIndex()
      || flag.isSet(UPDATE_PERIODIC_BC)
      || (PeriodicBcDataList.size() > 0 && manualPeriodicBC.size() == 0)) {
      manualPeriodicBC.clear();
      std::vector<PeriodicBcData>::iterator it;
      for (it = PeriodicBcDataList.begin(); it != PeriodicBcDataList.end(); it++)
	it->addToList(getFeSpace(it->row), manualPeriodicBC);
    }

    // apply periodic boundary conditions
    for (size_t k = 0; k < manualPeriodicBC.size(); k++)
      applyPeriodicBC(manualPeriodicBC[k], asmMatrix, asmVector);
    
    // update dirichlet data
    if (oldMeshChangeIdx != getMesh()->getChangeIndex()
      || flag.isSet(UPDATE_DIRICHLET_BC)
      || (DirichletBcDataList.size() > 0 && singularDirichletBC.size() == 0)) {
      singularDirichletBC.clear();
      std::vector<DirichletBcDataBase*>::iterator it;
      for (it = DirichletBcDataList.begin(); it != DirichletBcDataList.end(); it++) {
	(*it)->addToList(getFeSpace((*it)->row), singularDirichletBC);
      }
    }
    
    // apply dirichlet boundary conditions
    for (size_t k = 0; k < singularDirichletBC.size(); k++)
      applyDirichletBC(singularDirichletBC[k], asmMatrix, asmVector);
    
    // update solverMatrix
    if (asmMatrix && (singularDirichletBC.size() > 0 || manualPeriodicBC.size() > 0)) {
      solverMatrix.setMatrix(*getSystemMatrix());
    }
  }

  
  //////////////////////////////////////////////////////////////////////////////
  void solve(AdaptInfo *adaptInfo,
	      bool createMatrixData = true,
	      bool storeMatrixData = false)
  {
    ProblemStat_::solve(adaptInfo, createMatrixData, storeMatrixData);

    oldMeshChangeIdx = getMesh()->getChangeIndex();
  }

  
  /// Add arbitrary boundary condition to system
  //////////////////////////////////////////////////////////////////////////////
  void addBoundaryCondition(BoundaryCondition *bc, int row, int col)
  {
    boundaryConditionSet = true;

    if (systemMatrix && (*systemMatrix)[row][col])
      (*systemMatrix)[row][col]->getBoundaryManager()->addBoundaryCondition(bc);
    if (rhs)
      rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(bc);
  }

  
  //============================================================================
  
  
  //////////////////////////////////////////////////////////////////////////////
  template<typename ValueContainer>
  void addDirichletBC(BoundaryType type, int row, int col,
		      ValueContainer *values)
  {
    BoundaryTypeContainer *bound = new BoundaryTypeContainer(type);
    DirichletBcDataList.push_back(
      new DirichletBcData<BoundaryTypeContainer, ValueContainer>(
	row, col, *bound, *values));
  }
  
  // general templatized version of addDirichletBC
  
  template<typename PosType, typename ValueContainer>
  void addDirichletBC(PosType *pos, int row, int col,
		      ValueContainer *values)
  {
    DirichletBcDataList.push_back(
      new DirichletBcData<PosType, ValueContainer>(row, col, *pos, *values));
  }
  
  /**
    * Dirichlet boundary condition at DOF-Index of vertex near to coords 'pos'.
    * Value defined by AbstractFunction, that is evaluated at 'pos'
    **/
  template<typename ValueContainer>
  void addSingularDirichletBC(WorldVector<double> &pos, int row, int col,
			      ValueContainer &values)
  {
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(&pos, row, col, &values);
  }

  
  /**
    * Dirichlet boundary condition at DOF-Index 'idx'. Value defined by double.
    **/
  template<typename ValueContainer>
  void addDirichletBC(DofIndex* idx, int row, int col,
		      ValueContainer *values)
  {
    #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    if (MPI::COMM_WORLD.Get_rank() == 0)
    #endif
    DirichletBcDataList.push_back(
      new DirichletBcData<DofIndex, ValueContainer>(row, col, *idx, *values));
  }
  
  template<typename ValueContainer>
  void addSingularDirichletBC(DegreeOfFreedom idx, int row, int col,
			      ValueContainer &values)
  {
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    DofIndex* dofIndex = new DofIndex(idx);
    addDirichletBC(dofIndex, row, col, &values);
  }

  
  /**
    * set dirichlet boundary condition on implicitly defined boundary by zero levelset
    * of signed distance function.
    * An algebraic equation is forced in the region with dist<0
    * 
    **/
  template<typename ValueContainer>
  void addImplicitDirichletBC(DOFVector<double> &signedDist,
			      int row, int col,
			      ValueContainer &values)
  {
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(&signedDist, row, col, &values);
  }
  
  
  ///
  template<typename ValueContainer>
  void addImplicitDirichletBC(AbstractFunction<double, WorldVector<double> > &signedDist,
			      int row, int col,
			      ValueContainer &values)
  {
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(&signedDist, row, col, &values);
  }

  
  ///  
  template<typename ValueContainer>
  void addDirichletBC(BoundaryType nr, AbstractFunction<bool, WorldVector<double> >* meshIndicator,
		      int row, int col,
		      ValueContainer *values)
  {
    DirichletBcDataList.push_back(
      new DirichletBcData<AbstractFunction<bool, WorldVector<double> >, ValueContainer>(
	row, col, nr, *meshIndicator, *values));
  }
    
  template<typename ValueContainer>
  void addManualDirichletBC(BoundaryType nr, AbstractFunction<bool, WorldVector<double> >* meshIndicator,
			    int row, int col,
			    ValueContainer &values)
  {
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(nr, meshIndicator, row, col, &values);
  }
  

  ///
  template<typename ValueContainer>
  void addManualDirichletBC(AbstractFunction<bool, WorldVector<double> >* meshIndicator,
			    int row, int col,
			    ValueContainer &values)
  {
    WARNING("Depreciated! Use addDirichletBC instead!!!\n");
    addDirichletBC(meshIndicator, row, col, &values);
  }
  
  
  // ===========================================================================
  
  
  //////////////////////////////////////////////////////////////////////////////
  void addManualPeriodicBC(int row,
			    BoundaryType nr, AbstractFunction<bool, WorldVector<double> >* meshIndicator,
			    AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap)
  {
    PeriodicBcDataList.push_back(PeriodicBcData(row, nr, meshIndicator, periodicMap));
  }

  
  //////////////////////////////////////////////////////////////////////////////
  void addManualPeriodicBC(int row,
			    AbstractFunction<bool, WorldVector<double> >* meshIndicator,
			    AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap)
  {
    PeriodicBcDataList.push_back(PeriodicBcData(row, meshIndicator, periodicMap));
  }

  
  /// write Systemmatrix to file in matlab-format
  void writeMatrix(std::string filename)
  {
    mtl::io::matrix_market_ostream out(filename);
    SolverMatrix<Matrix<DOFMatrix*> > solverMatrix;
    solverMatrix.setMatrix(*getSystemMatrix());
    out << solverMatrix.getMatrix();
    out.close();
  }
  
  
  /// Returns the name of the problem
  inline string getComponentName(int comp = 0)
  {
    TEST_EXIT(comp < static_cast<int>(componentNames.size()) && comp >= 0)
      ("invalid component number\n");
    return componentNames[comp];
  }
    
protected:
  // traverse matrix rows and set unity row where dirichlet condition shall be applied.
  void applyDirichletBC(size_t row_, size_t col_, 
			DegreeOfFreedom idx, double value, 
			bool asmMatrix = true, bool asmVector = true)
  {
    using namespace mtl;
    typedef DOFMatrix::base_matrix_type Matrix;

    Matrix::size_type idx_= idx;
    bool value1set = false;

    if (asmMatrix) {
      typedef mtl::traits::range_generator<tag::row, Matrix>::type c_type;
      typedef mtl::traits::range_generator<tag::nz, c_type>::type  ic_type;
      
      // if matrix-block for the identity row is NULL, create new DOFMatrix
      if (getSystemMatrix(row_, col_) == NULL) {
	TEST_EXIT(row_ != col_)("should have been created already\n");
	(*systemMatrix)[row_][col_] = new DOFMatrix(getFeSpace(row_), getFeSpace(col_), "");
	(*systemMatrix)[row_][col_]->setCoupleMatrix(true);
	(*systemMatrix)[row_][col_]->getBoundaryManager()->
	  setBoundaryConditionMap((*systemMatrix)[row_][row_]->getBoundaryManager()->
				  getBoundaryConditionMap());
      }

      size_t numCols = static_cast<size_t>(getNumComponents());
      for (size_t col = 0; col < numCols; col++) {
	if (getSystemMatrix(row_, col) == NULL)
	  continue;

	// set Dirichlet-row in matrix
	Matrix &m = getSystemMatrix(row_, col)->getBaseMatrix();

	mtl::traits::row<Matrix>::type r(m);
	mtl::traits::col<Matrix>::type c(m);
	mtl::traits::value<Matrix>::type v(m);

	c_type cursor(begin<tag::row>(m)+idx_);
	for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor) {
	  value1set = value1set || (r(*icursor) == c(*icursor) && col == col_);
	  v(*icursor, (r(*icursor) == c(*icursor) && col == col_ ? 1.0 : 0.0));
	}
      }

      if (!value1set) {
	matrix::inserter<Matrix, update_plus<double> > ins(getSystemMatrix(row_, col_)->getBaseMatrix());
	ins[idx_][idx_] << 1.0;
      }
    }
    if (asmVector)
      (*getRhsVector(row_))[idx] = value; // set Dirichlet-Value at rhs-vector
      
    (*solution->getDOFVector(col_))[idx] = value; // set Dirichlet-value for solution component
  }
  
  
  void applyDirichletBC(SingularDirichletBC &sbc, bool asmMatrix = true, bool asmVector = true)
  {
    applyDirichletBC(sbc.row, sbc.col, sbc.idx, sbc.value, asmMatrix, asmVector);
  }
  

  void applyPeriodicBC(size_t row, std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &indices, bool asmMatrix = true, bool asmVector = true)
  {
    using namespace mtl;
    typedef DOFMatrix::base_matrix_type Matrix;

    typedef mtl::traits::range_generator<tag::row, Matrix>::type c_type;
    typedef mtl::traits::range_generator<tag::nz, c_type>::type  ic_type;

    size_t numCols = static_cast<size_t>(getNumComponents());
    for (size_t col = 0; col < numCols; col++) {
      TEST_EXIT(getSystemMatrix(row, col) != NULL)
	("SystemMatrix block (%d,%d) must not be NULL! Insert a Simple_ZOT(0.0) as Workaround.\n",row,col);

      Matrix &m = getSystemMatrix(row, col)->getBaseMatrix();

      mtl::traits::row<Matrix>::type r(m);
      mtl::traits::col<Matrix>::type c(m);
      mtl::traits::value<Matrix>::type v(m);
      
      std::vector<std::vector<std::pair<DegreeOfFreedom, double> > > row_values;
      row_values.resize(indices.size());

      // erase the rows for all first indices
      for (size_t i = 0; i < indices.size(); i++) {
	if (indices[i].first == indices[i].second)
	  continue;
	
	if (asmMatrix) {
	  c_type cursor(begin<tag::row>(m)+indices[i].first);
	  for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor) {
	    row_values[i].push_back(std::make_pair(c(*icursor), v(*icursor)));
	    v(*icursor, 0.0);
	  }
	}
	if (asmVector) {
	  (*(getRhsVector(row)))[indices[i].second] += (*(getRhsVector(row)))[indices[i].first];
	  (*(getRhsVector(row)))[indices[i].first] = 0.0;
	}
      }

      // add periodic associations of first and second indices, but only in the diagonal blocks
      if (asmMatrix) {
	matrix::inserter<Matrix, update_plus<double> > ins(m);
	if (row == col) {
	  for (size_t i = 0; i < indices.size(); i++) {
	    if (indices[i].first == indices[i].second)
	      continue;
	    ins[indices[i].first][indices[i].first] << 1.0;
	    ins[indices[i].first][indices[i].second] << -1.0;
	  }
	}
	for (size_t i = 0; i < indices.size(); i++) {
	  if (indices[i].first == indices[i].second)
	    continue;
	  for (size_t j = 0; j < row_values[i].size(); j++) {
	    ins[indices[i].second][row_values[i][j].first] << row_values[i][j].second;
	  }
	}
      }
    }
  }
  

  void applyPeriodicBC(ManualPeriodicBC &pbc, bool asmMatrix = true, bool asmVector = true)
  {
    applyPeriodicBC(pbc.row, pbc.indices, asmMatrix, asmVector);
  }

  
  void applyAlgebraicEquation(size_t row,
			      std::vector<DegreeOfFreedom> &row_idx,
			      std::vector<std::vector<double> > &coefficients,
			      std::vector<double> &rhs)
  {
    using namespace mtl;
    typedef DOFMatrix::base_matrix_type Matrix;

    typedef mtl::traits::range_generator<tag::row, Matrix>::type c_type;
    typedef mtl::traits::range_generator<tag::nz, c_type>::type  ic_type;

    size_t numCols = static_cast<size_t>(getNumComponents());
    TEST_EXIT(row_idx.size() == coefficients.size() && row_idx.size() == rhs.size() && rhs.size()>0)
      ("rhs_idx, coefficients and rhs must have the same size and size >! 0\n");
    TEST_EXIT(coefficients[0].size() == numCols)
      ("You have to give coefficients for all variables\n");

    for (size_t col = 0; col < numCols; col++) {
      TEST_EXIT(getSystemMatrix(row, col) != NULL)
	("SystemMatrix block (%d,%d) must not be NULL! Insert a Simple_ZOT(0.0) as Workaround.\n",row,col);

      Matrix &m = getSystemMatrix(row, col)->getBaseMatrix();

      mtl::traits::row<Matrix>::type r(m);
      mtl::traits::col<Matrix>::type c(m);
      mtl::traits::value<Matrix>::type v(m);

      // erase the rows for all row-indices and set rhs values
      for (size_t i = 0; i < coefficients.size(); i++) {
	c_type cursor(begin<tag::row>(m)+row_idx[i]);
	for (ic_type icursor(begin<tag::nz>(cursor)), icend(end<tag::nz>(cursor)); icursor != icend; ++icursor) {
	  v(*icursor, (r(*icursor) == c(*icursor) ? coefficients[i][col] : 0.0));
	}
	(*(getRhsVector(row)))[row_idx[i]] = rhs[i];
      }
    }
  }
  

  void applyAlgebraicEquation(size_t row,
			      DegreeOfFreedom &row_idx0,
			      std::vector<double> &coefficients0,
			      double rhs0)
  {
    std::vector<DegreeOfFreedom> row_idx; row_idx.push_back(row_idx0);
    std::vector<std::vector<double> > coefficients; coefficients.push_back(coefficients0);
    std::vector<double> rhs; rhs.push_back(rhs0);

    applyAlgebraicEquation(row, row_idx, coefficients, rhs);
  }

private:
  
  int oldMeshChangeIdx;
  std::vector<DOFVector<double>*> exactSolutions;

  // data for periodic boundary conditions
  std::vector<ManualPeriodicBC> manualPeriodicBC;
  std::vector<PeriodicBcData> PeriodicBcDataList;

  // data for dirichlet boundary conditions
  std::vector<SingularDirichletBC> singularDirichletBC;
  std::vector<DirichletBcDataBase*> DirichletBcDataList;

  std::map<const FiniteElemSpace*, bool> feSpaceVisited;
  
};
#endif // EXTENSIONS_EXTENDED_PROBLEM_STAT_H