POperators.cc 79.8 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
26
27
28
Phase_SOT::Phase_SOT(DOFVectorBase<double>* phaseDV_, double fac_)
: SecondOrderTerm(6), 
  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
94
95
96
Phase_FOT::Phase_FOT(DOFVectorBase<double>* phaseDV_, double fac_)
: FirstOrderTerm(6), 
  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
150
151
152
Phase_ZOT::Phase_ZOT(DOFVectorBase<double>* phaseDV_, double fac_)
: ZeroOrderTerm(6),
  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
200
201
202
203
PhaseInverse_ZOT::PhaseInverse_ZOT(DOFVectorBase<double>* phaseDV_, double fac_)
: ZeroOrderTerm(6),
  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
213
214
215
216
217

PhaseInverse_ZOT::PhaseInverse_ZOT(DOFVectorBase<double>* phaseDV_, double fac_, WorldVector<double> *facVec_, int component_)
: ZeroOrderTerm(6),
  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
266
267
268
Pow3_ZOT::Pow3_ZOT(DOFVectorBase<double> *rhoDV_, double fac_)
: ZeroOrderTerm(3), 
  rhoDV(rhoDV_), 
  fac(fac_)
269
{
270
271
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
272
}
273
274
275
276
277

void Pow3_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
278
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
279
}
280

281
282
283
284
285
void Pow3_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
286
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
287
}
288

289
290
void Pow3_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{ 
291
292
  for (int iq = 0; iq < nPoints; iq++)
    C[iq] += f(iq) * fac;
293
}
294

295
296
297
298
299
300
301
void Pow3_ZOT::eval(int nPoints,
	const mtl::dense_vector<double>& uhAtQP,
	const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
	const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
	mtl::dense_vector<double>& result,
	double opFactor)
{
302
303
304
  double factor = opFactor * fac;
  for (int iq = 0; iq < nPoints; iq++)
    result[iq] += factor * f(iq) * uhAtQP[iq];
305
}
306

307
308
double Pow3_ZOT::f(const int iq) const
{
309
  return rho[iq]*sqr(rho[iq]);
310
}
311
312
313

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

314
315
316
Pow3Phase_ZOT::Pow3Phase_ZOT(DOFVectorBase<double>* phaseDV, DOFVectorBase<double> *rhoDV_, double fac_)
: Phase_ZOT(phaseDV, fac_),
  rhoDV(rhoDV_)
317
{
318
319
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
320
}
321
322
323
324
325

void Pow3Phase_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
326
327
  Phase_ZOT::initElement(elInfo,subAssembler,quad);
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
328
}
329

330
331
332
333
334
void Pow3Phase_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
335
336
  Phase_ZOT::initElement(largeElInfo,smallElInfo,subAssembler,quad);
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
337
}
338

339
340
double Pow3Phase_ZOT::f(const int iq) const
{
341
  return rho[iq]*sqr(rho[iq]);
342
}
343
344
345

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

346
347
348
349
Pow2_ZOT::Pow2_ZOT(DOFVectorBase<double> *rhoDV_, double fac_)
: ZeroOrderTerm(2), 
  rhoDV(rhoDV_), 
  fac(fac_)
350
{
351
352
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
353
}
354
355
356
357
358

void Pow2_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
359
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
360
}
361

362
363
364
365
366
void Pow2_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
367
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
368
}
369

370
371
void Pow2_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{ 
372
373
  for (int iq = 0; iq < nPoints; iq++)
    C[iq] += f(iq) * fac;
374
}
375

376
377
378
379
380
381
382
void Pow2_ZOT::eval(int nPoints,
	const mtl::dense_vector<double>& uhAtQP,
	const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
	const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
	mtl::dense_vector<double>& result,
	double opFactor)
{
383
384
385
  double factor = opFactor * fac;
  for (int iq = 0; iq < nPoints; iq++)
    result[iq] += factor * f(iq) * uhAtQP[iq];
386
}
387

388
389
double Pow2_ZOT::f(const int iq) const
{
390
  return sqr(rho[iq]);
391
}
392
393
394

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

395
396
397
Pow2Phase_ZOT::Pow2Phase_ZOT(DOFVectorBase<double>* phaseDV, DOFVectorBase<double> *rhoDV_, double fac_)
: Phase_ZOT(phaseDV, fac_), 
  rhoDV(rhoDV_)
398
{
399
400
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
401
}
402
403
404
405
406

void Pow2Phase_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
407
408
  Phase_ZOT::initElement(elInfo,subAssembler,quad);
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
409
}
410

411
412
413
414
415
void Pow2Phase_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
416
417
  Phase_ZOT::initElement(largeElInfo,smallElInfo,subAssembler,quad);
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
418
}
419

420
421
double Pow2Phase_ZOT::f(const int iq) const
{
422
  return sqr(rho[iq]);
423
}
424

425
426
427
428
429
/* ----------------------------------------------------------- */

ConstrainedFrac_ZOT::ConstrainedFrac_ZOT(DOFVectorBase<double> *phaseDV_, DOFVectorBase<double> *rhoDV_, double fac_)
: Phase_ZOT(phaseDV_,fac_), 
  rhoDV(rhoDV_)
430
{
431
432
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
433
}
434
435
436
437
438

void ConstrainedFrac_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
439
440
  Phase_ZOT::initElement(elInfo, subAssembler,quad);
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
441
}
442

443
444
445
446
447
void ConstrainedFrac_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
448
449
  Phase_ZOT::initElement(largeElInfo, smallElInfo, subAssembler,quad);
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
450
}
451

452
453
double ConstrainedFrac_ZOT::f(const int iq) const
{
454
  return 1.0/std::max(1.e-6, rho[iq]+0.9);
455
}
456

457
/* ----------------------------------------------------------- */
458

459
460
461
ConstrainedFracSqr_ZOT::ConstrainedFracSqr_ZOT(DOFVectorBase<double> *phaseDV_, DOFVectorBase<double> *rhoDV_, double factor_)
: Phase_ZOT(phaseDV_,factor_), 
  rhoDV(rhoDV_)
462
{
463
464
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
465
}
466
467
468
469
470

void ConstrainedFracSqr_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
471
472
  Phase_ZOT::initElement(elInfo, subAssembler,quad);
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
473
}
474

475
476
477
478
479
void ConstrainedFracSqr_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
480
481
  Phase_ZOT::initElement(largeElInfo, smallElInfo, subAssembler,quad);
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
482
}
483

484
485
double ConstrainedFracSqr_ZOT::f(const int iq) const
{
486
  return -1.0/std::max(1.e-6, sqr(rho[iq]+0.9));
487
}
488

489
/* ----------------------------------------------------------- */
490

491
492
493
494
495
496
ConstrainedPowImpl_ZOT::ConstrainedPowImpl_ZOT(DOFVectorBase<double> *rhoDV_, double p_, double rhoMin_, double fac_)
: ZeroOrderTerm(2), 
  rhoDV(rhoDV_), 
  p(p_), 
  rhoMin(rhoMin_), 
  fac(fac_)
497
{
498
499
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
500
}
501
502
503
504
505

void ConstrainedPowImpl_ZOT::initElement(const ElInfo* elInfo, 
            SubAssembler* subAssembler,
            Quadrature *quad) 
{
506
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
507
}
508

509
510
511
512
513
void ConstrainedPowImpl_ZOT::initElement(const ElInfo* largeElInfo,
            const ElInfo* smallElInfo,
            SubAssembler* subAssembler,
            Quadrature *quad)
{
514
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
515
}
516

517
518
void ConstrainedPowImpl_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{ 
519
520
  for (int iq = 0; iq < nPoints; iq++)
      C[iq] += f(iq)*fac;
521
}
522

523
524
525
526
527
528
529
void ConstrainedPowImpl_ZOT::eval(int nPoints,
    const mtl::dense_vector<double>& uhAtQP,
    const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
    const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
    mtl::dense_vector<double>& result,
    double opFactor)
{
530
531
532
  double factor = opFactor*fac;
  for (int iq = 0; iq < nPoints; iq++)
    result[iq] += factor * f(iq) * uhAtQP[iq];
533
}
534

535
536
double ConstrainedPowImpl_ZOT::f(const int iq) const
{
537
  return 2.0*p*Helpers::toInt(0.0<=rhoMin-rho[iq]);
538
}
539

540
/* ----------------------------------------------------------- */
541

542
543
544
545
546
547
ConstrainedPowExpl_ZOT::ConstrainedPowExpl_ZOT(DOFVectorBase<double> *rhoDV_, double p_, double rhoMin_, double fac_) 
: ZeroOrderTerm(2), 
  rhoDV(rhoDV_), 
  p(p_), 
  rhoMin(rhoMin_), 
  fac(fac_)
548
{
549
550
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
551
}
552
553
554
555
556

void ConstrainedPowExpl_ZOT::initElement(const ElInfo* elInfo, 
            SubAssembler* subAssembler,
            Quadrature *quad) 
{
557
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
558
}
559

560
561
562
563
564
void ConstrainedPowExpl_ZOT::initElement(const ElInfo* largeElInfo,
            const ElInfo* smallElInfo,
            SubAssembler* subAssembler,
            Quadrature *quad)
{
565
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
566
}
567

568
569
void ConstrainedPowExpl_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{ 
570
571
  for (int iq = 0; iq < nPoints; iq++)
      C[iq] += f(iq)*fac;
572
}
573

574
575
576
577
578
579
580
void ConstrainedPowExpl_ZOT::eval(int nPoints,
    const mtl::dense_vector<double>& uhAtQP,
    const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
    const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
    mtl::dense_vector<double>& result,
    double opFactor)
{
581
582
583
  double factor = opFactor*fac;
  for (int iq = 0; iq < nPoints; iq++)
      result[iq] += factor * f(iq) * uhAtQP[iq];
584
}
585

586
587
double ConstrainedPowExpl_ZOT::f(const int iq) const
{
588
  return -2.0*p*std::max(0.0, rhoMin-rho[iq]);
589
}
590

591
592
593
594
595
596
597
598
/* ----------------------------------------------------------- */

ConstrainedPow_ZOT::ConstrainedPow_ZOT(DOFVectorBase<double> *rhoDV_, double p_, double rhoMin_, double fac_) 
: ZeroOrderTerm(2), 
  rhoDV(rhoDV_), 
  p(p_), 
  rhoMin(rhoMin_), 
  fac(fac_)
599
{
600
601
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
602
}
603
604
605
606
607

void ConstrainedPow_ZOT::initElement(const ElInfo* elInfo, 
            SubAssembler* subAssembler,
            Quadrature *quad) 
{
608
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
609
}
610

611
612
613
614
615
void ConstrainedPow_ZOT::initElement(const ElInfo* largeElInfo,
            const ElInfo* smallElInfo,
            SubAssembler* subAssembler,
            Quadrature *quad)
{
616
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
617
}
618

619
620
void ConstrainedPow_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{ 
621
622
  for (int iq = 0; iq < nPoints; iq++)
      C[iq] += f(iq)*fac;
623
}
624

625
626
627
628
629
630
631
void ConstrainedPow_ZOT::eval(int nPoints,
    const mtl::dense_vector<double>& uhAtQP,
    const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
    const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
    mtl::dense_vector<double>& result,
    double opFactor)
{
632
633
634
  double factor = opFactor*fac;
  for (int iq = 0; iq < nPoints; iq++)
      result[iq] += factor * f(iq) * uhAtQP[iq];
635
}
636

637
638
double ConstrainedPow_ZOT::f(const int iq) const
{
639
  return -2.0*p*std::max(0.0, rhoMin-rho[iq]);
640
}
641

642
/* ----------------------------------------------------------- */
643

644
645
646
FactorPhase_ZOT::FactorPhase_ZOT(DOFVectorBase<double> *phaseDV_, DOFVectorBase<double> *rhoDV_, double fac_)
: Phase_ZOT(phaseDV_,fac_), 
  rhoDV(rhoDV_)
647
{
648
649
  TEST_EXIT(rhoDV_)("No vector rho!\n");
  auxFeSpaces.insert(rhoDV_->getFeSpace());
650
}
651
652
653
654
655

void FactorPhase_ZOT::initElement(const ElInfo* elInfo, 
			SubAssembler* subAssembler,
			Quadrature *quad) 
{
656
657
  Phase_ZOT::initElement(elInfo, subAssembler,quad);
  getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
658
}
659

660
661
662
663
664
void FactorPhase_ZOT::initElement(const ElInfo* largeElInfo,
			const ElInfo* smallElInfo,
			SubAssembler* subAssembler,
			Quadrature *quad)
{
665
666
  Phase_ZOT::initElement(largeElInfo, smallElInfo, subAssembler,quad);
  getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
667
}
668

669
670
double FactorPhase_ZOT::f(const int iq) const
{
671
  return rho[iq];
672
}
673
674
675
676

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

WorldVecAndVecFct_FOT::WorldVecAndVecFct_FOT(WorldVector<DOFVectorBase<double>*> vecs_, DOFVectorBase<double> *phaseDV_, double fac_)
677
678
679
: FirstOrderTerm(2), 
  phaseDV(phaseDV_), 
  fac(fac_)
680
{
681
682
683
684
685
686
687
688
689
690
691
692
693
  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];
694
}
695

696
697
698
699
void WorldVecAndVecFct_FOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
700
701
702
703
  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);
704
}
705

706
707
708
709
void WorldVecAndVecFct_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
710
711
712
713
  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);
714
}
715

716
717
718
void WorldVecAndVecFct_FOT::getLb(const ElInfo *elInfo, 
			vector<mtl::dense_vector<double> >& result) const
{
719
720
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
721
	
722
723
724
725
726
727
728
  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);
  }
729
}
730

731
732
733
734
735
736
737
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)
{	
738
739
740
741
742
743
744
745
746
747
  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;
  }
748
}
749
750
751
752

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

WorldVec_FOT::WorldVec_FOT(WorldVector<DOFVector<double>*> vecs_, double fac_)
753
754
: FirstOrderTerm(2), 
  fac(fac_)
755
{
756
757
758
759
760
761
762
763
764
765
766
  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];
767
}
768

769
770
771
772
void WorldVec_FOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
773
774
775
  getVectorAtQPs(vec0DV, elInfo, subAssembler, quad, vec0);
  getVectorAtQPs(vec1DV, elInfo, subAssembler, quad, vec1);
  if(numVecs>=3) getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
776
}
777

778
779
780
781
void WorldVec_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
					SubAssembler* subAssembler,
					Quadrature *quad) 
{
782
783
784
  getVectorAtQPs(vec0DV, smallElInfo, largeElInfo, subAssembler, quad, vec0);
  getVectorAtQPs(vec1DV, smallElInfo, largeElInfo, subAssembler, quad, vec1);
  if(numVecs>=3) getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
785
}
786

787
788
789
void WorldVec_FOT::getLb(const ElInfo *elInfo, 
			vector<mtl::dense_vector<double> >& result) const
{
790
791
  const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
  const int nPoints = static_cast<int>(result.size());
792
	
793
794
795
796
797
798
799
  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);
  }
800
}
801

802
803
804
805
806
807
808
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)
{	
809
810
811
812
813
814
815
816
817
818
  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;
  }
819
}
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874

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

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

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

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

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

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

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

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

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

WorldVecPhase_FOT::WorldVecPhase_FOT(DOFVectorBase<double> *phaseDV_, WorldVector<DOFVector<double>*> vecs_,double fac_)
875
876
877
: FirstOrderTerm(2), 
  phaseDV(phaseDV_), 
  fac(fac_)
878
{
879
880
881
882
883
884
885
886
887
888
889
890
891
892
  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];
893
}
894

895
896
897
898
void WorldVecPhase_FOT::initElement(const ElInfo* elInfo, 
                    SubAssembler* subAssembler,
                    Quadrature *quad) 
{
899
900
901
902
  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);
903
}
904

905
906
907
908
void WorldVecPhase_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
                    SubAssembler* subAssembler,
                    Quadrature *quad) 
{
909
910
911
912
  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);
913
}
914

915
916
917
void WorldVecPhase_FOT::getLb(const ElInfo *elInfo, 
            vector<mtl::dense_vector<double> >& result) const
{
918
919
920
921
922
923
924
925
926
927
  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);
  }
928
}
929

930
931
932
933
934
935
936
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)
{   
937
938
939
940
941
942
943
944
945
946
  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;
  }
947
}
948

Praetorius, Simon's avatar
Praetorius, Simon committed
949
950
/* ----------------------------------------------------------- */

951
952
953
954
955
VecAndWorldVec_FOT::VecAndWorldVec_FOT(DOFVector<double>* vec0_, WorldVector<DOFVector<double>*> vecs_, AbstractFunction<double, double>* fct_, double fac_)
: FirstOrderTerm(2), 
  vDV(vec0_), 
  fct(fct_), 
  fac(fac_)
Praetorius, Simon's avatar
Praetorius, Simon committed
956
957
958
959
960
961
{
  numVecs = vecs_.size();
  TEST_EXIT(numVecs == 2 || numVecs == 3)
	("Only Dim=2 or Dim=3 possible\n");
  
  for (int i = 0; i < numVecs; i++) {
962
963
    TEST_EXIT(vecs_[i])("One vector is NULL!\n");	
    auxFeSpaces.insert(vecs_[i]->getFeSpace());
Praetorius, Simon's avatar
Praetorius, Simon committed
964
965
966
967
968
  }
  
  vec0DV = vecs_[0];
  vec1DV = vecs_[1];
  if (numVecs >= 3)
969
    vec2DV = vecs_[2];
970
}
971

Praetorius, Simon's avatar
Praetorius, Simon committed
972
973
974
975
976
977
978
979
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)
980
    getVectorAtQPs(vec2DV, elInfo, subAssembler, quad, vec2);
981
}
982

Praetorius, Simon's avatar
Praetorius, Simon committed
983
984
985
986
987
988
989
990
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)
991
    getVectorAtQPs(vec2DV, smallElInfo, largeElInfo, subAssembler, quad, vec2);
992
}
993

Praetorius, Simon's avatar
Praetorius, Simon committed
994
995
996
997
998
999
1000
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++) {
1001
1002
1003
1004
1005
1006
1007
    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
1008
  }
1009
}
1010

Praetorius, Simon's avatar
Praetorius, Simon committed
1011
1012
1013
1014
1015
1016
1017
1018
1019
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++) {
1020
1021
1022
1023
1024
1025
1026
1027
    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
1028
  }
1029
}
Praetorius, Simon's avatar
Praetorius, Simon committed
1030

1031
1032
1033
1034
1035
/* ----------------------------------------------------------- */

// vec1*vec2*d/dxi(...)
VecAndPartialDerivative_FOT::VecAndPartialDerivative_FOT(DOFVectorBase<double> *dv1,
      int component_, double fac_)
1036
1037
1038
1039
: FirstOrderTerm(3), 
  vec1DV(dv1), 
  fac(fac_), 
  component(component_)
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
{
  auxFeSpaces.insert(vec1DV->getFeSpace());

  setB(component);
}

void VecAndPartialDerivative_FOT::initElement(const ElInfo* elInfo, 
        SubAssembler* subAssembler,
        Quadrature *quad) 
{