Assembler.cc 9.95 KB
Newer Older
1
2
3
4
5
6
7
#include "Assembler.h"
#include "Operator.h"
#include "Element.h"
#include "QPsiPhi.h"
#include "ElementMatrix.h"
#include "ElementVector.h"
#include "DOFVector.h"
8
#include "OpenMP.h"
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#include <vector>
#include <algorithm>

namespace AMDiS {

  Assembler::Assembler(Operator  *op,
		       const FiniteElemSpace *rowFESpace_,
		       const FiniteElemSpace *colFESpace_) 
    : operat(op),
      rowFESpace(rowFESpace_),
      colFESpace(colFESpace_ ? colFESpace_ : rowFESpace_),
      nRow(rowFESpace->getBasisFcts()->getNumber()),
      nCol(colFESpace->getBasisFcts()->getNumber()),
      remember(true),
      rememberElMat(false),
      rememberElVec(false),
      elementMatrix(NULL),
      elementVector(NULL),
      lastMatEl(NULL),
      lastVecEl(NULL),
      lastTraverseId(-1)
  
Thomas Witkowski's avatar
Thomas Witkowski committed
31
  {}
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

  void Assembler::calculateElementMatrix(const ElInfo *elInfo, 
					 ElementMatrix *userMat,
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

    if (remember && ((factor != 1.0) || (operat->uhOld))) {
      rememberElMat = true;
    }
  
    if (rememberElMat && !elementMatrix)
      elementMatrix = NEW ElementMatrix(nRow, nCol);

    Element *el = elInfo->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
48
    
49
50
51
    checkForNewTraverse();

    checkQuadratures();
Thomas Witkowski's avatar
Thomas Witkowski committed
52
    
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
      initElement(elInfo);
    }

    if (el != lastMatEl || !operat->isOptimized()) {
      if (rememberElMat) {
	elementMatrix->set(0.0);
      }
      lastMatEl = el;
    } else {
      if (rememberElMat) {
	axpy(factor, *elementMatrix, *userMat);
	return;
      }
    }
  
    ElementMatrix *mat = rememberElMat ? elementMatrix : userMat;

    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
80
    if (rememberElMat && userMat) {
81
      axpy(factor, *elementMatrix, *userMat);
82
    }      
83
84
  }

85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
  void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
					 const ElInfo *colElInfo,
					 ElementMatrix *userMat,
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

    if (remember && ((factor != 1.0) || (operat->uhOld))) {
      rememberElMat = true;
    }
  
    if (rememberElMat && !elementMatrix)
      elementMatrix = NEW ElementMatrix(nRow, nCol);

    Element *el = rowElInfo->getElement();

    
    //    checkForNewTraverse();
    lastVecEl = lastMatEl = NULL;

    checkQuadratures();
    
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
      initElement(rowElInfo);
    }

    if (el != lastMatEl || !operat->isOptimized()) {
      if (rememberElMat) {
	elementMatrix->set(0.0);
      }
      lastMatEl = el;
    } else {
      if (rememberElMat) {
	axpy(factor, *elementMatrix, *userMat);
	return;
      }
    }
  
    ElementMatrix *mat = rememberElMat ? elementMatrix : userMat;

    if (secondOrderAssembler)
      secondOrderAssembler->calculateElementMatrix(rowElInfo, mat);
    if (firstOrderAssemblerGrdPsi)
      firstOrderAssemblerGrdPsi->calculateElementMatrix(rowElInfo, mat);
    if (firstOrderAssemblerGrdPhi)
      firstOrderAssemblerGrdPhi->calculateElementMatrix(rowElInfo, mat);
    if (zeroOrderAssembler) {
      zeroOrderAssembler->calculateElementMatrix(rowElInfo, mat);
    }

    if (rememberElMat && userMat) {
      axpy(factor, *elementMatrix, *userMat);
    }      
  }

140
141
142
143
144
145
  void Assembler::calculateElementVector(const ElInfo *elInfo, 
					 ElementVector *userVec,
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

Thomas Witkowski's avatar
Thomas Witkowski committed
146
    if (remember && factor != 1.0) {
147
148
149
      rememberElVec = true;
    }

150
    if (rememberElVec && !elementVector) {
151
      elementVector = NEW ElementVector(nRow);
152
    }
153
154
155
156
157

    Element *el = elInfo->getElement();
    checkForNewTraverse();
    checkQuadratures();

Thomas Witkowski's avatar
Thomas Witkowski committed
158
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
159
160
      initElement(elInfo);
    }
161
    
Thomas Witkowski's avatar
Thomas Witkowski committed
162
163
    if (el != lastVecEl || !operat->isOptimized()) {
      if (rememberElVec) {
164
165
166
167
	elementVector->set(0.0);
      }
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
168
      if (rememberElVec) {
169
170
171
172
173
	axpy(factor, *elementVector, *userVec);
	return;
      }
    }
    ElementVector *vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
174
    if (operat->uhOld && remember) {
175
      matVecAssemble(elInfo, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
176
      if (rememberElVec) {
177
178
179
180
	axpy(factor, *elementVector, *userVec);
      }
      return;
    } 
181
    if (firstOrderAssemblerGrdPsi) {
182
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
183
184
    }
    if (zeroOrderAssembler) {
185
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
186
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
187
    if (rememberElVec) {
188
      axpy(factor, *elementVector, *userVec);
189
    }      
190
191
192
193
194
195
196
197
  }

  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector *vec)
  {
    FUNCNAME("Assembler::matVecAssemble()");

    Element *el = elInfo->getElement(); 
    const BasisFunction *basFcts = rowFESpace->getBasisFcts();
198
199
    int nBasFcts = basFcts->getNumber();
    double *uhOldLoc = new double[nBasFcts];
200

201
202
    operat->uhOld->getLocalVector(el, uhOldLoc);
    
Thomas Witkowski's avatar
Thomas Witkowski committed
203
    if (el != lastMatEl) {
204
205
      calculateElementMatrix(elInfo, NULL);
    }
206
207
208
209
210
    
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
	val += (*elementMatrix)[i][j] * uhOldLoc[j];
211
212
213
      }
      (*vec)[i] += val;
    }
214
215
216
    

    delete [] uhOldLoc;       
217
218
219
220
221
222
  }

  void Assembler::initElement(const ElInfo *elInfo, Quadrature *quad)
  {
    checkQuadratures();

Thomas Witkowski's avatar
Thomas Witkowski committed
223
    if (secondOrderAssembler) 
224
      secondOrderAssembler->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
225
    if (firstOrderAssemblerGrdPsi)
226
      firstOrderAssemblerGrdPsi->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
227
    if (firstOrderAssemblerGrdPhi)
228
      firstOrderAssemblerGrdPhi->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
229
    if (zeroOrderAssembler)
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
      zeroOrderAssembler->initElement(elInfo, quad);
  }

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

    // create sub assemblers
Thomas Witkowski's avatar
Thomas Witkowski committed
245
    secondOrderAssembler = 
246
      SecondOrderAssembler::getSubAssembler(op, this, quad2, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
247
    firstOrderAssemblerGrdPsi = 
248
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
249
    firstOrderAssemblerGrdPhi = 
250
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
251
    zeroOrderAssembler = 
252
253
254
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, opt);
  }

255
  StandardAssembler::StandardAssembler(Operator *op,
256
257
258
259
260
261
262
263
264
265
266
				       Quadrature *quad2,
				       Quadrature *quad1GrdPsi,
				       Quadrature *quad1GrdPhi,
				       Quadrature *quad0,
				       const FiniteElemSpace *rowFESpace_,
				       const FiniteElemSpace *colFESpace_) 
    : Assembler(op, rowFESpace_, colFESpace_)
  {
    remember = false;

    // create sub assemblers
Thomas Witkowski's avatar
Thomas Witkowski committed
267
    secondOrderAssembler = 
268
      SecondOrderAssembler::getSubAssembler(op, this, quad2, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
269
    firstOrderAssemblerGrdPsi = 
270
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
271
    firstOrderAssemblerGrdPhi = 
272
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
273
    zeroOrderAssembler = 
274
275
276
277
278
279
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);
  }

  ElementMatrix *Assembler::initElementMatrix(ElementMatrix *elMat, 
					      const ElInfo *elInfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
280
    if (!elMat) {
281
282
283
284
      elMat = NEW ElementMatrix(nRow, nCol);
    }

    elMat->set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
285
286
         
    rowFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
Thomas Witkowski's avatar
Thomas Witkowski committed
287
						   rowFESpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
288
						   &(elMat->rowIndices));
289

290

Thomas Witkowski's avatar
Thomas Witkowski committed
291
292
293
    if (rowFESpace == colFESpace) {
      elMat->colIndices = elMat->rowIndices;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
294
      colFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
Thomas Witkowski's avatar
Thomas Witkowski committed
295
296
297
						     colFESpace->getAdmin(),
						     &(elMat->colIndices));
    }
298
    
299
300
301
302
303
304
    return elMat;
  }

  ElementVector *Assembler::initElementVector(ElementVector *elVec, 
					      const ElInfo *elInfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
305
    if (!elVec) {
306
307
308
309
310
311
312
313
314
      elVec = NEW ElementVector(nRow);
    }

    elVec->set(0.0);
  
    DOFAdmin *rowAdmin = rowFESpace->getAdmin();

    Element *element = elInfo->getElement();

315
    rowFESpace->getBasisFcts()->getLocalIndicesVec(element, rowAdmin, &(elVec->dofIndices));
316
317
318
319
320

    return elVec;
  }

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

}