POperators.cc 87.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
#include "POperators.h"
#include "Helpers.h"

using namespace std;
using namespace AMDiS;

Phase_SOT::Phase_SOT(DOFVectorBase<double>* phaseDV_,double fac_) :
	SecondOrderTerm(6), phaseDV(phaseDV_), fac(fac_)
{
	TEST_EXIT(phaseDV_)("no vector phase!\n");
	auxFeSpaces.insert(phaseDV_->getFeSpace());
12
}
13
14
15
16
void Phase_SOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) {
	getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
17
}
18
19
20
21
void Phase_SOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad) {
	getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
22
}
23
24
25
26
27
28
29
30
void Phase_SOT::getLALt(const ElInfo *elInfo,
			vector<mtl::dense2D<double> > &LALt) const
{
	const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(LALt.size());
    
	for (int iq = 0; iq < nPoints; iq++)
		l1lt(Lambda, LALt[iq], f(iq) * phase[iq] * fac);
31
}
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
void Phase_SOT::eval(int nPoints,
			const mtl::dense_vector<double>& uhAtQP,
			const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			mtl::dense_vector<double>& result,
			double opFactor)
{	
	if (num_rows(D2UhAtQP) > 0) {
		for (int iq = 0; iq < nPoints; iq++) {
			double feval = f(iq) * phase[iq] * opFactor * fac;
			double resultQP = 0.0;
			for (int i = 0; i < dimOfWorld; i++)
				resultQP += D2UhAtQP[iq][i][i];	
			result[iq] += feval * resultQP;
		}
	}
48
}
49
50
51
52
53
54
55
56
void Phase_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
			std::vector<WorldVector<double> > &result)
{
	int nPoints = grdUhAtQP.size();
	for (int iq = 0; iq < nPoints; iq++) {
		double factor = f(iq) * phase[iq] * fac;
		axpy(factor, grdUhAtQP[iq], result[iq]);
	}
57
}
58
59
60
double Phase_SOT::f(const int iq) const
{
	return 1.0;
61
}
62
63
64
65
66
67
68
69

/* ----------------------------------------------------------- */

Phase_FOT::Phase_FOT(DOFVectorBase<double>* phaseDV_,double fac_) :
  FirstOrderTerm(6), phaseDV(phaseDV_), fac(fac_)
{
  TEST_EXIT(phaseDV_)("no vector phase!\n");
  auxFeSpaces.insert(phaseDV_->getFeSpace());
70
}
71
72
73
74
void Phase_FOT::initElement(const ElInfo* elInfo, 
      SubAssembler* subAssembler,
      Quadrature *quad) {
  getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
75
}
76
77
78
79
void Phase_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
      SubAssembler* subAssembler,
      Quadrature *quad) {
  getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
80
}
81
82
83
84
85
86
87
88
89
void Phase_FOT::getLb(const ElInfo *elInfo, 
            vector<mtl::dense_vector<double> >& result) const
{
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(result.size());
    
    for (int iq = 0; iq < nPoints; iq++) {
        lb(Lambda, f(iq), result[iq], phase[iq]*fac);
    }
90
}
91
92
93
94
95
96
97
98
99
100
101
102
void Phase_FOT::eval(int nPoints,
            const mtl::dense_vector<double>& uhAtQP,
            const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
            const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
            mtl::dense_vector<double>& result,
            double opFactor)
{   
    double factor = fac * opFactor;
    for (int iq = 0; iq < nPoints; iq++) {
        double resultQP = grdUhAtQP[iq] * f(iq);
        result[iq] += phase[iq] * factor * resultQP;
    }
103
}
104
105
106
107
108
WorldVector<double> Phase_FOT::f(const int iq) const
{
  WorldVector<double> result;
  result.set(1.0);
  return result;
109
}
110
111
112
113
114
115

/* ----------------------------------------------------------- */

Phase_ZOT::Phase_ZOT(DOFVectorBase<double>* phaseDV_, double fac_) : ZeroOrderTerm(6),
	phaseDV(phaseDV_), fac(fac_)
{	
Praetorius, Simon's avatar
Praetorius, Simon committed
116
  TEST_EXIT(phaseDV_)("No vector Phase!\n");
117
	
Praetorius, Simon's avatar
Praetorius, Simon committed
118
  auxFeSpaces.insert(phaseDV_->getFeSpace());
119
}
120
121
122
123
void Phase_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
Praetorius, Simon's avatar
Praetorius, Simon committed
124
  getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
125
}
126
127
128
129
130
void Phase_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
Praetorius, Simon's avatar
Praetorius, Simon committed
131
  getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
132
}
133
134
void Phase_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{ 
Praetorius, Simon's avatar
Praetorius, Simon committed
135
136
  for (int iq = 0; iq < nPoints; iq++)
    C[iq] += f(iq) * phase[iq] * fac;
137
}
138
139
140
141
142
143
144
void Phase_ZOT::eval(int nPoints,
	const mtl::dense_vector<double>& uhAtQP,
	const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
	const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
	mtl::dense_vector<double>& result,
	double opFactor)
{
Praetorius, Simon's avatar
Praetorius, Simon committed
145
146
147
  double factor = opFactor * fac;
  for (int iq = 0; iq < nPoints; iq++)
    result[iq] += factor * f(iq) * phase[iq] * uhAtQP[iq];
148
}
149
150
double Phase_ZOT::f(const int iq) const
{
Praetorius, Simon's avatar
Praetorius, Simon committed
151
  return 1.0;
152
}
153
154
155
156
157
158
159
160
161
162
163

/* ----------------------------------------------------------- */

PhaseInverse_ZOT::PhaseInverse_ZOT(DOFVectorBase<double>* phaseDV_, double fac_) : ZeroOrderTerm(6),
	phaseDV(phaseDV_), fac(fac_), component(0)
{	
	TEST_EXIT(phaseDV_)("No vector Phase!\n");
	facVec = new WorldVector<double>();
	facVec->set(1.0);

	auxFeSpaces.insert(phaseDV_->getFeSpace());
164
}
165
166
167
168
169
170
PhaseInverse_ZOT::PhaseInverse_ZOT(DOFVectorBase<double>* phaseDV_, double fac_, WorldVector<double> *facVec_, int component_) : ZeroOrderTerm(6),
	phaseDV(phaseDV_), fac(fac_), component(component_), facVec(facVec_)
{	
	TEST_EXIT(phaseDV_)("No vector Phase!\n");
	
	auxFeSpaces.insert(phaseDV_->getFeSpace());
171
}
172
173
174
175
176
void PhaseInverse_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
	getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
177
}
178
179
180
181
182
183
void PhaseInverse_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
	getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
184
}
185
186
187
188
189
void PhaseInverse_ZOT::getC(const ElInfo *, int nPoints, ElementVector &C)
{ 
	double factor = fac * (*facVec)[component];
	for (int iq = 0; iq < nPoints; iq++)
		C[iq] += f(iq) * factor;
190
}
191
192
193
194
195
196
197
198
199
200
void PhaseInverse_ZOT::eval(int nPoints,
				const mtl::dense_vector<double>& uhAtQP,
				const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
				const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
				mtl::dense_vector<double>& result,
				double opFactor) 
{
	double factor = opFactor * fac * (*facVec)[component];
	for (int iq = 0; iq < nPoints; iq++)
		result[iq] += factor * f(iq) * uhAtQP[iq];
201
}
202
203
204
double PhaseInverse_ZOT::f(const int iq) const
{
	return std::max(0.0, std::min( 1.0, (1.0-phase[iq]) ) );
205
}
206
207
208
209
210
211
212
213

/* ----------------------------------------------------------- */

Pow3_ZOT::Pow3_ZOT(DOFVectorBase<double> *rhoDV_, double fac_) :
	ZeroOrderTerm(3), rhoDV(rhoDV_), fac(fac_)
{
	TEST_EXIT(rhoDV_)("No vector rho!\n");
	auxFeSpaces.insert(rhoDV_->getFeSpace());
214
}
215
216
217
218
219
220

void Pow3_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
	getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
221
}
222
223
224
225
226
227
void Pow3_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
	getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
228
}
229
230
231
232
void Pow3_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{ 
	for (int iq = 0; iq < nPoints; iq++)
		C[iq] += f(iq) * fac;
233
}
234
235
236
237
238
239
240
241
242
243
void Pow3_ZOT::eval(int nPoints,
	const mtl::dense_vector<double>& uhAtQP,
	const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
	const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
	mtl::dense_vector<double>& result,
	double opFactor)
{
	double factor = opFactor * fac;
	for (int iq = 0; iq < nPoints; iq++)
		result[iq] += factor * f(iq) * uhAtQP[iq];
244
}
245
246
247
double Pow3_ZOT::f(const int iq) const
{
	return rho[iq]*sqr(rho[iq]);
248
}
249
250
251
252
253
254
255
256

/* ----------------------------------------------------------- */

Pow3Phase_ZOT::Pow3Phase_ZOT(DOFVectorBase<double>* phaseDV, DOFVectorBase<double> *rhoDV_, double fac_) :
	Phase_ZOT(phaseDV, fac_), rhoDV(rhoDV_)
{
	TEST_EXIT(rhoDV_)("No vector rho!\n");
	auxFeSpaces.insert(rhoDV_->getFeSpace());
257
}
258
259
260
261
262
263
264

void Pow3Phase_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
	Phase_ZOT::initElement(elInfo,subAssembler,quad);
	getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
265
}
266
267
268
269
270
271
272
void Pow3Phase_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
	Phase_ZOT::initElement(largeElInfo,smallElInfo,subAssembler,quad);
	getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
273
}
274
275
276
double Pow3Phase_ZOT::f(const int iq) const
{
	return rho[iq]*sqr(rho[iq]);
277
}
278
279
280
281
282
283
284
285

/* ----------------------------------------------------------- */

Pow2_ZOT::Pow2_ZOT(DOFVectorBase<double> *rhoDV_, double fac_) :
	ZeroOrderTerm(2), rhoDV(rhoDV_), fac(fac_)
{
	TEST_EXIT(rhoDV_)("No vector rho!\n");
	auxFeSpaces.insert(rhoDV_->getFeSpace());
286
}
287
288
289
290
291
292

void Pow2_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
	getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
293
}
294
295
296
297
298
299
void Pow2_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
	getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
300
}
301
302
303
304
void Pow2_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{ 
	for (int iq = 0; iq < nPoints; iq++)
		C[iq] += f(iq) * fac;
305
}
306
307
308
309
310
311
312
313
314
315
void Pow2_ZOT::eval(int nPoints,
	const mtl::dense_vector<double>& uhAtQP,
	const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
	const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
	mtl::dense_vector<double>& result,
	double opFactor)
{
	double factor = opFactor * fac;
	for (int iq = 0; iq < nPoints; iq++)
		result[iq] += factor * f(iq) * uhAtQP[iq];
316
}
317
318
319
double Pow2_ZOT::f(const int iq) const
{
	return sqr(rho[iq]);
320
}
321
322
323
324
325
326
327
328

/* ----------------------------------------------------------- */

Pow2Phase_ZOT::Pow2Phase_ZOT(DOFVectorBase<double>* phaseDV, DOFVectorBase<double> *rhoDV_, double fac_) :
	Phase_ZOT(phaseDV, fac_), rhoDV(rhoDV_)
{
	TEST_EXIT(rhoDV_)("No vector rho!\n");
	auxFeSpaces.insert(rhoDV_->getFeSpace());
329
}
330
331
332
333
334
335
336

void Pow2Phase_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
	Phase_ZOT::initElement(elInfo,subAssembler,quad);
	getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
337
}
338
339
340
341
342
343
344
void Pow2Phase_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
	Phase_ZOT::initElement(largeElInfo,smallElInfo,subAssembler,quad);
	getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
345
}
346
347
348
double Pow2Phase_ZOT::f(const int iq) const
{
	return sqr(rho[iq]);
349
}
350
351
352
353
354
355
356
357
358
/* ----------------------------------------------------------- 
   1/(rho+0.9)
*/

ConstrainedFrac_ZOT::ConstrainedFrac_ZOT(DOFVectorBase<double> *phaseDV_, DOFVectorBase<double> *rhoDV_, double fac_) :
	Phase_ZOT(phaseDV_,fac_), rhoDV(rhoDV_)
{
	TEST_EXIT(rhoDV_)("No vector rho!\n");
	auxFeSpaces.insert(rhoDV_->getFeSpace());
359
}
360
361
362
363
364
365
366

void ConstrainedFrac_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
	Phase_ZOT::initElement(elInfo, subAssembler,quad);
	getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
367
}
368
369
370
371
372
373
374
void ConstrainedFrac_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
	Phase_ZOT::initElement(largeElInfo, smallElInfo, subAssembler,quad);
	getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
375
}
376
377
378
double ConstrainedFrac_ZOT::f(const int iq) const
{
	return 1.0/std::max(1.e-6, rho[iq]+0.9);
379
}
380
381
382
383
384
385
386
387
388
389

/* ----------------------------------------------------------- 
   -1/(rho+0.9)^2
*/

ConstrainedFracSqr_ZOT::ConstrainedFracSqr_ZOT(DOFVectorBase<double> *phaseDV_, DOFVectorBase<double> *rhoDV_, double factor_) :
	Phase_ZOT(phaseDV_,factor_), rhoDV(rhoDV_)
{
	TEST_EXIT(rhoDV_)("No vector rho!\n");
	auxFeSpaces.insert(rhoDV_->getFeSpace());
390
}
391
392
393
394
395
396
397

void ConstrainedFracSqr_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
	Phase_ZOT::initElement(elInfo, subAssembler,quad);
	getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
398
}
399
400
401
402
403
404
405
void ConstrainedFracSqr_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
	Phase_ZOT::initElement(largeElInfo, smallElInfo, subAssembler,quad);
	getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
406
}
407
408
409
double ConstrainedFracSqr_ZOT::f(const int iq) const
{
	return -1.0/std::max(1.e-6, sqr(rho[iq]+0.9));
410
}
411
412
413
414
415
416
417


ConstrainedPowImpl_ZOT::ConstrainedPowImpl_ZOT(DOFVectorBase<double> *rhoDV_, double p_, double rhoMin_, double fac_) :
    ZeroOrderTerm(2), rhoDV(rhoDV_), p(p_), rhoMin(rhoMin_), fac(fac_)
{
    TEST_EXIT(rhoDV_)("No vector rho!\n");
    auxFeSpaces.insert(rhoDV_->getFeSpace());
418
}
419
420
421
422
423
424

void ConstrainedPowImpl_ZOT::initElement(const ElInfo* elInfo, 
            SubAssembler* subAssembler,
            Quadrature *quad) 
{
    getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
425
}
426
427
428
429
430
431
void ConstrainedPowImpl_ZOT::initElement(const ElInfo* largeElInfo,
            const ElInfo* smallElInfo,
            SubAssembler* subAssembler,
            Quadrature *quad)
{
    getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
432
}
433
434
435
436
void ConstrainedPowImpl_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{ 
    for (int iq = 0; iq < nPoints; iq++)
        C[iq] += f(iq)*fac;
437
}
438
439
440
441
442
443
444
445
446
447
void ConstrainedPowImpl_ZOT::eval(int nPoints,
    const mtl::dense_vector<double>& uhAtQP,
    const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
    const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
    mtl::dense_vector<double>& result,
    double opFactor)
{
    double factor = opFactor*fac;
    for (int iq = 0; iq < nPoints; iq++)
        result[iq] += factor * f(iq) * uhAtQP[iq];
448
}
449
450
451
double ConstrainedPowImpl_ZOT::f(const int iq) const
{
    return 2.0*p*Helpers::toInt(0.0<=rhoMin-rho[iq]);
452
}
453
454
455
456
457
458
459


ConstrainedPowExpl_ZOT::ConstrainedPowExpl_ZOT(DOFVectorBase<double> *rhoDV_, double p_, double rhoMin_, double fac_) :
    ZeroOrderTerm(2), rhoDV(rhoDV_), p(p_), rhoMin(rhoMin_), fac(fac_)
{
    TEST_EXIT(rhoDV_)("No vector rho!\n");
    auxFeSpaces.insert(rhoDV_->getFeSpace());
460
}
461
462
463
464
465
466

void ConstrainedPowExpl_ZOT::initElement(const ElInfo* elInfo, 
            SubAssembler* subAssembler,
            Quadrature *quad) 
{
    getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
467
}
468
469
470
471
472
473
void ConstrainedPowExpl_ZOT::initElement(const ElInfo* largeElInfo,
            const ElInfo* smallElInfo,
            SubAssembler* subAssembler,
            Quadrature *quad)
{
    getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
474
}
475
476
477
478
void ConstrainedPowExpl_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{ 
    for (int iq = 0; iq < nPoints; iq++)
        C[iq] += f(iq)*fac;
479
}
480
481
482
483
484
485
486
487
488
489
void ConstrainedPowExpl_ZOT::eval(int nPoints,
    const mtl::dense_vector<double>& uhAtQP,
    const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
    const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
    mtl::dense_vector<double>& result,
    double opFactor)
{
    double factor = opFactor*fac;
    for (int iq = 0; iq < nPoints; iq++)
        result[iq] += factor * f(iq) * uhAtQP[iq];
490
}
491
492
493
double ConstrainedPowExpl_ZOT::f(const int iq) const
{
    return -2.0*p*std::max(0.0, rhoMin-rho[iq]);
494
}
495
496
497
498
499
500

ConstrainedPow_ZOT::ConstrainedPow_ZOT(DOFVectorBase<double> *rhoDV_, double p_, double rhoMin_, double fac_) :
    ZeroOrderTerm(2), rhoDV(rhoDV_), p(p_), rhoMin(rhoMin_), fac(fac_)
{
    TEST_EXIT(rhoDV_)("No vector rho!\n");
    auxFeSpaces.insert(rhoDV_->getFeSpace());
501
}
502
503
504
505
506
507

void ConstrainedPow_ZOT::initElement(const ElInfo* elInfo, 
            SubAssembler* subAssembler,
            Quadrature *quad) 
{
    getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
508
}
509
510
511
512
513
514
void ConstrainedPow_ZOT::initElement(const ElInfo* largeElInfo,
            const ElInfo* smallElInfo,
            SubAssembler* subAssembler,
            Quadrature *quad)
{
    getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
515
}
516
517
518
519
void ConstrainedPow_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{ 
    for (int iq = 0; iq < nPoints; iq++)
        C[iq] += f(iq)*fac;
520
}
521
522
523
524
525
526
527
528
529
530
void ConstrainedPow_ZOT::eval(int nPoints,
    const mtl::dense_vector<double>& uhAtQP,
    const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
    const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
    mtl::dense_vector<double>& result,
    double opFactor)
{
    double factor = opFactor*fac;
    for (int iq = 0; iq < nPoints; iq++)
        result[iq] += factor * f(iq) * uhAtQP[iq];
531
}
532
533
534
double ConstrainedPow_ZOT::f(const int iq) const
{
    return -2.0*p*std::max(0.0, rhoMin-rho[iq]);
535
}
536
537
538
539
540
541
542
543
544
545

/* -----------------------------------------------------------
   phi*rho
   --------- */

FactorPhase_ZOT::FactorPhase_ZOT(DOFVectorBase<double> *phaseDV_, DOFVectorBase<double> *rhoDV_, double fac_) :
	Phase_ZOT(phaseDV_,fac_), rhoDV(rhoDV_)
{
	TEST_EXIT(rhoDV_)("No vector rho!\n");
	auxFeSpaces.insert(rhoDV_->getFeSpace());
546
}
547
548
549
550
551
552
553

void FactorPhase_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
	Phase_ZOT::initElement(elInfo, subAssembler,quad);
	getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
554
}
555
556
557
558
559
560
561
void FactorPhase_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
	Phase_ZOT::initElement(largeElInfo, smallElInfo, subAssembler,quad);
	getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
562
}
563
564
565
double FactorPhase_ZOT::f(const int iq) const
{
	return rho[iq];
566
}
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585

/* ----------------------------------------------------------- */

WorldVecAndVecFct_FOT::WorldVecAndVecFct_FOT(WorldVector<DOFVectorBase<double>*> vecs_, DOFVectorBase<double> *phaseDV_, double fac_)
	: FirstOrderTerm(2), phaseDV(phaseDV_), fac(fac_)
{
	numVecs=vecs_.size();
	TEST_EXIT(numVecs==2 || numVecs==3)("Only Dim=2 or Dim=3 possible\n");
	
	TEST_EXIT(phaseDV_)("phaseDV is NULL!\n");
	for (int i = 0; i < numVecs; i++) {
		TEST_EXIT(vecs_[i])("One vector is NULL!\n");
	
		auxFeSpaces.insert(vecs_[i]->getFeSpace());
	}
	
	vec0DV=vecs_[0];
	vec1DV=vecs_[1];
	if(numVecs>=3) vec2DV=vecs_[2];
586
}
587
588
589
590
591
592
593
594
void WorldVecAndVecFct_FOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
	getVectorAtQPs(vec0DV, elInfo, subAssembler, quad, vec0);
	getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
	if(numVecs>=3) getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
	getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
595
}
596
597
598
599
600
601
602
603
void WorldVecAndVecFct_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
	getVectorAtQPs(vec0DV, smallElInfo, largeElInfo, subAssembler, quad, vec0);
	getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
	if(numVecs>=3) getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
	getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
604
}
605
606
607
608
609
610
611
612
613
614
615
616
617
void WorldVecAndVecFct_FOT::getLb(const ElInfo *elInfo, 
			vector<mtl::dense_vector<double> >& result) const
{
	const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(result.size());
	
	for (int iq = 0; iq < nPoints; iq++) {		
		WorldVector<double> vec;
		vec[0]=vec0[iq];
		vec[1]=vec1[iq];
		if(numVecs>=3) vec[2]=vec2[iq];
		lb(Lambda, vec, result[iq], phase[iq]*fac);
	}
618
}
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
void WorldVecAndVecFct_FOT::eval(int nPoints,
			const mtl::dense_vector<double>& uhAtQP,
			const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			mtl::dense_vector<double>& result,
			double opFactor)
{	
	double factor = fac * opFactor;
	for (int iq = 0; iq < nPoints; iq++) {
		double resultQP = 0.0;
				
		resultQP += grdUhAtQP[iq][0] * vec0[iq];
		resultQP += grdUhAtQP[iq][1] * vec1[iq];
		if(numVecs>=3) resultQP += grdUhAtQP[iq][2] * vec2[iq];
		
		result[iq] += phase[iq] * factor * resultQP;
	}
636
}
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653

/* ----------------------------------------------------------- */

WorldVec_FOT::WorldVec_FOT(WorldVector<DOFVector<double>*> vecs_, double fac_)
	: FirstOrderTerm(2), fac(fac_)
{
	numVecs=vecs_.size();
	TEST_EXIT(numVecs==2 || numVecs==3)("Only Dim=2 or Dim=3 possible\n");
	
	for (int i = 0; i < numVecs; i++) {
		TEST_EXIT(vecs_[i])("One vector is NULL!\n");	
		auxFeSpaces.insert(vecs_[i]->getFeSpace());
	}
	
	vec0DV=vecs_[0];
	vec1DV=vecs_[1];
	if(numVecs>=3) vec2DV=vecs_[2];
654
}
655
656
657
658
659
660
661
void WorldVec_FOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
	getVectorAtQPs(vec0DV, elInfo, subAssembler, quad, vec0);
	getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
	if(numVecs>=3) getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
662
}
663
664
665
666
667
668
669
void WorldVec_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
	getVectorAtQPs(vec0DV, smallElInfo, largeElInfo, subAssembler, quad, vec0);
	getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
	if(numVecs>=3) getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
670
}
671
672
673
674
675
676
677
678
679
680
681
682
683
void WorldVec_FOT::getLb(const ElInfo *elInfo, 
			vector<mtl::dense_vector<double> >& result) const
{
	const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(result.size());
	
	for (int iq = 0; iq < nPoints; iq++) {
		WorldVector<double> vec;
		vec[0]=vec0[iq];
		vec[1]=vec1[iq];
		if(numVecs>=3) vec[2]=vec2[iq];
		lb(Lambda, vec, result[iq], fac);
	}
684
}
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
void WorldVec_FOT::eval(int nPoints,
			const mtl::dense_vector<double>& uhAtQP,
			const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			mtl::dense_vector<double>& result,
			double opFactor)
{	
	double factor = fac * opFactor;
	for (int iq = 0; iq < nPoints; iq++) {
		double resultQP = 0.0;
				
		resultQP += grdUhAtQP[iq][0] * vec0[iq];
		resultQP += grdUhAtQP[iq][1] * vec1[iq];
		if(numVecs>=3) resultQP += grdUhAtQP[iq][2] * vec2[iq];
		
		result[iq] += factor * resultQP;
	}
702
}
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777

/* ----------------------------------------------------------- */

WorldVector_FOT::WorldVector_FOT(DOFVector<WorldVector<double> >* dv, double fac_)
: FirstOrderTerm(1),
  vec(dv),
  fac(fac_)
{
  TEST_EXIT(dv)("No vector!\n");

  auxFeSpaces.insert(dv->getFeSpace());
}


void WorldVector_FOT::initElement(const ElInfo* elInfo,
				  SubAssembler* subAssembler,
				  Quadrature *quad)
{
  getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
}


void WorldVector_FOT::initElement(const ElInfo* largeElInfo,
				  const ElInfo* smallElInfo,
				  SubAssembler* subAssembler,
				  Quadrature *quad)
{
  getVectorAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad, vecAtQPs);
}


void WorldVector_FOT::getLb(const ElInfo *elInfo,
			    vector<mtl::dense_vector<double> >& Lb) const
{
  const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(Lb.size());

  for (int iq = 0; iq < nPoints; iq++)
    lb(grdLambda, vecAtQPs[iq], Lb[iq], fac);
}


void WorldVector_FOT::eval(int nPoints,
			   const mtl::dense_vector<double>&,
			   const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			   const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			   mtl::dense_vector<double>& result,
		 	   double opFactor)
{
  if (num_rows(grdUhAtQP) > 0) {
    double fac_ = fac * opFactor;
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += (vecAtQPs[iq] * grdUhAtQP[iq]) * fac_;
  }
}

/* ----------------------------------------------------------- */

WorldVecPhase_FOT::WorldVecPhase_FOT(DOFVectorBase<double> *phaseDV_, WorldVector<DOFVector<double>*> vecs_,double fac_)
    : FirstOrderTerm(2), phaseDV(phaseDV_), fac(fac_)
{
    numVecs=vecs_.size();
    TEST_EXIT(numVecs==2 || numVecs==3)("Only Dim=2 or Dim=3 possible\n");
    
    TEST_EXIT(phaseDV_)("phaseDV is NULL!\n");
    for (int i = 0; i < numVecs; i++) {
        TEST_EXIT(vecs_[i])("One vector is NULL!\n");
    
        auxFeSpaces.insert(vecs_[i]->getFeSpace());
    }
    auxFeSpaces.insert(phaseDV_->getFeSpace());
    
    vec0DV=vecs_[0];
    vec1DV=vecs_[1];
    if(numVecs>=3) vec2DV=vecs_[2];
778
}
779
780
781
782
783
784
785
786
void WorldVecPhase_FOT::initElement(const ElInfo* elInfo, 
                    SubAssembler* subAssembler,
                    Quadrature *quad) 
{
    getVectorAtQPs(vec0DV, elInfo, subAssembler, quad, vec0);
    getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
    if(numVecs>=3) getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
    getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
787
}
788
789
790
791
792
793
794
795
void WorldVecPhase_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
                    SubAssembler* subAssembler,
                    Quadrature *quad) 
{
    getVectorAtQPs(vec0DV, smallElInfo, largeElInfo, subAssembler, quad, vec0);
    getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
    if(numVecs>=3) getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
    getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
796
}
797
798
799
800
801
802
803
804
805
806
807
808
809
void WorldVecPhase_FOT::getLb(const ElInfo *elInfo, 
            vector<mtl::dense_vector<double> >& result) const
{
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(result.size());
    
    for (int iq = 0; iq < nPoints; iq++) {      
        WorldVector<double> vec;
        vec[0]=vec0[iq];
        vec[1]=vec1[iq];
        if(numVecs>=3) vec[2]=vec2[iq];
        lb(Lambda, vec, result[iq], phase[iq]*fac);
    }
810
}
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
void WorldVecPhase_FOT::eval(int nPoints,
            const mtl::dense_vector<double>& uhAtQP,
            const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
            const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
            mtl::dense_vector<double>& result,
            double opFactor)
{   
    double factor = fac * opFactor;
    for (int iq = 0; iq < nPoints; iq++) {
        double resultQP = 0.0;
                
        resultQP += grdUhAtQP[iq][0] * vec0[iq];
        resultQP += grdUhAtQP[iq][1] * vec1[iq];
        if(numVecs>=3) resultQP += grdUhAtQP[iq][2] * vec2[iq];
        
        result[iq] += phase[iq] * factor * resultQP;
    }
828
}
829

Praetorius, Simon's avatar
Praetorius, Simon committed
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850

/* ----------------------------------------------------------- */

VecAndWorldVec_FOT::VecAndWorldVec_FOT(DOFVector<double>* vec0_, 
									   WorldVector<DOFVector<double>*> vecs_, 
									   AbstractFunction<double, double>* fct_, double fac_)
	: FirstOrderTerm(2), vDV(vec0_), fct(fct_), fac(fac_)
{
  numVecs = vecs_.size();
  TEST_EXIT(numVecs == 2 || numVecs == 3)
	("Only Dim=2 or Dim=3 possible\n");
  
  for (int i = 0; i < numVecs; i++) {
	TEST_EXIT(vecs_[i])("One vector is NULL!\n");	
	auxFeSpaces.insert(vecs_[i]->getFeSpace());
  }
  
  vec0DV = vecs_[0];
  vec1DV = vecs_[1];
  if (numVecs >= 3)
	vec2DV = vecs_[2];
851
}
Praetorius, Simon's avatar
Praetorius, Simon committed
852
853
854
855
856
857
858
859
860
void VecAndWorldVec_FOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
  getVectorAtQPs(vDV, elInfo, subAssembler, quad, v);
  getVectorAtQPs(vec0DV, elInfo, subAssembler, quad, vec0);
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
  if (numVecs >= 3)
	getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
861
}
Praetorius, Simon's avatar
Praetorius, Simon committed
862
863
864
865
866
867
868
869
870
void VecAndWorldVec_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
  getVectorAtQPs(vDV, smallElInfo, largeElInfo, subAssembler, quad, v);
  getVectorAtQPs(vec0DV, smallElInfo, largeElInfo, subAssembler, quad, vec0);
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
  if (numVecs >= 3)
	getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
871
}
Praetorius, Simon's avatar
Praetorius, Simon committed
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
void VecAndWorldVec_FOT::getLb(const ElInfo *elInfo, 
			vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
  
  for (int iq = 0; iq < nPoints; iq++) {
	double f = (*fct)(v[iq]);
	WorldVector<double> vec;
	vec[0] = vec0[iq];
	vec[1] = vec1[iq];
	if (numVecs >= 3)
	  vec[2] = vec2[iq];
	lb(Lambda, f*vec, result[iq], fac);
  }
887
}
Praetorius, Simon's avatar
Praetorius, Simon committed
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
void VecAndWorldVec_FOT::eval(int nPoints,
			const mtl::dense_vector<double>& uhAtQP,
			const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			mtl::dense_vector<double>& result,
			double opFactor)
{	
  double factor = fac * opFactor;
  for (int iq = 0; iq < nPoints; iq++) {
	double resultQP = 0.0;
			
	resultQP += grdUhAtQP[iq][0] * vec0[iq];
	resultQP += grdUhAtQP[iq][1] * vec1[iq];
	if (numVecs >= 3)
	  resultQP += grdUhAtQP[iq][2] * vec2[iq];
	
	result[iq] += factor * resultQP * (*fct)(v[iq]);
  }
906
}
Praetorius, Simon's avatar
Praetorius, Simon committed
907
908
909



910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
/* ----------------------------------------------------------- */

// vec1*vec2*d/dxi(...)
VecAndPartialDerivative_FOT::VecAndPartialDerivative_FOT(DOFVectorBase<double> *dv1,
      int component_, double fac_)
: FirstOrderTerm(3), vec1DV(dv1), fac(fac_), component(component_)
{
  auxFeSpaces.insert(vec1DV->getFeSpace());

  setB(component);
}

void VecAndPartialDerivative_FOT::initElement(const ElInfo* elInfo, 
        SubAssembler* subAssembler,
        Quadrature *quad) 
{ 
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
927
}
928
929
930
931
932
void VecAndPartialDerivative_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
  SubAssembler* subAssembler,
  Quadrature *quad)
{ 
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
933
}
934
935
936
937
938
939
940
941
942
void VecAndPartialDerivative_FOT::getLb(const ElInfo *elInfo,
      vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(result.size());
    
  for(int iq = 0; iq < nPoints; iq++) {
    lb_one(Lambda, result[iq], vec1[iq] * fac);
  }
943
}
944
945
946
947
948
949
950
951
952
void VecAndPartialDerivative_FOT::eval(int nPoints,
      const mtl::dense_vector<double>& uhAtQP,
      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
      mtl::dense_vector<double>& result,
      double opFactor)
{
  for(int iq = 0; iq < nPoints; iq++) 
    result[iq] += opFactor * vec1[iq] * grdUhAtQP[iq][component] * fac;
953
}
954

Praetorius, Simon's avatar
Praetorius, Simon committed
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
/* ----------------------------------------------------------- */

// vec1*vec2*d/dxi(...)
VecAndPartialDerivativeIJ_FOT::VecAndPartialDerivativeIJ_FOT(DOFVectorBase<double> *dv1, 
															 DOFVectorBase<double> *dv2,
															 int i_,
															 int j_,
															 double fac_)
: FirstOrderTerm(4), vec1DV(dv1), vec2DV(dv2), fct(NULL), fac(fac_), xi(i_), xj(j_)
{
  auxFeSpaces.insert(vec1DV->getFeSpace());
  auxFeSpaces.insert(vec2DV->getFeSpace());

  setB(xi);
}

VecAndPartialDerivativeIJ_FOT::VecAndPartialDerivativeIJ_FOT(DOFVectorBase<double> *dv1, 
							  DOFVectorBase<double> *dv2,
							  int i_,
							  int j_,
							  BinaryAbstractFunction<double, double, double>* fct_, // = f(v1, d_j(v2))
							  double fac_)
: FirstOrderTerm(4), vec1DV(dv1), vec2DV(dv2), fct(fct_), fac(fac_), xi(i_), xj(j_)
{
  auxFeSpaces.insert(vec1DV->getFeSpace());
  auxFeSpaces.insert(vec2DV->getFeSpace());

  setB(xi);
}

void VecAndPartialDerivativeIJ_FOT::initElement(const ElInfo* elInfo, 
        SubAssembler* subAssembler,
        Quadrature *quad) 
{ 
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
  getGradientsAtQPs(vec2DV, elInfo, subAssembler, quad, grad);
}

void VecAndPartialDerivativeIJ_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
  SubAssembler* subAssembler,
  Quadrature *quad)
{ 
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
  getGradientsAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, grad);
}

void VecAndPartialDerivativeIJ_FOT::getLb(const ElInfo *elInfo,
      vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
  
  if (fct == NULL) {
	for(int iq = 0; iq < nPoints; iq++) {
	  lb_one(Lambda, result[iq], vec1[iq] * grad[iq][xj] * fac);
	}
  } else {
	for(int iq = 0; iq < nPoints; iq++) {
	  lb_one(Lambda, result[iq], (*fct)(vec1[iq], grad[iq][xj]) * fac);
	}
  }
}

void VecAndPartialDerivativeIJ_FOT::eval(int nPoints,
      const mtl::dense_vector<double>& uhAtQP,
      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
      mtl::dense_vector<double>& result,
      double opFactor)
{
  if (fct == NULL) {
	for(int iq = 0; iq < nPoints; iq++) 
	  result[iq] += opFactor * vec1[iq] * grad[iq][xj] * grdUhAtQP[iq][xi] * fac;
  } else {
	for(int iq = 0; iq < nPoints; iq++) 
	  result[iq] += opFactor * (*fct)(vec1[iq], grad[iq][xj]) * grdUhAtQP[iq][xi] * fac;
  }
}


/* ----------------------------------------------------------- */

// vec1*vec2*d/dxi(...)
Vec2AndPartialDerivative_FOT::Vec2AndPartialDerivative_FOT(
			       DOFVectorBase<double> *dv1, 
			       DOFVectorBase<double> *dv2,
			       int component_, 
			       double fac_)
: FirstOrderTerm(3), vec1DV(dv1), vec2DV(dv2), fct(NULL), fac(fac_), component(component_)
{
  auxFeSpaces.insert(vec1DV->getFeSpace());
  auxFeSpaces.insert(vec2DV->getFeSpace());

  setB(component);
}


Vec2AndPartialDerivative_FOT::Vec2AndPartialDerivative_FOT(
			       DOFVectorBase<double> *dv1, 
			       DOFVectorBase<double> *dv2,
			       int component_,
			       BinaryAbstractFunction<double, double, double>* fct_, // = f(v1, v2)
			       double fac_)
: FirstOrderTerm(3), vec1DV(dv1), vec2DV(dv2), fct(fct_), fac(fac_), component(component_)
{
  auxFeSpaces.insert(vec1DV->getFeSpace());
  auxFeSpaces.insert(vec2DV->getFeSpace());

  setB(component);
}

void Vec2AndPartialDerivative_FOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad) 
{ 
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
  getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
}

void Vec2AndPartialDerivative_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
	SubAssembler* subAssembler,
	Quadrature *quad)
{ 
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
  getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
}

void Vec2AndPartialDerivative_FOT::getLb(const ElInfo *elInfo,
			vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
  
  if (fct == NULL) {
    for(int iq = 0; iq < nPoints; iq++)
      lb_one(Lambda, result[iq], vec1[iq] * vec2[iq] * fac);
  } else {
    for(int iq = 0; iq < nPoints; iq++)
      lb_one(Lambda, result[iq], (*fct)(vec1[iq], vec2[iq]) * fac);
  }
1095
}
Praetorius, Simon's avatar
Praetorius, Simon committed
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
void Vec2AndPartialDerivative_FOT::eval(int nPoints,
			const mtl::dense_vector<double>& uhAtQP,
			const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			mtl::dense_vector<double>& result,
			double opFactor)
{
  if (fct == NULL) {
    for(int iq = 0; iq < nPoints; iq++) 
      result[iq] += opFactor * vec1[iq] * vec2[iq] * grdUhAtQP[iq][component] * fac;
  } else {
    for(int iq = 0; iq < nPoints; iq++) 
      result[iq] += opFactor * (*fct)(vec1[iq], vec2[iq]) * grdUhAtQP[iq][component] * fac;
  }
1110
}
Praetorius, Simon's avatar
Praetorius, Simon committed
1111
1112


1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
/* ----------------------------------------------------------- */

// vec1*vec2*d/dxi(...)
PhaseAndViscosityPartialDerivativeBoundary_FOT::PhaseAndViscosityPartialDerivativeBoundary_FOT(DOFVectorBase<double> *phaseDV_, DOFVectorBase<double> *chDV_,
      int component_, int normalComponent_, double viscosity1_, double viscosity2_)
: FirstOrderTerm(3), phaseDV(phaseDV_), chDV(chDV_), component(component_), normalComponent(normalComponent_), v1(viscosity1_), v2(viscosity2_)
{
  auxFeSpaces.insert(phaseDV->getFeSpace());
  auxFeSpaces.insert(chDV->getFeSpace());

  setB(component);
}

void PhaseAndViscosityPartialDerivativeBoundary_FOT::initElement(const ElInfo* elInfo, 
        SubAssembler* subAssembler,
        Quadrature *quad) 
{ 
  getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
  getVectorAtQPs(chDV, elInfo, subAssembler, quad, ch);

  WorldVector<double> normal;
  int side= -1;
  for (int i= 0; i < elInfo->getMesh()->getDim()+1; ++i) {
    if(elInfo->getBoundary(i) > 2) {
      side= i;
      break;
    }
  }
//   elInfo->getNormal(side, normal);
  normalJ = 1.0; //normal[normalComponent];
1143
}
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
void PhaseAndViscosityPartialDerivativeBoundary_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
  SubAssembler* subAssembler,
  Quadrature *quad)
{ 
  getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
  getVectorAtQPs(chDV, smallElInfo, largeElInfo, subAssembler, quad, ch);

  WorldVector<double> normal;
  int side= -1;
  for (int i= 0; i < largeElInfo->getMesh()->getDim()+1; ++i) {
    if(largeElInfo->getBoundary(i) > 2) {
      side= i;
      break;
    }
  }
//   largeElInfo->getNormal(side, normal);
  normalJ = 1.0; //normal[normalComponent];
1161
}
1162
1163
1164
1165
1166
1167
1168
1169
1170
void PhaseAndViscosityPartialDerivativeBoundary_FOT::getLb(const ElInfo *elInfo,
      vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
    
  for(int iq = 0; iq < nPoints; iq++) {
    lb_one(Lambda, result[iq], phase[iq] * normalJ * f(iq));
  }
1171
}
1172
1173
1174
1175
1176
1177
1178
1179
1180
void PhaseAndViscosityPartialDerivativeBoundary_FOT::eval(int nPoints,
      const mtl::dense_vector<double>& uhAtQP,
      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
      mtl::dense_vector<double>& result,
      double opFactor)
{
  for(int iq = 0; iq < nPoints; iq++) 
    result[iq] += opFactor * phase[iq] * grdUhAtQP[iq][component] * normalJ * f(iq);
1181
}
1182
1183
1184
1185
1186
1187

double PhaseAndViscosityPartialDerivativeBoundary_FOT::f(const int iq) const
{
    double ph= std::max(-1.0,std::min(1.0,ch[iq]));
    double viscosity= v1*0.5*(ph+1.0)-v2*0.5*(ph-1.0);
    return viscosity;
1188
}
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217

/* ----------------------------------------------------------- */

// vec1*vec2*d/dxi(...)
ViscosityPartialDerivativeBoundary_FOT::ViscosityPartialDerivativeBoundary_FOT(DOFVectorBase<double> *chDV_,
      int component_, int normalComponent_, double viscosity1_, double viscosity2_)
: FirstOrderTerm(3), chDV(chDV_), component(component_), normalComponent(normalComponent_), v1(viscosity1_), v2(viscosity2_)
{
  auxFeSpaces.insert(chDV->getFeSpace());

  setB(component);
}

void ViscosityPartialDerivativeBoundary_FOT::initElement(const ElInfo* elInfo, 
        SubAssembler* subAssembler,
        Quadrature *quad) 
{ 
  getVectorAtQPs(chDV, elInfo, subAssembler, quad, ch);

  WorldVector<double> normal;
  int side= -1;
  for (int i= 0; i < elInfo->getMesh()->getDim()+1; ++i) {
    if(elInfo->getBoundary(i) > 2) {
      side= i;
      break;
    }
  }
//   elInfo->getNormal(side, normal);
  normalJ = 1.0; //normal[normalComponent];
1218
}
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
void ViscosityPartialDerivativeBoundary_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
  SubAssembler* subAssembler,
  Quadrature *quad)
{ 
  getVectorAtQPs(chDV, smallElInfo, largeElInfo, subAssembler, quad, ch);

  WorldVector<double> normal;
  int side= -1;
  for (int i= 0; i < largeElInfo->getMesh()->getDim()+1; ++i) {
    if(largeElInfo->getBoundary(i) > 2) {
      side= i;
      break;
    }
  }
//   largeElInfo->getNormal(side, normal);
  normalJ = 1.0; //normal[normalComponent];
1235
}
1236
1237
1238
1239
1240
1241
1242
1243
1244
void ViscosityPartialDerivativeBoundary_FOT::getLb(const ElInfo *elInfo,
      vector<mtl::dense_vector<double> >& result) const
{
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
    
  for(int iq = 0; iq < nPoints; iq++) {
    lb_one(Lambda, result[iq], normalJ * f(iq));
  }
1245
}
1246
1247
1248
1249
1250
1251
1252
1253
1254
void ViscosityPartialDerivativeBoundary_FOT::eval(int nPoints,
      const mtl::dense_vector<double>& uhAtQP,
      const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
      const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
      mtl::dense_vector<double>& result,
      double opFactor)
{
  for(int iq = 0; iq < nPoints; iq++) 
    result[iq] += opFactor * grdUhAtQP[iq][component] * normalJ * f(iq);
1255
}
1256
1257
1258
1259
1260
1261

double ViscosityPartialDerivativeBoundary_FOT::f(const int iq) const
{
    double ph= std::max(-1.0,std::min(1.0,ch[iq]));
    double viscosity= v1*0.5*(ph+1.0)-v2*0.5*(ph-1.0);
    return viscosity;
1262
}
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
/* ----------------------------------------------------------- */

// vec1*vec2*d/dxi(...)
Vec2ProductPartial_FOT::Vec2ProductPartial_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
			int component_, double fac_)
: FirstOrderTerm(3), vec1DV(dv1), vec2DV(dv2), fac(fac_), component(component_)
{
	auxFeSpaces.insert(vec1DV->getFeSpace());
	auxFeSpaces.insert(vec2DV->getFeSpace());

	setB(component);
}

void Vec2ProductPartial_FOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad) 
{ 
	getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
	getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
1282
}
1283
1284
1285
1286
1287
1288
void Vec2ProductPartial_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
	SubAssembler* subAssembler,
	Quadrature *quad)
{ 
	getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
	getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
1289
}
1290
1291
1292
1293
1294
1295
1296
1297
1298
void Vec2ProductPartial_FOT::getLb(const ElInfo *elInfo,
			vector<mtl::dense_vector<double> >& result) const
{
	const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(result.size());
    
	for(int iq = 0; iq < nPoints; iq++) {
		lb_one(Lambda, result[iq], vec1[iq] * vec2[iq] * fac);
	}
1299
}
1300
1301
1302
1303
1304
1305
1306
1307
1308
void Vec2ProductPartial_FOT::eval(int nPoints,
			const mtl::dense_vector<double>& uhAtQP,
			const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			mtl::dense_vector<double>& result,
			double opFactor)
{
	for(int iq = 0; iq < nPoints; iq++) 
		result[iq] += opFactor * vec1[iq] * vec2[iq] * grdUhAtQP[iq][component] * fac;
1309
}
1310
1311
1312
1313
1314
1315
1316

/* ----------------------------------------------------------- */

PartialDerivative_FOT::PartialDerivative_FOT(int component_, double fac_)
	: FirstOrderTerm(1), component(component_), fac(fac_)
{
	setB(component);
1317
}
1318
1319
void PartialDerivative_FOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
1320
					Quadrature *quad) {}
1321
1322
void PartialDerivative_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
1323
					Quadrature *quad) {}
1324
1325
1326
1327
1328
1329
1330
1331
1332
void PartialDerivative_FOT::getLb(const ElInfo *elInfo, 
			vector<mtl::dense_vector<double> >& result) const
{
	const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
    const int nPoints = static_cast<int>(result.size());
	
	for (int iq = 0; iq < nPoints; iq++) {	
		lb_one(Lambda, result[iq], fac);
	}
1333
}
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
void PartialDerivative_FOT::eval(int nPoints,
			const mtl::dense_vector<double>& uhAtQP,
			const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
			const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
			mtl::dense_vector<double>& result,
			double opFactor)
{	
	for (int iq = 0; iq < nPoints; iq++) {		
		result[iq] += opFactor * grdUhAtQP[iq][component] * fac;
	}
1344
}
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355

/* ----------------------------------------------------------- */

PartialDerivative_ZOT::PartialDerivative_ZOT(DOFVectorBase<double> *vecDV_, int component_, double fac_)
	: ZeroOrderTerm(1), vecDV(vecDV_), fac(fac_), component(component_)
{
	TEST_EXIT(vecDV_)("No first vector!\n");
	auxFeSpaces.insert(vecDV_->getFeSpace());

	facVec = new WorldVector<double>(); facVec->set(1.0);
	component2 = 0;
1356
}
1357
1358
1359
1360
1361
PartialDerivative_ZOT::PartialDerivative_ZOT(DOFVectorBase<double> *vecDV_, int component_, double fac_, WorldVector<double> *facVec_, int component2_)
	: ZeroOrderTerm(1), vecDV(vecDV_), fac(fac_), facVec(facVec_), component(component_), component2(component2_)
{
	TEST_EXIT(vecDV_)("No first vector!\n");
	auxFeSpaces.insert(vecDV_->getFeSpace());
1362
}
1363
1364
1365
1366
1367
void PartialDerivative_ZOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad)
{
  getGradientsAtQPs(vecDV, elInfo, subAssembler, quad, grad);
1368
}
1369
1370
1371
1372
1373
void PartialDerivative_ZOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad)
{
  getGradientsAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, grad);
1374
}
1375
1376
1377
1378
1379
void PartialDerivative_ZOT::getC(const ElInfo *, int nPoints, ElementVector &C)
{ 	
	double factor = fac * (*facVec)[component2];
	for (int iq = 0; iq < nPoints; iq++)
		C[iq] += grad[iq][component]*factor;
1380
}
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
void PartialDerivative_ZOT::eval(int nPoints,
				const mtl::dense_vector<double>& uhAtQP,
				const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
				const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
				mtl::dense_vector<double>& result,
				double opFactor)
{	
	double factor = opFactor * fac * (*facVec)[component2];
	for (int iq = 0; iq < nPoints; iq++) {		
		result[iq] += grad[iq][component] * uhAtQP[iq] * factor;
	}
1392
}
1393
1394
1395
1396
1397
1398
1399
1400
1401


/* ----------------------------------------------------------- */

VecAndPartialDerivative_ZOT::VecAndPartialDerivative_ZOT(DOFVectorBase<double> *vecDV_, DOFVectorBase<double> *gradDV_, int component_,  double fac_)
: ZeroOrderTerm(2),
  vecDV(vecDV_),
  gradDV(gradDV_),
  component(component_),
Praetorius, Simon's avatar
Praetorius, Simon committed
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
  fct(NULL),
  fac(fac_)
{
  TEST_EXIT(vecDV_)("No value vector!\n");
  TEST_EXIT(gradDV_)("No gradient vector!\n");

  auxFeSpaces.insert(vecDV_->getFeSpace());
  auxFeSpaces.insert(gradDV_->getFeSpace());
}


VecAndPartialDerivative_ZOT::VecAndPartialDerivative_ZOT(DOFVectorBase<double> *vecDV_, 
							DOFVectorBase<double> *gradDV_, 
							int component_,
							AbstractFunction<double, double>* fct_,
							double fac_)
: ZeroOrderTerm(2),
  vecDV(vecDV_),
  gradDV(gradDV_),
  component(component_),
  fct(fct_),
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
  fac(fac_)
{
  TEST_EXIT(vecDV_)("No value vector!\n");
  TEST_EXIT(gradDV_)("No gradient vector!\n");

  auxFeSpaces.insert(vecDV_->getFeSpace());
  auxFeSpaces.insert(gradDV_->getFeSpace());
}


void VecAndPartialDerivative_ZOT::initElement(const ElInfo* elInfo,
					      SubAssembler* subAssembler,
					      Quadrature *quad)
{
  getVectorAtQPs(vecDV, elInfo, subAssembler, quad, vec);
  getGradientsAtQPs(gradDV, elInfo, subAssembler, quad, grad);
}


void VecAndPartialDerivative_ZOT::initElement(const ElInfo* largeElInfo,
					      const ElInfo* smallElInfo,
					      SubAssembler* subAssembler,
					      Quadrature *quad)
{
  getVectorAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, vec);
  getGradientsAtQPs(gradDV, smallElInfo, largeElInfo, subAssembler, quad, grad);
}


void VecAndPartialDerivative_ZOT::getC(const ElInfo *, int nPoints, ElementVector &C)
{
Praetorius, Simon's avatar
Praetorius, Simon committed
1454
1455
1456
1457
1458
1459
1460
1461
  if (fct == NULL) {
	for (int iq = 0; iq < nPoints; iq++)
	  C[iq] += vec[iq] * grad[iq][component] * fac;
  } else {
	for (int iq = 0; iq < nPoints; iq++)
	  C[iq] += (*fct)(vec[iq]) * grad[iq][component] * fac;
  }
	
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
}


void VecAndPartialDerivative_ZOT::eval(int nPoints,
				       const mtl::dense_vector<double>& uhAtQP,
				       const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
				       const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
				       mtl::dense_vector<double>& result,
				       double opFactor)
{
  double fac_ = fac * opFactor;
Praetorius, Simon's avatar
Praetorius, Simon committed
1473
1474
1475
1476
1477
1478
1479
  if (fct == NULL) {
	for (int iq = 0; iq < nPoints; iq++)
	  result[iq] += vec[iq] * grad[iq][component] * uhAtQP[iq] * fac_;
  } else {
	for (int iq = 0; iq < nPoints; iq++)
	  result[iq] += (*fct)(vec[iq]) * grad[iq][component] * uhAtQP[iq] * fac_;
  }
1480
1481
1482
1483
}

/* ----------------------------------------------------------- */

Praetorius, Simon's avatar
Praetorius, Simon committed
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
Vec2AndPartialDerivative_ZOT::Vec2AndPartialDerivative_ZOT(DOFVectorBase<double> *vec1DV_, 
							    DOFVectorBase<double> *vec2DV_, 
							    DOFVectorBase<double> *gradDV_,
							    int component_,  
							    double fac_)
	: ZeroOrderTerm(3), vec1DV(vec1DV_), vec2DV(vec2DV_), gradDV(gradDV_), fct(NULL), fac(fac_), component(component_)
{
  TEST_EXIT(vec1DV_)("No first vector!\n");
  TEST_EXIT(vec2DV_)("No second vector!\n");
  TEST_EXIT(gradDV_)("No gradient vector!\n");
  auxFeSpaces.insert(vec1DV_->getFeSpace());
  auxFeSpaces.insert(vec2DV_->getFeSpace());
  auxFeSpaces.insert(gradDV_->getFeSpace());
1497
}
Praetorius, Simon's avatar
Praetorius, Simon committed
1498
1499
1500
1501
1502
1503
1504
1505

Vec2AndPartialDerivative_ZOT::Vec2AndPartialDerivative_ZOT(DOFVectorBase<double> *vec1DV_, 
							  DOFVectorBase<double> *vec2DV_, 
							  DOFVectorBase<double> *gradDV_,
							  int component_,  
							  BinaryAbstractFunction<double, double, double>* fct_,
							  double fac_)
	: ZeroOrderTerm(3), vec1DV(vec1DV_), vec2DV(vec2DV_), gradDV(gradDV_), fct(fct_), fac(fac_), component(component_)
1506
{
Praetorius, Simon's avatar
Praetorius, Simon committed
1507
1508
1509
1510
1511
1512
  TEST_EXIT(vec1DV_)("No first vector!\n");
  TEST_EXIT(vec2DV_)("No second vector!\n");
  TEST_EXIT(gradDV_)("No gradient vector!\n");
  auxFeSpaces.insert(vec1DV_->getFeSpace());
  auxFeSpaces.insert(vec2DV_->getFeSpace());
  auxFeSpaces.insert(gradDV_->getFeSpace());
1513
}
Praetorius, Simon's avatar
Praetorius, Simon committed
1514

1515
1516
1517
1518
void Vec2AndPartialDerivative_ZOT::initElement(const ElInfo* elInfo,
					SubAssembler* subAssembler,
					Quadrature *quad)
{
Praetorius, Simon's avatar