PeriodicBC.cc 8.28 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
20
21


22
23
24
25
26
27
28
29
30
31
32
#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 {

33
  std::vector<PeriodicDOFMapping*> PeriodicDOFMapping::mappings;
34

35

36
  PeriodicDOFMapping*
37
  PeriodicDOFMapping::providePeriodicDOFMapping(const BasisFunction *fcts)
38
  {
39
    for (std::vector<PeriodicDOFMapping*>::iterator it = mappings.begin();
40
41
	 it != mappings.end(); ++it)
      if ((*it)->basFcts == fcts)
42
	return *it;
43

44
45
    PeriodicDOFMapping *newMapping = new PeriodicDOFMapping(fcts);
    mappings.push_back(newMapping);
46
47
48
49
    return newMapping;
  }


50
51
  PeriodicDOFMapping::PeriodicDOFMapping(const BasisFunction *fcts)
    : basFcts(fcts)
52
53
  {
    FUNCNAME("PeriodicDOFMapping::PeriodicDOFMapping()");
54
    TEST_EXIT(basFcts->getDim() > 1)("dim == 1\n");
55
    int nBasFcts = basFcts->getNumber();
56
    DimVec<double> *lambda;
57
58
59
    for (int i = 0; i < nBasFcts; i++) {
      lambda = basFcts->getCoords(i);
      indexOfCoords[*lambda] = i;
60
61
62
    }
  }

63

64
65
  PeriodicDOFMapping::~PeriodicDOFMapping()
  {
66
    std::map<DimVec<int>, DegreeOfFreedom*, DimVecLess<int> >::iterator it;
67
    for (it = dofPermutation.begin(); it != dofPermutation.end(); ++it)
68
69
      if (it->second)
	delete [] it->second;
70
71
  }

72

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

77
    if (dofPermutation[vertexPermutation] == NULL) {
78
79
      int dim = basFcts->getDim();
      int nBasFcts = basFcts->getNumber();
80
      int sum = 0;
81
      for (int i = 0; i < dim + 1; i++) {
82
	sum += i - vertexPermutation[i];
83
	TEST_EXIT(vertexPermutation[i] < dim + 1)
84
85
	  ("invalid vertexPermuation\n");
      }
86
      TEST_EXIT(sum == 0)("invalid vertexPermutation\n");
87
88
89
90
91

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

92
      DegreeOfFreedom *mapping = new DegreeOfFreedom[nBasFcts];
93

94
95
      for (int i = 0; i < nBasFcts; i++) {
	lambda = basFcts->getCoords(i);
96
	for (int j = 0; j < dim + 1; j++)
97
	  newLambda[vertexPermutation[j]] = (*lambda)[j];
98

99
	mapping[i] = indexOfCoords[newLambda];
100
101
      }

102
      dofPermutation[vertexPermutation] = mapping;
103
104
    }

105
    return dofPermutation[vertexPermutation];
106
107
108
  }


109
  PeriodicBC::PeriodicBC(BoundaryType type, const FiniteElemSpace *rowSpace, bool diagonal)
110
    : BoundaryCondition(type, rowSpace, NULL),
111
112
      masterMatrix(NULL),
      isDiagonal(diagonal)
113
  {
114
    if (rowFeSpace->getMesh()->getDim() > 1)
115
      periodicDOFMapping =
116
        PeriodicDOFMapping::providePeriodicDOFMapping(rowFeSpace->getBasisFcts());
117
    else
118
      periodicDOFMapping = NULL;
119
120
  }

121

122
123
124
  PeriodicBC::~PeriodicBC()
  {}

125

126
127
128
129
  void PeriodicBC::initMatrix(DOFMatrix* matrix)
  {
    FUNCNAME("PeriodicBC::initMatrix()");

130
131
    if (!masterMatrix) {
      masterMatrix = matrix;
132
      Mesh *mesh = matrix->getRowFeSpace()->getMesh();
133
      associated = &(mesh->getPeriodicAssociations(boundaryType, rowFeSpace->getAdmin()));
134

Thomas Witkowski's avatar
Thomas Witkowski committed
135
      TEST_EXIT(associated)
136
        ("No associations for periodic boundary condition %d!\n", boundaryType);
137
138
139
    }
  }

140

141
142
  void PeriodicBC::fillBoundaryCondition(DOFMatrix *matrix,
					 ElInfo *elInfo,
143
					 const DegreeOfFreedom *dofIndices,
144
145
					 const BoundaryType *localBound,
					 int nBasFcts)
146
  {
147
    FUNCNAME_DBG("PeriodicBC::fillBoundaryCondition()");
148

Thomas Witkowski's avatar
Thomas Witkowski committed
149
150
    if (matrix != masterMatrix)
      return;
151

Thomas Witkowski's avatar
Thomas Witkowski committed
152
153
154
    int dim = rowFeSpace->getMesh()->getDim();
    if (dim == 1)
      return;
155

Thomas Witkowski's avatar
Thomas Witkowski committed
156
    DOFAdmin *admin = rowFeSpace->getAdmin();
157
    int iadmin = rowFeSpace->getMesh()->getAdminIndex(admin);
Thomas Witkowski's avatar
Thomas Witkowski committed
158
159
160
161
162
163
164
165
166
    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);
167
168
169

    for (int side = 0; side <= dim; side++) {
      if (elInfo->getBoundary(sideGeoIndex, side) == boundaryType) {
Thomas Witkowski's avatar
Thomas Witkowski committed
170
171
172
173
	for (int vertex = 0; vertex < dim; vertex++) {
	  int index = element->getVertexOfPosition(sideGeoIndex, side, vertex);
	  periodicDOFs[vertex] = (*associated)[dofIndices[index]];
	}
174

Thomas Witkowski's avatar
Thomas Witkowski committed
175
176
	Element *neigh = elInfo->getNeighbour(side);
	TEST_EXIT_DBG(neigh)("Wrong neighbour information at side %d!\n", side);
177

Thomas Witkowski's avatar
Thomas Witkowski committed
178
	basFcts->getLocalIndices(neigh, admin, neighIndices);
179

Thomas Witkowski's avatar
Thomas Witkowski committed
180
181
182
183
184
185
	int oppVertex = 0;
	for (int i = 0; i < dim + 1; i++) {
	  // get vertex permutation
	  if (i == side) {
	    vertexPermutation[i] = 0;
	  } else {
186
	    DegreeOfFreedom periodicDOF =
Thomas Witkowski's avatar
Thomas Witkowski committed
187
	      periodicDOFs[element->getPositionOfVertex(side, i)];
188

Thomas Witkowski's avatar
Thomas Witkowski committed
189
190
	    int j = 0;
	    for (; j < dim + 1; j++)
191
          if (neigh->getDof(j, iadmin) == periodicDOF)
Thomas Witkowski's avatar
Thomas Witkowski committed
192
		break;
193

Thomas Witkowski's avatar
Thomas Witkowski committed
194
	    vertexPermutation[i] = j;
195
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
196
	  oppVertex += i - vertexPermutation[i];
197
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
198
	vertexPermutation[side] = oppVertex;
199

Thomas Witkowski's avatar
Thomas Witkowski committed
200
	// get DOF permutation
201
	const DegreeOfFreedom *dofPermutation =
Thomas Witkowski's avatar
Thomas Witkowski committed
202
	  periodicDOFMapping->getDOFPermutation(vertexPermutation);
203

Thomas Witkowski's avatar
Thomas Witkowski committed
204
205
206
207
	// set associated dofs
	for (int i = 0; i < num; i++)
	  if ((*(basFcts->getCoords(i)))[side] == 0)
	    (*associated)[dofIndices[i]] = neighIndices[dofPermutation[i]];
208
209
210
211
      }
    }
  }

212

213
214
215
216
  void PeriodicBC::exitMatrix(DOFMatrix* matrix)
  {
    FUNCNAME("PeriodicBC::exitMatrix()");

Thomas Witkowski's avatar
Thomas Witkowski committed
217
218
    TEST_EXIT(matrix)("No matrix\n");
    TEST_EXIT(associated)("No associated vector!\n");
219

220
    if (matrix == masterMatrix)
221
      masterMatrix = NULL;
222
223

    using namespace mtl;
224
225
226
    typedef DOFMatrix::base_matrix_type Matrix;
    typedef mtl::traits::range_generator<mtl::tag::row, Matrix>::type c_type;
    typedef mtl::traits::range_generator<mtl::tag::nz, c_type>::type  ic_type;
227

228
    DOFAdmin* admin = rowFeSpace->getAdmin();
229
230
231
232
233
234
235
236
    int adminSize = admin->getUsedSize();
    Matrix &A= matrix->getBaseMatrix();

    mtl::traits::col<Matrix>::type c(A);
    mtl::traits::value<Matrix>::type v(A);

    std::vector<std::vector<std::pair<DegreeOfFreedom, double> > > row_values;
    row_values.resize(adminSize);
237

238
239
240
    for (DegreeOfFreedom i = 0; i < adminSize; i++) {
      if (i < (*associated)[i]) {
        c_type cursor(begin<mtl::tag::row>(A) + i);
241

242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
        for (ic_type icursor(begin<mtl::tag::nz>(cursor)), icend(end<mtl::tag::nz>(cursor));
            icursor != icend; ++icursor) {
          row_values[i].push_back(std::make_pair(c(*icursor), v(*icursor)));
          v(*icursor, 0.0);
        }
      }
    }
    matrix::inserter<Matrix, update_plus<double> > ins(A);
    if (isDiagonal) {
      for (DegreeOfFreedom i = 0; i < adminSize; i++) {
        if (i < (*associated)[i]) {
          ins[i][i] << 1.0;
          ins[i][(*associated)[i]] << -1.0;
        }
      }
    }
    for (DegreeOfFreedom i = 0; i < adminSize; i++) {
      if (i < (*associated)[i]) {
        for (size_t j = 0; j < row_values[i].size(); j++)
          ins[(*associated)[i]][row_values[i][j].first] << row_values[i][j].second;
      }
    }
264
265
  }

266

267
268
  void PeriodicBC::exitVector(DOFVectorBase<double>* vector)
  {
269
270
271
    FUNCNAME("PeriodicBC::exitVector()");
    TEST_EXIT(rowFeSpace == vector->getFeSpace())("Should not happen.\n");

272
    DOFIterator<double> vecIt(vector, USED_DOFS);
273
    Mesh *mesh = vector->getFeSpace()->getMesh();
274
    VertexVector *associated = &(mesh->getPeriodicAssociations(boundaryType, rowFeSpace->getAdmin()));
275

276
    for (vecIt.reset(); !vecIt.end(); ++vecIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
277
278
      DegreeOfFreedom dof = vecIt.getDOFIndex();
      DegreeOfFreedom newDOF = (*associated)[dof];
279

280
      if (dof < newDOF) {
281
282
        (*vector)[newDOF] += (*vector)[dof];
        (*vector)[dof] = 0.0;
283
284
285
286
287
      }
    }
  }

}