MeshStructure.cc

#include "MeshStructure.h"
#include "MeshStructure_ED.h"
#include "PartitionElementData.h"
#include "Mesh.h"
#include "Element.h"
#include "Traverse.h"
#include "ElInfo.h"
#include "RefinementManager.h"

namespace AMDiS {

  const int MeshStructure::unsignedLongSize = sizeof(unsigned long int) * 8;


  void MeshStructure::insertElement(bool isLeaf) 
  {
    // overflow? -> next index
    if (pos >= unsignedLongSize) {
      code.push_back(currentCode);
      pos = 0;
      currentCode = 0;
    }

    // insert element in binary code
    if (!isLeaf) {
      unsigned long int one = 1;
      currentCode += (one << pos);
    } 

    pos++;
    nElements++;
  }


  void MeshStructure::clear()
  {
    currentCode = 0;
    code.resize(0);
    pos = 0;
    nElements = 0;
    currentElement = 0;
  }


  void MeshStructure::init(Mesh *mesh) 
  {
    clear();

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      insertElement(elInfo->getElement()->isLeaf());
      elInfo = stack.traverseNext(elInfo);
    }
  
    commit();
  }


  void MeshStructure::init(BoundaryObject &bound)
  {
    FUNCNAME("MeshStructure::init()");

    Element *el = bound.el;

    TEST_EXIT_DBG(el)("No element!\n");

    clear();

    int s1 = el->getSubObjOfChild(0, bound.subObj, bound.ithObj, bound.elType);
    int s2 = el->getSubObjOfChild(1, bound.subObj, bound.ithObj, bound.elType);
    
    TEST_EXIT(s1 != -1 || s2 != -1)("This should not happen!\n");

    if (debugMode) {
      MSG("addAlondSide(%d, %d, %d, %d)\n",
	  bound.elIndex, bound.ithObj, bound.elType, bound.reverseMode);
      MSG("Element is leaf: %d\n", el->isLeaf());
      MSG("s1 = %d    s2 = %d\n", s1, s2);
    }
    
    if (!el->isLeaf()) {
      if (s1 == -1)
	addAlongSide(el->getSecondChild(), bound.subObj, s2, 
		     el->getChildType(bound.elType), bound.reverseMode);
      else if (s2 == -1)
	addAlongSide(el->getFirstChild(), bound.subObj, s1, 
		     el->getChildType(bound.elType), bound.reverseMode);
      else
	addAlongSide(el, bound.subObj, bound.ithObj, bound.elType, bound.reverseMode);
    }

    commit();    
  }


  void MeshStructure::addAlongSide(Element *el, GeoIndex subObj, int ithObj, 
				   int elType, bool reverseOrder)
  {
    FUNCNAME("MeshStructure::addAlongSide()");

    if (debugMode) {
      MSG("addAlondSide(%d, %d, %d, %d)\n",
	  el->getIndex(), ithObj, elType, reverseOrder);
      MSG("Element is leaf: %d\n", el->isLeaf());
    }
    
    insertElement(el->isLeaf());
    
    if (!el->isLeaf()) {
      int s1 = el->getSubObjOfChild(0, subObj, ithObj, elType);
      int s2 = el->getSubObjOfChild(1, subObj, ithObj, elType);

      if (debugMode) {
	MSG("Child index %d  %d\n", 
	    el->getFirstChild()->getIndex(),
	    el->getSecondChild()->getIndex());
	MSG("s1 = %d    s2 = %d\n", s1, s2);
	MSG("   \n");
      }

      if (!reverseOrder) {
	if (s1 != -1) 
	  addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), reverseOrder);
	if (s2 != -1)
	  addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), reverseOrder);
      } else {
	if (s2 != -1)
	  addAlongSide(el->getSecondChild(), subObj, s2, el->getChildType(elType), reverseOrder);
	if (s1 != -1) 
	  addAlongSide(el->getFirstChild(), subObj, s1, el->getChildType(elType), reverseOrder);
      }
    }
  }


  void MeshStructure::reset() 
  {
    currentIndex = 0;
    pos = 0;
    currentElement = 0;

    if (code.size() > 0)
      currentCode = code[0];
    else 
      currentCode = 0;
  }


  bool MeshStructure::nextElement(MeshStructure *insert)
  {
    FUNCNAME("MeshStructure::nextElement()");

    if (insert)
      insert->insertElement(isLeafElement());

    pos++;
    currentElement++;

    if (currentElement >= nElements) 
      return false;

    if (pos >= unsignedLongSize) {
      currentIndex++;
      TEST_EXIT_DBG(currentIndex < static_cast<int>(code.size()))
	("End of structure reached!\n");
      pos = 0;
      currentCode = code[currentIndex];
    } else {
      currentCode >>= 1;
    }

    return true;
  }


  bool MeshStructure::skipBranch(MeshStructure *insert)
  {
    FUNCNAME("MeshStructure::skipBranch()");

    if (isLeafElement()) {
      return nextElement(insert);
    } else {
      bool cont = nextElement(insert);
      cont = skipBranch(insert); // left branch
      TEST_EXIT_DBG(cont)("Invalid structure!\n");
      cont = skipBranch(insert); // righ branch
      return cont;
    }
  }


  void MeshStructure::merge(MeshStructure *structure1,
			    MeshStructure *structure2,
			    MeshStructure *result)
  {
    FUNCNAME("MeshStructure::merge()");

    result->clear();
    structure1->reset();
    structure2->reset();

    bool cont = true;
    while (cont) {
      bool cont1, cont2;
      if (structure1->isLeafElement() == structure2->isLeafElement()) {
	cont1 = structure1->nextElement(result);
	cont2 = structure2->nextElement();
      } else {
	if (structure1->isLeafElement()) {
	  cont1 = structure1->nextElement();
	  cont2 = structure2->skipBranch(result);
	} else {
	  cont1 = structure1->skipBranch(result);
	  cont2 = structure2->nextElement();
	}
      }
      TEST_EXIT_DBG(cont1 == cont2)("Structures don't match!\n");
      cont = cont1;
    }

    result->commit();
  }


  void MeshStructure::fitMeshToStructure(Mesh *mesh,
					 RefinementManager *manager,
					 bool checkPartition,
					 bool debugMode,
					 int macroElIndex) 
  {
    FUNCNAME("MeshStructure::fitMeshToStructure()");

    bool cont = true;

    // decorate leaf data
    reset();
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      if (macroElIndex >= 0 && 
	  elInfo->getMacroElement()->getElement()->getIndex() != macroElIndex) {
	elInfo = stack.traverseNext(elInfo);
	continue;
      }

      Element *element = elInfo->getElement();

      TEST_EXIT(cont)("unexpected structure code end!\n");

      if (isLeafElement()) {
	TEST_EXIT(element->isLeaf())("mesh finer than code\n");
      }

      if (element->isLeaf() && !isLeafElement()) {
	MeshStructure *structure = new MeshStructure();
	cont = skipBranch(structure);
	structure->commit();

	bool decorate = true;
	if (checkPartition) {
	  PartitionElementData *partitionData = dynamic_cast<PartitionElementData*>
	    (element->getElementData(PARTITION_ED));
	  TEST_EXIT_DBG(partitionData)("no partition element data\n");
	  PartitionStatus status = partitionData->getPartitionStatus();
	  if (debugMode == false && (status == OUT || status == UNDEFINED))
	    decorate = false;
	}

	if (decorate) {
	  MeshStructure_ED *elData = new MeshStructure_ED(element->getElementData());
	  elData->setStructure(structure);
	  element->setElementData(elData);
	} else {
	  delete structure;
	}
      } else {
	cont = nextElement();
      }

      elInfo = stack.traverseNext(elInfo);
    }

    // refine mesh
    bool finished = true;

    do {
      finished = true;
      elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
      while (elInfo) {
	if (macroElIndex >= 0 && 
	    elInfo->getMacroElement()->getElement()->getIndex() != macroElIndex) {
	  elInfo = stack.traverseNext(elInfo);
	  continue;
	}

	Element *element = elInfo->getElement();
	if (element->getElementData(MESH_STRUCTURE) != NULL) {
	  element->setMark(1);
	  finished = false;
	} else {
	  element->setMark(0);
	}
	elInfo = stack.traverseNext(elInfo);
      }
      manager->refineMesh(mesh);
    } while (!finished);
  }


  bool MeshStructure::compare(MeshStructure &other)
  {
    return (other.getCode() == code);
  }


  void MeshStructure::writeMeshFile(Mesh *mesh, std::string filename)
  {
    FUNCNAME("MeshStructure::writeMeshFile()");

    int nMacroElements = 0;
    int macroElIndex = -1;
    std::ofstream file;
    file.open(filename.c_str(), std::ios::out | std::ios::binary | std::ios::trunc);

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
    while (elInfo) {
      nMacroElements++;
      elInfo = stack.traverseNext(elInfo);
    }

    file.write(reinterpret_cast<char*>(&nMacroElements), sizeof(nMacroElements));
    
    elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER);
    while (elInfo) {
      if (elInfo->getLevel() == 0) {
	if (macroElIndex != -1) {
	  commit();
	  writeMacroElement(file, macroElIndex);
	}

	clear();
	macroElIndex = elInfo->getElement()->getIndex();
      }

      insertElement(elInfo->getElement()->isLeaf());
      elInfo = stack.traverseNext(elInfo);
    }

    // And write the last macro element to file.
    TEST_EXIT_DBG(macroElIndex != -1)("Should not happen!\n");
    commit();
    writeMacroElement(file, macroElIndex);

    file.close();
  }


  void MeshStructure::readMeshFile(Mesh *mesh, 
				   RefinementManager *refManager,
				   std::string filename)
  {
    FUNCNAME("MeshStructure::readMeshFile()");

    std::ifstream file;
    file.open(filename.c_str(), std::ios::in | std::ios::binary);

    int nMacroElements = -1;
    file.read(reinterpret_cast<char*>(&nMacroElements), sizeof(nMacroElements));

    for (int i = 0; i < nMacroElements; i++) {
      int macroElIndex = -1;
      int codeSize = 0;
      code.clear();

      file.read(reinterpret_cast<char*>(&macroElIndex), sizeof(macroElIndex));
      file.read(reinterpret_cast<char*>(&nElements), sizeof(nElements));
      file.read(reinterpret_cast<char*>(&codeSize), sizeof(codeSize));

      for (int j = 0; j < codeSize; j++) {
	unsigned long int readCode = 0;
	file.read(reinterpret_cast<char*>(&readCode), sizeof(readCode));
	code.push_back(readCode);
      }
	
      TEST_EXIT_DBG(macroElIndex >= 0)("Should not happen!\n");

      reset();
      fitMeshToStructure(mesh, refManager, false, false, macroElIndex);
      
    }

    file.close();
  }


  void MeshStructure::writeMacroElement(std::ofstream& file, int macroElIndex)
  {
    unsigned int size = code.size();

    file.write(reinterpret_cast<char*>(&macroElIndex), sizeof(macroElIndex));
    file.write(reinterpret_cast<char*>(&nElements), sizeof(nElements));
    file.write(reinterpret_cast<char*>(&size), sizeof(size));

    for (unsigned int i = 0; i < code.size(); i++)
      file.write(reinterpret_cast<char*>(&(code[i])), sizeof(code[i]));
  }

}