PeriodicBC.cc

/******************************************************************************
 *
 * 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 Vey, Thomas Witkowski, Andreas Naumann, 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.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 *
 ******************************************************************************/


#include "PeriodicBC.h"
#include "BasisFunction.h"
#include "DOFVector.h"
#include "LeafData.h"
#include "DOFMatrix.h"
#include "Traverse.h"
#include "Boundary.h"
#include "VertexVector.h"

namespace AMDiS {

  std::vector<PeriodicDOFMapping*> PeriodicDOFMapping::mappings;


  PeriodicDOFMapping*
  PeriodicDOFMapping::providePeriodicDOFMapping(const BasisFunction *fcts)
  {
    for (std::vector<PeriodicDOFMapping*>::iterator it = mappings.begin();
	 it != mappings.end(); ++it)
      if ((*it)->basFcts == fcts)
	return *it;

    PeriodicDOFMapping *newMapping = new PeriodicDOFMapping(fcts);
    mappings.push_back(newMapping);
    return newMapping;
  }


  PeriodicDOFMapping::PeriodicDOFMapping(const BasisFunction *fcts)
    : basFcts(fcts)
  {
    FUNCNAME("PeriodicDOFMapping::PeriodicDOFMapping()");
    TEST_EXIT(basFcts->getDim() > 1)("dim == 1\n");
    int nBasFcts = basFcts->getNumber();
    DimVec<double> *lambda;
    for (int i = 0; i < nBasFcts; i++) {
      lambda = basFcts->getCoords(i);
      indexOfCoords[*lambda] = i;
    }
  }


  PeriodicDOFMapping::~PeriodicDOFMapping()
  {
    std::map<DimVec<int>, DegreeOfFreedom*, DimVecLess<int> >::iterator it;
    for (it = dofPermutation.begin(); it != dofPermutation.end(); ++it)
      if (it->second)
	delete [] it->second;
  }


  const DegreeOfFreedom *PeriodicDOFMapping::getDOFPermutation(const DimVec<int> &vertexPermutation)
  {
    FUNCNAME("PeriodicDOFMapping::getDOFPermutation()");

    if (dofPermutation[vertexPermutation] == NULL) {
      int dim = basFcts->getDim();
      int nBasFcts = basFcts->getNumber();
      int sum = 0;
      for (int i = 0; i < dim + 1; i++) {
	sum += i - vertexPermutation[i];
	TEST_EXIT(vertexPermutation[i] < dim + 1)
	  ("invalid vertexPermuation\n");
      }
      TEST_EXIT(sum == 0)("invalid vertexPermutation\n");

      // create dof permutation
      DimVec<double> *lambda;
      DimVec<double> newLambda(dim, NO_INIT);

      DegreeOfFreedom *mapping = new DegreeOfFreedom[nBasFcts];

      for (int i = 0; i < nBasFcts; i++) {
	lambda = basFcts->getCoords(i);
	for (int j = 0; j < dim + 1; j++)
	  newLambda[vertexPermutation[j]] = (*lambda)[j];

	mapping[i] = indexOfCoords[newLambda];
      }

      dofPermutation[vertexPermutation] = mapping;
    }

    return dofPermutation[vertexPermutation];
  }


  PeriodicBC::PeriodicBC(BoundaryType type, const FiniteElemSpace *rowSpace, bool diagonal)
    : BoundaryCondition(type, rowSpace, NULL),
      masterMatrix(NULL),
      isDiagonal(diagonal)
  {
    if (rowFeSpace->getMesh()->getDim() > 1)
      periodicDOFMapping =
        PeriodicDOFMapping::providePeriodicDOFMapping(rowFeSpace->getBasisFcts());
    else
      periodicDOFMapping = NULL;
  }


  PeriodicBC::~PeriodicBC()
  {}


  void PeriodicBC::initMatrix(DOFMatrix* matrix)
  {
    FUNCNAME("PeriodicBC::initMatrix()");

    if (!masterMatrix) {
      masterMatrix = matrix;
      Mesh *mesh = matrix->getRowFeSpace()->getMesh();
      associated = &(mesh->getPeriodicAssociations(boundaryType, rowFeSpace->getAdmin()));

      TEST_EXIT(associated)
        ("No associations for periodic boundary condition %d!\n", boundaryType);
    }
  }


  void PeriodicBC::fillBoundaryCondition(DOFMatrix *matrix,
					 ElInfo *elInfo,
					 const DegreeOfFreedom *dofIndices,
					 const BoundaryType *localBound,
					 int nBasFcts)
  {
    FUNCNAME_DBG("PeriodicBC::fillBoundaryCondition()");

    if (matrix != masterMatrix)
      return;

    int dim = rowFeSpace->getMesh()->getDim();
    if (dim == 1)
      return;

    DOFAdmin *admin = rowFeSpace->getAdmin();
    int iadmin = rowFeSpace->getMesh()->getAdminIndex(admin);
    FixVec<int, WORLD> elFace(dim, NO_INIT);
    FixVec<int, WORLD> neighFace(dim, NO_INIT);
    DimVec<int> vertexPermutation(dim, NO_INIT);
    const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
    int num = basFcts->getNumber();
    Element *element = elInfo->getElement();
    DimVec<DegreeOfFreedom> periodicDOFs(dim - 1, NO_INIT);
    GeoIndex sideGeoIndex = INDEX_OF_DIM(dim - 1, dim);
    std::vector<DegreeOfFreedom> neighIndices(num);

    for (int side = 0; side <= dim; side++) {
      if (elInfo->getBoundary(sideGeoIndex, side) == boundaryType) {
	for (int vertex = 0; vertex < dim; vertex++) {
	  int index = element->getVertexOfPosition(sideGeoIndex, side, vertex);
	  periodicDOFs[vertex] = (*associated)[dofIndices[index]];
	}

	Element *neigh = elInfo->getNeighbour(side);
	TEST_EXIT_DBG(neigh)("Wrong neighbour information at side %d!\n", side);

	basFcts->getLocalIndices(neigh, admin, neighIndices);

	int oppVertex = 0;
	for (int i = 0; i < dim + 1; i++) {
	  // get vertex permutation
	  if (i == side) {
	    vertexPermutation[i] = 0;
	  } else {
	    DegreeOfFreedom periodicDOF =
	      periodicDOFs[element->getPositionOfVertex(side, i)];

	    int j = 0;
	    for (; j < dim + 1; j++)
          if (neigh->getDof(j, iadmin) == periodicDOF)
		break;

	    vertexPermutation[i] = j;
	  }
	  oppVertex += i - vertexPermutation[i];
	}
	vertexPermutation[side] = oppVertex;

	// get DOF permutation
	const DegreeOfFreedom *dofPermutation =
	  periodicDOFMapping->getDOFPermutation(vertexPermutation);

	// set associated dofs
	for (int i = 0; i < num; i++)
	  if ((*(basFcts->getCoords(i)))[side] == 0)
	    (*associated)[dofIndices[i]] = neighIndices[dofPermutation[i]];
      }
    }
  }


  void PeriodicBC::exitMatrix(DOFMatrix* matrix)
  {
    FUNCNAME("PeriodicBC::exitMatrix()");

    TEST_EXIT(matrix)("No matrix\n");
    TEST_EXIT(associated)("No associated vector!\n");

    if (matrix == masterMatrix)
      masterMatrix = NULL;

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

    DOFAdmin* admin = rowFeSpace->getAdmin();
    int adminSize = admin->getUsedSize();
    Matrix &A= matrix->getBaseMatrix();

    mtl::traits::col<Matrix>::type c(A);
    mtl::traits::value<Matrix>::type v(A);

    std::vector<std::vector<std::pair<DegreeOfFreedom, double> > > row_values;
    row_values.resize(adminSize);

    for (DegreeOfFreedom i = 0; i < adminSize; i++) {
      if (i < (*associated)[i]) {
        c_type cursor(begin<mtl::tag::row>(A) + i);

        for (ic_type icursor(begin<mtl::tag::nz>(cursor)), icend(end<mtl::tag::nz>(cursor));
            icursor != icend; ++icursor) {
          row_values[i].push_back(std::make_pair(c(*icursor), v(*icursor)));
          v(*icursor, 0.0);
        }
      }
    }
    matrix::inserter<Matrix, update_plus<double> > ins(A);
    if (isDiagonal) {
      for (DegreeOfFreedom i = 0; i < adminSize; i++) {
        if (i < (*associated)[i]) {
          ins[i][i] << 1.0;
          ins[i][(*associated)[i]] << -1.0;
        }
      }
    }
    for (DegreeOfFreedom i = 0; i < adminSize; i++) {
      if (i < (*associated)[i]) {
        for (size_t j = 0; j < row_values[i].size(); j++)
          ins[(*associated)[i]][row_values[i][j].first] << row_values[i][j].second;
      }
    }
  }


  void PeriodicBC::exitVector(DOFVectorBase<double>* vector)
  {
    FUNCNAME("PeriodicBC::exitVector()");
    TEST_EXIT(rowFeSpace == vector->getFeSpace())("Should not happen.\n");

    DOFIterator<double> vecIt(vector, USED_DOFS);
    Mesh *mesh = vector->getFeSpace()->getMesh();
    VertexVector *associated = &(mesh->getPeriodicAssociations(boundaryType, rowFeSpace->getAdmin()));

    for (vecIt.reset(); !vecIt.end(); ++vecIt) {
      DegreeOfFreedom dof = vecIt.getDOFIndex();
      DegreeOfFreedom newDOF = (*associated)[dof];

      if (dof < newDOF) {
        (*vector)[newDOF] += (*vector)[dof];
        (*vector)[dof] = 0.0;
      }
    }
  }

}