POperators.cc 43.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
/******************************************************************************
 *
 * Extension of AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: Simon Praetorius et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/


19
20
21
22
23
24
#include "POperators.h"
#include "Helpers.h"

using namespace std;
using namespace AMDiS;

25
Phase_SOT::Phase_SOT(DOFVectorBase<double>* phaseDV_, double fac_)
26
: SecondOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()), 
27
28
  phaseDV(phaseDV_), 
  fac(fac_)
29
{
30
31
  TEST_EXIT(phaseDV_)("no vector phase!\n");
  auxFeSpaces.insert(phaseDV_->getFeSpace());
32
}
33

34
35
void Phase_SOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
36
37
38
			Quadrature *quad)
{
  getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
39
}
40

41
42
void Phase_SOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
43
44
45
			Quadrature *quad)
{
  getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
46
}
47

48
49
50
void Phase_SOT::getLALt(const ElInfo *elInfo,
			vector<mtl::dense2D<double> > &LALt) const
{
51
52
53
54
55
  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);
56
}
57

58
void Phase_SOT::eval(int nPoints,
59
60
61
62
63
		    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)
64
{	
65
66
67
68
69
70
71
72
73
  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;
    }
  }
74
}
75

76
77
78
void Phase_SOT::weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
			std::vector<WorldVector<double> > &result)
{
79
80
81
82
83
  int nPoints = grdUhAtQP.size();
  for (int iq = 0; iq < nPoints; iq++) {
    double factor = f(iq) * phase[iq] * fac;
    axpy(factor, grdUhAtQP[iq], result[iq]);
  }
84
}
85

86
87
double Phase_SOT::f(const int iq) const
{
88
  return 1.0;
89
}
90
91
92

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

93
Phase_FOT::Phase_FOT(DOFVectorBase<double>* phaseDV_, double fac_)
94
: FirstOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()), 
95
96
  phaseDV(phaseDV_), 
  fac(fac_)
97
98
99
{
  TEST_EXIT(phaseDV_)("no vector phase!\n");
  auxFeSpaces.insert(phaseDV_->getFeSpace());
100
}
101

102
103
void Phase_FOT::initElement(const ElInfo* elInfo, 
      SubAssembler* subAssembler,
104
105
      Quadrature *quad)
{
106
  getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
107
}
108

109
110
void Phase_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
      SubAssembler* subAssembler,
111
112
      Quadrature *quad)
{
113
  getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
114
}
115

116
117
118
void Phase_FOT::getLb(const ElInfo *elInfo, 
            vector<mtl::dense_vector<double> >& result) const
{
119
120
121
122
123
  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);
124
}
125

126
void Phase_FOT::eval(int nPoints,
127
128
129
130
131
		    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)
132
{   
133
134
135
136
137
  double factor = fac * opFactor;
  for (int iq = 0; iq < nPoints; iq++) {
    double resultQP = grdUhAtQP[iq] * f(iq);
    result[iq] += phase[iq] * factor * resultQP;
  }
138
}
139

140
141
142
143
144
WorldVector<double> Phase_FOT::f(const int iq) const
{
  WorldVector<double> result;
  result.set(1.0);
  return result;
145
}
146
147
148

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

149
Phase_ZOT::Phase_ZOT(DOFVectorBase<double>* phaseDV_, double fac_)
150
: ZeroOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()),
151
152
  phaseDV(phaseDV_), 
  fac(fac_)
153
{	
Praetorius, Simon's avatar
Praetorius, Simon committed
154
  TEST_EXIT(phaseDV_)("No vector Phase!\n");
155
	
Praetorius, Simon's avatar
Praetorius, Simon committed
156
  auxFeSpaces.insert(phaseDV_->getFeSpace());
157
}
158

159
160
161
162
void Phase_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
Praetorius, Simon's avatar
Praetorius, Simon committed
163
  getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
164
}
165

166
167
168
169
170
void Phase_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
Praetorius, Simon's avatar
Praetorius, Simon committed
171
  getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
172
}
173

174
175
void Phase_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{ 
Praetorius, Simon's avatar
Praetorius, Simon committed
176
177
  for (int iq = 0; iq < nPoints; iq++)
    C[iq] += f(iq) * phase[iq] * fac;
178
}
179

180
void Phase_ZOT::eval(int nPoints,
181
182
183
184
185
		    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)
186
{
Praetorius, Simon's avatar
Praetorius, Simon committed
187
188
189
  double factor = opFactor * fac;
  for (int iq = 0; iq < nPoints; iq++)
    result[iq] += factor * f(iq) * phase[iq] * uhAtQP[iq];
190
}
191

192
193
double Phase_ZOT::f(const int iq) const
{
Praetorius, Simon's avatar
Praetorius, Simon committed
194
  return 1.0;
195
}
196
197
198

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

199
PhaseInverse_ZOT::PhaseInverse_ZOT(DOFVectorBase<double>* phaseDV_, double fac_)
200
: ZeroOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()),
201
202
203
  phaseDV(phaseDV_), 
  fac(fac_), 
  component(0)
204
{	
205
206
207
  TEST_EXIT(phaseDV_)("No vector Phase!\n");
  facVec = new WorldVector<double>();
  facVec->set(1.0);
208

209
  auxFeSpaces.insert(phaseDV_->getFeSpace());
210
}
211
212

PhaseInverse_ZOT::PhaseInverse_ZOT(DOFVectorBase<double>* phaseDV_, double fac_, WorldVector<double> *facVec_, int component_)
213
: ZeroOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()),
214
215
216
217
  phaseDV(phaseDV_), 
  fac(fac_), 
  component(component_), 
  facVec(facVec_)
218
{	
219
220
221
  TEST_EXIT(phaseDV_)("No vector Phase!\n");
  
  auxFeSpaces.insert(phaseDV_->getFeSpace());
222
}
223

224
225
226
227
void PhaseInverse_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
228
  getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
229
}
230

231
232
233
234
235
void PhaseInverse_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
236
  getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
237
}
238

239
240
void PhaseInverse_ZOT::getC(const ElInfo *, int nPoints, ElementVector &C)
{ 
241
242
243
  double factor = fac * (*facVec)[component];
  for (int iq = 0; iq < nPoints; iq++)
    C[iq] += f(iq) * factor;
244
}
245

246
void PhaseInverse_ZOT::eval(int nPoints,
247
248
249
250
251
			    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) 
252
{
253
254
255
  double factor = opFactor * fac * (*facVec)[component];
  for (int iq = 0; iq < nPoints; iq++)
    result[iq] += factor * f(iq) * uhAtQP[iq];
256
}
257

258
259
double PhaseInverse_ZOT::f(const int iq) const
{
260
  return std::max(0.0, std::min( 1.0, (1.0-phase[iq]) ) );
261
}
262
263
264
265

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

WorldVecAndVecFct_FOT::WorldVecAndVecFct_FOT(WorldVector<DOFVectorBase<double>*> vecs_, DOFVectorBase<double> *phaseDV_, double fac_)
266
: FirstOrderTerm(vecs_[0]->getFeSpace()->getBasisFcts()->getDegree() + phaseDV_->getFeSpace()->getBasisFcts()->getDegree()), 
267
268
  phaseDV(phaseDV_), 
  fac(fac_)
269
{
270
271
272
273
274
275
276
277
278
279
280
281
282
  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];
283
}
284

285
286
287
288
void WorldVecAndVecFct_FOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
289
290
291
292
  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);
293
}
294

295
296
297
298
void WorldVecAndVecFct_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
299
300
301
302
  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);
303
}
304

305
306
307
void WorldVecAndVecFct_FOT::getLb(const ElInfo *elInfo, 
			vector<mtl::dense_vector<double> >& result) const
{
308
309
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
310
	
311
312
313
314
315
316
317
  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);
  }
318
}
319

320
321
322
323
324
325
326
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)
{	
327
328
329
330
331
332
333
334
335
336
  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;
  }
337
}
338
339
340
341

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

WorldVec_FOT::WorldVec_FOT(WorldVector<DOFVector<double>*> vecs_, double fac_)
342
: FirstOrderTerm(vecs_[0]->getFeSpace()->getBasisFcts()->getDegree()), 
343
  fac(fac_)
344
{
345
346
347
348
349
350
351
352
353
354
355
  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];
356
}
357

358
359
360
361
void WorldVec_FOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
362
363
364
  getVectorAtQPs(vec0DV, elInfo, subAssembler, quad, vec0);
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
  if(numVecs>=3) getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
365
}
366

367
368
369
370
void WorldVec_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
371
372
373
  getVectorAtQPs(vec0DV, smallElInfo, largeElInfo, subAssembler, quad, vec0);
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
  if(numVecs>=3) getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
374
}
375

376
377
378
void WorldVec_FOT::getLb(const ElInfo *elInfo, 
			vector<mtl::dense_vector<double> >& result) const
{
379
380
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
381
	
382
383
384
385
386
387
388
  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);
  }
389
}
390

391
392
393
394
395
396
397
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)
{	
398
399
400
401
402
403
404
405
406
407
  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;
  }
408
}
409
410
411
412

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

WorldVector_FOT::WorldVector_FOT(DOFVector<WorldVector<double> >* dv, double fac_)
413
: FirstOrderTerm(dv->getFeSpace()->getBasisFcts()->getDegree()),
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
  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_)
464
: FirstOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree() + vecs_[0]->getFeSpace()->getBasisFcts()->getDegree()), 
465
466
  phaseDV(phaseDV_), 
  fac(fac_)
467
{
468
469
470
471
472
473
474
475
476
477
478
479
480
481
  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];
482
}
483

484
485
486
487
void WorldVecPhase_FOT::initElement(const ElInfo* elInfo, 
                    SubAssembler* subAssembler,
                    Quadrature *quad) 
{
488
489
490
491
  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);
492
}
493

494
495
496
497
void WorldVecPhase_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
                    SubAssembler* subAssembler,
                    Quadrature *quad) 
{
498
499
500
501
  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);
502
}
503

504
505
506
void WorldVecPhase_FOT::getLb(const ElInfo *elInfo, 
            vector<mtl::dense_vector<double> >& result) const
{
507
508
509
510
511
512
513
514
515
516
  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);
  }
517
}
518

519
520
521
522
523
524
525
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)
{   
526
527
528
529
530
531
532
533
534
535
  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;
  }
536
}
537

Praetorius, Simon's avatar
Praetorius, Simon committed
538
539
/* ----------------------------------------------------------- */

540
VecAndWorldVec_FOT::VecAndWorldVec_FOT(DOFVector<double>* vec0_, WorldVector<DOFVector<double>*> vecs_, AbstractFunction<double, double>* fct_, double fac_)
541
: FirstOrderTerm(fct_->getDegree() + vecs_[0]->getFeSpace()->getBasisFcts()->getDegree()), 
542
543
544
  vDV(vec0_), 
  fct(fct_), 
  fac(fac_)
Praetorius, Simon's avatar
Praetorius, Simon committed
545
546
547
548
549
550
{
  numVecs = vecs_.size();
  TEST_EXIT(numVecs == 2 || numVecs == 3)
	("Only Dim=2 or Dim=3 possible\n");
  
  for (int i = 0; i < numVecs; i++) {
551
552
    TEST_EXIT(vecs_[i])("One vector is NULL!\n");	
    auxFeSpaces.insert(vecs_[i]->getFeSpace());
Praetorius, Simon's avatar
Praetorius, Simon committed
553
554
555
556
557
  }
  
  vec0DV = vecs_[0];
  vec1DV = vecs_[1];
  if (numVecs >= 3)
558
    vec2DV = vecs_[2];
559
}
560

Praetorius, Simon's avatar
Praetorius, Simon committed
561
562
563
564
565
566
567
568
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)
569
    getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
570
}
571

Praetorius, Simon's avatar
Praetorius, Simon committed
572
573
574
575
576
577
578
579
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)
580
    getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
581
}
582

Praetorius, Simon's avatar
Praetorius, Simon committed
583
584
585
586
587
588
589
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++) {
590
591
592
593
594
595
596
    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);
Praetorius, Simon's avatar
Praetorius, Simon committed
597
  }
598
}
599

Praetorius, Simon's avatar
Praetorius, Simon committed
600
601
602
603
604
605
606
607
608
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++) {
609
610
611
612
613
614
615
616
    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]);
Praetorius, Simon's avatar
Praetorius, Simon committed
617
  }
618
}
Praetorius, Simon's avatar
Praetorius, Simon committed
619

620
621
622
623
624
/* ----------------------------------------------------------- */

// vec1*vec2*d/dxi(...)
VecAndPartialDerivative_FOT::VecAndPartialDerivative_FOT(DOFVectorBase<double> *dv1,
      int component_, double fac_)
625
: FirstOrderTerm(dv1->getFeSpace()->getBasisFcts()->getDegree()), 
626
627
628
  vec1DV(dv1), 
  fac(fac_), 
  component(component_)
629
630
631
632
633
634
635
636
637
638
639
{
  auxFeSpaces.insert(vec1DV->getFeSpace());

  setB(component);
}

void VecAndPartialDerivative_FOT::initElement(const ElInfo* elInfo, 
        SubAssembler* subAssembler,
        Quadrature *quad) 
{ 
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
640
}
641

642
643
644
645
646
void VecAndPartialDerivative_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
  SubAssembler* subAssembler,
  Quadrature *quad)
{ 
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
647
}
648

649
650
651
652
653
654
655
656
657
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);
  }
658
}
659

660
661
662
663
664
665
666
667
668
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;
669
}
670

Praetorius, Simon's avatar
Praetorius, Simon committed
671
672
673
674
/* ----------------------------------------------------------- */

// vec1*vec2*d/dxi(...)
VecAndPartialDerivativeIJ_FOT::VecAndPartialDerivativeIJ_FOT(DOFVectorBase<double> *dv1, 
675
676
677
678
							    DOFVectorBase<double> *dv2,
							    int i_,
							    int j_,
							    double fac_)
679
: FirstOrderTerm(dv1->getFeSpace()->getBasisFcts()->getDegree() + dv2->getFeSpace()->getBasisFcts()->getDegree() - 1), 
680
681
682
683
684
685
  vec1DV(dv1), 
  vec2DV(dv2), 
  fct(NULL), 
  fac(fac_), 
  xi(i_), 
  xj(j_)
Praetorius, Simon's avatar
Praetorius, Simon committed
686
687
688
689
690
691
692
693
694
695
696
697
698
{
  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_)
699
: FirstOrderTerm(fct_->getDegree()), 
700
701
702
703
704
705
  vec1DV(dv1), 
  vec2DV(dv2), 
  fct(fct_), 
  fac(fac_), 
  xi(i_), 
  xj(j_)
Praetorius, Simon's avatar
Praetorius, Simon committed
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
{
  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) {
736
737
738
    for(int iq = 0; iq < nPoints; iq++) {
      lb_one(Lambda, result[iq], vec1[iq] * grad[iq][xj] * fac);
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
739
  } else {
740
741
742
    for(int iq = 0; iq < nPoints; iq++) {
      lb_one(Lambda, result[iq], (*fct)(vec1[iq], grad[iq][xj]) * fac);
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
743
744
745
746
747
748
749
750
751
752
753
  }
}

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) {
754
755
    for(int iq = 0; iq < nPoints; iq++) 
      result[iq] += opFactor * vec1[iq] * grad[iq][xj] * grdUhAtQP[iq][xi] * fac;
Praetorius, Simon's avatar
Praetorius, Simon committed
756
  } else {
757
758
    for(int iq = 0; iq < nPoints; iq++) 
      result[iq] += opFactor * (*fct)(vec1[iq], grad[iq][xj]) * grdUhAtQP[iq][xi] * fac;
Praetorius, Simon's avatar
Praetorius, Simon committed
759
760
761
762
763
764
765
766
767
768
769
  }
}

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

// vec1*vec2*d/dxi(...)
Vec2AndPartialDerivative_FOT::Vec2AndPartialDerivative_FOT(
			       DOFVectorBase<double> *dv1, 
			       DOFVectorBase<double> *dv2,
			       int component_, 
			       double fac_)
770
: FirstOrderTerm(dv1->getFeSpace()->getBasisFcts()->getDegree() + dv2->getFeSpace()->getBasisFcts()->getDegree()), 
771
772
773
774
775
  vec1DV(dv1), 
  vec2DV(dv2), 
  fct(NULL), 
  fac(fac_), 
  component(component_)
Praetorius, Simon's avatar
Praetorius, Simon committed
776
777
778
779
780
781
782
783
784
785
786
787
788
789
{
  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_)
790
: FirstOrderTerm(fct_->getDegree()), 
791
792
793
794
795
  vec1DV(dv1), 
  vec2DV(dv2), 
  fct(fct_), 
  fac(fac_), 
  component(component_)
Praetorius, Simon's avatar
Praetorius, Simon committed
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
{
  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);
  }
832
}
833

Praetorius, Simon's avatar
Praetorius, Simon committed
834
835
836
837
838
839
840
841
842
843
844
845
846
847
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;
  }
848
}
Praetorius, Simon's avatar
Praetorius, Simon committed
849

850
851
/* ----------------------------------------------------------- */

852
853
PartialDerivative_FOT::PartialDerivative_FOT(int component_, double fac_)
: FirstOrderTerm(0), 
854
  component(component_), 
855
  fac(fac_)
856
{
857
  setB(component);
858
}
859

860
861
void PartialDerivative_FOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
862
					Quadrature *quad) {}
863
					
864
865
void PartialDerivative_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
866
					Quadrature *quad) {}
867
					
868
869
870
void PartialDerivative_FOT::getLb(const ElInfo *elInfo, 
			vector<mtl::dense_vector<double> >& result) const
{
871
872
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
873
	
874
875
876
  for (int iq = 0; iq < nPoints; iq++) {	
    lb_one(Lambda, result[iq], fac);
  }
877
}
878

879
880
881
882
883
884
885
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)
{	
886
887
  for (int iq = 0; iq < nPoints; iq++)	
    result[iq] += opFactor * grdUhAtQP[iq][component] * fac;
888
}
889
890
891
892

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

PartialDerivative_ZOT::PartialDerivative_ZOT(DOFVectorBase<double> *vecDV_, int component_, double fac_)
893
: ZeroOrderTerm(vecDV_->getFeSpace()->getBasisFcts()->getDegree() - 1), 
894
895
896
  vecDV(vecDV_), 
  fac(fac_), 
  component(component_)
897
{
898
899
  TEST_EXIT(vecDV_)("No first vector!\n");
  auxFeSpaces.insert(vecDV_->getFeSpace());
900

901
902
  facVec = new WorldVector<double>(); facVec->set(1.0);
  component2 = 0;
903
}
904

905
PartialDerivative_ZOT::PartialDerivative_ZOT(DOFVectorBase<double> *vecDV_, int component_, double fac_, WorldVector<double> *facVec_, int component2_)
906
: ZeroOrderTerm(vecDV_->getFeSpace()->getBasisFcts()->getDegree() - 1), 
907
908
909
910
911
  vecDV(vecDV_), 
  fac(fac_), 
  facVec(facVec_), 
  component(component_), 
  component2(component2_)
912
{
913
914
  TEST_EXIT(vecDV_)("No first vector!\n");
  auxFeSpaces.insert(vecDV_->getFeSpace());
915
}
916

917
918
919
920
921
void PartialDerivative_ZOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad)
{
  getGradientsAtQPs(vecDV, elInfo, subAssembler, quad, grad);
922
}
923

924
925
926
927
928
void PartialDerivative_ZOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad)
{
  getGradientsAtQPs(vecDV, smallElInfo, largeElInfo, subAssembler, quad, grad);
929
}
930

931
932
933
934
935
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;
936
}
937

938
939
940
941
942
943
944
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)
{	
945
946
947
  double factor = opFactor * fac * (*facVec)[component2];
  for (int iq = 0; iq < nPoints; iq++)		
    result[iq] += grad[iq][component] * uhAtQP[iq] * factor;
948
}
949
950
951
952

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

VecAndPartialDerivative_ZOT::VecAndPartialDerivative_ZOT(DOFVectorBase<double> *vecDV_, DOFVectorBase<double> *gradDV_, int component_,  double fac_)
953
: ZeroOrderTerm(vecDV_->getFeSpace()->getBasisFcts()->getDegree() + gradDV_->getFeSpace()->getBasisFcts()->getDegree() - 1),
954
955
956
  vecDV(vecDV_),
  gradDV(gradDV_),
  component(component_),
Praetorius, Simon's avatar
Praetorius, Simon committed
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
  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_)
972
: ZeroOrderTerm(fct_->getDegree() + gradDV_->getFeSpace()->getBasisFcts()->getDegree() - 1),
Praetorius, Simon's avatar
Praetorius, Simon committed
973
974
975
976
  vecDV(vecDV_),
  gradDV(gradDV_),
  component(component_),
  fct(fct_),
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
  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
1005
  if (fct == NULL) {
1006
1007
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += vec[iq] * grad[iq][component] * fac;
Praetorius, Simon's avatar
Praetorius, Simon committed
1008
  } else {
1009
1010
1011
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*fct)(vec[iq]) * grad[iq][component] * fac;
  }	
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
}

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
1022
  if (fct == NULL) {
1023
1024
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += vec[iq] * grad[iq][component] * uhAtQP[iq] * fac_;
Praetorius, Simon's avatar
Praetorius, Simon committed
1025
  } else {
1026
1027
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += (*fct)(vec[iq]) * grad[iq][component] * uhAtQP[iq] * fac_;
Praetorius, Simon's avatar
Praetorius, Simon committed
1028
  }
1029
1030
1031
1032
}

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

Praetorius, Simon's avatar
Praetorius, Simon committed
1033
1034
1035
1036
1037
Vec2AndPartialDerivative_ZOT::Vec2AndPartialDerivative_ZOT(DOFVectorBase<double> *vec1DV_, 
							    DOFVectorBase<double> *vec2DV_, 
							    DOFVectorBase<double> *gradDV_,
							    int component_,  
							    double fac_)
1038
: ZeroOrderTerm(vec1DV_->getFeSpace()->getBasisFcts()->getDegree() + vec2DV_->getFeSpace()->getBasisFcts()->getDegree() + gradDV_->getFeSpace()->getBasisFcts()->getDegree() - 1), 
1039
1040
1041
1042
1043
1044
  vec1DV(vec1DV_), 
  vec2DV(vec2DV_), 
  gradDV(gradDV_), 
  fct(NULL), 
  fac(fac_), 
  component(component_)
Praetorius, Simon's avatar
Praetorius, Simon committed
1045
1046
1047
1048
1049
1050
1051
{
  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());
1052
}
Praetorius, Simon's avatar
Praetorius, Simon committed
1053
1054
1055
1056
1057
1058
1059

Vec2AndPartialDerivative_ZOT::Vec2AndPartialDerivative_ZOT(DOFVectorBase<double> *vec1DV_, 
							  DOFVectorBase<double> *vec2DV_, 
							  DOFVectorBase<double> *gradDV_,
							  int component_,  
							  BinaryAbstractFunction<double, double, double>* fct_,
							  double fac_)
1060
: ZeroOrderTerm(fct_->getDegree() + gradDV_->getFeSpace()->getBasisFcts()->getDegree() - 1), 
1061
1062
1063
1064
1065
1066
  vec1DV(vec1DV_), 
  vec2DV(vec2DV_), 
  gradDV(gradDV_), 
  fct(fct_), 
  fac(fac_), 
  component(component_)
1067
{
Praetorius, Simon's avatar
Praetorius, Simon committed
1068
1069
1070
1071
1072
1073
  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());
1074
}
Praetorius, Simon's avatar
Praetorius, Simon committed
1075

1076
1077
1078
1079
void Vec2AndPartialDerivative_ZOT::initElement(const ElInfo* elInfo,
					SubAssembler* subAssembler,
					Quadrature *quad)
{
Praetorius, Simon's avatar
Praetorius, Simon committed
1080
1081
1082
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
  getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
  getGradientsAtQPs(gradDV, elInfo, subAssembler, quad, grad);
1083
}
1084

1085
1086
1087
1088
void Vec2AndPartialDerivative_ZOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad)
{
Praetorius, Simon's avatar
Praetorius, Simon committed
1089
1090
1091
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
  getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
  getGradientsAtQPs(gradDV, smallElInfo, largeElInfo, subAssembler, quad, grad);
1092
}
1093

1094
1095
void Vec2AndPartialDerivative_ZOT::getC(const ElInfo *, int nPoints, ElementVector &C)
{
Praetorius, Simon's avatar
Praetorius, Simon committed
1096
1097
1098
1099
1100
1101
1102
  if (fct == NULL) {
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += vec1[iq]*vec2[iq]*grad[iq][component]*fac;
  } else {
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*fct)(vec1[iq],vec2[iq])*grad[iq][component]*fac;
  }
1103
}
1104

1105
1106
1107
1108
1109
1110
1111
void Vec2AndPartialDerivative_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
1112
1113
1114
1115
1116
1117
1118
1119
1120
  if (fct == NULL) {
    for (int iq = 0; iq < nPoints; iq++) {
      result[iq] += (vec1[iq] * vec2[iq]) * grad[iq][component] * uhAtQP[iq] * opFactor * fac;
    }
  } else {
    for (int iq = 0; iq < nPoints; iq++) {
      result[iq] += (*fct)(vec1[iq], vec2[iq]) * grad[iq][component] * uhAtQP[iq] * opFactor * fac;
    }
  }
1121
}
1122
1123
1124
1125

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

MatrixPhase_SOT::MatrixPhase_SOT(DOFVectorBase<double>* phaseDV_, int comp_, double fac_) :
1126
    SecondOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()), phaseDV(phaseDV_), comp(comp_), fac(fac_), symmetric(true)
1127
1128
1129
{
    TEST_EXIT(phaseDV_)("no vector phase!\n");
    auxFeSpaces.insert(phaseDV_->getFeSpace());
1130
}
1131
1132
1133
1134
void MatrixPhase_SOT::initElement(const ElInfo* elInfo, 
            SubAssembler* subAssembler,
            Quadrature *quad) {
    getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
1135
}
1136
1137
1138
1139
void MatrixPhase_SOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
            SubAssembler* subAssembler,
            Quadrature *quad) {
    getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
1140
}
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164

/// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to the diagonal matrix.
void MatrixPhase_SOT::lDlt(const DimVec<WorldVector<double> >& Lambda, 
            const WorldVector<double>& vec,
            mtl::dense2D<double>& LALt,
            double factor) const
{
    const int dim = num_rows(LALt);
    
    for (int i = 0; i < dim; i++) {
        double val = 0.0;
        for (int k = 0; k < dimOfWorld; k++)
            val += vec[k] * Lambda[i][k