MeshManipulation.cc

//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


#include "parallel/MeshManipulation.h"
#include "Mesh.h"
#include "MeshStructure.h"
#include "BasisFunction.h"
#include "Traverse.h"
#include "Debug.h"

namespace AMDiS {

  MeshManipulation::MeshManipulation(FiniteElemSpace *space)
    : feSpace(space),
      mesh(space->getMesh())
  {
    switch (mesh->getDim()) {
    case 2:
      refineManager = new RefinementManager2d();
      break;
    case 3:
      refineManager = new RefinementManager3d();
      break;
    default:
      ERROR_EXIT("invalid dim!\n");
    }
  }


  MeshManipulation::~MeshManipulation()
  {
    delete refineManager;
  }


  void MeshManipulation::deleteDoubleDofs(std::set<MacroElement*>& newMacroEl, 
					  ElementObjects &objects)
  {
    FUNCNAME("MeshManipulation::deleteDoubleDofs()");

    // Create data structure that maps macro element indices to the 
    // corresponding pointers.
    map<int, MacroElement*> macroIndexMap;
    for (std::set<MacroElement*>::iterator it = newMacroEl.begin();
	 it != newMacroEl.end(); ++it)
      macroIndexMap[(*it)->getIndex()] = *it;

    std::set<int> macrosProcessed;
    map<const DegreeOfFreedom*, const DegreeOfFreedom*> mapDelDofs;

    
    // === Traverse mesh and put all "old" macro element to macrosProcessed  ===
    // === that stores all macro elements which are really connected to the  ===
    // === overall mesh structure.                                           ===

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
      if (newMacroEl.count(elInfo->getMacroElement()) == 0) {
	int index = elInfo->getMacroElement()->getIndex();

	macrosProcessed.insert(index);
	macroIndexMap[index] = elInfo->getMacroElement();
      }

      elInfo = stack.traverseNext(elInfo);
    }


    // === Traverse all new elements and connect them to the overall mesh     ===
    // === structure by deleting all double DOFs on macro element's vertices, ===
    // === edges and faces.                                                   ===

    for (std::set<MacroElement*>::iterator it = newMacroEl.begin(); 
	 it != newMacroEl.end(); ++it) {

      // === Traverse all vertices of the new element. ===

      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
	vector<ElementObjectData> &vertexEl = 
	  objects.getElementsVertex((*it)->getIndex(), i);

	for (vector<ElementObjectData>::iterator elIt = vertexEl.begin();
	     elIt != vertexEl.end(); ++elIt) {
	  if (elIt->elIndex == (*it)->getIndex())
	    continue;

	  if (macrosProcessed.count(elIt->elIndex) == 1) {
	    TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex) == 1)
	      ("Should not happen!\n");

	    Element *el0 = (*it)->getElement();
	    Element *el1 = macroIndexMap[elIt->elIndex]->getElement();

	    const DegreeOfFreedom *dof0 = el0->getDof(i);
	    const DegreeOfFreedom *dof1 = el1->getDof(elIt->ithObject);

	    mapDelDofs[dof0] = dof1;
    
	    break;
	  } 
	}
      }

      for (int i = 0; i < mesh->getGeo(EDGE); i++) {
	ElementObjectData elObj((*it)->getIndex(), i);

	vector<ElementObjectData> &edgeEl = 
	  objects.getElementsEdge((*it)->getIndex(), i);

	for (vector<ElementObjectData>::iterator elIt = edgeEl.begin();
	     elIt != edgeEl.end(); ++elIt) {
	  if (elIt->elIndex == (*it)->getIndex())
	    continue;

	  if (macrosProcessed.count(elIt->elIndex) == 1) {
	    TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex))("Should not happen!\n");

	    Element *el0 = (*it)->getElement();	    
	    Element *el1 = macroIndexMap[elIt->elIndex]->getElement();

	    bool reverseMode = objects.getEdgeReverseMode(elObj, *elIt);

	    BoundaryObject b0(el0, 0, EDGE, i, reverseMode);
	    BoundaryObject b1(el1, 0, EDGE, elIt->ithObject, false);

	    DofContainer dofs0, dofs1;
	    
	    el0->getVertexDofs(feSpace, b0, dofs0);
	    el0->getNonVertexDofs(feSpace, b0, dofs0);
	    el1->getVertexDofs(feSpace, b1, dofs1);
	    el1->getNonVertexDofs(feSpace, b1, dofs1);

#if (DEBUG != 0)
	    //	    debug::testDofsByCoords(feSpace, dofs0, dofs1);
#endif

 	    for (unsigned int i = 0; i < dofs0.size(); i++)
 	      mapDelDofs[dofs0[i]] = dofs1[i];

	    break;
	  }
	}
      }


      for (int i = 0; i < mesh->getGeo(FACE); i++) {
	ElementObjectData elObj((*it)->getIndex(), i);

	vector<ElementObjectData> &faceEl = 
	  objects.getElementsFace((*it)->getIndex(), i);

	for (vector<ElementObjectData>::iterator elIt = faceEl.begin();
	     elIt != faceEl.end(); ++elIt) {
	  if (elIt->elIndex == (*it)->getIndex())
	    continue;

	  if (macrosProcessed.count(elIt->elIndex) == 1) {
	    TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex))("Should not happen!\n");

	    Element *el0 = (*it)->getElement();	    
	    Element *el1 = macroIndexMap[elIt->elIndex]->getElement();

	    bool reverseMode = objects.getFaceReverseMode(elObj, *elIt);

	    BoundaryObject b0(el0, 0, FACE, i, reverseMode);
	    BoundaryObject b1(el1, 0, FACE, elIt->ithObject, false);

	    DofContainer dofs0, dofs1;
	    
	    el0->getVertexDofs(feSpace, b0, dofs0);
	    el0->getNonVertexDofs(feSpace, b0, dofs0);
	    el1->getVertexDofs(feSpace, b1, dofs1);
	    el1->getNonVertexDofs(feSpace, b1, dofs1);

#if (DEBUG != 0)
	    //	    debug::testDofsByCoords(feSpace, dofs0, dofs1);
#endif

 	    for (unsigned int i = 0; i < dofs0.size(); i++)
 	      mapDelDofs[dofs0[i]] = dofs1[i];

	    break;
	  }
	}
      }

      macrosProcessed.insert((*it)->getIndex());
    }


    // === Remove all DOF replacments of the form A -> B, B -> C by A -> C. ===

    bool changed = false;
    do {
      changed = false;
      for (map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin();
	   it != mapDelDofs.end(); ++it) {
	map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator findIt = mapDelDofs.find(it->second);
	if (findIt != mapDelDofs.end()) {
	  TEST_EXIT_DBG(it->first != findIt->second)
	    ("Cycle %d -> %d and %d -> %d! Should not happen!\n",
	     *(it->first), *(it->second), *(findIt->first), *(findIt->second));
	  it->second = findIt->second;   
	  changed = true;
	}
      } 
    } while (changed);


    // === Set new DOF pointers in all elements of the mesh. ===

    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      for (int i = 0; i < mesh->getGeo(VERTEX); i++)
	if (mapDelDofs.count(elInfo->getElement()->getDof(i)) == 1)
	  elInfo->getElement()->setDof(i, const_cast<DegreeOfFreedom*>(mapDelDofs[elInfo->getElement()->getDof(i)]));	

      elInfo = stack.traverseNext(elInfo);
    }


    // === And delete all double DOFs. ===

    for (map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin();
	 it != mapDelDofs.end(); ++it)
      mesh->freeDof(const_cast<DegreeOfFreedom*>(it->first), VERTEX);
  }


  bool MeshManipulation::fitElementToMeshCode(MeshStructure &code,
					      BoundaryObject &boundEl)
  {
    FUNCNAME("MeshManipulation::fitElementToMeshCode()");

    TEST_EXIT_DBG(boundEl.el)("No element given!\n");

    // If the code is empty, the element does not matter and the function can
    // return without chaning the mesh.
    if (code.empty())
      return false;

    // s0 and s1 are the number of the edge/face in both child of the element,
    // which contain the edge/face the function has to traverse through. If the
    // edge/face is not contained in one of the children, s0 or s1 is -1.
    int s0 = 
      boundEl.el->getSubObjOfChild(0, boundEl.subObj, boundEl.ithObj, boundEl.elType);
    int s1 = 
      boundEl.el->getSubObjOfChild(1, boundEl.subObj, boundEl.ithObj, boundEl.elType);

    TEST_EXIT_DBG(s0 != -1 || s1 != -1)("This should not happen!\n");

    bool meshChanged = false;
    Flag traverseFlag = 
      Mesh::CALL_EVERY_EL_PREORDER | Mesh::FILL_NEIGH | Mesh::FILL_BOUND;

    // Test for reverse mode, in which the left and right children of elements
    // are flipped.
    if (boundEl.reverseMode)
      traverseFlag |= Mesh::CALL_REVERSE_MODE;    


    // === If the edge/face is contained in both children. ===

    Mesh *mesh = boundEl.el->getMesh();

    if (s0 != -1 && s1 != -1) {
      // Create traverse stack and traverse within the mesh until the element,
      // which should be fitted to the mesh structure code, is reached.
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
      while (elInfo && elInfo->getElement() != boundEl.el)
	elInfo = stack.traverseNext(elInfo);      

      TEST_EXIT_DBG(elInfo->getElement() == boundEl.el)("This should not happen!\n");

      return fitElementToMeshCode(code, stack, boundEl.subObj, 
				  boundEl.ithObj, boundEl.reverseMode);
    }


    // === The edge/face is contained in only on of the both children. ===

    if (boundEl.el->isLeaf()) {

      // If element is leaf and code contains only one leaf element, we are finished.
      if (code.getNumElements() == 1 && code.isLeafElement())
	return false;     

      // Create traverse stack and traverse the mesh to the element.
      TraverseStack stack;
      ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);
      while (elInfo && elInfo->getElement() != boundEl.el)
	elInfo = stack.traverseNext(elInfo);      

      TEST_EXIT_DBG(elInfo)("This should not happen!\n");

      // Code is not leaf, therefore refine the element.
      boundEl.el->setMark(1);
      refineManager->setMesh(boundEl.el->getMesh());
      refineManager->setStack(&stack);
      refineManager->refineFunction(elInfo);
      meshChanged = true;
    }

    Element *child0 = boundEl.el->getFirstChild();
    Element *child1 = boundEl.el->getSecondChild();
    if (boundEl.reverseMode) {
      swap(s0, s1);
      swap(child0, child1);    
    }

    // === We know that the edge/face is contained in only one of the children. ===
    // === Therefore, traverse the mesh to this children and fit this element   ===
    // === To the mesh structure code.                                          ===

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, traverseFlag);

    if (s0 != -1) {
      while (elInfo && elInfo->getElement() != child0)
	elInfo = stack.traverseNext(elInfo);     

      meshChanged |= 
	fitElementToMeshCode(code, stack, boundEl.subObj, s0, boundEl.reverseMode);
    } else {
      while (elInfo && elInfo->getElement() != child1) 
	elInfo = stack.traverseNext(elInfo);      

      meshChanged |= 
	fitElementToMeshCode(code, stack, boundEl.subObj, s1, boundEl.reverseMode);
    }


    return meshChanged;
  }


  bool MeshManipulation::fitElementToMeshCode(MeshStructure &code, 
					      TraverseStack &stack,
					      GeoIndex subObj,
					      int ithObj, 
					      bool reverseMode)
  {
    FUNCNAME("MeshManipulation::fitElementToMeshCode()");


    // === Test if there are more elements in stack to check with the code. ===

    ElInfo *elInfo = stack.getElInfo();
    if (!elInfo)
      return false;


    // === Test if code contains a leaf element. If this is the case, the ===
    // === current element is finished. Traverse the mesh to the next     ===
    // === coarser element.                                               ===

    if (code.isLeafElement()) {
      int level = elInfo->getLevel();

      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);

      return false;
    }


    bool meshChanged = false;
    Element *el = elInfo->getElement();


    // === If element is leaf (and the code is not), refine the element. ===

    if (el->isLeaf()) {
      TEST_EXIT_DBG(elInfo->getLevel() < 255)("This should not happen!\n");

      el->setMark(1);
      refineManager->setMesh(el->getMesh());
      refineManager->setStack(&stack);
      refineManager->refineFunction(elInfo);
      meshChanged = true;
    }


    // === Continue fitting the mesh structure code to the children of the ===
    // === current element.                                                ===

    int s0 = el->getSubObjOfChild(0, subObj, ithObj, elInfo->getType());
    int s1 = el->getSubObjOfChild(1, subObj, ithObj, elInfo->getType());
    Element *child0 = el->getFirstChild();
    Element *child1 = el->getSecondChild();
    if (reverseMode) {
      swap(s0, s1);
      swap(child0, child1);
    }

    
    // === Traverse left child. ===

    if (s0 != -1) {
      // The edge/face is contained in the left child, therefore fit this
      // child to the mesh structure code.

      stack.traverseNext(elInfo);
      code.nextElement();
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s0, reverseMode);
      elInfo = stack.getElInfo();
    } else {
      // The edge/face is not contained in the left child. Hence we need
      // to traverse through all subelements of the left child until we get
      // the second child of the current element.

      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getElement() != child1); 

      TEST_EXIT_DBG(elInfo != NULL)("This should not happen!\n");
    }  

    TEST_EXIT_DBG(elInfo->getElement() == child1)
      ("Got wrong child with idx = %d! Searched for child idx = %d\n",
       elInfo->getElement()->getIndex(), child1->getIndex());


    // === Traverse right child. ===

    if (s1 != -1) {
      // The edge/face is contained in the right child, therefore fit this
      // child to the mesh structure code.

      code.nextElement();
      meshChanged |= fitElementToMeshCode(code, stack, subObj, s1, reverseMode);
    } else {
      // The edge/face is not contained in the right child. Hence we need
      // to traverse through all subelements of the right child until we have
      // finished traversing the current element with all its subelements.

      int level = elInfo->getLevel();

      do {
	elInfo = stack.traverseNext(elInfo);
      } while (elInfo && elInfo->getLevel() > level);
    }


    return meshChanged;
  }
  
}