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
#include "POperators.h"
#include "Helpers.h"

22
23
namespace AMDiS { namespace extensions {
  
24
25
using namespace std;

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

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

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

49
50
51
void Phase_SOT::getLALt(const ElInfo *elInfo,
			vector<mtl::dense2D<double> > &LALt) const
{
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++)
Praetorius, Simon's avatar
Praetorius, Simon committed
56
    l1lt(Lambda, LALt[iq], f(iq) * fac);
57
}
58

59
void Phase_SOT::eval(int nPoints,
60
61
62
63
64
		    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)
65
{	
66
67
  if (num_rows(D2UhAtQP) > 0) {
    for (int iq = 0; iq < nPoints; iq++) {
Praetorius, Simon's avatar
Praetorius, Simon committed
68
      double feval = f(iq) * opFactor * fac;
69
70
71
72
73
74
      double resultQP = 0.0;
      for (int i = 0; i < dimOfWorld; i++)
	resultQP += D2UhAtQP[iq][i][i];	
      result[iq] += feval * resultQP;
    }
  }
75
}
76

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

87
88
double Phase_SOT::f(const int iq) const
{
Praetorius, Simon's avatar
Praetorius, Simon committed
89
  return std::max(0.0, std::min( 1.0, phase[iq] ) );
90
}
91
92
93

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

WorldVec_FOT::WorldVec_FOT(WorldVector<DOFVector<double>*> vecs_, double fac_)
343
: FirstOrderTerm(vecs_[0]->getFeSpace()->getBasisFcts()->getDegree()), 
344
  fac(fac_)
345
{
346
  numVecs=vecs_.getSize();
347
348
349
350
351
352
353
  
  for (int i = 0; i < numVecs; i++) {
    TEST_EXIT(vecs_[i])("One vector is NULL!\n");	
    auxFeSpaces.insert(vecs_[i]->getFeSpace());
  }
  
  vec0DV=vecs_[0];
Praetorius, Simon's avatar
Praetorius, Simon committed
354
  if(numVecs>=2) vec1DV=vecs_[1];
355
  if(numVecs>=3) vec2DV=vecs_[2];
356
}
357

358
359
360
361
void WorldVec_FOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
362
  getVectorAtQPs(vec0DV, elInfo, subAssembler, quad, vec0);
Praetorius, Simon's avatar
Praetorius, Simon committed
363
  if(numVecs>=2) getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
364
  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
  getVectorAtQPs(vec0DV, smallElInfo, largeElInfo, subAssembler, quad, vec0);
Praetorius, Simon's avatar
Praetorius, Simon committed
372
  if(numVecs>=2) getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
373
  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
  for (int iq = 0; iq < nPoints; iq++) {
    WorldVector<double> vec;
    vec[0]=vec0[iq];
Praetorius, Simon's avatar
Praetorius, Simon committed
385
    if(numVecs>=2) vec[1]=vec1[iq];
386
387
388
    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
  double factor = fac * opFactor;
  for (int iq = 0; iq < nPoints; iq++) {
    double resultQP = 0.0;
		    
    resultQP += grdUhAtQP[iq][0] * vec0[iq];
Praetorius, Simon's avatar
Praetorius, Simon committed
403
    if(numVecs>=2) resultQP += grdUhAtQP[iq][1] * vec1[iq];
404
405
406
407
    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
  numVecs=vecs_.getSize();
469
470
471
472
473
474
475
476
477
478
479
480
481
  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
  numVecs = vecs_.getSize();
Praetorius, Simon's avatar
Praetorius, Simon committed
547
548
549
550
  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
  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);
For faster browsing, not all history is shown. View entire blame