ZeroOrderAssembler.cc 9.29 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
#include <vector>
14
#include <boost/numeric/mtl/mtl.hpp>
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#include "Assembler.h"
#include "ZeroOrderAssembler.h"
#include "Operator.h"
#include "QPsiPhi.h"
#include "FiniteElemSpace.h"
#include "Quadrature.h"
#include "DOFVector.h"

namespace AMDiS {

  std::vector<SubAssembler*> ZeroOrderAssembler::optimizedSubAssemblers;
  std::vector<SubAssembler*> ZeroOrderAssembler::standardSubAssemblers;

  ZeroOrderAssembler::ZeroOrderAssembler(Operator *op,
					 Assembler *assembler,
					 Quadrature *quad,
					 bool optimized)
    : SubAssembler(op, assembler, quad, 0, optimized)
  {}

35

36
37
38
39
40
  ZeroOrderAssembler* ZeroOrderAssembler::getSubAssembler(Operator* op,
							  Assembler *assembler,
							  Quadrature *quad,
							  bool optimized)
  {
41
    // check if an assembler is needed at all
42
    if (op->zeroOrder.size() == 0)
43
      return NULL;   
44
45
46
47

    ZeroOrderAssembler *newAssembler;

    std::vector<SubAssembler*> *subAssemblers =
48
      optimized ? &optimizedSubAssemblers : &standardSubAssemblers;
49

50
    std::vector<OperatorTerm*> opTerms = op->zeroOrder;
51
52
53
54
55

    sort(opTerms.begin(), opTerms.end());

    // check if a new assembler is needed
    if (quad) {
Thomas Witkowski's avatar
Thomas Witkowski committed
56
      for (int i = 0; i < static_cast<int>(subAssemblers->size()); i++) {
57
58
59
60
	std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());

	sort(assTerms.begin(), assTerms.end());

61
	if (opTerms == assTerms && (*subAssemblers)[i]->getQuadrature() == quad)
62
63
64
65
66
67
	  return dynamic_cast<ZeroOrderAssembler*>((*subAssemblers)[i]);
      }
    }
 
    // check if all terms are pw_const
    bool pwConst = true;
68
    for (unsigned int i = 0; i < op->zeroOrder.size(); i++) {
69
      if (!op->zeroOrder[i]->isPWConst()) {
70
71
72
73
74
75
	pwConst = false;
	break;
      }
    }  

    // create new assembler
76
    if (!optimized) {
77
      newAssembler = new StandardZOA(op, assembler, quad);
78
    } else {
79
      if (pwConst)
80
	newAssembler = new PrecalcZOA(op, assembler, quad);
81
      else
82
     	newAssembler = new FastQuadZOA(op, assembler, quad);      
83
    }
84
85
86
87
88

    subAssemblers->push_back(newAssembler);
    return newAssembler;
  }

89

90
  StandardZOA::StandardZOA(Operator *op, Assembler *assembler, Quadrature *quad)
91
92
93
94
95
    : ZeroOrderAssembler(op, assembler, quad, false)      
  {
    name = "standard zero order assembler";
  }

96

97
98
  void StandardZOA::calculateElementMatrix(const ElInfo *elInfo, 
					   ElementMatrix& mat)
99
  {
100
101
    const BasisFunction *psi = rowFeSpace->getBasisFcts();
    const BasisFunction *phi = colFeSpace->getBasisFcts();
102
103

    int nPoints = quadrature->getNumPoints();
104
    ElementVector c(nPoints, 0.0);
105
    std::vector<double> phival(nCol);
106

107
108
    for (std::vector<OperatorTerm*>::iterator termIt = terms.begin(); 
	 termIt != terms.end(); ++termIt)
109
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
110

111
    if (symmetric) {
112
113
      TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
            
114
115
116
117
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

	// calculate phi at QPs only once!
118
	for (int i = 0; i < nCol; i++)
119
120
121
122
	  phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));

	for (int i = 0; i < nRow; i++) {
	  double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
123
	  mat[i][i] += quadrature->getWeight(iq) * c[iq] * psival * phival[i];
124
125
	  for (int j = i + 1; j < nCol; j++) {
	    double val = quadrature->getWeight(iq) * c[iq] * psival * phival[j];
126
127
	    mat[i][j] += val;
	    mat[j][i] += val;
128
129
130
131
132
133
134
135
	  }
	}
      }
    } else {      //  non symmetric assembling 
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

	// calculate phi at QPs only once!
136
	for (int i = 0; i < nCol; i++)
137
138
139
140
	  phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));

	for (int i = 0; i < nRow; i++) {
	  double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
141
	  for (int j = 0; j < nCol; j++)
142
	    mat[i][j] += quadrature->getWeight(iq) * c[iq] * psival * phival[j];	  
143
144
145
146
147
	}
      }
    }
  }

148

149
  void StandardZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
150
151
  {
    int nPoints = quadrature->getNumPoints();
152
    ElementVector c(nPoints, 0.0);
153

154
    std::vector<OperatorTerm*>::iterator termIt;
155
    for (termIt = terms.begin(); termIt != terms.end(); ++termIt)
156
157
158
159
160
161
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);

    for (int iq = 0; iq < nPoints; iq++) {
      c[iq] *= elInfo->getDet();

      for (int i = 0; i < nRow; i++) {
162
	double psi = (*(rowFeSpace->getBasisFcts()->getPhi(i)))
163
	  (quadrature->getLambda(iq));
164
	vec[i] += quadrature->getWeight(iq) * c[iq] * psi;
165
166
167
168
      }
    }
  }

169

170
171
  FastQuadZOA::FastQuadZOA(Operator *op, Assembler *assembler, Quadrature *quad)
    : ZeroOrderAssembler(op, assembler, quad, true)
172
  {
173
    name = "fast quadrature zero order assembler";
174
  }
175

176

177
  void FastQuadZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
178
179
180
181
  {
    int nPoints = quadrature->getNumPoints();

    if (firstCall) {
182
183
184
185
186
      const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
      basFcts = colFeSpace->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
      firstCall = false;      
187
188
    }

189
    if (num_rows(c) < static_cast<unsigned int>(nPoints))
190
191
      c.change_dim(nPoints);
    c = 0.0;
192

193
194
    for (std::vector<OperatorTerm*>::iterator termIt = terms.begin(); 
	 termIt != terms.end(); ++termIt)
195
196
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);

197
198
199
    const mtl::dense2D<double>& psi = psiFast->getPhi();
    const mtl::dense2D<double>& phi = phiFast->getPhi();

200
    if (symmetric) {
201
202
      TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");

203
      for (int iq = 0; iq < nPoints; iq++) {
204
205
	c[iq] *= elInfo->getDet() * quadrature->getWeight(iq);

206
	for (int i = 0; i < nRow; i++) {
207
	  mat[i][i] += c[iq] * psi[iq][i] * phi[iq][i];
208
	  for (int j = i + 1; j < nCol; j++) {
209
	    double val = c[iq] * psi[iq][i] * phi[iq][j];
210
211
	    mat[i][j] += val;
	    mat[j][i] += val;
212
213
214
215
	  }
	}
      }
    } else {      /*  non symmetric assembling   */
216
      for (int iq = 0; iq < nPoints; iq++)
217
	c[iq] *= elInfo->getDet() * quadrature->getWeight(iq);
218

219
220
221
222
223
224
225
      for (int i = 0; i < nRow; i++)
	for (int j = 0; j < nCol; j++) {
	  double v = 0.0;
	  for (int iq = 0; iq < nPoints; iq++)
	    v += c[iq] * psi[iq][i] * phi[iq][j];
	  mat[i][j] += v;
	}
226
227
    }
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
228

229

230
  void FastQuadZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
231
  {
232
233
    int nPoints = quadrature->getNumPoints();

234
    if (firstCall) {
235
236
237
238
239
      const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
      basFcts = colFeSpace->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
      firstCall = false;      
240
241
    }

242
    ElementVector c(nPoints, 0.0);
243
    std::vector<OperatorTerm*>::iterator termIt;
244
    for (termIt = terms.begin(); termIt != terms.end(); ++termIt)
245
246
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);

247
248
    const mtl::dense2D<double>& psi = psiFast->getPhi();

249
250
251
    for (int iq = 0; iq < nPoints; iq++) {
      c[iq] *= elInfo->getDet();

252
      for (int i = 0; i < nRow; i++)
253
	vec[i] += quadrature->getWeight(iq) * c[iq] * psi[iq][i];
254
255
256
    }
  }

257

258
259
260
  PrecalcZOA::PrecalcZOA(Operator *op, Assembler *assembler, Quadrature *quad) 
    : ZeroOrderAssembler(op, assembler, quad, true)
  {
261
    name = "precalculated zero order assembler";
262
263
  }

264

265
266
  void PrecalcZOA::calculateElementMatrix(const ElInfo *elInfo, 
					  ElementMatrix& mat)
267
  {
268
269
    FUNCNAME("PrecalcZOA::calculateElementMatrix()");

270
    if (firstCall) {
271
272
273
274
275
      q00 = Q00PsiPhi::provideQ00PsiPhi(rowFeSpace->getBasisFcts(), 
					colFeSpace->getBasisFcts(), 
					quadrature);
      q0 = Q0Psi::provideQ0Psi(rowFeSpace->getBasisFcts(), quadrature);
      firstCall = false;      
276
277
    }

278
    ElementVector c(1, 0.0);
279
    int size = static_cast<int>(terms.size());
280

281
    for (int i = 0; i < size; i++)
282
      (static_cast<ZeroOrderTerm*>((terms[i])))->getC(elInfo, 1, c);    
283

284
    c[0] *= elInfo->getDet();
285
    
286
287
    if (symmetric) {
      for (int i = 0; i < nRow; i++) {
288
289
290
291
292
293
    	mat[i][i] += c[0] * q00->getValue(i,i);
    	for (int j = i + 1; j < nCol; j++) {
    	  double val = c[0] * q00->getValue(i, j);
    	  mat[i][j] += val;
    	  mat[j][i] += val;
    	}
294
295
      }
    } else {
296
297
      for (int i = 0; i < nRow; i++)
	for (int j = 0; j < nCol; j++)
298
	  mat[i][j] += c[0] * q00->getValue(i, j);
299
300
301
    }
  }

302

303
304
  void PrecalcZOA::calculateElementVector(const ElInfo *elInfo, 
					  ElementVector& vec)
305
  {
306
307
    FUNCNAME("PrecalcZOA::calculateElementVector()");

308
    if (firstCall) {
309
310
311
312
313
      q00 = Q00PsiPhi::provideQ00PsiPhi(rowFeSpace->getBasisFcts(), 
					colFeSpace->getBasisFcts(), 
					quadrature);
      q0 = Q0Psi::provideQ0Psi(rowFeSpace->getBasisFcts(), quadrature);	
      firstCall = false;      
314
315
316
317
    }

    std::vector<OperatorTerm*>::iterator termIt;

318
    ElementVector c(1, 0.0);
319
    for (termIt = terms.begin(); termIt != terms.end(); ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
320
      (static_cast<ZeroOrderTerm*>(*termIt))->getC(elInfo, 1, c);
321

322
    c[0] *= elInfo->getDet();
323
324

    for (int i = 0; i < nRow; i++)
325
      vec[i] += c[0] * q0->getValue(i);
326
327
328
  }

}