MeshManipulation.cc 5.5 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
    FUNCNAME("MeshManipulation::deleteDoubleDofs()");

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

28
29
30

    // Create data structure that maps macro element indices to the 
    // corresponding pointers.
31
32
33
34
35
36
37
38
    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;

39
40
41
42
    
    // === Traverse mesh and put all "old" macro element to macrosProcessed  ===
    // === that stores all macro elements which are really conected to the   ===
    // === overall mesh structure.                                           ===
43

44
    TraverseStack stack;
45
    ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL);
46
    while (elInfo) {
47
48
49
50
51
52
53
      if (newMacroEl.count(elInfo->getMacroElement()) == 0) {
	int index = elInfo->getMacroElement()->getIndex();

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

54
55
56
57
      elInfo = stack.traverseNext(elInfo);
    }


58
59
60
61
    // === 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.                                                   ===

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

65
66
67
68
      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);
69

70
71
72
73
	for (std::vector<ElementObjectData>::iterator elIt = vertexEl.begin();
	     elIt != vertexEl.end(); ++elIt) {
	  if (elIt->elIndex == (*it)->getIndex())
	    continue;
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
119
120
121
122
123
124
125
126
127
128
129
	  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;
130
131
132
133
	  }
	}
      }

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

136
137

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

140

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

143
144
145
146
147
148
149
    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()) {
150
151
152
	  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));
153
154
155
156
157
158
159
	  it->second = findIt->second;   
	  changed = true;
	}
      } 
    } while (changed);


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

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

      elInfo = stack.traverseNext(elInfo);
    }

171

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

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

179
  
180
}