FirstOrderTerm.cc 28 KB
Newer Older
1 2 3 4 5 6 7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9 10 11 12 13 14 15 16 17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, 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.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
20 21


22 23 24
#include "FirstOrderTerm.h"
#include "DOFVector.h"

25 26
using namespace std;

27 28 29 30 31 32 33
namespace AMDiS {

  // =========== VecAtQP_FOT ==========

  VecAtQP_FOT::VecAtQP_FOT(DOFVectorBase<double> *dv,
			   AbstractFunction<double, double> *af,
			   WorldVector<double> *wv)
34
    : FirstOrderTerm(af ? af->getDegree() : dv->getFeSpace()->getBasisFcts()->getDegree()), vec(dv), f(af), b(wv)
35 36 37
  {
    TEST_EXIT(dv)("No vector!\n");

38
    auxFeSpaces.insert(dv->getFeSpace());
39 40 41 42 43
  }

  VecAtQP_FOT::VecAtQP_FOT(DOFVectorBase<double> *dv,
			   AbstractFunction<double, double> *af,
			   int bIdx)
44
    : FirstOrderTerm(af ? af->getDegree() : dv->getFeSpace()->getBasisFcts()->getDegree()), vec(dv), f(af)
45 46 47 48
  {
    TEST_EXIT(dv)("No vector!\n");

    bOne = bIdx;
49
    auxFeSpaces.insert(dv->getFeSpace());
50 51
  }

52
  void VecAtQP_FOT::initElement(const ElInfo* elInfo,
53
				SubAssembler* subAssembler,
54
				Quadrature *quad)
55
  {
56
    getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
57 58
  }

59 60 61 62 63
  void VecAtQP_FOT::initElement(const ElInfo* smallElInfo,
				const ElInfo* largeElInfo,
				SubAssembler* subAssembler,
				Quadrature *quad)
  {
64
    getVectorAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad, vecAtQPs);
65 66
  }

67
  void VecAtQP_FOT::getLb(const ElInfo *elInfo,
68
			  vector<mtl::dense_vector<double> >& Lb) const
69
  {
70
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
71
    const int nPoints = static_cast<int>(Lb.size());
72 73

    if (bOne > -1) {
74 75 76 77 78 79 80
      if (f) {
	for (int iq = 0; iq < nPoints; iq++)
	  lb_one(grdLambda, Lb[iq], (*f)(vecAtQPs[iq]));
      } else {
	for (int iq = 0; iq < nPoints; iq++)
	  lb_one(grdLambda, Lb[iq], vecAtQPs[iq]);
      }
81
    } else if (b) {
82 83 84 85 86 87 88
      if (f) {
	for (int iq = 0; iq < nPoints; iq++)
	  lb(grdLambda, *b, Lb[iq], (*f)(vecAtQPs[iq]));
      } else {
	for (int iq = 0; iq < nPoints; iq++)
	  lb(grdLambda, *b, Lb[iq], vecAtQPs[iq]);
      }
89
    } else {
90 91
      if (f) {
	for (int iq = 0; iq < nPoints; iq++)
92
	  l1(grdLambda, Lb[iq], (*f)(vecAtQPs[iq]));
93 94
      } else {
	for (int iq = 0; iq < nPoints; iq++)
95
	  l1(grdLambda, Lb[iq], vecAtQPs[iq]);
96
      }
97 98 99
    }
  }

100 101 102 103 104
  void VecAtQP_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,
105
			  double fac)
106 107
  {
    int dow = Global::getGeo(WORLD);
108

109
    if (num_rows(grdUhAtQP) > 0) {
110
      for (int iq = 0; iq < nPoints; iq++) {
111
	double factor = (f ? (*f)(vecAtQPs[iq]) : vecAtQPs[iq]);
112 113
	double resultQP = 0.0;
	for (int i = 0; i < dow; i++)
114
	  resultQP += grdUhAtQP[iq][i];
115

116
	result[iq] += fac * factor * resultQP;
117
      }
118 119 120
    }
  }

121

122 123
  // =========== CoordsAtQP_FOT ===========

124
  void CoordsAtQP_FOT::initElement(const ElInfo* elInfo,
125
				   SubAssembler* subAssembler,
126
				   Quadrature *quad)
127
  {
128
    subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs);
129 130
  }

131
  void CoordsAtQP_FOT::getLb(const ElInfo *elInfo,
132
			     vector<mtl::dense_vector<double> >& Lb) const
133
  {
134
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
135 136
    const int nPoints = static_cast<int>(Lb.size());

137
    for (int iq = 0; iq < nPoints; iq++)
138
      l1(grdLambda, Lb[iq], (*g)(coordsAtQPs[iq]));
139 140 141
  }

  void CoordsAtQP_FOT::eval(int nPoints,
142 143 144 145
			    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,
146
			    double f)
147 148
  {
    int dow = Global::getGeo(WORLD);
149

150
    if (num_rows(grdUhAtQP) > 0) {
151 152 153 154
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*g)(coordsAtQPs[iq]);
	double resultQP = 0.0;
	for (int i = 0; i < dow; i++)
155
	  resultQP += grdUhAtQP[iq][i];
156

157 158 159 160 161 162 163 164
	result[iq] += f * factor * resultQP;
      }
    }
  }


  // ========== VecCoordsAtQP_FOT ==========

165
  void VecCoordsAtQP_FOT::initElement(const ElInfo* elInfo,
166
				      SubAssembler* subAssembler,
167
				      Quadrature *quad)
168
  {
169
    subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs);
170 171
  }

172 173
  void VecCoordsAtQP_FOT::getLb(const ElInfo *elInfo,
				vector<mtl::dense_vector<double> >& Lb) const
174
  {
175
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
176 177
    const int nPoints = static_cast<int>(Lb.size());

178
    for (int iq = 0; iq < nPoints; iq++)
179
      lb(grdLambda, *b, Lb[iq], (*g)(coordsAtQPs[iq]));
180 181
  }

182 183 184 185 186
  void VecCoordsAtQP_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,
187
			        double f)
188 189
  {
    int dow = Global::getGeo(WORLD);
190

191
    if (num_rows(grdUhAtQP) > 0) {
192 193 194 195 196
      for (int iq = 0; iq < nPoints; iq++) {
	double factor = (*g)(coordsAtQPs[iq]);
	double resultQP = 0.0;
	for (int i = 0; i < dow; i++)
	  resultQP += grdUhAtQP[iq][i];
197

198 199 200 201 202 203 204 205 206 207
	result[iq] += f * factor * resultQP;
      }
    }
  }


  // ========== VectorGradient_FOT ==========

  VectorGradient_FOT::VectorGradient_FOT(DOFVectorBase<double> *dv,
					 AbstractFunction<WorldVector<double>, WorldVector<double> > *af)
208
    : FirstOrderTerm(af ? af->getDegree() : dv->getFeSpace()->getBasisFcts()->getDegree()),
209
      vec(dv), f(af)
210
  {
211 212
    FUNCNAME("VectorGradient_FOT::VectorGradient_FOT()");

213 214
    TEST_EXIT(dv)("No vector!\n");

215
    auxFeSpaces.insert(dv->getFeSpace());
216 217
  }

218
  void VectorGradient_FOT::initElement(const ElInfo* elInfo,
219
				       SubAssembler* subAssembler,
220
				       Quadrature *quad)
221
  {
222
    getGradientsAtQPs(vec, elInfo, subAssembler, quad, gradAtQPs);
223 224
  }

225
  void VectorGradient_FOT::getLb(const ElInfo *elInfo,
226
				 vector<mtl::dense_vector<double> >& Lb) const
227
  {
228
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
229 230 231
    const int nPoints = static_cast<int>(Lb.size());

    if (f)
232
      for (int iq = 0; iq < nPoints; iq++)
233
        lb(grdLambda, (*f)(gradAtQPs[iq]), Lb[iq], 1.0);
234
    else
235
      for (int iq = 0; iq < nPoints; iq++)
236
        lb(grdLambda, gradAtQPs[iq], Lb[iq], 1.0);
237 238 239
  }

  void VectorGradient_FOT::eval(int nPoints,
240 241 242 243
				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,
244 245
				double factor)
  {
246
    if (num_rows(grdUhAtQP) > 0) {
247 248 249 250 251 252 253
      if (f) {
	for (int iq = 0; iq < nPoints; iq++) {
	  WorldVector<double> b = (*f)(gradAtQPs[iq]);
	  result[iq] += b * grdUhAtQP[iq] * factor;
	}
      } else {
	for (int iq = 0; iq < nPoints; iq++)
254
	  result[iq] += gradAtQPs[iq] * grdUhAtQP[iq] * factor;
255 256 257 258 259 260 261 262 263
      }
    }
  }


  // =========== VectorFct_FOT ==========

  VectorFct_FOT::VectorFct_FOT(DOFVectorBase<double> *dv,
			       AbstractFunction<WorldVector<double>, double> *fct)
264
    : FirstOrderTerm(fct->getDegree()), vec(dv), vecFct(fct)
265 266 267
  {
    TEST_EXIT(dv)("No vector!\n");

268
    auxFeSpaces.insert(dv->getFeSpace());
269 270
  }

271
  void VectorFct_FOT::initElement(const ElInfo* elInfo,
272
				  SubAssembler* subAssembler,
273
				  Quadrature *quad)
274
  {
275
    getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
276 277
  }

278
  void VectorFct_FOT::getLb(const ElInfo *elInfo,
279
			    vector<mtl::dense_vector<double> >& Lb) const
280
  {
281
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
282 283
    const int nPoints = static_cast<int>(Lb.size());

284
    for (int iq = 0; iq < nPoints; iq++)
285
      lb(grdLambda, (*vecFct)(vecAtQPs[iq]), Lb[iq], 1.0);
286 287
  }

288 289 290 291 292
  void VectorFct_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,
293
			    double factor)
294
  {
295
    if (num_rows(grdUhAtQP) > 0) {
296 297 298 299 300 301 302 303 304
      for (int iq = 0; iq < nPoints; iq++) {
	WorldVector<double> b = (*vecFct)(vecAtQPs[iq]);
	result[iq] += b * grdUhAtQP[iq] * factor;
      }
    }
  }


  // =========== VecFctAtQP_FOT ==========
305 306

  void VecFctAtQP_FOT::initElement(const ElInfo* elInfo,
307
				   SubAssembler* subAssembler,
308
				   Quadrature *quad)
309
  {
310
    subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs);
311 312 313
  }


314
  void VecFctAtQP_FOT::getLb(const ElInfo *elInfo,
315
			     vector<mtl::dense_vector<double> >& Lb) const
316
  {
317
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
318 319
    const int nPoints = static_cast<int>(Lb.size());

320
    for (int iq = 0; iq < nPoints; iq++)
321
      lb(grdLambda, (*g)(coordsAtQPs[iq]), Lb[iq], 1.0);
322 323 324
  }

  void VecFctAtQP_FOT::eval(int nPoints,
325 326 327 328
			    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,
329
			    double f)
330 331
  {
    int dow = Global::getGeo(WORLD);
332

333
    if (num_rows(grdUhAtQP) > 0) {
334 335 336 337
      for (int iq = 0; iq < nPoints; iq++) {
	double resultQP = 0.0;
	const WorldVector<double> &b = (*g)(coordsAtQPs[iq]);
	for (int i = 0; i < dow; i++)
338
	  resultQP += b[i] * grdUhAtQP[iq][i];
339 340 341 342 343 344 345 346 347 348
	result[iq] += f * resultQP;
      }
    }
  }


  // =========== VecGrad_FOT ===========

  VecGrad_FOT::VecGrad_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
			   BinaryAbstractFunction<WorldVector<double>, double, WorldVector<double> > *fct)
349
    : FirstOrderTerm(fct->getDegree()), vec1(dv1), vec2(dv2), vecFct(fct)
350 351 352 353
  {
    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");

354
    auxFeSpaces.insert(dv1->getFeSpace());
355
    auxFeSpaces.insert(dv2->getFeSpace());
356 357
  }

358
  void VecGrad_FOT::initElement(const ElInfo* elInfo,
359
				SubAssembler* subAssembler,
360
				Quadrature *quad)
361
  {
362
    getVectorAtQPs(vec1, elInfo, subAssembler, quad, vecAtQPs);
363
    getGradientsAtQPs(vec2, elInfo, subAssembler, quad, gradAtQPs);
364 365 366 367 368 369 370
  }

  void VecGrad_FOT::initElement(const ElInfo* smallElInfo,
				const ElInfo* largeElInfo,
				SubAssembler* subAssembler,
				Quadrature *quad)
  {
371
    getVectorAtQPs(vec1, smallElInfo, largeElInfo, subAssembler, quad, vecAtQPs);
372
    getGradientsAtQPs(vec2, smallElInfo, largeElInfo, subAssembler, quad, gradAtQPs);
373 374
  }

375
  void VecGrad_FOT::getLb(const ElInfo *elInfo,
376
			  vector<mtl::dense_vector<double> >& Lb) const
377
  {
378
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
379 380
    const int nPoints = static_cast<int>(Lb.size());

381
    for (int iq = 0; iq < nPoints; iq++)
382
      lb(grdLambda, (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]), Lb[iq], 1.0);
383
  }
384

385 386 387 388 389
  void VecGrad_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,
390
			  double factor)
391 392
  {
    if (num_rows(grdUhAtQP) > 0) {
393 394 395 396 397 398 399 400 401 402
      for (int iq = 0; iq < nPoints; iq++) {
	WorldVector<double> b = (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]);
	result[iq] += b * grdUhAtQP[iq] * factor;
      }
    }
  }


  // =========== FctGrad2_FOT ===========

403
  void FctGrad2_FOT::initElement(const ElInfo* elInfo,
404
				SubAssembler* subAssembler,
405
				Quadrature *quad)
406
  {
407 408
    getGradientsAtQPs(vec1, elInfo, subAssembler, quad, grad1AtQPs);
    getGradientsAtQPs(vec2, elInfo, subAssembler, quad, grad2AtQPs);
409 410
  }

411
  void FctGrad2_FOT::getLb(const ElInfo *elInfo,
412
			   vector<mtl::dense_vector<double> >& Lb) const
413
  {
414
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
415 416
    const int nPoints = static_cast<int>(Lb.size());

417
    for (int iq = 0; iq < nPoints; iq++)
418
      lb(grdLambda, (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]), Lb[iq], 1.0);
419
  }
420

421
  void FctGrad2_FOT::eval(int nPoints,
422 423 424 425
			  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,
426
			  double factor)
427
  {
428
    if (num_rows(grdUhAtQP) > 0) {
429 430 431 432 433 434 435 436 437 438
      for (int iq = 0; iq < nPoints; iq++) {
	WorldVector<double> b = (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]);
	result[iq] += b * grdUhAtQP[iq] * factor;
      }
    }
  }


  // ========== Vec2AndGradAtQP_FOT ==========

439
  Vec2AndGradAtQP_FOT::Vec2AndGradAtQP_FOT(DOFVectorBase<double> *dv1,
440 441
					   DOFVectorBase<double> *dv2,
					   TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f_, WorldVector<double> *b_)
442 443 444 445
    : FirstOrderTerm(f_->getDegree()),
      vec1(dv1),
      vec2(dv2),
      f(f_),
446
      b(b_)
447 448
  {
    auxFeSpaces.insert(dv1->getFeSpace());
449
    auxFeSpaces.insert(dv2->getFeSpace());
450 451
  }

452
  void Vec2AndGradAtQP_FOT::initElement(const ElInfo* elInfo,
453
					SubAssembler* subAssembler,
454 455
					Quadrature *quad)
  {
456 457 458
    getVectorAtQPs(vec1, elInfo, subAssembler, quad, vec1AtQPs);
    getVectorAtQPs(vec2, elInfo, subAssembler, quad, vec2AtQPs);
    getGradientsAtQPs(vec1, elInfo, subAssembler, quad, gradAtQPs1);
459
  }
460 461 462

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

467 468
    for (int iq = 0; iq < nPoints; iq++) {
      if (b)
469
	lb(grdLambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq],gradAtQPs1[iq]));
470
      else
471
	l1(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], gradAtQPs1[iq]));
472 473
    }
  }
474

475 476 477 478 479
  void Vec2AndGradAtQP_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,
480 481
				  double fac)
  {
482
    if (num_rows(grdUhAtQP) > 0)
483
      for (int iq = 0; iq < nPoints; iq++)
484 485
	result[iq] += fac *
	  (*f)(vec1AtQPs[iq], vec2AtQPs[iq], gradAtQPs1[iq]) *
486 487 488 489 490 491 492 493 494 495 496 497
	  ((*b) * grdUhAtQP[iq]);
  }


  // ========== FctVecAtQP_FOT ==========


  FctVecAtQP_FOT::FctVecAtQP_FOT(DOFVectorBase<double> *dv,
				 BinaryAbstractFunction<double, WorldVector<double>, double> *f_,
				 WorldVector<double> *b_)
    : FirstOrderTerm(f_->getDegree()), vec(dv), f(f_), b(b_)
  {
498
    auxFeSpaces.insert(dv->getFeSpace());
499 500
  }

501
  void FctVecAtQP_FOT::initElement(const ElInfo* elInfo,
502
				   SubAssembler* subAssembler,
503 504
				   Quadrature *quad)
  {
505
    getVectorAtQPs(vec, elInfo, subAssembler, quad, vecAtQPs);
506
    subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs);
507 508
  }

509
  void FctVecAtQP_FOT::getLb(const ElInfo *elInfo,
510
			     vector<mtl::dense_vector<double> >& Lb) const
511
  {
512
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
513 514
    const int nPoints = static_cast<int>(Lb.size());

515 516
    for (int iq = 0; iq < nPoints; iq++) {
      if (b)
517
	lb(grdLambda, *b, Lb[iq], (*f)(coordsAtQPs[iq], vecAtQPs[iq]));
518
      else
519
	l1(grdLambda, Lb[iq], (*f)(coordsAtQPs[iq], vecAtQPs[iq]));
520 521
    }
  }
522

523
  void FctVecAtQP_FOT::eval(int nPoints,
524 525 526 527
			    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,
528
			    double fac)
529
  {
530
    if (num_rows(grdUhAtQP) > 0)
531
      for (int iq = 0; iq < nPoints; iq++)
532
	result[iq] +=
533 534 535 536 537 538 539 540 541
	  fac * (*f)(coordsAtQPs[iq], vecAtQPs[iq]) * ((*b) * grdUhAtQP[iq]);
  }


  // ========== Vec2AtQP_FOT ==========

  Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
			     BinaryAbstractFunction<double, double, double> *af,
			     WorldVector<double> *b_)
542 543
    : FirstOrderTerm(af ?
		     af->getDegree() :
544
		     dv1->getFeSpace()->getBasisFcts()->getDegree()
545
		     + dv2->getFeSpace()->getBasisFcts()->getDegree()),
546
      vec1(dv1), vec2(dv2), f(af), b(b_)
547
  {
548 549 550
    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");

551
    auxFeSpaces.insert(dv1->getFeSpace());
552
    auxFeSpaces.insert(dv2->getFeSpace());
553 554 555 556 557
  }

  Vec2AtQP_FOT::Vec2AtQP_FOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
			     BinaryAbstractFunction<double, double, double> *af,
			     int bIdx)
558 559
    : FirstOrderTerm(af ?
		     af->getDegree() :
560
		     dv1->getFeSpace()->getBasisFcts()->getDegree()
561
		     + dv2->getFeSpace()->getBasisFcts()->getDegree()),
562
      vec1(dv1), vec2(dv2), f(af)
563
  {
564 565 566 567
    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");

    bOne = bIdx;
568
    auxFeSpaces.insert(dv1->getFeSpace());
569
    auxFeSpaces.insert(dv2->getFeSpace());
570 571
  }

572
  void Vec2AtQP_FOT::initElement(const ElInfo* elInfo,
573
				 SubAssembler* subAssembler,
574 575
				 Quadrature *quad)
  {
576 577
    getVectorAtQPs(vec1, elInfo, subAssembler, quad, vec1AtQPs);
    getVectorAtQPs(vec2, elInfo, subAssembler, quad, vec2AtQPs);
578
  }
579

580
  void Vec2AtQP_FOT::getLb(const ElInfo *elInfo,
581
			   vector<mtl::dense_vector<double> >& Lb) const
582
  {
583
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
584
    const int nPoints = static_cast<int>(Lb.size());
585 586

    if (bOne > -1) {
587 588 589 590 591 592 593
      if (f) {
	for (int iq = 0; iq < nPoints; iq++)
	  lb_one(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
      } else {
	for (int iq = 0; iq < nPoints; iq++)
	  lb_one(grdLambda, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
      }
594
    } else if (b) {
595 596 597 598 599 600 601
      if (f) {
	for (int iq = 0; iq < nPoints; iq++)
	  lb(grdLambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
      } else {
	for (int iq = 0; iq < nPoints; iq++)
	  lb(grdLambda, *b, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
      }
602
    } else {
603 604 605 606 607 608 609
      if (f) {
	for (int iq = 0; iq < nPoints; iq++)
	  l1(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq]));
      } else {
	for (int iq = 0; iq < nPoints; iq++)
	  l1(grdLambda, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq]);
      }
610 611
    }
  }
612

613
  void Vec2AtQP_FOT::eval(int nPoints,
614 615 616 617
			  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,
618
			  double fac)
619
  {
620 621
    if (num_rows(grdUhAtQP) > 0) {
      if (f) {
622
	for (int iq = 0; iq < nPoints; iq++)
623 624
	  result[iq] += fac * (*f)(vec1AtQPs[iq], vec2AtQPs[iq]) * ((*b) * grdUhAtQP[iq]);
      } else {
625
	for (int iq = 0; iq < nPoints; iq++)
626 627 628
	  result[iq] += fac * vec1AtQPs[iq] * vec2AtQPs[iq] * ((*b) * grdUhAtQP[iq]);
      }
    }
629 630
  }

631

632
  // ========== Vec3AtQP_FOT ==========
633

634 635
  Vec3AtQP_FOT::Vec3AtQP_FOT(DOFVectorBase<double> *dv1,
				   DOFVectorBase<double> *dv2,
636 637 638
				   DOFVectorBase<double> *dv3,
				   TertiaryAbstractFunction<double, double, double, double> *f_,
				   WorldVector<double> *bvec)
639 640
    : FirstOrderTerm(f_ ?
		     f_->getDegree() :
641 642
		     dv1->getFeSpace()->getBasisFcts()->getDegree()
		     + dv2->getFeSpace()->getBasisFcts()->getDegree()
643 644 645 646
		     + dv3->getFeSpace()->getBasisFcts()->getDegree()),
      vec1(dv1),
      vec2(dv2),
      vec3(dv3),
647 648
      f(f_),
      b(bvec)
649 650 651 652
  {
    auxFeSpaces.insert(dv1->getFeSpace());
    auxFeSpaces.insert(dv2->getFeSpace());
    auxFeSpaces.insert(dv3->getFeSpace());
653 654
  }

655 656
  Vec3AtQP_FOT::Vec3AtQP_FOT(DOFVectorBase<double> *dv1,
				   DOFVectorBase<double> *dv2,
657 658 659
				   DOFVectorBase<double> *dv3,
				   TertiaryAbstractFunction<double, double, double, double> *f_,
				   int b)
660 661
    : FirstOrderTerm(f_ ?
		     f_->getDegree() :
662 663
		     dv1->getFeSpace()->getBasisFcts()->getDegree()
		     + dv2->getFeSpace()->getBasisFcts()->getDegree()
664 665 666 667
		     + dv3->getFeSpace()->getBasisFcts()->getDegree()),
      vec1(dv1),
      vec2(dv2),
      vec3(dv3),
668
      f(f_)
669
  {
670
    bOne = b;
671 672 673
    auxFeSpaces.insert(dv1->getFeSpace());
    auxFeSpaces.insert(dv2->getFeSpace());
    auxFeSpaces.insert(dv3->getFeSpace());
674 675
  }

676
  void Vec3AtQP_FOT::initElement(const ElInfo* elInfo,
677
				    SubAssembler* subAssembler,
678 679
				    Quadrature *quad)
  {
680 681 682
    getVectorAtQPs(vec1, elInfo, subAssembler, quad, vec1AtQPs);
    getVectorAtQPs(vec2, elInfo, subAssembler, quad, vec2AtQPs);
    getVectorAtQPs(vec3, elInfo, subAssembler, quad, vec3AtQPs);
683
  }
684

685
  void Vec3AtQP_FOT::getLb(const ElInfo *elInfo,
686
			      vector<mtl::dense_vector<double> >& Lb) const
687
  {
688
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
689 690
    const int nPoints = static_cast<int>(Lb.size());

691
    if (bOne > -1) {
692 693 694 695 696 697 698
      if (f) {
	for (int iq = 0; iq < nPoints; iq++)
	  lb_one(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
      } else {
	for (int iq = 0; iq < nPoints; iq++)
	  lb_one(grdLambda, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq] * vec3AtQPs[iq]);
      }
699
    } else {
700 701 702 703 704 705 706 707 708 709 710 711 712 713
      if (f) {
	for (int iq = 0; iq < nPoints; iq++) {
	  if (b)
	    lb(grdLambda, *b, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
	  else
	    l1(grdLambda, Lb[iq], (*f)(vec1AtQPs[iq], vec2AtQPs[iq], vec3AtQPs[iq]));
	}
      } else {
	for (int iq = 0; iq < nPoints; iq++) {
	  if (b)
	    lb(grdLambda, *b, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq] * vec3AtQPs[iq]);
	  else
	    l1(grdLambda, Lb[iq], vec1AtQPs[iq] * vec2AtQPs[iq] * vec3AtQPs[iq]);
	}
714 715 716 717
      }
    }
  }

718
  void Vec3AtQP_FOT::eval(int nPoints,
719 720 721 722
			     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,
723
			     double fac)
724
  {
725 726 727 728 729 730
    if (num_rows(grdUhAtQP) == 0)
      return;

    if (bOne > -1) {
      ERROR_EXIT("Not yet implemented!\n");
    } else {
731 732
      if (f) {
	for (int iq = 0; iq < nPoints; iq++)
733
	  result[iq] += fac * (*f)(vec1AtQPs[iq], vec2AtQPs[iq] ,vec3AtQPs[iq]) *
734 735 736
	    ((*b) * grdUhAtQP[iq]);
      } else {
	for (int iq = 0; iq < nPoints; iq++)
737
	  result[iq] += fac * vec1AtQPs[iq] * vec2AtQPs[iq] * vec3AtQPs[iq] *
738 739
	    ((*b) * grdUhAtQP[iq]);
      }
740
    }
741 742 743 744 745
  }


  // =========== General_FOT ==========

746 747
  General_FOT::General_FOT(vector<DOFVectorBase<double>*> vecs,
			   vector<DOFVectorBase<double>*> grads,
748
			   TertiaryAbstractFunction<WorldVector<double>,
749
			   WorldVector<double>,
750
			   vector<double>,
751
			   vector<WorldVector<double> > > *af)
752 753
    : FirstOrderTerm(af->getDegree()), vecs_(vecs), grads_(grads), f_(af)
  {
754
    vecsAtQPs.resize(vecs_.size());
755 756 757
    gradsAtQPs_.resize(grads_.size());

    for (int i = 0; i < static_cast<int>(vecs.size()); i++) {
758
      TEST_EXIT(vecs[i])("One vector is NULL!\n");
759

760
      auxFeSpaces.insert(vecs[i]->getFeSpace());
761
    }
762 763

    for (int i = 0; i < static_cast<int>(grads.size()); i++) {
764
      TEST_EXIT(grads[i])("One gradient vector is NULL!\n");
765

766
      auxFeSpaces.insert(grads[i]->getFeSpace());
767
    }
768 769
  }

770

771
  void General_FOT::initElement(const ElInfo* elInfo,
772 773 774 775 776 777
				SubAssembler* subAssembler,
				Quadrature *quad)
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

778
    subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs);
779 780

    for (int i = 0; i < nVecs; i++)
781
      getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad, vecsAtQPs[i]);
782
    for (int i = 0; i < nGrads; i++)
783
      getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad, gradsAtQPs_[i]);
784 785
  }

786

787
  void General_FOT::getLb(const ElInfo *elInfo,
788
			  vector<mtl::dense_vector<double> >& Lb) const
789 790 791 792
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

793 794
    vector<double> vecsArg(nVecs);
    vector<WorldVector<double> > gradsArg(nGrads);
795

796
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
797
    const int nPoints = static_cast<int>(Lb.size());
798 799 800

    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++)
801 802 803
	vecsArg[i] = vecsAtQPs[i][iq];
      for (int i = 0; i < nGrads; i++)
	gradsArg[i] = gradsAtQPs_[i][iq];
804

805
      lb(grdLambda, (*f_)(coordsAtQPs[iq], vecsArg, gradsArg), Lb[iq], 1.0);
806 807 808
    }
  }

809

810 811 812 813 814
  void General_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,
815
			  double factor)
816
  {
817
    int dow = Global::getGeo(WORLD);
818 819 820
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

821 822
    vector<double> vecsArg(nVecs);
    vector<WorldVector<double> > gradsArg(nGrads);
823

824
    if (num_rows(grdUhAtQP) > 0) {
825 826 827 828
      for (int iq = 0; iq < nPoints; iq++) {
	double resultQP = 0.0;

	for (int i = 0; i < nVecs; i++)
829
	  vecsArg[i] = vecsAtQPs[i][iq];
830
	for (int i = 0; i < nGrads; i++)
831
	  gradsArg[i] = gradsAtQPs_[i][iq];
832

833
	const WorldVector<double> &b = (*f_)(coordsAtQPs[iq], vecsArg, gradsArg);
834

835
	for (int i = 0; i < dow; i++)
836
	  resultQP += grdUhAtQP[iq][i] * b[i];
837

838 839 840 841 842 843 844 845
	result[iq] += factor * resultQP;
      }
    }
  }


  // =========== GeneralParametric_FOT ==========

846 847
  GeneralParametric_FOT::GeneralParametric_FOT(vector<DOFVectorBase<double>*> vecs,
					       vector<DOFVectorBase<double>*> grads,
848
					       QuartAbstractFunction<WorldVector<double>,
849 850
					       WorldVector<double>,
					       WorldVector<double>,
851
					       vector<double>,
852
					       vector<WorldVector<double> > > *af)
853 854
    : FirstOrderTerm(af->getDegree()), vecs_(vecs), grads_(grads), f_(af)
  {
855
    vecsAtQPs.resize(vecs_.size());
856 857 858
    gradsAtQPs_.resize(grads_.size());

    for (int i = 0; i < static_cast<int>(vecs.size()); i++) {
859
      TEST_EXIT(vecs[i])("One vector is NULL!\n");
860

861
      auxFeSpaces.insert(vecs[i]->getFeSpace());
862
    }
863 864

    for (int i = 0; i < static_cast<int>(grads.size()); i++) {
865
      TEST_EXIT(grads[i])("One gradient vector is NULL!\n");
866

867
      auxFeSpaces.insert(grads[i]->getFeSpace());
868
    }
869 870 871
  }


872
  void GeneralParametric_FOT::initElement(const ElInfo* elInfo,
873 874 875 876 877 878 879
					  SubAssembler* subAssembler,
					  Quadrature *quad)
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

    elInfo->getElementNormal(elementNormal_);
880
    subAssembler->getCoordsAtQPs(elInfo, quad, coordsAtQPs);
881 882

    for (int i = 0; i < nVecs; i++)
883
      getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad, vecsAtQPs[i]);
884
    for (int i = 0; i < nGrads; i++)
885
      getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad, gradsAtQPs_[i]);
886 887
  }

888

889
  void GeneralParametric_FOT::getLb(const ElInfo *elInfo,
890
				    vector<mtl::dense_vector<double> >& Lb) const
891 892 893 894
  {
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

895 896
    vector<double> vecsArg(nVecs);
    vector<WorldVector<double> > gradsArg(nGrads);
897

898
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
899 900
    const int nPoints = static_cast<int>(Lb.size());

901 902
    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++)
903
	vecsArg[i] = vecsAtQPs[i][iq];
904 905
      for (int i = 0; i < nGrads; i++)
	gradsArg[i] = gradsAtQPs_[i][iq];
906 907

      lb(grdLambda, (*f_)(coordsAtQPs[iq],  elementNormal_, vecsArg, gradsArg),
908
	 Lb[iq], 1.0);
909 910 911
    }
  }

912

913 914 915 916 917
  void GeneralParametric_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,
918
				   double factor)
919 920 921 922 923
  {
    int dow = Global::getGeo(WORLD);
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());

924 925
    vector<double> vecsArg(nVecs);
    vector<WorldVector<double> > gradsArg(nGrads);
926

927
    if (num_rows(grdUhAtQP) > 0) {
928 929 930 931
      for (int iq = 0; iq < nPoints; iq++) {
	double resultQP = 0.0;

	for (int i = 0; i < nVecs; i++)
932
	  vecsArg[i] = vecsAtQPs[i][iq];
933
	for (int i = 0; i < nGrads; i++)
934
	  gradsArg[i] = gradsAtQPs_[i][iq];
935

936
	const WorldVector<double> &b =
937
	  (*f_)(coordsAtQPs[iq],  elementNormal_, vecsArg, gradsArg);
938 939 940

	for (int i = 0; i < dow; i++)
	  resultQP += grdUhAtQP[iq][i] * b[i];
941

942 943 944 945 946 947 948 949
	result[iq] += factor * resultQP;
      }
    }
  }



}