ZeroOrderTerm.cc 31.1 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
14
15
16
17
18
19
20
21
#include "ZeroOrderTerm.h"
#include "DOFVector.h"

namespace AMDiS {

  // ========== VecAtQP_ZOT ==========

  VecAtQP_ZOT::VecAtQP_ZOT(DOFVectorBase<double> *dv,
			   AbstractFunction<double, double> *af)
22
23
    : ZeroOrderTerm(af ? af->getDegree() : dv->getFeSpace()->getBasisFcts()->getDegree()), 
      vec(dv), f(af)
24
  {
25
26
    FUNCNAME("VecAtQP_ZOT::VecAtQP_ZOT()");

27
28
    TEST_EXIT(dv)("No vector!\n");

29
    auxFeSpaces.insert(dv->getFeSpace());
30
31
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
32

33
34
35
36
  void VecAtQP_ZOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad) 
  {
37
    getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
38
39
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
40

41
42
43
44
45
  void VecAtQP_ZOT::initElement(const ElInfo* smallElInfo,
				const ElInfo* largeElInfo,
				SubAssembler* subAssembler,
				Quadrature *quad)
  {
46
    getVectorAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad, vecAtQPs);
47
48
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
49

50
  void VecAtQP_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector& C)  
51
52
53
  { 
    if (f) {
      for (int iq = 0; iq < nPoints; iq++)
54
        C[iq] += (*f)(vecAtQPs[iq]);
55
56
    } else {
      for (int iq = 0; iq < nPoints; iq++)
57
        C[iq] += vecAtQPs[iq];
58
59
60
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
61

62
  void VecAtQP_ZOT::eval(int nPoints,
63
			 const ElementVector& uhAtQP,
64
65
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
66
			 ElementVector& result,
67
68
69
70
			 double fac) 
  {
    if (f) {
      for (int iq = 0; iq < nPoints; iq++)
71
        result[iq] += fac * (*f)(vecAtQPs[iq]) * uhAtQP[iq];
72
73
    } else {
      for (int iq = 0; iq < nPoints; iq++)
74
        result[iq] += fac * vecAtQPs[iq] * uhAtQP[iq];
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
    }
  }

  
  // ========== MultVecAtQP_ZOT ==========

  MultVecAtQP_ZOT::MultVecAtQP_ZOT(DOFVectorBase<double> *dv1,
				   DOFVectorBase<double> *dv2,
				   AbstractFunction<double, double> *af1,
				   AbstractFunction<double, double> *af2)
    : ZeroOrderTerm(af1->getDegree() + af2->getDegree()), 
      vec1(dv1), vec2(dv2), f1(af1), f2(af2)
  {
    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");

91
92
    auxFeSpaces.insert(dv1->getFeSpace());
    auxFeSpaces.insert(dv2->getFeSpace());
93
94
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
95

96
97
98
99
  void MultVecAtQP_ZOT::initElement(const ElInfo* elInfo, 
				    SubAssembler* subAssembler,
				    Quadrature *quad) 
  {
100
101
    getVectorAtQPs(vec1, elInfo, subAssembler, quad, vecAtQPs1);
    getVectorAtQPs(vec2, elInfo, subAssembler, quad, vecAtQPs2);
102
103
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
104

105
  void MultVecAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)  
106
107
108
109
110
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f1)(vecAtQPs1[iq]) * (*f2)(vecAtQPs2[iq]);    
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
111

112
  void MultVecAtQP_ZOT::eval(int nPoints,
113
			     const ElementVector& uhAtQP,
114
115
			     const WorldVector<double> *grdUhAtQP,
			     const WorldMatrix<double> *D2UhAtQP,
116
			     ElementVector& result,
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
			     double fac) 
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += fac * (*f1)(vecAtQPs1[iq]) * (*f2)(vecAtQPs2[iq]) * uhAtQP[iq];
  }

  
  // ========== Vec2AtQP_ZOT ==========

  Vec2AtQP_ZOT::Vec2AtQP_ZOT(DOFVectorBase<double> *dv1,
			     DOFVectorBase<double> *dv2,
			     BinaryAbstractFunction<double, double, double> *af)
    : ZeroOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), f(af)
  {
    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");

134
135
    auxFeSpaces.insert(dv1->getFeSpace());
    auxFeSpaces.insert(dv2->getFeSpace());
136
137
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
138

139
140
141
142
  void Vec2AtQP_ZOT::initElement(const ElInfo* elInfo, 
				 SubAssembler* subAssembler,
				 Quadrature *quad) 
  {
143
144
    getVectorAtQPs(vec1, elInfo, subAssembler, quad, vecAtQPs1);
    getVectorAtQPs(vec2, elInfo, subAssembler, quad, vecAtQPs2);
145
146
147
148
149
150
151
  }

  void Vec2AtQP_ZOT::initElement(const ElInfo* smallElInfo, 
				 const ElInfo* largeElInfo,
				 SubAssembler* subAssembler,
				 Quadrature *quad)
  {
152
    TEST_EXIT(vec1->getFeSpace() == vec2->getFeSpace())("Not yet implemented!\n");
153

154
155
    getVectorAtQPs(vec1, smallElInfo, largeElInfo, subAssembler, quad, vecAtQPs1);
    getVectorAtQPs(vec2, smallElInfo, largeElInfo, subAssembler, quad, vecAtQPs2);
156
157
  }

158
  void Vec2AtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)  
159
160
161
162
163
164
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);    
  }

  void Vec2AtQP_ZOT::eval(int nPoints,
165
			  const ElementVector& uhAtQP,
166
167
			  const WorldVector<double> *grdUhAtQP,
			  const WorldMatrix<double> *D2UhAtQP,
168
			  ElementVector& result,
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
			  double fac) 
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq]) * uhAtQP[iq];    
  }


  // ========== Vec3AtQP_ZOT ==========

  Vec3AtQP_ZOT::Vec3AtQP_ZOT(DOFVectorBase<double> *dv1,
			     DOFVectorBase<double> *dv2,
			     DOFVectorBase<double> *dv3,
			     TertiaryAbstractFunction<double, double, double, double> *af)
    : ZeroOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), vec3(dv3), f(af)
  {
    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");
    TEST_EXIT(dv3)("No thierd vector!\n");

188
189
190
    auxFeSpaces.insert(dv1->getFeSpace());
    auxFeSpaces.insert(dv2->getFeSpace());
    auxFeSpaces.insert(dv3->getFeSpace());
191
192
193
194
195
196
  }

  void Vec3AtQP_ZOT::initElement(const ElInfo* elInfo, 
				 SubAssembler* subAssembler,
				 Quadrature *quad) 
  {
197
198
199
    getVectorAtQPs(vec1, elInfo, subAssembler, quad, vecAtQPs1);
    getVectorAtQPs(vec2, elInfo, subAssembler, quad, vecAtQPs2);
    getVectorAtQPs(vec3, elInfo, subAssembler, quad, vecAtQPs3);
200
201
  }

202
  void Vec3AtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)  
203
204
205
206
207
208
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]);    
  }

  void Vec3AtQP_ZOT::eval(int nPoints,
209
			  const ElementVector& uhAtQP,
210
211
			  const WorldVector<double> *grdUhAtQP,
			  const WorldMatrix<double> *D2UhAtQP,
212
			  ElementVector& result,
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
			  double fac) 
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += 
	fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]) * uhAtQP[iq];    
  }


  // ========== FctGradientCoords_ZOT ==========

  FctGradientCoords_ZOT::FctGradientCoords_ZOT(DOFVectorBase<double> *dv,
					       BinaryAbstractFunction<double, WorldVector<double>, WorldVector<double> > *af)
    : ZeroOrderTerm(af->getDegree()), vec(dv), f(af)
  {
    TEST_EXIT(dv)("No vector!\n");

229
    auxFeSpaces.insert(dv->getFeSpace());
230
231
232
233
234
235
236
237
238
239
  }

  void FctGradientCoords_ZOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad) 
  {
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

240
  void FctGradientCoords_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)  
241
242
243
244
245
246
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(gradAtQPs[iq], coordsAtQPs[iq]);    
  }
 
  void FctGradientCoords_ZOT::eval(int nPoints,
247
				   const ElementVector& uhAtQP,
248
249
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
250
				   ElementVector& result,
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
				   double fac) 
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += 
	fac * (*f)(gradAtQPs[iq], coordsAtQPs[iq]) * uhAtQP[iq];    
  }

  
  // ========== VecGradCoordsAtQP_ZOT ==========
 
  VecGradCoordsAtQP_ZOT::VecGradCoordsAtQP_ZOT(DOFVectorBase<double> *dv,
					       TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *af)
    : ZeroOrderTerm(af->getDegree()), vec(dv), f(af) 
  {
    TEST_EXIT(dv)("No vector!\n");

267
    auxFeSpaces.insert(dv->getFeSpace());
268
269
270
271
272
273
  }

  void VecGradCoordsAtQP_ZOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad) 
  {
274
    getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
275
276
277
278
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

279
  void VecGradCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)  
280
281
282
283
284
285
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]);
  }
 
  void VecGradCoordsAtQP_ZOT::eval(int nPoints,
286
				   const ElementVector& uhAtQP,
287
288
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
289
				   ElementVector& result,
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
				   double fac) 
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += 
	fac * 
	(*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]) * 
	uhAtQP[iq];
  }


  // ========== VecAndCoordsAtQP_ZOT ==========

  VecAndCoordsAtQP_ZOT::VecAndCoordsAtQP_ZOT(DOFVectorBase<double> *dv,
					     BinaryAbstractFunction<double, double, WorldVector<double> > *af)
    : ZeroOrderTerm(af->getDegree()), vec(dv), f(af) 
  {
    TEST_EXIT(dv)("No vector!\n");

308
    auxFeSpaces.insert(dv->getFeSpace());
309
310
311
312
313
314
  }

  void VecAndCoordsAtQP_ZOT::initElement(const ElInfo* elInfo, 
					 SubAssembler* subAssembler,
					 Quadrature *quad) 
  {
315
    getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
316
317
318
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

319
  void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)  
320
321
322
323
324
325
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(vecAtQPs[iq], coordsAtQPs[iq]);    
  }

  void VecAndCoordsAtQP_ZOT::eval(int nPoints,
326
				  const ElementVector& uhAtQP,
327
328
				  const WorldVector<double> *grdUhAtQP,
				  const WorldMatrix<double> *D2UhAtQP,
329
				  ElementVector& result,
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
				  double fac) 
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += fac * (*f)(vecAtQPs[iq], coordsAtQPs[iq]) * uhAtQP[iq];    
  }

 
  // ========== Vec2AndGradAtQP_ZOT ==========

  Vec2AndGradAtQP_ZOT::Vec2AndGradAtQP_ZOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
					   TertiaryAbstractFunction<double, double, WorldVector<double>, double > *af)
    : ZeroOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), f(af) 
  {
    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");

346
347
    auxFeSpaces.insert(dv1->getFeSpace());
    auxFeSpaces.insert(dv2->getFeSpace());
348
349
350
351
352
353
  }

  void Vec2AndGradAtQP_ZOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) 
  {
354
355
    getVectorAtQPs(vec1, elInfo, subAssembler, quad, vecAtQPs1);
    getVectorAtQPs(vec2, elInfo, subAssembler, quad, vecAtQPs2);
356
357
358
359
    gradAtQPs = getGradientsAtQPs(vec1, elInfo, subAssembler, quad);
  }

  void Vec2AndGradAtQP_ZOT::eval(int nPoints,
360
				 const ElementVector& uhAtQP,
361
362
				 const WorldVector<double> *grdUhAtQP,
				 const WorldMatrix<double> *D2UhAtQP,
363
				 ElementVector& result,
364
365
366
367
368
369
370
371
372
				 double fac) 
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += 
	fac * 
	(*f)(vecAtQPs1[iq], gradAtQPs[iq], vecAtQPs2[iq]) * 
	uhAtQP[iq];
  }

373
  void Vec2AndGradAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)  
374
375
376
377
378
379
380
381
382
383
384
385
386
387
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(vecAtQPs1[iq], gradAtQPs[iq], vecAtQPs2[iq]);
  }


  // ========== FctGradient_ZOT ==========

  FctGradient_ZOT::FctGradient_ZOT(DOFVectorBase<double> *dv,
				   AbstractFunction<double, WorldVector<double> > *af)
    : ZeroOrderTerm(af->getDegree()), vec(dv), f(af)
  {
    TEST_EXIT(dv)("No vector!\n");

388
    auxFeSpaces.insert(dv->getFeSpace());
389
390
391
392
393
394
395
396
397
  }

  void FctGradient_ZOT::initElement(const ElInfo* elInfo, 
				    SubAssembler* subAssembler,
				    Quadrature *quad) 
  {
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
  }

398
  void FctGradient_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C) 
399
400
401
402
403
404
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(gradAtQPs[iq]);
  }
 
  void FctGradient_ZOT::eval(int nPoints,
405
			     const ElementVector& uhAtQP,
406
407
			     const WorldVector<double> *grdUhAtQP,
			     const WorldMatrix<double> *D2UhAtQP,
408
			     ElementVector& result,
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
			     double fac) 
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += fac * (*f)(gradAtQPs[iq]) * uhAtQP[iq];    
  }


  // ========== VecAndGradAtQP_ZOT ==========
 
  VecAndGradAtQP_ZOT::VecAndGradAtQP_ZOT(DOFVectorBase<double> *dv,
					 BinaryAbstractFunction<double, double, WorldVector<double> > *af)
    : ZeroOrderTerm(af->getDegree()), vec(dv), f(af) 
  {
    TEST_EXIT(dv)("No vector!\n");

424
    auxFeSpaces.insert(dv->getFeSpace());
425
426
427
428
429
430
  }

  void VecAndGradAtQP_ZOT::initElement(const ElInfo* elInfo, 
				       SubAssembler* subAssembler,
				       Quadrature *quad) 
  {
431
    getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
432
433
434
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
  }

435
  void VecAndGradAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C) 
436
437
438
439
440
441
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
  }
 
  void VecAndGradAtQP_ZOT::eval(int nPoints,
442
				const ElementVector& uhAtQP,
443
444
				const WorldVector<double> *grdUhAtQP,
				const WorldMatrix<double> *D2UhAtQP,
445
				ElementVector& result,
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
				double fac) 
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += fac * (*f)(vecAtQPs[iq], gradAtQPs[iq]) * uhAtQP[iq];   
  }


  // ========== VecAndGradVecAtQP_ZOT ==========
 
  VecAndGradVecAtQP_ZOT::VecAndGradVecAtQP_ZOT(DOFVectorBase<double> *dv, 
					       DOFVectorBase<double> *dGrd,
					       BinaryAbstractFunction<double, double, WorldVector<double> > *af) 
    : ZeroOrderTerm(af->getDegree()), vec(dv), vecGrd(dGrd), f(af) 
  {
    TEST_EXIT(dv)("No vector!\n");
    TEST_EXIT(dGrd)("No gradient vector!\n");

463
464
    auxFeSpaces.insert(dv->getFeSpace());
    auxFeSpaces.insert(dGrd->getFeSpace());
465
466
467
468
469
470
  }

  void VecAndGradVecAtQP_ZOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad) 
  {
471
    getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
472
473
474
    gradAtQPs = getGradientsAtQPs(vecGrd, elInfo, subAssembler, quad);
  }

475
  void VecAndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C) 
476
477
478
479
480
481
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
  }
 
  void VecAndGradVecAtQP_ZOT::eval(int nPoints,
482
				   const ElementVector& uhAtQP,
483
484
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
485
				   ElementVector& result,
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
				   double fac) 
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += fac * (*f)(vecAtQPs[iq], gradAtQPs[iq]) * uhAtQP[iq];
  }


  // ========== VecAndGradVec2AtQP_ZOT ==========

  VecAndGradVec2AtQP_ZOT::VecAndGradVec2AtQP_ZOT(DOFVectorBase<double> *dv, 
						 DOFVectorBase<double> *dGrd1, 
						 DOFVectorBase<double> *dGrd2, 
						 TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *af) 
    : ZeroOrderTerm(af->getDegree()), vec(dv), vecGrd1(dGrd1), vecGrd2(dGrd2), f(af) 
  {
    TEST_EXIT(dv)("No vector!\n");
    TEST_EXIT(dGrd1)("No first gradient vector!\n");
    TEST_EXIT(dGrd2)("No second gradient vector!\n");

505
506
507
    auxFeSpaces.insert(dv->getFeSpace());
    auxFeSpaces.insert(dGrd1->getFeSpace());
    auxFeSpaces.insert(dGrd2->getFeSpace());
508
509
510
511
512
513
  }

  void VecAndGradVec2AtQP_ZOT::initElement(const ElInfo* elInfo, 
					   SubAssembler* subAssembler,
					   Quadrature *quad) 
  {
514
    getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
515
516
517
518
    grad1AtQPs = getGradientsAtQPs(vecGrd1, elInfo, subAssembler, quad);
    grad2AtQPs = getGradientsAtQPs(vecGrd2, elInfo, subAssembler, quad);
  }

519
  void VecAndGradVec2AtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
520
521
522
523
524
525
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(vecAtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]);
  }
 
  void VecAndGradVec2AtQP_ZOT::eval(int nPoints,
526
				    const ElementVector& uhAtQP,
527
528
				    const WorldVector<double> *grdUhAtQP,
				    const WorldMatrix<double> *D2UhAtQP,
529
				    ElementVector& result,
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
				    double fac) 
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += 
	fac * (*f)(vecAtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]) * uhAtQP[iq];  
  }


  // ========== VecOfDOFVecsAtQP_ZOT ==========
 
  VecOfDOFVecsAtQP_ZOT::VecOfDOFVecsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
					     AbstractFunction<double, std::vector<double> > *af)
    : ZeroOrderTerm(af->getDegree()), vecs(dv), f(af) 
  {
    vecsAtQPs.resize(vecs.size());

Thomas Witkowski's avatar
Thomas Witkowski committed
546
    for (unsigned int i = 0; i < dv.size(); i++) {
547
548
      TEST_EXIT(dv[i])("One vector is NULL!\n");

549
      auxFeSpaces.insert(dv[i]->getFeSpace());
550
551
552
    }
  } 

Thomas Witkowski's avatar
Thomas Witkowski committed
553

554
555
556
557
  void VecOfDOFVecsAtQP_ZOT::initElement(const ElInfo* elInfo, 
					 SubAssembler* subAssembler,
					 Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
558
559
    unsigned int size = vecs.size();
    for (unsigned int i = 0; i < size; i++)
560
      getVectorAtQPs(vecs[i], elInfo, subAssembler, quad, vecsAtQPs[i]);
561
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
562

563
 
564
  void VecOfDOFVecsAtQP_ZOT::getC(const ElInfo *, int nPoints,  ElementVector& C)
565
  { 
Thomas Witkowski's avatar
Thomas Witkowski committed
566
    unsigned int size = vecs.size();
567
568
569
    std::vector<double> arg(size);

    for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
570
      for (unsigned int i = 0; i < size; i++)
571
572
573
574
575
576
	arg[i] = vecsAtQPs[i][iq];
      
      C[iq] += (*f)(arg);
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
577

578
  void VecOfDOFVecsAtQP_ZOT::eval(int nPoints,
579
				  const ElementVector& uhAtQP,
580
581
				  const WorldVector<double> *grdUhAtQP,
				  const WorldMatrix<double> *D2UhAtQP,
582
				  ElementVector& result,
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
				  double fac) 
  {
    int size = static_cast<int>(vecs.size());
    std::vector<double> arg(size);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < size; i++)
	arg[i] = vecsAtQPs[i][iq];
      
      result[iq] += fac * (*f)(arg) * uhAtQP[iq];
    }
  }


  // ========== VecOfGradientsAtQP_ZOT ==========

  VecOfGradientsAtQP_ZOT::VecOfGradientsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv,
						 AbstractFunction<double, std::vector<WorldVector<double>*> > *af) 
    : ZeroOrderTerm(af->getDegree()), vecs(dv), f(af) 
  {
    gradsAtQPs.resize(vecs.size());

    for (int i = 0; i < static_cast<int>(dv.size()); i++) {
      TEST_EXIT(dv[i])("One vector is NULL!\n");

608
      auxFeSpaces.insert(dv[i]->getFeSpace());
609
610
611
612
613
614
615
616
617
618
619
620
    }
  } 

  void VecOfGradientsAtQP_ZOT::initElement(const ElInfo* elInfo, 
					   SubAssembler* subAssembler,
					   Quadrature *quad) 
  {
    int size = static_cast<int>(vecs.size());
    for (int i = 0; i < size; i++)
      gradsAtQPs[i] = getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad);    
  }

621
  void VecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
622
623
624
625
626
627
628
629
630
631
632
633
634
  { 
    int size = static_cast<int>(vecs.size());
    std::vector<WorldVector<double>*> arg(size);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < size; i++)
	arg[i] = &(gradsAtQPs[i][iq]);
      
      C[iq] += (*f)(arg);
    }
  }

  void VecOfGradientsAtQP_ZOT::eval(int nPoints,
635
				    const ElementVector& uhAtQP,
636
637
				    const WorldVector<double> *grdUhAtQP,
				    const WorldMatrix<double> *D2UhAtQP,
638
				    ElementVector& result,
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
				    double fac) 
  {
    int size = static_cast<int>(vecs.size());
    std::vector<WorldVector<double>*> arg(size);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < size; i++)
	arg[i] = &(gradsAtQPs[i][iq]);
      
      result[iq] += fac * (*f)(arg) * uhAtQP[iq];
    }
  }


  // ========== VecDivergence_ZOT ==========
 
  VecDivergence_ZOT::VecDivergence_ZOT(int nComponents,
				       DOFVectorBase<double> *vec0,
				       DOFVectorBase<double> *vec1,
				       DOFVectorBase<double> *vec2)
    : ZeroOrderTerm(0)
  {
    vecs.resize(nComponents);
    gradsAtQPs.resize(nComponents);
    vecs[0] = vec0;
    vecs[1] = vec1;
    vecs[2] = vec2;

667
    auxFeSpaces.insert(vec0->getFeSpace());
668
    if (vec1) 
669
      auxFeSpaces.insert(vec1->getFeSpace());
670
    if (vec2) 
671
      auxFeSpaces.insert(vec2->getFeSpace());
672
673
674
675
676
677
678
679
680
681
682
683
  }

  void VecDivergence_ZOT::initElement(const ElInfo* elInfo, 
				      SubAssembler* subAssembler,
				      Quadrature *quad) 
  {
    int size = static_cast<int>(vecs.size());
    for (int i = 0; i < size; i++)
      gradsAtQPs[i] = getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad);    
  }


684
  void VecDivergence_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
685
686
687
688
689
690
691
692
693
  { 
    int size = static_cast<int>(vecs.size());

    for (int iq = 0; iq < nPoints; iq++)
      for (int i = 0; i < size; i++)
	C[iq] += gradsAtQPs[i][iq][i];    
  }

  void VecDivergence_ZOT::eval(int nPoints,
694
			       const ElementVector& uhAtQP,
695
696
			       const WorldVector<double> *grdUhAtQP,
			       const WorldMatrix<double> *D2UhAtQP,
697
			       ElementVector& result,
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
			       double fac) 
  {
    int size = static_cast<int>(vecs.size());

    for (int iq = 0; iq < nPoints; iq++) {
      double d = 0.0;
      for (int i = 0; i < size; i++)
	d += gradsAtQPs[i][iq][i];
      
      result[iq] += d * uhAtQP[iq] * fac;
    }
  }


  // ========== VecAndVecOfGradientsAtQP_ZOT =========

  VecAndVecOfGradientsAtQP_ZOT::VecAndVecOfGradientsAtQP_ZOT(DOFVector<double> *v,
							     const std::vector<DOFVector<double>*>& dv,
							     BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *af)
    : ZeroOrderTerm(af->getDegree()), vec(v), vecs(dv), f(af)
  {
    gradsAtQPs.resize(vecs.size());

    TEST_EXIT(v)("No vector!\n");

723
    auxFeSpaces.insert(v->getFeSpace());
724
725
726
    for (int i = 0; i < static_cast<int>(dv.size()); i++) {
      TEST_EXIT(dv[i])("One gradient vector is NULL!\n");

727
      auxFeSpaces.insert(dv[i]->getFeSpace());
728
729
730
731
732
733
734
735
736
737
738
    }
  }

  void VecAndVecOfGradientsAtQP_ZOT::initElement(const ElInfo* elInfo,
						 SubAssembler* subAssembler,
						 Quadrature *quad)
  {
    int size = static_cast<int>(vecs.size());
    for (int i = 0; i < size; i++)
      gradsAtQPs[i] = getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad);

739
    getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
740
741
742
  }

  void VecAndVecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints,
743
					  ElementVector& C)
744
745
746
747
748
749
750
751
752
753
754
755
756
  {
    int size = static_cast<int>(vecs.size());
    std::vector<WorldVector<double>*> arg(size);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < size; i++)
	arg[i] = &(gradsAtQPs[i][iq]);

      C[iq] += (*f)(vecAtQPs[iq], arg);
    }
  }

  void VecAndVecOfGradientsAtQP_ZOT::eval(int nPoints,
757
					  const ElementVector& uhAtQP,
758
759
					  const WorldVector<double> *grdUhAtQP,
					  const WorldMatrix<double> *D2UhAtQP,
760
					  ElementVector& result,
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
					  double fac) 
  {
    int size = static_cast<int>(vecs.size());
    std::vector<WorldVector<double>*> arg(size);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < size; i++)
	arg[i] = &(gradsAtQPs[i][iq]);

      result[iq] += fac * (*f)(vecAtQPs[iq], arg) * uhAtQP[iq];
    }
  }


  // ========== Vec2AndGrad2AtQP_ZOT ==========

  Vec2AndGrad2AtQP_ZOT::Vec2AndGrad2AtQP_ZOT(DOFVectorBase<double> *dv1, 
					     DOFVectorBase<double> *dv2,
					     QuartAbstractFunction<double, double, double, WorldVector<double>,WorldVector<double> > *af)
    : ZeroOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), f(af) 
  {
    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");
    
785
786
    auxFeSpaces.insert(dv1->getFeSpace());
    auxFeSpaces.insert(dv2->getFeSpace());
787
788
789
790
791
792
  }

  void Vec2AndGrad2AtQP_ZOT::initElement(const ElInfo* elInfo, 
					 SubAssembler* subAssembler,
					 Quadrature *quad) 
  {
793
794
    getVectorAtQPs(vec1, elInfo, subAssembler, quad, vecAtQPs1);
    getVectorAtQPs(vec2, elInfo, subAssembler, quad, vecAtQPs2);
795
796
797
798
799
800
    gradAtQPs1 = getGradientsAtQPs(vec1, elInfo, subAssembler, quad);
    gradAtQPs2 = getGradientsAtQPs(vec2, elInfo, subAssembler, quad);
  }
  
 
  void Vec2AndGrad2AtQP_ZOT::eval(int nPoints,
801
				  const ElementVector& uhAtQP,
802
803
				  const WorldVector<double> *grdUhAtQP,
				  const WorldMatrix<double> *D2UhAtQP,
804
				  ElementVector& result,
805
806
807
808
809
810
811
812
				  double fac) 
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += 
	fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq], gradAtQPs1[iq], gradAtQPs2[iq]) * 
	uhAtQP[iq];
  }

813
  void Vec2AndGrad2AtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)  
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], gradAtQPs1[iq], gradAtQPs2[iq]);
  }

  
  // ========== Vec2AndGradVecAtQP_ZOT ==========

  Vec2AndGradVecAtQP_ZOT::Vec2AndGradVecAtQP_ZOT(DOFVectorBase<double> *dv1, 
						 DOFVectorBase<double> *dv2, 
						 DOFVectorBase<double> *dGrd,
						 TertiaryAbstractFunction<double, double,double, WorldVector<double> > *af) 
    : ZeroOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), vecGrd(dGrd), f(af) 
  {
    TEST_EXIT(dv1)("No vector!\n");
    TEST_EXIT(dv2)("No vector!\n");
    TEST_EXIT(dGrd)("No gradient vector!\n");
    
832
833
834
    auxFeSpaces.insert(dv1->getFeSpace());
    auxFeSpaces.insert(dv2->getFeSpace());
    auxFeSpaces.insert(dGrd->getFeSpace());
835
836
837
838
839
840
  }
  
  void Vec2AndGradVecAtQP_ZOT::initElement(const ElInfo* elInfo, 
					   SubAssembler* subAssembler,
					   Quadrature *quad) 
  {
841
842
    getVectorAtQPs(vec1, elInfo, subAssembler, quad, vec1AtQPs);
    getVectorAtQPs(vec2, elInfo, subAssembler, quad, vec2AtQPs);
843
844
845
    gradAtQPs = getGradientsAtQPs(vecGrd, elInfo, subAssembler, quad);
  }
  
846
  void Vec2AndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)  
847
848
849
850
851
852
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*f)(vec1AtQPs[iq], vec2AtQPs[iq], gradAtQPs[iq]);
  }
  
  void Vec2AndGradVecAtQP_ZOT::eval(int nPoints,
853
				   const ElementVector& uhAtQP,
854
855
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
856
				   ElementVector& result,
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
				   double fac) 
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += 
	fac * (*f)(vec1AtQPs[iq], vec2AtQPs[iq], gradAtQPs[iq]) * uhAtQP[iq];
  }


  // =========== General_ZOT ==========
  
  General_ZOT::General_ZOT(std::vector<DOFVectorBase<double>*> vecs,
			   std::vector<DOFVectorBase<double>*> grads,
			   TertiaryAbstractFunction<double, WorldVector<double>, std::vector<double>, std::vector<WorldVector<double> > > *af)
    : ZeroOrderTerm(af->getDegree()), vecs_(vecs), grads_(grads), f_(af)
  {
872
    vecsAtQPs.resize(vecs_.size());
873
874
875
876
877
    gradsAtQPs_.resize(grads_.size());

    for (int i = 0; i < static_cast<int>(vecs.size()); i++) {
      TEST_EXIT(vecs[i])("One vector is NULL!\n");

878
      auxFeSpaces.insert(vecs[i]->getFeSpace());
879
880
881
882
883
    }   

    for (int i = 0; i < static_cast<int>(grads.size()); i++) {
      TEST_EXIT(grads[i])("One gradient vector is NULL!\n");

884
      auxFeSpaces.insert(grads[i]->getFeSpace());
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
    }   

    vecsArg.resize(vecs_.size());
    gradsArg.resize(grads_.size());
  }

  void General_ZOT::initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
				Quadrature *quad)
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);

    for (int i = 0; i < nVecs; i++)
900
      getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad, vecsAtQPs[i]);
901
902
903
904
    for (int i = 0; i < nGrads; i++)
      gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
  }

905
  void General_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
906
907
908
909
910
911
  {
    unsigned int nVecs = vecs_.size();
    unsigned int nGrads = grads_.size();

    for (int iq = 0; iq < nPoints; iq++) {
      for (unsigned int i = 0; i < nVecs; i++)
912
	vecsArg[i] = vecsAtQPs[i][iq];
913
914
915
916
917
918
919
920
      for (unsigned int i = 0; i < nGrads; i++)
	gradsArg[i] = gradsAtQPs_[i][iq];

      C[iq] += (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
    }
  }

  void General_ZOT::eval(int nPoints,
921
			 const ElementVector& uhAtQP,
922
923
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
924
			 ElementVector& result,
925
926
927
928
929
930
931
			 double fac) 
  {
    unsigned int nVecs = vecs_.size();
    unsigned int nGrads = grads_.size();

    for (int iq = 0; iq < nPoints; iq++) {
      for (unsigned int i = 0; i < nVecs; i++)
932
	vecsArg[i] = vecsAtQPs[i][iq];      
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
      for (unsigned int i = 0; i < nGrads; i++) 
	gradsArg[i] = gradsAtQPs_[i][iq];      

      result[iq] += fac * (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg) * uhAtQP[iq];
    }
  }


  // ========== GeneralParametric_ZOT ==========

  GeneralParametric_ZOT::GeneralParametric_ZOT(std::vector<DOFVectorBase<double>*> vecs,
					       std::vector<DOFVectorBase<double>*> grads,
					       QuartAbstractFunction<double, 
					       WorldVector<double>,
					       WorldVector<double>,
					       std::vector<double>, 
					       std::vector<WorldVector<double> > > *af)
    : ZeroOrderTerm(af->getDegree()), vecs_(vecs), grads_(grads), f_(af)
  {
952
953
    vecsAtQPs.resize(vecs_.size());
    gradsAtQPs.resize(grads_.size());
954
955
956
957

    for (int i = 0; i < static_cast<int>(vecs.size()); i++) {
      TEST_EXIT(vecs[i])("One vector is NULL!\n");

958
      auxFeSpaces.insert(vecs[i]->getFeSpace());
959
960
961
962
963
    }   

    for (int i = 0; i < static_cast<int>(grads.size()); i++) {
      TEST_EXIT(grads[i])("One gradient vector is NULL!\n");

964
      auxFeSpaces.insert(grads[i]->getFeSpace());
965
966
967
968
969
970
971
972
973
974
975
976
977
    }   
  }

  void GeneralParametric_ZOT::initElement(const ElInfo* elInfo, SubAssembler* 
					  subAssembler, Quadrature *quad)
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    elInfo->getElementNormal(elementNormal_);
    coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);

    for (int i = 0; i < nVecs; i++)
978
      getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad, vecsAtQPs[i]);
979
    for (int i = 0; i < nGrads; i++)
980
      gradsAtQPs[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
981
982
  }

983
  void GeneralParametric_ZOT::getC(const ElInfo *, int nPoints, ElementVector& C)
984
985
986
987
988
989
990
991
992
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++)
993
	vecsArg[i] = vecsAtQPs[i][iq];      
994
      for (int i = 0; i < nGrads; i++)
995
	gradsArg[i] = gradsAtQPs[i][iq];
996
997
998
999
1000
1001
      
      C[iq] += (*f_)(coordsAtQPs_[iq],  elementNormal_, vecsArg, gradsArg);
    }
  }

  void GeneralParametric_ZOT::eval(int nPoints,
1002
				   const ElementVector& uhAtQP,
1003
1004
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
1005
				   ElementVector& result,
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
				   double fac) 
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++)
1016
	vecsArg[i] = vecsAtQPs[i][iq];     
1017
      for (int i = 0; i < nGrads; i++)
1018
	gradsArg[i] = gradsAtQPs[i][iq];      
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
      result[iq] += 
	fac * (*f_)(coordsAtQPs_[iq],  elementNormal_, vecsArg, gradsArg) * uhAtQP[iq];
    }
  }


  // ========== CoordsAtQP_ZOT ==========

  void CoordsAtQP_ZOT::initElement(const ElInfo* elInfo, 
				   SubAssembler* subAssembler,
				   Quadrature *quad) 
  {
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

1034
  void CoordsAtQP_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector& C)
1035
1036
1037
1038
1039
1040
  { 
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*g)(coordsAtQPs[iq]);    
  }

  void CoordsAtQP_ZOT::eval(int nPoints,
1041
			    const ElementVector& uhAtQP,
1042
1043
			    const WorldVector<double> *grdUhAtQP,
			    const WorldMatrix<double> *D2UhAtQP,
1044
			    ElementVector& result,
1045
1046
1047
1048
1049
1050
1051
1052
			    double fac) 
  {
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += fac * (*g)(coordsAtQPs[iq]) * uhAtQP[iq];    
  }


}