PeriodicBC.cc 6.69 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
16
17
18
19
20
21
22
23
#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 {

24
  std::vector<PeriodicDOFMapping*> PeriodicDOFMapping::mappings;
25

26

27
  PeriodicDOFMapping* 
28
  PeriodicDOFMapping::providePeriodicDOFMapping(const BasisFunction *fcts)
29
  {
30
31
32
    for (std::vector<PeriodicDOFMapping*>::iterator it = mappings.begin(); 
	 it != mappings.end(); ++it)
      if ((*it)->basFcts == fcts)
33
	return *it;
34

35
36
    PeriodicDOFMapping *newMapping = new PeriodicDOFMapping(fcts);
    mappings.push_back(newMapping);
37
38
39
40
    return newMapping;
  }


41
42
  PeriodicDOFMapping::PeriodicDOFMapping(const BasisFunction *fcts) 
    : basFcts(fcts) 
43
44
  {
    FUNCNAME("PeriodicDOFMapping::PeriodicDOFMapping()");
45
46
    TEST_EXIT(basFcts->getDim() > 1)("dim == 1\n");  
    int nBasFcts = basFcts->getNumber();
47
    DimVec<double> *lambda;
48
49
50
    for (int i = 0; i < nBasFcts; i++) {
      lambda = basFcts->getCoords(i);
      indexOfCoords[*lambda] = i;
51
52
53
    }
  }

54

55
56
  PeriodicDOFMapping::~PeriodicDOFMapping()
  {
57
    std::map<DimVec<int>, DegreeOfFreedom*, DimVecLess<int> >::iterator it;
58
    for (it = dofPermutation.begin(); it != dofPermutation.end(); ++it)
59
60
      if (it->second)
	delete [] it->second;
61
62
  }

63

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

68
69
70
    if (dofPermutation[vertexPermutation] == NULL) {
      int dim = basFcts->getDim();
      int nBasFcts = basFcts->getNumber();
71
      int sum = 0;
72
      for (int i = 0; i < dim + 1; i++) {
73
	sum += i - vertexPermutation[i];
74
	TEST_EXIT(vertexPermutation[i] < dim + 1)
75
76
	  ("invalid vertexPermuation\n");
      }
77
      TEST_EXIT(sum == 0)("invalid vertexPermutation\n");
78
79
80
81
82

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

83
      DegreeOfFreedom *mapping = new DegreeOfFreedom[nBasFcts]; 
84

85
86
      for (int i = 0; i < nBasFcts; i++) {
	lambda = basFcts->getCoords(i);
87
	for (int j = 0; j < dim + 1; j++) 
88
	  newLambda[vertexPermutation[j]] = (*lambda)[j];
89

90
	mapping[i] = indexOfCoords[newLambda];
91
92
      }

93
      dofPermutation[vertexPermutation] = mapping;
94
95
    }

96
    return dofPermutation[vertexPermutation];
97
98
99
  }


100
101
102
  PeriodicBC::PeriodicBC(BoundaryType type, FiniteElemSpace *rowSpace) 
    : BoundaryCondition(type, rowSpace, NULL),
      masterMatrix(NULL)
103
  {
104
    if (rowFeSpace->getMesh()->getDim() > 1)
105
      periodicDOFMapping = 
106
	PeriodicDOFMapping::providePeriodicDOFMapping(rowFeSpace->getBasisFcts());
107
108
    else
      periodicDOFMapping = NULL;    
109
110
  }

111

112
113
114
  PeriodicBC::~PeriodicBC()
  {}

115

116
117
118
119
  void PeriodicBC::initMatrix(DOFMatrix* matrix)
  {
    FUNCNAME("PeriodicBC::initMatrix()");

120
121
    if (!masterMatrix) {
      masterMatrix = matrix;
122
      Mesh *mesh = matrix->getRowFeSpace()->getMesh();
123
      associated = mesh->getPeriodicAssociations()[boundaryType];
124
      
Thomas Witkowski's avatar
Thomas Witkowski committed
125
126
      TEST_EXIT(associated)
	("No associations for periodic boundary condition %d!\n", boundaryType);
127
128
129
    }
  }

130

131
132
  void PeriodicBC::fillBoundaryCondition(DOFMatrix *matrix,
					 ElInfo *elInfo,
133
					 const DegreeOfFreedom *dofIndices,
134
135
					 const BoundaryType *localBound,
					 int nBasFcts)
136
  {
137
138
    FUNCNAME("PeriodicBC::fillBoundaryCondition()");

Thomas Witkowski's avatar
Thomas Witkowski committed
139
140
    if (matrix != masterMatrix)
      return;
141
  
Thomas Witkowski's avatar
Thomas Witkowski committed
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
    int dim = rowFeSpace->getMesh()->getDim();
    if (dim == 1)
      return;
    
    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);
    GeoIndex sideGeoIndex = INDEX_OF_DIM(dim - 1, dim);
    std::vector<DegreeOfFreedom> neighIndices(num);
    
    for (int side = 0; side <= dim; side++) {      
      if (elInfo->getBoundary(sideGeoIndex, side) == boundaryType) {	
	for (int vertex = 0; vertex < dim; vertex++) {
	  int index = element->getVertexOfPosition(sideGeoIndex, side, vertex);
	  periodicDOFs[vertex] = (*associated)[dofIndices[index]];
	}
	
	Element *neigh = elInfo->getNeighbour(side);
	TEST_EXIT_DBG(neigh)("Wrong neighbour information at side %d!\n", side);
	
	basFcts->getLocalIndices(neigh, admin, neighIndices);
168
	
Thomas Witkowski's avatar
Thomas Witkowski committed
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
	int oppVertex = 0;
	for (int i = 0; i < dim + 1; i++) {
	  // get vertex permutation
	  if (i == side) {
	    vertexPermutation[i] = 0;
	  } else {
	    DegreeOfFreedom periodicDOF = 
	      periodicDOFs[element->getPositionOfVertex(side, i)];
	    
	    int j = 0;
	    for (; j < dim + 1; j++)
	      if (neigh->getDof(j, 0) == periodicDOF) 
		break;
	    
	    vertexPermutation[i] = j;
184
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
185
	  oppVertex += i - vertexPermutation[i];
186
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
187
188
189
190
191
192
193
194
195
196
	vertexPermutation[side] = oppVertex;
	
	// get DOF permutation
	const DegreeOfFreedom *dofPermutation = 
	  periodicDOFMapping->getDOFPermutation(vertexPermutation);
	
	// set associated dofs
	for (int i = 0; i < num; i++)
	  if ((*(basFcts->getCoords(i)))[side] == 0)
	    (*associated)[dofIndices[i]] = neighIndices[dofPermutation[i]];
197
198
199
200
      }
    }
  }

201

202
203
204
205
  void PeriodicBC::exitMatrix(DOFMatrix* matrix)
  {
    FUNCNAME("PeriodicBC::exitMatrix()");

Thomas Witkowski's avatar
Thomas Witkowski committed
206
207
208
    TEST_EXIT(matrix)("No matrix\n");

    TEST_EXIT(associated)("No associated vector!\n");
209

210
211
    if (matrix == masterMatrix)
      masterMatrix = NULL;    
212
213
214

    using namespace mtl;

215
    DOFAdmin* admin = rowFeSpace->getAdmin();
Thomas Witkowski's avatar
Thomas Witkowski committed
216
    std::vector<int> dofMap(admin->getUsedSize());
217
    for (int i = 0; i < admin->getUsedSize(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
218
      dofMap[i] = (*associated)[i];
219

220
    // Compute reorder matrix (newRow and newCol yields transposed!!!)
Thomas Witkowski's avatar
Thomas Witkowski committed
221
    matrix::traits::reorder<>::type R= matrix::reorder(dofMap);
222
223
224
225
    DOFMatrix::base_matrix_type &A= matrix->getBaseMatrix(), C;

    C = R * A * trans(R) + A;
    A = 0.5 * C;
226
227
  }

228

229
230
231
232
233
  void PeriodicBC::exitVector(DOFVectorBase<double>* vector)
  {
    FUNCNAME("PeriodicBC::exitVector()");

    DOFIterator<double> vecIt(vector, USED_DOFS);
234
    Mesh *mesh = vector->getFeSpace()->getMesh();
235
    VertexVector *associated = mesh->getPeriodicAssociations()[boundaryType];
236

237
    for (vecIt.reset(); !vecIt.end(); ++vecIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
238
239
      DegreeOfFreedom dof = vecIt.getDOFIndex();
      DegreeOfFreedom newDOF = (*associated)[dof];
240

241
      if (dof < newDOF) {
Thomas Witkowski's avatar
Thomas Witkowski committed
242
	double entry = ((*vector)[dof] + (*vector)[newDOF]) * 0.5;
243
244
245
246
247
248
249
	(*vector)[dof] = entry;
	(*vector)[newDOF] = entry;
      }
    }
  }

}