MeshManipulation.cc 4.93 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
15
#include "parallel/MeshManipulation.h"
#include "Mesh.h"
#include "Traverse.h"
16
#include "Debug.h"
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

namespace AMDiS {

  void MeshManipulation::deleteDoubleDofs(std::set<MacroElement*>& newMacroEl)
  {
    std::map<int, MacroElement*> leafInMacroEl;
    
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL);
    while (elInfo) {
      leafInMacroEl[elInfo->getElement()->getIndex()] = elInfo->getMacroElement();
      elInfo = stack.traverseNext(elInfo);
    }

    deleteDoubleDofs(leafInMacroEl, newMacroEl, 0);
    deleteDoubleDofs(leafInMacroEl, newMacroEl, 1);
  }


  void MeshManipulation::deleteDoubleDofs(std::map<int, MacroElement*>& leafInMacroEl,
					  std::set<MacroElement*>& newMacroEl,
					  int mode)
  {
    FUNCNAME("MeshManipulation::deleteDoubleDofs()");

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

    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH);
    while (elInfo) {
      Element *el = elInfo->getElement();

      for (int i = 0; i < mesh->getGeo(NEIGH); i++) {
51

52
53
54
	Element *neigh = elInfo->getNeighbour(i);

	if (!neigh)
55
56
57
58
59
60
61
	  continue;

	TEST_EXIT_DBG(leafInMacroEl.count(el->getIndex()) == 1)
	  ("Should not happen!\n");
	TEST_EXIT_DBG(leafInMacroEl.count(neigh->getIndex()) == 1)
	  ("No macro element for el %d, that is %d-ith neigbour of element %d!\n",
	   neigh->getIndex(), i, el->getIndex());
62
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

	int elMacroIndex = leafInMacroEl[el->getIndex()]->getIndex();
	int neighMacroIndex = leafInMacroEl[neigh->getIndex()]->getIndex();

	if (elMacroIndex == neighMacroIndex)
	  continue;

	if ((mode == 0 && 
	     newMacroEl.count(leafInMacroEl[el->getIndex()]) == 1 &&
	     newMacroEl.count(leafInMacroEl[neigh->getIndex()]) == 1 &&
	     elMacroIndex > neighMacroIndex) ||
	    (mode == 1 &&
	     newMacroEl.count(leafInMacroEl[el->getIndex()]) == 0 &&
	     newMacroEl.count(leafInMacroEl[neigh->getIndex()]) == 1)) {

	  if (el->getEdge(i) != neigh->getEdge(elInfo->getOppVertex(i))) {
	    std::vector<int> elEdgeIndex(2);
	    elEdgeIndex[0] = el->getVertexOfEdge(i, 0);
	    elEdgeIndex[1] = el->getVertexOfEdge(i, 1);

	    std::vector<int> neighEdgeIndex(2);
	    neighEdgeIndex[0] = neigh->getVertexOfEdge(elInfo->getOppVertex(i), 0);
	    neighEdgeIndex[1] = neigh->getVertexOfEdge(elInfo->getOppVertex(i), 1);

	    std::vector<DegreeOfFreedom> elEdgeDof(2);
	    elEdgeDof[0] = el->getDof(elEdgeIndex[0], 0);
	    elEdgeDof[1] = el->getDof(elEdgeIndex[1], 0);

	    std::vector<DegreeOfFreedom> neighEdgeDof(2);
	    neighEdgeDof[0] = neigh->getDof(neighEdgeIndex[0], 0);
	    neighEdgeDof[1] = neigh->getDof(neighEdgeIndex[1], 0);

	    if ((elEdgeDof[0] < elEdgeDof[1] && neighEdgeDof[0] > neighEdgeDof[1]) ||
		(elEdgeDof[0] > elEdgeDof[1] && neighEdgeDof[0] < neighEdgeDof[1])) {
	      std::swap(neighEdgeIndex[0], neighEdgeIndex[1]);
	      std::swap(neighEdgeDof[0], neighEdgeDof[1]);
	    }

100

101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
	    for (int j = 0; j < 2; j++) {

	      if (neighEdgeDof[j] != elEdgeDof[j]) {
		const DegreeOfFreedom *delDof = neigh->getDof(neighEdgeIndex[j]);
		const DegreeOfFreedom *replaceDof = el->getDof(elEdgeIndex[j]);

		if (mapDelDofs.count(delDof) == 1 && mapDelDofs[delDof] != replaceDof) {
		  if (mapDofsMacro[mapDelDofs[delDof]] > elMacroIndex) {
		    mapDelDofs[replaceDof] = mapDelDofs[delDof];
		  } else {
		    mapDelDofs[mapDelDofs[delDof]] = replaceDof;
		    mapDelDofs[delDof] = replaceDof;
		    mapDofsMacro[replaceDof] = elMacroIndex;
		  }
		} else {
		  mapDelDofs[delDof] = replaceDof;
		  mapDofsMacro[replaceDof] = elMacroIndex;
		}
	      }
120
	    }	     
121
122
123
124
	  }
	}
      }

125

126
127
128
      elInfo = stack.traverseNext(elInfo);
    }

129

130
131
132
133
134
135
136
137
138
139
140
141
142
143
    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()) {
	  it->second = findIt->second;   
	  changed = true;
	}
      } 
    } while (changed);


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

      elInfo = stack.traverseNext(elInfo);
    }

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

}