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
  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
55
    l1lt(Lambda, LALt[iq], f(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
  if (num_rows(D2UhAtQP) > 0) {
    for (int iq = 0; iq < nPoints; iq++) {
Praetorius, Simon's avatar
Praetorius, Simon committed
67
      double feval = f(iq) * opFactor * fac;
68
69
70
71
72
73
      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
  int nPoints = grdUhAtQP.size();
  for (int iq = 0; iq < nPoints; iq++) {
Praetorius, Simon's avatar
Praetorius, Simon committed
81
    double factor = f(iq) * fac;
82
83
    axpy(factor, grdUhAtQP[iq], result[iq]);
  }
84
}
85

86
87
double Phase_SOT::f(const int iq) const
{
Praetorius, Simon's avatar
Praetorius, Simon committed
88
  return std::max(0.0, std::min( 1.0, phase[iq] ) );
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
  numVecs=vecs_.size();
  
  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
353
  if(numVecs>=2) vec1DV=vecs_[1];
354
  if(numVecs>=3) vec2DV=vecs_[2];
355
}
356

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

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

375
376
377
void WorldVec_FOT::getLb(const ElInfo *elInfo, 
			vector<mtl::dense_vector<double> >& result) const
{
378
379
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
380
	
381
382
383
  for (int iq = 0; iq < nPoints; iq++) {
    WorldVector<double> vec;
    vec[0]=vec0[iq];
Praetorius, Simon's avatar
Praetorius, Simon committed
384
    if(numVecs>=2) vec[1]=vec1[iq];
385
386
387
    if(numVecs>=3) vec[2]=vec2[iq];
    lb(Lambda, vec, result[iq], fac);
  }
388
}
389

390
391
392
393
394
395
396
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)
{	
397
398
399
400
401
  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
402
    if(numVecs>=2) resultQP += grdUhAtQP[iq][1] * vec1[iq];
403
404
405
406
    if(numVecs>=3) resultQP += grdUhAtQP[iq][2] * vec2[iq];
    
    result[iq] += factor * resultQP;
  }
407
}
408
409
410
411

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

WorldVector_FOT::WorldVector_FOT(DOFVector<WorldVector<double> >* dv, double fac_)
412
: FirstOrderTerm(dv->getFeSpace()->getBasisFcts()->getDegree()),
413
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
  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_)
463
: FirstOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree() + vecs_[0]->getFeSpace()->getBasisFcts()->getDegree()), 
464
465
  phaseDV(phaseDV_), 
  fac(fac_)
466
{
467
468
469
470
471
472
473
474
475
476
477
478
479
480
  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];
481
}
482

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

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

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

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

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

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

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

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

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

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

619
620
621
622
623
/* ----------------------------------------------------------- */

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

  setB(component);
}

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

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

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

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

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

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

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) {
753
754
    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
755
  } else {
756
757
    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
758
759
760
761
762
763
764
765
766
767
768
  }
}

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

// vec1*vec2*d/dxi(...)
Vec2AndPartialDerivative_FOT::Vec2AndPartialDerivative_FOT(
			       DOFVectorBase<double> *dv1, 
			       DOFVectorBase<double> *dv2,
			       int component_, 
			       double fac_)
769
: FirstOrderTerm(dv1->getFeSpace()->getBasisFcts()->getDegree() + dv2->getFeSpace()->getBasisFcts()->getDegree()), 
770
771
772
773
774
  vec1DV(dv1), 
  vec2DV(dv2), 
  fct(NULL), 
  fac(fac_), 
  component(component_)
Praetorius, Simon's avatar
Praetorius, Simon committed
775
776
777
778
779
780
781
782
783
784
785
786
787
788
{
  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_)
789
: FirstOrderTerm(fct_->getDegree()), 
790
791
792
793
794
  vec1DV(dv1), 
  vec2DV(dv2), 
  fct(fct_), 
  fac(fac_), 
  component(component_)
Praetorius, Simon's avatar
Praetorius, Simon committed
795
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
{
  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);
  }
831
}
832

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

849
850
/* ----------------------------------------------------------- */

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

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

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

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

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

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

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

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

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

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

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

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

VecAndPartialDerivative_ZOT::VecAndPartialDerivative_ZOT(DOFVectorBase<double> *vecDV_, DOFVectorBase<double> *gradDV_, int component_,  double fac_)
952
: ZeroOrderTerm(vecDV_->getFeSpace()->getBasisFcts()->getDegree() + gradDV_->getFeSpace()->getBasisFcts()->getDegree() - 1),
953
954
955
  vecDV(vecDV_),
  gradDV(gradDV_),
  component(component_),
Praetorius, Simon's avatar
Praetorius, Simon committed
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
  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_)
971
: ZeroOrderTerm(fct_->getDegree() + gradDV_->getFeSpace()->getBasisFcts()->getDegree() - 1),
Praetorius, Simon's avatar
Praetorius, Simon committed
972
973
974
975
  vecDV(vecDV_),
  gradDV(gradDV_),
  component(component_),
  fct(fct_),
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
  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