MeshManipulation.cc 4.73 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// 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.


13
14
#include "parallel/MeshManipulation.h"
#include "Mesh.h"
15
#include "BasisFunction.h"
16
#include "Traverse.h"
17
#include "Debug.h"
18
19
20

namespace AMDiS {

21
22
  void MeshManipulation::deleteDoubleDofs(std::set<MacroElement*>& newMacroEl, 
					  ElementObjects &objects)
23
  {
24
25
26
27
28
29
30
31
32
33
34
35
36
    FUNCNAME("MeshManipulation::deleteDoubleDofs()");

    TEST_EXIT(mesh->getDim() == 2)("Not yet supported for dim != 2!\n");

    std::map<int, MacroElement*> macroIndexMap;
    for (std::set<MacroElement*>::iterator it = newMacroEl.begin();
	 it != newMacroEl.end(); ++it)
      macroIndexMap[(*it)->getIndex()] = *it;

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


37
    TraverseStack stack;
38
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
39
    while (elInfo) {
40
41
42
43
44
45
46
      if (newMacroEl.count(elInfo->getMacroElement()) == 0) {
	int index = elInfo->getMacroElement()->getIndex();

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

47
48
49
50
      elInfo = stack.traverseNext(elInfo);
    }


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

54
55
56
57
      for (int i = 0; i < mesh->getGeo(VERTEX); i++) {
	ElementObjectData elObj((*it)->getIndex(), i);
	DegreeOfFreedom vertex = objects.getVertexLocalMap(elObj);
	std::vector<ElementObjectData> &vertexEl = objects.getElements(vertex);
58

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

64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
	  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);
	DofEdge edge = objects.getEdgeLocalMap(elObj);
	std::vector<ElementObjectData> &edgeEl = objects.getElements(edge);

	for (std::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) == 1)
	      ("Should not happen!\n");

	    TEST_EXIT(mesh->getDim() == 2)("In 3D set correct reverse mode!\n");

	    TEST_EXIT(feSpace->getBasisFcts()->getDegree() == 1)
	      ("Not yet implemented!\n");

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

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

	    DofContainer dofs0, dofs1;
	    
	    el0->getVertexDofs(feSpace, b0, dofs0);
	    el1->getVertexDofs(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;
119
120
121
122
	  }
	}
      }

123
      TEST_EXIT(mesh->getDim() == 2)("Add face traverse here!\n");
124

125
126

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

129

130
131
132
133
134
135
136
    bool changed = false;
    do {
      changed = false;
      for (std::map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin();
	   it != mapDelDofs.end(); ++it) {
	std::map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator findIt = mapDelDofs.find(it->second);
	if (findIt != mapDelDofs.end()) {
137
138
139
	  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));
140
141
142
143
144
145
146
	  it->second = findIt->second;   
	  changed = true;
	}
      } 
    } while (changed);


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

      elInfo = stack.traverseNext(elInfo);
    }

156

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

162
  
163
}