Assembler.cc 10.3 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
  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)
126
      secondOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo, mat);
127
    if (firstOrderAssemblerGrdPsi)
128
      firstOrderAssemblerGrdPsi->calculateElementMatrix(rowElInfo, colElInfo, mat);
129
    if (firstOrderAssemblerGrdPhi)
130
      firstOrderAssemblerGrdPhi->calculateElementMatrix(rowElInfo, colElInfo, mat);
131
    if (zeroOrderAssembler) {
132
      zeroOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo, mat);
133
    }
134
 
135
136
137
138
139
    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
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);
  }

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
291
292
293
    if (rowFESpace == colFESpace) {
      elMat->colIndices = elMat->rowIndices;
    } else {
294
295
296
297
298
299
300
301
302
303
      if (colElInfo) {
	colFESpace->getBasisFcts()->getLocalIndicesVec(colElInfo->getElement(),
						       colFESpace->getAdmin(),
						       &(elMat->colIndices));
      } else {
	// If there is no colElInfo pointer, use rowElInfo the get the indices.
	colFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
						       colFESpace->getAdmin(),
						       &(elMat->colIndices));
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
304
    }
305
    
306
307
308
309
310
311
    return elMat;
  }

  ElementVector *Assembler::initElementVector(ElementVector *elVec, 
					      const ElInfo *elInfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
312
    if (!elVec) {
313
314
315
316
317
318
319
320
      elVec = NEW ElementVector(nRow);
    }

    elVec->set(0.0);
  
    DOFAdmin *rowAdmin = rowFESpace->getAdmin();
    Element *element = elInfo->getElement();

321
    rowFESpace->getBasisFcts()->getLocalIndicesVec(element, rowAdmin, &(elVec->dofIndices));
322
323
324
325
326

    return elVec;
  }

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

}