Element.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 "Element.h"
#include "DOFAdmin.h"
#include "Mesh.h"
#include "CoarseningManager.h"
#include "FixVec.h"
#include "ElementRegion_ED.h"
#include "Serializer.h"
#include "MeshStructure.h"
#include "BasisFunction.h"

namespace AMDiS {

  std::map<DegreeOfFreedom*, bool> Element::deletedDOFs;


  Element::Element(Mesh *aMesh)
  {
    mesh = aMesh;
    index = mesh ? mesh->getNextElementIndex() : -1; 
    child[0] = NULL;
    child[1] = NULL;
    newCoord = NULL;
    elementData = NULL;
    mark = 0;

    if (mesh)
      createNewDofPtrs();
    else
      dof = NULL;    
  }


  // call destructor through Mesh::freeElement !!!
  Element::~Element()
  {  
    if (child[0])
      delete child[0];
    
    if (child[1])
      delete child[1];   

    if (newCoord)
      delete newCoord;    

    if (elementData) {
      elementData->deleteDecorated();
      delete elementData;
    }
  }


  int Element::getRegion() const 
  {
    if (!elementData) 
      return -1;

    ElementRegion_ED* red = 
      dynamic_cast<ElementRegion_ED*>(elementData->getElementData(ELEMENT_REGION));

    if (red)
      return red->getRegion();    
    
    return -1;
  }


  void Element::createNewDofPtrs()
  {
    FUNCNAME("Element::setDofPtrs()");

    TEST_EXIT_DBG(mesh)("no mesh!\n");

    dof = mesh->createDofPtrs();  
  }


  bool Element::deleteElementData(int typeID)
  {
    FUNCNAME("Element::deleteElementData()");

    if (elementData) {
      if (elementData->isOfType(typeID)) {
	ElementData *tmp = elementData;
	elementData = elementData->getDecorated();
	delete tmp;
	tmp = NULL;
	return true;
      } else {
	return elementData->deleteDecorated(typeID);
      }
    }
    return false;
  }


  void Element::deleteElementDOFs()
  {
    int dim = mesh->getDim();
    int j = 0;

    for (int pos = 0; pos <= dim; pos++) {
      GeoIndex position = INDEX_OF_DIM(pos, dim);
      int ndof = 0;
     
      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
	ndof += mesh->getDofAdmin(i).getNumberOfDofs(position);

      if (ndof > 0) {
	for (int i = 0; i < mesh->getGeo(position); i++) {
	  if (dof[j]) {
	    if (deletedDOFs.count(dof[j]) == 0) {
	      deletedDOFs[dof[j]] = true;
	      delete [] dof[j];
	    }
	  }  
	  j++;
	}
      }
    }

    delete [] dof;
    
    if (child[0])
      child[0]->deleteElementDOFs();
    if (child[1])
      child[1]->deleteElementDOFs();
  }


  Element* Element::cloneWithDOFs()
  {
    Element *el;
    
    if (isLine()) {
      el = new Line(NULL);
    } else if (isTriangle()) {
      el = new Triangle(NULL);
    } else {
      el = new Tetrahedron(NULL);
    }

    el->mesh = mesh;
    el->index = index;
    el->mark = mark;
    if (newCoord) {
      WorldVector<double> *nc = new WorldVector<double>();
      *nc = *newCoord;
      el->newCoord = nc;
    }
    
    /* =========== And here we clone the DOFs =========== */
   
    el->dof = new DegreeOfFreedom*[mesh->getNumberOfNodes()];

    int dim = mesh->getDim();
    int j = 0;

    for (int pos = 0; pos <= dim; pos++) {
      GeoIndex position = INDEX_OF_DIM(pos, dim);
      int ndof = 0;

      for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
	ndof += mesh->getDofAdmin(i).getNumberOfDofs(position);

      if (ndof > 0) {
	for (int i = 0; i < mesh->getGeo(position); i++) {
	  if (dof[j] != NULL) {
	    std::pair<DegreeOfFreedom, int> idx = std::make_pair(dof[j][0], pos);

	    if (Mesh::serializedDOFs[idx] == NULL) {
	      el->dof[j] = new DegreeOfFreedom[ndof];
	      for (int k = 0; k < ndof; k++)
		el->dof[j][k] = dof[j][k];

	      Mesh::serializedDOFs[idx] = el->dof[j];
	    } else {
	      el->dof[j] = Mesh::serializedDOFs[idx];
	    }
	  } else {
	    el->dof[j] = NULL;
	  }
	  j++;
	}
      }
    }
    
    /* =========== And clone the children ============= */

    if (child[0]) 
      el->child[0] = child[0]->cloneWithDOFs();
    if (child[1])
      el->child[1] = child[1]->cloneWithDOFs();

    return el;
  }


  /****************************************************************************/
  /*  ATTENTION:                                                              */
  /*  new_dof_fct() destroys new_dof !!!!!!!!!!                               */
  /*  should be used only at the end of dof_compress()!!!!!                   */
  /****************************************************************************/

  void Element::newDofFct1(const DOFAdmin* admin, std::vector<int>& newDofIndex)
  {
    int n0, nd, nd0;

    if ((nd = admin->getNumberOfDofs(VERTEX)))  {
      int vertices = mesh->getGeo(VERTEX);
      nd0 = admin->getNumberOfPreDofs(VERTEX);
      n0 = admin->getMesh()->getNode(VERTEX);
      for (int i = 0; i < vertices; i++)
	changeDofs1(admin, newDofIndex, n0, nd0, nd, i);
    }

    if (mesh->getDim() > 1) {
      if ((nd = admin->getNumberOfDofs(EDGE)))  {
	int edges = mesh->getGeo(EDGE); 
	nd0 = admin->getNumberOfPreDofs(EDGE);
	n0 = admin->getMesh()->getNode(EDGE);
	for (int i = 0; i < edges; i++)
	  changeDofs1(admin, newDofIndex, n0, nd0, nd, i);
      }
    }

    if (mesh->getDim() == 3) {
      if ((nd = admin->getNumberOfDofs(FACE)))  {
	int faces = mesh->getGeo(FACE);
	nd0 = admin->getNumberOfPreDofs(FACE);
	n0 = admin->getMesh()->getNode(FACE);
	for (int i = 0; i < faces; i++)
	  changeDofs1(admin, newDofIndex, n0, nd0, nd, i);
      }
    }

    if ((nd = admin->getNumberOfDofs(CENTER)))  {
      nd0 = admin->getNumberOfPreDofs(CENTER);
      n0 = admin->getMesh()->getNode(CENTER);      
      changeDofs1(admin, newDofIndex, n0, nd0, nd, 0);
    }
  }


  void Element::newDofFct2(const DOFAdmin* admin)
  {
    int n0, nd0;

    int nd = admin->getNumberOfDofs(VERTEX);
    if (nd) {
      int vertices = mesh->getGeo(VERTEX);
      nd0 = admin->getNumberOfPreDofs(VERTEX);
      n0 = admin->getMesh()->getNode(VERTEX);
      for (int i = 0; i < vertices; i++)
	changeDofs2(n0, nd0, nd, i);
    }

    if (mesh->getDim() > 1) {
      nd = admin->getNumberOfDofs(EDGE);
      if (nd) {
	int edges = mesh->getGeo(EDGE); 
	nd0 = admin->getNumberOfPreDofs(EDGE);
	n0 = admin->getMesh()->getNode(EDGE);
	for (int i = 0; i < edges; i++)
	  changeDofs2(n0, nd0, nd, i);
      }
    }

    if (mesh->getDim() == 3) {
      nd = admin->getNumberOfDofs(FACE);
      if (nd) {
	int faces = mesh->getGeo(FACE);
	nd0 = admin->getNumberOfPreDofs(FACE);
	n0 = admin->getMesh()->getNode(FACE);
	for (int i = 0; i < faces; i++)
	  changeDofs2(n0, nd0, nd, i);
      }
    }

    nd = admin->getNumberOfDofs(CENTER);
    if (nd) {
      nd0 = admin->getNumberOfPreDofs(CENTER);
      n0 = admin->getMesh()->getNode(CENTER);
      // only one center
      changeDofs2(n0, nd0, nd, 0);	  
    }
  }


  void Element::changeDofs1(const DOFAdmin* admin, std::vector<int>& newDofIndex,
			    int n0, int nd0, int nd, int pos)
  {
    FUNCNAME("Element::changeDofs1()");

    DegreeOfFreedom *ldof = dof[n0 + pos] + nd0;
    for (int j = 0; j < nd; j++) {
      int k = ldof[j];
      if (k >= 0)
	ldof[j] = -newDofIndex[k] - 1;
    }
  }
  
  
  void Element::changeDofs2(int n0, int nd0, int nd, int pos)
  {
    DegreeOfFreedom *ldof = dof[n0 + pos] + nd0;
    for (int j = 0; j < nd; j++) {
      int k = ldof[j];
      if (k < 0)
	ldof[j] = -k - 1;      
    }
  }


  /****************************************************************************/
  /* opp_vertex checks whether the face with vertices dof[0],..,dof[DIM-1] is */
  /* part of mel's boundary. returns the opposite vertex if true, -1 else     */
  /****************************************************************************/

  int Element::oppVertex(FixVec<DegreeOfFreedom*, DIMEN> pdof) const
  {
    int nv = 0;
    int ov = 0;
    int vertices = mesh->getGeo(VERTEX);
    int dim = mesh->getDim();

    for (int i = 0; i < vertices; i++) {
      if (nv < i - 1)  
	return(-1);

      for (int j = 0; j < dim; j++) {
	if (dof[i] == pdof[j]) {
	  // i is a common vertex 
	  ov += i;
	  nv++;
	  break;
	}
      }

    }
    
    if (nv != mesh->getDim()) 
      return(-1);

    // The opposite vertex is 3(6) - (sum of indices of common vertices) in 2D(3D)

    switch(mesh->getDim()) {
    case 1:
      return ov;
      break;
    case 2:
      return 3 - ov;
      break;
    case 3:
      return 6 - ov;
      break;
    default:
      ERROR_EXIT("invalid dim\n");
      return 0;
    }
  }


  void Element::eraseNewCoord() 
  {
    if (newCoord != NULL) {
      delete newCoord;
      newCoord = NULL;
    }
  }
 

  void Element::serialize(std::ostream &out) 
  {
    // write children
    if (child[0]) {
      out << child[0]->getTypeName() << "\n";
      child[0]->serialize(out);
      child[1]->serialize(out);
    } else {
      out << "NULL\n";
    }

    // write dofs
    int dim = mesh->getDim();
    int nodes = mesh->getNumberOfNodes();
    int j = 0;
    SerUtil::serialize(out, nodes);
    int dofValid = (dof != NULL ? 1 : 0);
    SerUtil::serialize(out, dofValid);
   
    if (dof != NULL) {
      for (int pos = 0; pos <= dim; pos++) {
	GeoIndex position = INDEX_OF_DIM(pos, dim);
	int ndof = 0;
	
	for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++)
	  ndof += mesh->getDofAdmin(i).getNumberOfDofs(position);	      
	
	if (ndof > 0) {
	  for (int i = 0; i < mesh->getGeo(position); i++) {
	    if (dof[j] != NULL) {
	      // Create index to check if the dofs were already written.
	      std::pair<DegreeOfFreedom, int> idx = std::make_pair(dof[j][0], pos);
	      
	      if (Mesh::serializedDOFs[idx] == NULL) {
		Mesh::serializedDOFs[idx] = dof[j];
		SerUtil::serialize(out, ndof);
		SerUtil::serialize(out, pos);
		out.write(reinterpret_cast<const char*>(dof[j]), 
			  ndof * sizeof(DegreeOfFreedom));
	      } else {
		int minusOne = -1;
		SerUtil::serialize(out, minusOne);
		SerUtil::serialize(out, pos);
		out.write(reinterpret_cast<const char*>(&(dof[j][0])), 
			  sizeof(DegreeOfFreedom));
	      }
	    } else {
	      int zero = 0;
	      SerUtil::serialize(out, zero);
	      SerUtil::serialize(out, pos);
	    }
	    j++;
	  }
	}
      }
    }

    // write index
    SerUtil::serialize(out, index);

    // write mark
    SerUtil::serialize(out, mark);

    // write newCoord
    if (newCoord) {
      out << "WorldVector\n";
      newCoord->serialize(out);
    } else {
      out << "NULL\n";
    }

    // write element data
    if (elementData) {
      out << elementData->getTypeName() << "\n";
      elementData->serialize(out);
    } else {
      out << "NULL\n";
    }
  }


  void Element::deserialize(std::istream &in)
  {
    FUNCNAME("Element::deserialize()");

    std::string typeName = "";

    // read children
    in >> typeName;
    in.get();

    if (typeName != "NULL") {
      if (typeName == "Line") {
	child[0] = new Line(NULL);
	child[1] = new Line(NULL);      
      }else  if (typeName == "Triangle") {
	child[0] = new Triangle(NULL);
	child[1] = new Triangle(NULL);      
      } else  if (typeName == "Tetrahedron") {
	child[0] = new Tetrahedron(NULL);
	child[1] = new Tetrahedron(NULL);      
      } else {
	ERROR_EXIT("Wrong element type!\n");
      }

      child[0]->deserialize(in);
      child[1]->deserialize(in);
    } else {
      child[0] = child[1] = NULL;
    }

    // read dofs
    int nodes;
    SerUtil::deserialize(in, nodes);
    int dofValid = 0;
    SerUtil::deserialize(in, dofValid);

    TEST_EXIT_DBG(nodes > 0)("Should not happen!\n");    
    dof = (dofValid ? new DegreeOfFreedom*[nodes] : NULL);
     
    if (dofValid) {
      for (int i = 0; i < nodes; i++) {
	int nDofs, pos;
	SerUtil::deserialize(in, nDofs);
	SerUtil::deserialize(in, pos);
	
	if (nDofs) {
	  if (nDofs != -1) {
	    dof[i] = new DegreeOfFreedom[nDofs];
	    in.read(reinterpret_cast<char*>(dof[i]), nDofs * sizeof(DegreeOfFreedom));
	    
	    // Create index to check if the dofs were alread read from file.
	    std::pair<DegreeOfFreedom, int> idx = std::make_pair(dof[i][0], pos);
	    
	    if (Mesh::serializedDOFs[idx] != NULL) {
	      DegreeOfFreedom *dofPtr = Mesh::serializedDOFs[idx];
	      delete [] dof[i];
	      dof[i] = dofPtr;
	    } else {
	      Mesh::serializedDOFs[idx] = dof[i];
	    }
	    
	  } else {
	    DegreeOfFreedom index;
	    SerUtil::deserialize(in, index);
	    
	    std::pair<DegreeOfFreedom, int> idx = std::make_pair(index, pos);
	    TEST_EXIT(Mesh::serializedDOFs.find(idx) !=  Mesh::serializedDOFs.end())
	      ("This should never happen!\n");
	    dof[i] = Mesh::serializedDOFs[idx];
	  }
	} else {
	  dof[i] = NULL;
	}
      }
    }
      
    // read index
    SerUtil::deserialize(in, index);

    // read mark
    SerUtil::deserialize(in, mark);

    // read newCoord
    in >> typeName;
    in.get();

    if (typeName != "NULL") {
      if (typeName == "WorldVector") {
	newCoord = new WorldVector<double>;
	newCoord->deserialize(in);
      } else {
	ERROR_EXIT("unexpected type name\n");
      }
    } else {
      newCoord = NULL;
    }

    // read element data
    in >> typeName;
    in.get();

    if (typeName != "NULL") {
      elementData = 
	CreatorMap<ElementData>::getCreator(typeName, "deserializetion from file")->create();

      if (elementData)
	elementData->deserialize(in);
      else
	ERROR_EXIT("unexpected type name\n");      
    } else {
      elementData = NULL;
    }
  }


  void writeDotFile(Element *el, std::string filename, int maxLevels)
  {
    std::vector<int> graphNodes;
    std::vector<std::pair<int, int> > graphEdges;

    std::vector<Element*> traverseElements, nextTraverseElements;
    traverseElements.push_back(el);

    int nLevels = 0;

    while (!traverseElements.empty() && (maxLevels == -1 || nLevels < maxLevels)) {
      nextTraverseElements.clear();

      for (unsigned int i = 0; i < traverseElements.size(); i++) {
	graphNodes.push_back(traverseElements[i]->getIndex());

	if (!traverseElements[i]->isLeaf() && 
	    (maxLevels == -1 || nLevels + 1 < maxLevels)) {
	  graphEdges.push_back(std::make_pair(traverseElements[i]->getIndex(),
					      traverseElements[i]->getChild(0)->getIndex()));
	  graphEdges.push_back(std::make_pair(traverseElements[i]->getIndex(),
					      traverseElements[i]->getChild(1)->getIndex()));
	  nextTraverseElements.push_back(traverseElements[i]->getChild(0));
	  nextTraverseElements.push_back(traverseElements[i]->getChild(1));
	}
      }

      traverseElements = nextTraverseElements;
      nLevels++;
    }

    std::ofstream file;
    file.open(filename.c_str());

    file << "digraph G\n";
    file << "{\n";

    for (unsigned int i = 0; i < graphNodes.size(); i++)
      file << "   node" << graphNodes[i] 
	   << " [ label = \"" << graphNodes[i] << "\"];\n";    

    for (unsigned int i = 0; i < graphEdges.size(); i++) 
      file << "   \"node" << graphEdges[i].first << "\" -> \"node"
	   << graphEdges[i].second << "\";\n";

    file << "}\n";

    file.close();
  }


  void Element::getAllDofs(const FiniteElemSpace* feSpace, 
			   BoundaryObject bound, 
			   DofContainer& dofs,
			   bool baseDofPtr)
  {
    getNodeDofs(feSpace, bound, dofs, baseDofPtr);

    if (feSpace->getBasisFcts()->getDegree() > 1)
      getHigherOrderDofs(feSpace, bound, dofs, baseDofPtr);
  }

}