Assembler.cc 11.7 KB
Newer Older
1
2
3
#include <vector>
#include <algorithm>
#include <boost/numeric/mtl/mtl.hpp>
4
5
6
7
8
#include "Assembler.h"
#include "Operator.h"
#include "Element.h"
#include "QPsiPhi.h"
#include "DOFVector.h"
9
#include "OpenMP.h"
10
11
12

namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
13
14
15
  Assembler::Assembler(Operator *op,
		       const FiniteElemSpace *row,
		       const FiniteElemSpace *col) 
16
    : operat(op),
Thomas Witkowski's avatar
Thomas Witkowski committed
17
18
      rowFESpace(row),
      colFESpace(col ? col : row),
19
20
21
22
23
      nRow(rowFESpace->getBasisFcts()->getNumber()),
      nCol(colFESpace->getBasisFcts()->getNumber()),
      remember(true),
      rememberElMat(false),
      rememberElVec(false),
24
25
      elementMatrix(nRow, nCol),
      elementVector(nRow),
26
27
28
29
      lastMatEl(NULL),
      lastVecEl(NULL),
      lastTraverseId(-1)
  
30
  {}
Thomas Witkowski's avatar
Thomas Witkowski committed
31
32

  Assembler::~Assembler()
33
  {}
34
35

  void Assembler::calculateElementMatrix(const ElInfo *elInfo, 
36
					 ElementMatrix& userMat,
37
38
39
40
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

41
    if (remember && (factor != 1.0 || operat->uhOld))
42
      rememberElMat = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
43

44
    Element *el = elInfo->getElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
45

46
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
47
48
49
      initElement(elInfo);

    if (el != lastMatEl || !operat->isOptimized()) {
50
51
52
      if (rememberElMat)
	set_to_zero(elementMatrix);

53
54
55
      lastMatEl = el;
    } else {
      if (rememberElMat) {
56
	userMat += factor * elementMatrix;
57
58
59
	return;
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
60
 
61
    ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
62
63
64
65
66
67
68
69
70
71

    if (secondOrderAssembler)
      secondOrderAssembler->calculateElementMatrix(elInfo, mat);
    if (firstOrderAssemblerGrdPsi)
      firstOrderAssemblerGrdPsi->calculateElementMatrix(elInfo, mat);
    if (firstOrderAssemblerGrdPhi)
      firstOrderAssemblerGrdPhi->calculateElementMatrix(elInfo, mat);
    if (zeroOrderAssembler)
      zeroOrderAssembler->calculateElementMatrix(elInfo, mat);

Thomas Witkowski's avatar
Thomas Witkowski committed
72
    if (rememberElMat && &userMat != &elementMatrix)
73
      userMat += factor * elementMatrix;
74
75
  }

76

77
78
  void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
					 const ElInfo *colElInfo,
79
80
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
81
					 ElementMatrix& userMat,
82
83
84
85
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

86
    if (remember && (factor != 1.0 || operat->uhOld))
87
88
      rememberElMat = true;
  
Thomas Witkowski's avatar
Thomas Witkowski committed
89
    Element *el = smallElInfo->getElement();   
90
    lastVecEl = lastMatEl = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
91
   
92
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
Thomas Witkowski's avatar
Thomas Witkowski committed
93
      initElement(smallElInfo, largeElInfo);
94
95

    if (el != lastMatEl || !operat->isOptimized()) {
96
97
98
      if (rememberElMat)
	set_to_zero(elementMatrix);

99
100
101
      lastMatEl = el;
    } else {
      if (rememberElMat) {
102
	userMat += factor * elementMatrix;
103
104
105
	return;
      }
    }
106
 
107
    ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
108

109
    if (secondOrderAssembler) {
110
      ERROR_EXIT("Da muss i noch ma testen!\n");
111
      secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
112
113
114

      ElementMatrix m(nRow, nRow);
      smallElInfo->getSubElemCoordsMat_so(m, rowFESpace->getBasisFcts()->getDegree());
115

116
117
      ElementMatrix tmpMat(nRow, nRow);
      tmpMat = m * mat;
118
      mat = tmpMat;      
119
120
121
    }

    if (firstOrderAssemblerGrdPsi) {
122
      ERROR_EXIT("Not yet implemented for first order assembler!\n");
123
      firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
124
125
126
    }

    if (firstOrderAssemblerGrdPhi) {
127
      ERROR_EXIT("Not yet implemented for first order assembler!\n");
128
      firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
129
    }
130

131
132
    if (zeroOrderAssembler) {
      zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
133

134
135
      ElementMatrix m(nRow, nRow);
      smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
136

137
      ElementMatrix tmpMat(nRow, nRow);
138
139
140
141

      if (smallElInfo == colElInfo) 
	tmpMat = m * mat;	
      else
142
143
	tmpMat = mat * trans(m);

144
      mat = tmpMat;     
145
    }
146

147
148
    if (rememberElMat && &userMat != &elementMatrix)
      userMat += factor * elementMatrix;   
149
150
  }

151

152
  void Assembler::calculateElementVector(const ElInfo *elInfo, 
153
					 ElementVector& userVec,
154
155
156
157
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

158
    if (remember && factor != 1.0)
159
160
161
162
      rememberElVec = true;

    Element *el = elInfo->getElement();

163
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
164
      initElement(elInfo);
165
    
Thomas Witkowski's avatar
Thomas Witkowski committed
166
    if (el != lastVecEl || !operat->isOptimized()) {
167
168
169
      if (rememberElVec)
	set_to_zero(elementVector);
	
170
171
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
172
      if (rememberElVec) {
173
	userVec += factor * elementVector;
174
175
176
	return;
      }
    }
177
178

    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
179
    if (operat->uhOld && remember) {
180
      matVecAssemble(elInfo, vec);
181
      if (rememberElVec)
182
	userVec += factor * elementVector;      
183

184
185
      return;
    } 
186
187

    if (firstOrderAssemblerGrdPsi)
188
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
189
    if (zeroOrderAssembler)
190
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
191
      
192
    if (rememberElVec)
193
      userVec += factor * elementVector;    
194
195
  }

196

Thomas Witkowski's avatar
Thomas Witkowski committed
197
198
199
200
  void Assembler::calculateElementVector(const ElInfo *mainElInfo, 
					 const ElInfo *auxElInfo,
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
201
					 ElementVector& userVec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
202
203
204
205
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

206
    if (remember && factor != 1.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
207
208
209
210
      rememberElVec = true;

    Element *el = mainElInfo->getElement();

211
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
212
213
      initElement(smallElInfo, largeElInfo);
   
Thomas Witkowski's avatar
Thomas Witkowski committed
214
    if (el != lastVecEl || !operat->isOptimized()) {
215
216
217
      if (rememberElVec)
	set_to_zero(elementVector);

Thomas Witkowski's avatar
Thomas Witkowski committed
218
219
220
      lastVecEl = el;
    } else {
      if (rememberElVec) {
221
	userVec += factor * elementVector;
Thomas Witkowski's avatar
Thomas Witkowski committed
222
223
224
	return;
      }
    }
225
    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
226
227

    if (operat->uhOld && remember) {
228
      if (smallElInfo->getLevel() == largeElInfo->getLevel())
Thomas Witkowski's avatar
Thomas Witkowski committed
229
	matVecAssemble(auxElInfo, vec);
230
231
      else
	matVecAssemble(mainElInfo, auxElInfo, smallElInfo, largeElInfo, vec);      
Thomas Witkowski's avatar
Thomas Witkowski committed
232

233
      if (rememberElVec)
234
	userVec += factor * elementVector;      
235

Thomas Witkowski's avatar
Thomas Witkowski committed
236
237
238
239
240
      return;
    } 

    ERROR_EXIT("Not yet implemented!\n");

241
242
243
244
245
246
247
248
#if 0
    if (firstOrderAssemblerGrdPsi)
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);    
    if (zeroOrderAssembler)
      zeroOrderAssembler->calculateElementVector(elInfo, vec);    
    if (rememberElVec)
      axpy(factor, *elementVector, *userVec);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
249
250
  }

251

252
  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector& vec)
253
254
255
256
  {
    FUNCNAME("Assembler::matVecAssemble()");

    Element *el = elInfo->getElement(); 
257
    double *uhOldLoc = new double[nRow];
258

259
260
    operat->uhOld->getLocalVector(el, uhOldLoc);
    
261
262
    if (el != lastMatEl) {
      set_to_zero(elementMatrix);
263
      calculateElementMatrix(elInfo, elementMatrix);
264
265
    }

266
    for (int i = 0; i < nRow; i++) {
267
      double val = 0.0;
268
      for (int j = 0; j < nRow; j++)
269
	val += elementMatrix[i][j] * uhOldLoc[j];
270
      
271
      vec[i] += val;
272
    }   
273

Thomas Witkowski's avatar
Thomas Witkowski committed
274
275
276
    delete [] uhOldLoc;
  }

277

Thomas Witkowski's avatar
Thomas Witkowski committed
278
279
  void Assembler::matVecAssemble(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
				 const ElInfo *smallElInfo, const ElInfo *largeElInfo,
280
				 ElementVector& vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
281
282
283
284
285
286
287
288
289
290
291
292
293
294
  {
    FUNCNAME("Assembler::matVecAssemble()");

    TEST_EXIT(rowFESpace->getBasisFcts() == colFESpace->getBasisFcts())
      ("Works only for equal basis functions for different components!\n");

    TEST_EXIT(operat->uhOld->getFESpace()->getMesh() == auxElInfo->getMesh())
      ("Da stimmt was nicht!\n");

    Element *mainEl = mainElInfo->getElement(); 
    Element *auxEl = auxElInfo->getElement();

    const BasisFunction *basFcts = rowFESpace->getBasisFcts();
    int nBasFcts = basFcts->getNumber();
295
    std::vector<double> uhOldLoc(nBasFcts);
Thomas Witkowski's avatar
Thomas Witkowski committed
296

297
    operat->uhOld->getLocalVector(auxEl, &(uhOldLoc[0]));
Thomas Witkowski's avatar
Thomas Witkowski committed
298
299

    if (mainEl != lastMatEl) {
300
      set_to_zero(elementMatrix);
301
      calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, 
302
 			     elementMatrix);    
Thomas Witkowski's avatar
Thomas Witkowski committed
303
    }
304

Thomas Witkowski's avatar
Thomas Witkowski committed
305
306
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
307
      for (int j = 0; j < nBasFcts; j++)
308
 	val += elementMatrix[i][j] * uhOldLoc[j];
309
      vec[i] += val;
310
    }   
311
312
  }

313

Thomas Witkowski's avatar
Thomas Witkowski committed
314
315
316
  void Assembler::initElement(const ElInfo *smallElInfo, 
			      const ElInfo *largeElInfo,
			      Quadrature *quad)
317
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
318
    if (secondOrderAssembler) 
Thomas Witkowski's avatar
Thomas Witkowski committed
319
      secondOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
320
    if (firstOrderAssemblerGrdPsi)
Thomas Witkowski's avatar
Thomas Witkowski committed
321
      firstOrderAssemblerGrdPsi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
322
    if (firstOrderAssemblerGrdPhi)
Thomas Witkowski's avatar
Thomas Witkowski committed
323
      firstOrderAssemblerGrdPhi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
324
    if (zeroOrderAssembler)
Thomas Witkowski's avatar
Thomas Witkowski committed
325
      zeroOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
326
327
  }

328

329
  void Assembler::checkQuadratures()
Thomas Witkowski's avatar
Thomas Witkowski committed
330
331
  { 
    if (secondOrderAssembler) {
332
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
333
334
      if (!secondOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
335
336
337
338
339
	int degree = operat->getQuadratureDegree(2);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	secondOrderAssembler->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
340
    if (firstOrderAssemblerGrdPsi) {
341
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
342
343
      if (!firstOrderAssemblerGrdPsi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
344
345
346
347
348
	int degree = operat->getQuadratureDegree(1, GRD_PSI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPsi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
349
    if (firstOrderAssemblerGrdPhi) {
350
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
      if (!firstOrderAssemblerGrdPhi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
353
354
355
356
357
	int degree = operat->getQuadratureDegree(1, GRD_PHI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPhi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
358
    if (zeroOrderAssembler) {
359
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
360
361
      if (!zeroOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
362
363
364
365
366
367
368
	int degree = operat->getQuadratureDegree(0);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	zeroOrderAssembler->setQuadrature(quadrature);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
369
370
371
372
373
  void Assembler::finishAssembling()
  {
    lastVecEl = NULL;
    lastMatEl = NULL;
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422

  OptimizedAssembler::OptimizedAssembler(Operator  *op,
					 Quadrature *quad2,
					 Quadrature *quad1GrdPsi,
					 Quadrature *quad1GrdPhi,
					 Quadrature *quad0,
					 const FiniteElemSpace *row,
					 const FiniteElemSpace *col) 
    : Assembler(op, row, col)
  {
    bool opt = (row == col);

    // create sub assemblers
    secondOrderAssembler = 
      SecondOrderAssembler::getSubAssembler(op, this, quad2, opt);
    firstOrderAssemblerGrdPsi = 
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, opt);
    firstOrderAssemblerGrdPhi = 
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, opt);
    zeroOrderAssembler = 
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, opt);

    checkQuadratures();
  }

  StandardAssembler::StandardAssembler(Operator *op,
				       Quadrature *quad2,
				       Quadrature *quad1GrdPsi,
				       Quadrature *quad1GrdPhi,
				       Quadrature *quad0,
				       const FiniteElemSpace *row,
				       const FiniteElemSpace *col) 
    : Assembler(op, row, col)
  {
    remember = false;

    // create sub assemblers
    secondOrderAssembler = 
      SecondOrderAssembler::getSubAssembler(op, this, quad2, false);
    firstOrderAssemblerGrdPsi = 
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, false);
    firstOrderAssemblerGrdPhi = 
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, false);
    zeroOrderAssembler = 
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);

    checkQuadratures();
  }

423
}