PeriodicBC.cc 7.46 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include "PeriodicBC.h"
#include "BasisFunction.h"
#include "DOFVector.h"
#include "LeafData.h"
#include "DOFMatrix.h"
#include "Traverse.h"
#include "Boundary.h"
#include "VertexVector.h"

namespace AMDiS {

  ::std::vector<PeriodicDOFMapping*> PeriodicDOFMapping::mappings_;

  PeriodicDOFMapping* 
  PeriodicDOFMapping::providePeriodicDOFMapping(const BasisFunction *basFcts)
  {
    ::std::vector<PeriodicDOFMapping*>::iterator it;
    ::std::vector<PeriodicDOFMapping*>::iterator end = mappings_.end();
19
20
    for (it = mappings_.begin(); it != end; ++it) {
      if ((*it)->basFcts_ == basFcts) {
21
22
23
24
25
26
27
28
29
30
31
32
33
34
	return *it;
      }
    }
    PeriodicDOFMapping *newMapping = NEW PeriodicDOFMapping(basFcts);
    mappings_.push_back(newMapping);
    return newMapping;
  }


  PeriodicDOFMapping::PeriodicDOFMapping(const BasisFunction *basFcts) 
    : basFcts_(basFcts) 
  {
    FUNCNAME("PeriodicDOFMapping::PeriodicDOFMapping()");
    TEST_EXIT(basFcts_->getDim() > 1)("dim == 1\n");  
35
    int num = basFcts_->getNumber();
36
    DimVec<double> *lambda;
37
    for (int i = 0; i < num; i++) {
38
39
40
41
42
43
44
45
      lambda = basFcts_->getCoords(i);
      indexOfCoords_[*lambda] = i;
    }
  }

  PeriodicDOFMapping::~PeriodicDOFMapping()
  {
    ::std::map<DimVec<int>, DegreeOfFreedom*, DimVecLess<int> >::iterator it;
46
47
    for (it = dofPermutation_.begin(); it != dofPermutation_.end(); ++it) {
      if (it->second) {
48
49
50
51
52
53
54
55
56
57
	FREE_MEMORY(it->second, DegreeOfFreedom, basFcts_->getNumber());
      }
    }
  }

  const DegreeOfFreedom 
  *PeriodicDOFMapping::getDOFPermutation(const DimVec<int> &vertexPermutation)
  {
    FUNCNAME("PeriodicDOFMapping::getDOFPermutation()");

58
    if (dofPermutation_[vertexPermutation] == NULL) {
59
60
61
62
      int dim = basFcts_->getDim();
      int num = basFcts_->getNumber();

      int sum = 0;
63
      for (int i = 0; i < dim + 1; i++) {
64
	sum += i - vertexPermutation[i];
65
	TEST_EXIT_DBG(vertexPermutation[i] < dim + 1)
66
67
	  ("invalid vertexPermuation\n");
      }
68
      TEST_EXIT_DBG(sum == 0)("invalid vertexPermutation\n");
69
70
71
72
73
74
75

      // create dof permutation
      DimVec<double> *lambda;
      DimVec<double> newLambda(dim, NO_INIT);

      DegreeOfFreedom *mapping = GET_MEMORY(DegreeOfFreedom, num); 

76
      for (int i = 0; i < num; i++) {
77
	lambda = basFcts_->getCoords(i);
78
	for (int j = 0; j < dim + 1; j++) {
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
	  newLambda[vertexPermutation[j]] = (*lambda)[j];
	}
	mapping[i] = indexOfCoords_[newLambda];
      }

      dofPermutation_[vertexPermutation] = mapping;
    }

    return dofPermutation_[vertexPermutation];
  }


  PeriodicBC::PeriodicBC(BoundaryType type, FiniteElemSpace *rowFESpace_) 
    : BoundaryCondition(type, rowFESpace_, NULL),
      masterMatrix_(NULL)
  {
95
    if (rowFESpace->getMesh()->getDim() > 1) {
96
97
98
99
100
101
102
103
104
105
106
107
108
109
      periodicDOFMapping_ = 
	PeriodicDOFMapping::providePeriodicDOFMapping(rowFESpace_->getBasisFcts());
    } else {
      periodicDOFMapping_ = NULL;
    }
  }

  PeriodicBC::~PeriodicBC()
  {}

  void PeriodicBC::initMatrix(DOFMatrix* matrix)
  {
    FUNCNAME("PeriodicBC::initMatrix()");

110
    if (!masterMatrix_) {
111
112
113
114
115
116
      masterMatrix_ = matrix;

      Mesh *mesh = matrix->getRowFESpace()->getMesh();

      associated_ = mesh->getPeriodicAssociations()[boundaryType];
      
117
118
      TEST_EXIT_DBG(associated_)("no associations for periodic boundary condition %d\n",
				 boundaryType);
119
120
121
122
123
124
125
126
127
128
129
130
131
132

      const BasisFunction *basFcts = rowFESpace->getBasisFcts();
      int num = basFcts->getNumber();

      neighIndices_ = GET_MEMORY(DegreeOfFreedom, num);
    }
  }

  void PeriodicBC::fillBoundaryCondition(DOFMatrix             *matrix,
					 ElInfo                *elInfo,
					 const DegreeOfFreedom *dofIndices,
					 const BoundaryType    *localBound,
					 int                    nBasFcts)
  {
133
    if (matrix == masterMatrix_) {
134
135
  
      int dim = rowFESpace->getMesh()->getDim();
136
      if (dim > 1) {
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
	DOFAdmin *admin = rowFESpace->getAdmin();

	FixVec<int, WORLD> elFace(dim, NO_INIT);
	FixVec<int, WORLD> neighFace(dim, NO_INIT);
	DimVec<int> vertexPermutation(dim, NO_INIT);

	const BasisFunction *basFcts = rowFESpace->getBasisFcts();
	int num = basFcts->getNumber();

	Element *element = elInfo->getElement();

	DimVec<DegreeOfFreedom> periodicDOFs(dim-1, NO_INIT);

	int vertex, index, side;

	GeoIndex sideGeoIndex = INDEX_OF_DIM(dim-1, dim);

154
	for (side = 0; side < dim + 1; side++) {
155

156
	  if (elInfo->getBoundary(sideGeoIndex, side) == boundaryType) {
157

158
	    for (vertex = 0; vertex < dim; vertex++) {
159
160
161
162
163
164
165
166
167
168
169
	      index = element->getVertexOfPosition(sideGeoIndex,
						   side,
						   vertex);
	      periodicDOFs[vertex] = (*associated_)[dofIndices[index]];
	    }

	    Element *neigh = elInfo->getNeighbour(side);

	    basFcts->getLocalIndices(neigh, admin, neighIndices_);

	    int oppVertex = 0;
170
	    for (int i = 0; i < dim + 1; i++) {
171
	      // get vertex permutation
172
	      if (i == side) {
173
174
175
176
177
		vertexPermutation[i] = 0;
	      } else {
		DegreeOfFreedom periodicDOF = 
		  periodicDOFs[element->getPositionOfVertex(side, i)];

178
179
180
181
		int j;
		for (j = 0; j < dim + 1; j++) {
		  if (neigh->getDOF(j, 0) == periodicDOF) 
		    break;
182
183
184
185
186
187
188
189
190
191
192
193
		}
		vertexPermutation[i] = j;
	      }
	      oppVertex += i - vertexPermutation[i];
	    }
	    vertexPermutation[side] = oppVertex;

	    // get DOF permutation
	    const DegreeOfFreedom *dofPermutation = 
	      periodicDOFMapping_->getDOFPermutation(vertexPermutation);
	
	    // set associated dofs
194
195
	    for (int i = 0; i < num; i++) {
	      if ((*(basFcts->getCoords(i)))[side] == 0) {
196
197
198
199
200
201
202
203
204
205
206
207
208
		(*associated_)[dofIndices[i]] = neighIndices_[dofPermutation[i]];
	      }
	    }
	  }
	}
      }
    }
  }

  void PeriodicBC::exitMatrix(DOFMatrix* matrix)
  {
    FUNCNAME("PeriodicBC::exitMatrix()");

209
    TEST_EXIT_DBG(matrix)("no matrix\n");
210

211
    if (matrix == masterMatrix_) {
212
213
      FREE_MEMORY(neighIndices_, DegreeOfFreedom, 
		  rowFESpace->getBasisFcts()->getNumber());
214
215
216
217
218
219
220
221
222
      masterMatrix_ = NULL;
    }

    ::std::vector< ::std::vector<MatEntry> >::iterator rowIt;
    ::std::vector< ::std::vector<MatEntry> >::iterator rowEnd = matrix->end();
    int colIndex, rowSize;
    int row, col, newRow, newCol;
    double entry, *newEntryPtr;

223
    for (rowIt = matrix->begin(), row = 0; rowIt != rowEnd; ++rowIt, ++row) {
224
225
      rowSize = static_cast<int>(rowIt->size());
      newRow = (*associated_)[row];
226
      for (colIndex = 0; colIndex < rowSize; colIndex++) {
227
	col = (*rowIt)[colIndex].col;
228
229
230
231
	if (col == DOFMatrix::UNUSED_ENTRY) 
	  continue;
	if (col == DOFMatrix::NO_MORE_ENTRIES) 
	  break;
232
233
234
235
	newCol = (*associated_)[col];

	newEntryPtr = matrix->hasSparseDOFEntry(newRow, newCol);

236
237
238
	if ((row < newRow) || 
	    ((row == newRow) && (col < newCol)) ||
	    !newEntryPtr) 
239
	  {
240
	    if (!newEntryPtr) {
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
	      entry = 0.5 * (*rowIt)[colIndex].entry;
	      (*rowIt)[colIndex].entry = entry;
	      matrix->addSparseDOFEntry(1.0, newRow, newCol, entry, false);
	    } else {
	      entry = 0.5 * ((*rowIt)[colIndex].entry + *newEntryPtr);
	      (*rowIt)[colIndex].entry = *newEntryPtr = entry;
	    }
	  }
      }
    }
  }

  void PeriodicBC::exitVector(DOFVectorBase<double>* vector)
  {
    FUNCNAME("PeriodicBC::exitVector()");

    double entry;

    DegreeOfFreedom dof;
    DegreeOfFreedom newDOF;

    DOFIterator<double> vecIt(vector, USED_DOFS);

    Mesh *mesh = vector->getFESpace()->getMesh();
    VertexVector *associated_ = mesh->getPeriodicAssociations()[boundaryType];

267
    for (vecIt.reset(); !vecIt.end(); ++vecIt) {
268
269
270
271
      dof = vecIt.getDOFIndex();

      newDOF = (*associated_)[dof];

272
      if (dof < newDOF) {
273
274
275
276
277
278
279
280
	entry = ((*vector)[dof] + (*vector)[newDOF]) * 0.5;
	(*vector)[dof] = entry;
	(*vector)[newDOF] = entry;
      }
    }
  }

}