SecondOrderTerm.h 38.5 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 25 26 27 28 29 30 31

/** \file SecondOrderTerm.h */

#ifndef AMDIS_SECOND_ORDER_TERM_H
#define AMDIS_SECOND_ORDER_TERM_H

#include "AMDiS_fwd.h"
#include "OperatorTerm.h"
#include "ElInfo.h"
#include "AbstractFunction.h"
32
#include "MatrixVectorOperations.h"
33 34 35 36 37

namespace AMDiS {

  /**
   * \ingroup Assembler
38
   *
39 40 41 42 43 44 45
   * \brief
   * Describes the second order terms: \f$ \nabla \cdot A \nabla u(\vec{x}) \f$
   */
  class SecondOrderTerm : public OperatorTerm
  {
  public:
    /// Constructor.
46 47
    SecondOrderTerm(int deg)
      : OperatorTerm(deg)
48 49 50
    {}

    /// Destructor.
51
    virtual ~SecondOrderTerm()
52 53 54
    {}

    /// Evaluation of \f$ \Lambda A \Lambda^t \f$ at all quadrature points.
55
    virtual void getLALt(const ElInfo *elInfo,
56
			 std::vector<mtl::dense2D<double> > &result) const = 0;
57 58 59 60

    /// Evaluation of \f$ A \nabla u(\vec{x}) \f$ at all quadrature points.
    virtual void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
			  std::vector<WorldVector<double> > &result) = 0;
61 62 63 64 65

  protected:
    /// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$.
    void lalt(const DimVec<WorldVector<double> >& Lambda,
	      const WorldMatrix<double>& matrix,
66
	      mtl::dense2D<double>& LALt,
67 68 69
	      bool symm,
	      double factor) const;

70 71 72
    /// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for \f$ A \f$
    /// the matrix having a ONE in the position \f$ (K,L) \f$
    /// and ZEROS in all other positions.
73 74
    static void lalt_kl(const DimVec<WorldVector<double> >& Lambda,
			int k, int l,
75
			mtl::dense2D<double>& LALt,
76 77
			double factor);

78
    /// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to
79
    /// the identity.
80
    inline void l1lt(const DimVec<WorldVector<double> >& Lambda,
81
		     mtl::dense2D<double>& LALt,
82 83
		     double factor) const
    {
84
      const int dim = num_rows(LALt);
85

86 87 88 89 90 91 92 93 94 95 96 97 98
      for (int i = 0; i < dim; i++) {
	double val = 0.0;
	for (int k = 0; k < dimOfWorld; k++)
	  val += Lambda[i][k] * Lambda[i][k];
	val *= factor;
	LALt[i][i] += val;
	for (int j = i + 1; j < dim; j++) {
	  val = 0.0;
	  for (int k = 0; k < dimOfWorld; k++)
	    val += Lambda[i][k] * Lambda[j][k];
	  val *= factor;
	  LALt[i][j] += val;
	  LALt[j][i] += val;
99 100
	}
      }
101 102
    }

103 104 105 106 107
  };


  /**
   * \ingroup Assembler
108
   *
109 110 111 112
   * \brief
   * Implements the laplace operator multiplied with a scalar factor:
   * \f$ f \cdot \Delta u(\vec{x}) \f$
   */
113
  class Simple_SOT : public SecondOrderTerm
114 115 116
  {
  public:
    /// Constructor.
117 118
    Simple_SOT(double f = 1.0)
      : SecondOrderTerm(0)
119 120 121 122 123 124 125 126
    {
      factor = new double;
      *factor = f;

      setSymmetric(true);
    }

    /// Constructor.
127 128
    Simple_SOT(double *fptr)
      : SecondOrderTerm(0),
129 130 131 132 133 134
	factor(fptr)
    {
      setSymmetric(true);
    }

    /// Implements SecondOrderTerm::getLALt().
135
    inline void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const
136
    {
137
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
138 139
      const int nPoints = static_cast<int>(LALt.size());

140
      for (int iq = 0; iq < nPoints; iq++)
141
	l1lt(grdLambda, LALt[iq], (*factor));
142 143
    }

144
    /// Implemetation of SecondOrderTerm::eval().
145
    void eval(int nPoints,
146 147 148 149
	      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,
150
	      double f)
151 152 153
    {
      int dow = Global::getGeo(WORLD);

154
      if (num_rows(D2UhAtQP) > 0) {
155 156 157 158 159 160 161 162 163 164
	for (int iq = 0; iq < nPoints; iq++) {
	  double resultQP = 0.0;
	  for (int i = 0; i < dow; i++) {
	    resultQP += D2UhAtQP[iq][i][i];
	  }
	  result[iq] += resultQP * f * (*factor);
	}
      }
    }

165
    /// Implemetation of SecondOrderTerm::weakEval().
166
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
167
		  std::vector<WorldVector<double> > &result)
168 169 170 171 172 173 174 175 176 177 178 179 180 181
    {
      int nPoints = grdUhAtQP.size();
      for (int iq = 0; iq < nPoints; iq++)
	axpy(*factor, grdUhAtQP[iq], result[iq]);
    }

  private:
    /// Pointer to the factor.
    double *factor;
  };


  /**
   * \ingroup Assembler
182
   *
183 184 185 186 187
   * \brief
   * SecondOrderTerm where A is a function which maps a DOFVector evaluated at
   * a given quadrature point to a WolrdMatrix:
   * \f$ \nabla \cdot A(v(\vec{x})) \nabla u(\vec{x}) \f$
   */
188
  class MatrixFct_SOT : public SecondOrderTerm
189 190 191
  {
  public:
    /// Constructor.
192
    MatrixFct_SOT(DOFVectorBase<double> *dv,
193 194 195 196 197 198
		  AbstractFunction<WorldMatrix<double>, double> *fct,
		  AbstractFunction<WorldVector<double>, WorldMatrix<double> > *div,
		  bool sym = false);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
199
		     Quadrature *quad = NULL);
200 201

    /// Implements SecondOrderTerm::getLALt().
202
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
203

204
    /// Implemetation of SecondOrderTerm::eval().
205
    void eval(int nPoints,
206 207 208 209
	      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,
210 211
	      double factor);

212
    /// Implemetation of SecondOrderTerm::weakEval().
213 214 215 216 217 218 219 220
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Pointer to the values of the DOFVector at quadrature points.
221
    mtl::dense_vector<double> vecAtQPs;
222 223 224 225 226 227 228 229 230 231 232 233

    /// Function for A.
    AbstractFunction<WorldMatrix<double>, double>* matrixFct;

    ///
    AbstractFunction<WorldVector<double>, WorldMatrix<double> >* divFct;

    /// True, if \ref matrixFct produces always symmetric matrices.
    bool symmetric;
  };


234
  /**
235 236 237 238 239 240 241 242 243
   * \ingroup Assembler
   *
   * \brief
   * SecondOrderTerm where A is a given fixed WorldMatrix<double>:
   * \f$ \nabla \cdot A \nabla u(\vec{x}) \f$
   */
  class Matrix_SOT : public SecondOrderTerm {
  public:
    /// Constructor
244
    Matrix_SOT(WorldMatrix<double> mat)
245 246 247 248 249 250 251
      : SecondOrderTerm(0), matrix(mat)
    {
      symmetric = matrix.isSymmetric();
      setSymmetric(symmetric);
    }

    /// Implements SecondOrderTerm::getLALt().
252
    inline void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const
253 254
    {
      const DimVec<WorldVector<double> >& Lambda = elInfo->getGrdLambda();
255 256
      const int nPoints = static_cast<int>(LALt.size());

257
      for (int iq = 0; iq < nPoints; iq++)
258
	lalt(Lambda, matrix, LALt[iq], symmetric, 1.0);
259
    }
260

261
    /// Implemetation of SecondOrderTerm::eval().
262
    void eval(int nPoints,
263 264 265 266
	      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,
267 268
	      double factor);

269
    /// Implemetation of SecondOrderTerm::weakEval().
270 271 272 273 274 275 276 277 278 279 280 281
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    /// Matrix stroring A.
    WorldMatrix<double> matrix;

    /// True, if \ref matrix is symmetric.
    bool symmetric;
  };


282
  /**
283 284 285 286 287 288 289
   * \ingroup Assembler
   *
   * \brief
   * SecondOrderTerm where A is a WorldMatrix<double> having a ONE in position IJ
   * and ZERO in all other positions
   * \f$ \nabla \cdot A \nabla u(\vec{x}) \f$
   */
290
  class FactorIJ_SOT : public SecondOrderTerm
291 292 293
  {
  public:
    /// Constructor.
294
    FactorIJ_SOT(int x_i, int x_j, double f)
295 296 297 298 299 300 301 302 303
      : SecondOrderTerm(0), xi(x_i), xj(x_j)
    {
      factor = new double;
      *factor = f;

      setSymmetric(xi == xj);
    }

    /// Constructor.
304
    FactorIJ_SOT(int x_i, int x_j, double *fptr)
305 306 307 308 309 310
      : SecondOrderTerm(0), xi(x_i), xj(x_j), factor(fptr)
    {
      setSymmetric(xi == xj);
    }

    /// Implements SecondOrderTerm::getLALt().
311
    inline void getLALt(const ElInfo *elInfo,
312
			std::vector<mtl::dense2D<double> > &LALt) const
313
    {
314
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
315 316
      const int nPoints = static_cast<int>(LALt.size());

317
      for (int iq = 0; iq < nPoints; iq++)
318
	lalt_kl(grdLambda, xi, xj, LALt[iq], (*factor));
319 320
    }

321
    /// Implemetation of SecondOrderTerm::eval().
322
    void eval(int nPoints,
323 324 325 326
	      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,
327 328
	      double fac)
    {
329
      if (num_rows(D2UhAtQP) > 0)
330 331 332 333
	for (int iq = 0; iq < nPoints; iq++)
	  result[iq] += (*factor) * D2UhAtQP[iq][xi][xj] * fac;
    }

334
    /// Implemetation of SecondOrderTerm::weakEval().
335
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
336
		  std::vector<WorldVector<double> > &result)
337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353
    {
      int nPoints = grdUhAtQP.size();
      for (int iq = 0; iq < nPoints; iq++)
	result[iq][xi] += (*factor) * grdUhAtQP[iq][xj];
    }

  private:
    /// Directions for the derivatives.
    int xi, xj;

    /// Pointer to the factor.
    double *factor;
  };


  /**
   * \ingroup Assembler
354
   *
355 356 357 358 359
   * \brief
   * Laplace operator multiplied with a function evaluated at the quadrature
   * points of a given DOFVector:
   * \f$ f(v(\vec{x})) \Delta u(\vec{x}) \f$
   */
360
  class VecAtQP_SOT : public SecondOrderTerm
361
  {
362 363
  public:
    /// Constructor.
364
    VecAtQP_SOT(DOFVectorBase<double> *dv, AbstractFunction<double, double> *af = NULL,
365
		double factor_ = 1.0);
366 367 368

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
369
		     Quadrature *quad = NULL);
370 371 372 373 374

    /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
    void initElement(const ElInfo* smallElInfo,
		     const ElInfo* largeElInfo,
		     SubAssembler* subAssembler,
375
		     Quadrature *quad = NULL);
376 377

    /// Implements SecondOrderTerm::getLALt().
378
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
379

380
    /// Implemetation of SecondOrderTerm::eval().
381
    void eval(int nPoints,
382 383 384 385
	      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,
386 387
	      double factor);

388
    /// Implemetation of SecondOrderTerm::weakEval().
389 390 391 392 393 394 395 396
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
397
    mtl::dense_vector<double> vecAtQPs;
398 399 400

    /// Function evaluated at \ref vecAtQPs.
    AbstractFunction<double, double> *f;
401

402
    double factor;
403 404 405 406
  };

  /**
   * \ingroup Assembler
407
   *
408 409 410 411 412
   * \brief
   * Laplace operator multiplied with a function evaluated at the quadrature
   * points of a given DOFVector:
   * \f$ f(v(\vec{x}), w(\vec{x})) \Delta u(\vec{x}) \f$
   */
413
  class Vec2AtQP_SOT : public SecondOrderTerm
414
  {
415 416
  public:
    /// Constructor.
417
    Vec2AtQP_SOT(DOFVectorBase<double> *dv1,
418
		 DOFVectorBase<double> *dv2,
419
		 BinaryAbstractFunction<double, double, double> *af = NULL,
420
		double factor_ = 1.0);
421 422 423

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
424
		     Quadrature *quad = NULL);
425 426

    /// Implements SecondOrderTerm::getLALt().
427 428
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;

429
    /// Implemetation of SecondOrderTerm::eval().
430
    void eval(int nPoints,
431 432 433 434
	      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,
435
	      double factor);
436

437
    /// Implemetation of SecondOrderTerm::weakEval().
438 439
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);
440

441 442
  protected:
    /// DOFVector to be evaluated at quadrature points.
443
    DOFVectorBase<double> *vec1, *vec2;
444

445
    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
446
    mtl::dense_vector<double> vecAtQPs1, vecAtQPs2;
447

448 449
    /// Function evaluated at \ref vecAtQPs.
    BinaryAbstractFunction<double, double, double> *f;
450

451
    double factor;
452 453 454 455 456
  };


  /**
   * \ingroup Assembler
457
   *
458 459 460 461
   * \brief
   * Laplace multiplied with a function evaluated at each quadrature point:
   * \f$ f(\vec{x}) \Delta u(\vec{x}) \f$
   */
462
  class CoordsAtQP_SOT : public SecondOrderTerm
463
  {
464 465 466 467 468 469 470 471 472 473
  public:
    /// Constructor.
    CoordsAtQP_SOT(AbstractFunction<double, WorldVector<double> > *af)
      : SecondOrderTerm(af->getDegree()), g(af)
    {
      setSymmetric(true);
    }

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
474
		     Quadrature *quad = NULL);
475 476

    /// Implements SecondOrderTerm::getLALt().
477
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
478

479
    /// Implemetation of SecondOrderTerm::eval().
480
    void eval(int nPoints,
481 482 483 484
	      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,
485 486
	      double factor);

487
    /// Implemetation of SecondOrderTerm::weakEval().
488 489 490 491 492 493
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);


  protected:
    /// Stores coordinates at quadrature points. Set in \ref initElement().
494
    mtl::dense_vector<WorldVector<double> > coordsAtQPs;
495 496 497 498 499 500 501 502 503 504

    /// Function evaluated at quadrature points.
    AbstractFunction<double, WorldVector<double> > *g;
  };


  /**
   * \ingroup Assembler
   *
   * \brief
505
   * SecondOrderTerm where A is a function which maps the gradient of a
506 507 508 509 510 511 512 513 514 515 516 517 518 519
   * DOFVector at each quadrature point to WorldMatrix<double>:
   * \f$ \nabla \cdot A(\nabla v(\vec{x})) \nabla u(\vec{x})\f$
   */
  class MatrixGradient_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    MatrixGradient_SOT(DOFVectorBase<double> *dv,
		       AbstractFunction<WorldMatrix<double>, WorldVector<double> > *af,
		       AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
		       bool symm = false);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
520
		     Quadrature *quad = NULL);
521 522

    /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
523
    void initElement(const ElInfo* smallElInfo,
524 525
		     const ElInfo* largeElInfo,
		     SubAssembler* subAssembler,
526
		     Quadrature *quad = NULL);
527 528

    /// Implements SecondOrderTerm::getLALt().
529
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
530

531
    /// Implemetation of SecondOrderTerm::eval().
532
    void eval(int nPoints,
533 534 535 536
	      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,
537 538
	      double factor);

539
    /// Implemetation of SecondOrderTerm::weakEval().
540 541 542 543 544 545 546 547 548 549 550 551 552 553 554
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    DOFVectorBase<double>* vec;

    /// Function which maps each entry in \ref gradAtQPs to a WorldMatrix<double>.
    AbstractFunction<WorldMatrix<double>, WorldVector<double> > *f;

    AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divFct;

    /** \brief
     * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
     * at each quadrature point.
     */
555
    mtl::dense_vector<WorldVector<double> > gradAtQPs;
556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577

    /// True, if \ref f provides always a symmetric WorldMatrix<double>.
    bool symmetric;

  private:
    /// Temporary matrix used in \ref MatrixGradient_SOT::weakEval.
    WorldMatrix<double> tmpMat;
  };


  /**
   * \ingroup Assembler
   *
   * \brief
   * Laplace multiplied with a function which maps the gradient of a DOFVector
   * at each quadrature point to a double:
   * \f$ f(\nabla v(\vec{x})) \Delta u(\vec{x}) \f$
   */
  class FctGradient_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
578
    FctGradient_SOT(DOFVectorBase<double> *dv,
579 580 581 582
		    AbstractFunction<double, WorldVector<double> > *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
583
		     Quadrature *quad = NULL);
584 585

    /// Implements SecondOrderTerm::getLALt().
586
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
587

588
    /// Implemetation of SecondOrderTerm::eval().
589
    void eval(int nPoints,
590 591 592 593
	      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,
594 595
	      double factor);

596
    /// Implemetation of SecondOrderTerm::weakEval().
597 598 599 600 601 602 603 604 605 606 607 608 609
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    DOFVectorBase<double>* vec;

    /// Function wich maps \ref gradAtQPs to a double.
    AbstractFunction<double, WorldVector<double> > *f;

    /** \brief
     * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
     * at each quadrature point.
     */
610
    mtl::dense_vector<WorldVector<double> > gradAtQPs;
611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632
  };


  /**
   * \ingroup Assembler
   *
   * \brief
   * Laplace multiplied with a function which maps the gradient of a DOFVector
   * at each quadrature point to a double:
   * \f$ \nabla \cdot A(v(\vec{x}), \nabla v(\vec{x})) \nabla u(\vec{x}) \f$
   */
  class VecMatrixGradientAtQP_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    VecMatrixGradientAtQP_SOT(DOFVectorBase<double> *dv,
			      BinaryAbstractFunction<WorldMatrix<double>, double, WorldVector<double> > *af,
			      AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
			      bool symm = false);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
633
		     Quadrature *quad = NULL);
634 635

    /// Implements SecondOrderTerm::getLALt().
636
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
637

638
    /// Implemetation of SecondOrderTerm::eval().
639
    void eval(int nPoints,
640 641 642 643
	      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,
644 645
	      double factor);

646
    /// Implemetation of SecondOrderTerm::weakEval().
647 648
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);
649

650 651
  protected:
    DOFVectorBase<double>* vec;
652

653
    /// Function wich maps \ref gradAtQPs to a double.
654
    BinaryAbstractFunction<WorldMatrix<double>,
655
			   double, WorldVector<double> > *f;
656

657
    AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divFct;
658

659 660 661 662
    /** \brief
     * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
     * at each quadrature point.
     */
663 664
    mtl::dense_vector<double> vecAtQPs;

665
    mtl::dense_vector<WorldVector<double> > gradAtQPs;
666

667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688
    bool symmetric;
  };


  /**
   * \ingroup Assembler
   *
   * \brief
   * Laplace multiplied with a function which maps the gradient of a DOFVector
   * at each quadrature point to a double:
   * \f$ f(\nabla v(\vec{x})) \Delta u(\vec{x}) \f$
   */
  class VecGradCoordsAtQP_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    VecGradCoordsAtQP_SOT(DOFVectorBase<double> *dv,
			  TertiaryAbstractFunction<double, double,
			  WorldVector<double>, WorldVector<double> > *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
689
		     Quadrature *quad = NULL);
690 691

    /// Implements SecondOrderTerm::getLALt().
692
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
693

694
    /// Implemetation of SecondOrderTerm::eval().
695
    void eval(int nPoints,
696 697 698 699
	      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,
700 701
	      double factor);

702
    /// Implemetation of SecondOrderTerm::weakEval().
703 704 705 706 707 708 709 710 711 712
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    DOFVectorBase<double>* vec;

    /// Function wich maps \ref gradAtQPs to a double.
    TertiaryAbstractFunction<double, double, WorldVector<double>,
			     WorldVector<double> > *f;

713 714
    /// Pointer to a WorldVector<double> array containing the gradients of the
    /// DOFVector at each quadrature point.
715 716
    mtl::dense_vector<double> vecAtQPs;

717
    mtl::dense_vector<WorldVector<double> > gradAtQPs;
718

719
    mtl::dense_vector<WorldVector<double> > coordsAtQPs;
720 721 722 723 724
  };


  /**
   * \ingroup Assembler
725
   *
726 727 728 729 730 731 732 733
   * \brief
   * Laplace operator multiplied with a function evaluated at the quadrature
   * points of a given DOFVector:
   * \f$ -f(v(\vec{x})) \Delta u(\vec{x}) \f$
   */
  class VecAndCoordsAtQP_SOT : public SecondOrderTerm {
  public:
    /// Constructor.
734
    VecAndCoordsAtQP_SOT(DOFVectorBase<double> *dv,
735 736 737 738
			 BinaryAbstractFunction<double, double, WorldVector<double> > *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
739
		     Quadrature *quad = NULL);
740 741

    /// Implements SecondOrderTerm::eval().
742
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
743 744

    void eval(int nPoints,
745 746 747 748
	      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,
749 750 751 752 753 754 755 756 757 758
	      double factor);

    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
759
    mtl::dense_vector<double> vecAtQPs;
760

761
    mtl::dense_vector<WorldVector<double> > coordsAtQPs;
762 763 764 765 766 767 768 769 770 771

    /// Function evaluated at \ref vecAtQPs.
    BinaryAbstractFunction<double, double,  WorldVector<double> > *f;
  };


  /**
   * \ingroup Assembler
   *
   * \brief
772
   * SecondOrderTerm where A is a function which maps the gradient of a
773 774 775 776 777 778 779 780 781
   * DOFVector at each quadrature point to WorldMatrix<double>:
   * \f$ \nabla \cdot A(\nabla v(\vec{x})) \nabla u(\vec{x})\f$
   */
  class MatrixGradientAndCoords_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    MatrixGradientAndCoords_SOT(DOFVectorBase<double> *dv,
				BinaryAbstractFunction<WorldMatrix<double>,
782
				WorldVector<double> ,
783 784 785 786 787 788 789
				WorldVector<double> > *af,
				AbstractFunction<WorldVector<double>,
				WorldMatrix<double> > *divAf,
				bool symm = false);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
790
		     Quadrature *quad = NULL);
791 792

    /// Implements SecondOrderTerm::getLALt().
793
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
794

795
    /// Implemetation of SecondOrderTerm::eval().
796
    void eval(int nPoints,
797 798 799 800
	      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,
801 802
	      double factor);

803
    /// Implemetation of SecondOrderTerm::weakEval().
804 805 806 807 808 809 810 811 812 813 814 815
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    DOFVectorBase<double>* vec;

    /// Function which maps each entry in \ref gradAtQPs to a WorldMatrix<double>.
    BinaryAbstractFunction<WorldMatrix<double>,
			   WorldVector<double>, WorldVector<double> > *f;

    AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divFct;

816 817
    /// Pointer to a WorldVector<double> array containing the gradients of the
    /// DOFVector at each quadrature point.
818
    mtl::dense_vector<WorldVector<double> > gradAtQPs;
819

820
    mtl::dense_vector<WorldVector<double> > coordsAtQPs;
821 822 823 824 825 826

    /// True, if \ref f provides always a symmetric WorldMatrix<double>.
    bool symmetric;
  };


827
  class MatrixVec2_SOT : public SecondOrderTerm
828 829 830
  {
  public:
    /// Constructor.
831
    MatrixVec2_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
832 833 834 835 836 837
		  BinaryAbstractFunction<double, double, double> *f,
		  WorldMatrix<double> Af,
		  bool sym = false);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
838
		     Quadrature *quad = NULL);
839 840

    /// Implements SecondOrderTerm::getLALt().
841
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
842

843
    /// Implemetation of SecondOrderTerm::eval().
844
    void eval(int nPoints,
845 846 847 848
	      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,
849 850
	      double factor);

851
    /// Implemetation of SecondOrderTerm::weakEval().
852 853 854 855 856
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    /// DOFVector to be evaluated at quadrature points.
857
    DOFVectorBase<double> *vec1, *vec2;
858 859

    /// Pointer to the values of the DOFVector at quadrature points.
860
    mtl::dense_vector<double> vec1AtQPs, vec2AtQPs;
861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878

    /// Function for A.
    BinaryAbstractFunction<double, double, double> * fct;

    ///
    WorldMatrix<double> A;

    /// True, if \ref matrixFct produces always symmetric matrices.
    bool symmetric;
  };


  class General_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    General_SOT(std::vector<DOFVectorBase<double>*> vecs,
		std::vector<DOFVectorBase<double>*> grads,
879
		TertiaryAbstractFunction<WorldMatrix<double>,
880
		  WorldVector<double>,
881
		  std::vector<double>,
882
		  std::vector<WorldVector<double> > > *f,
883
		AbstractFunction<WorldVector<double>,
884 885 886 887
		  WorldMatrix<double> > *divFct,
		bool symmetric);

    /// Implementation of \ref OperatorTerm::initElement().
888
    void initElement(const ElInfo*, SubAssembler*, Quadrature *quad= NULL);
889 890

    /// Implements SecondOrderTerm::getLALt().
891
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
892

893
    /// Implemetation of SecondOrderTerm::eval().
894
    void eval(int nPoints,
895 896 897 898
	      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,
899 900
	      double factor);

901
    /// Implemetation of SecondOrderTerm::weakEval().
902 903 904 905
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
906
    std::vector<DOFVectorBase<double>*> vecs_;
907 908 909

    std::vector<DOFVectorBase<double>*> grads_;

910
    TertiaryAbstractFunction<WorldMatrix<double>,
911
			     WorldVector<double>,
912
			     std::vector<double>,
913 914
			     std::vector<WorldVector<double> > > *f_;

915
    AbstractFunction<WorldVector<double>,
916 917
		     WorldMatrix<double> > *divFct_;

918
    mtl::dense_vector<WorldVector<double> > coordsAtQPs;
919

920
    std::vector<mtl::dense_vector<double> > vecsAtQPs;
921

922
    std::vector<mtl::dense_vector<WorldVector<double> > > gradsAtQPs_;
923 924 925 926 927 928 929 930 931 932 933

    bool symmetric_;
  };


  class GeneralParametric_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    GeneralParametric_SOT(std::vector<DOFVectorBase<double>*> vecs,
			  std::vector<DOFVectorBase<double>*> grads,
934
			  QuartAbstractFunction<WorldMatrix<double>,
935 936
			  WorldVector<double>,
			  WorldVector<double>,
937
			  std::vector<double>,
938
			  std::vector<WorldVector<double> > > *f,
939
			  AbstractFunction<WorldVector<double>,
940 941 942 943
			  WorldMatrix<double> > *divFct,
			  bool symmetric);

    /// Implementation of \ref OperatorTerm::initElement().
944
    void initElement(const ElInfo*,
945
		     SubAssembler* ,
946
		     Quadrature *quad= NULL);
947 948

    /// Implements SecondOrderTerm::getLALt().
949
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
950

951
    /// Implemetation of SecondOrderTerm::eval().
952
    void eval(int nPoints,
953 954 955 956
	      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,
957 958
	      double factor);

959
    /// Implemetation of SecondOrderTerm::weakEval().
960 961 962 963
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
964
    std::vector<DOFVectorBase<double>*> vecs_;
965 966 967

    std::vector<DOFVectorBase<double>*> grads_;

968
    QuartAbstractFunction<WorldMatrix<double>,
969 970
			  WorldVector<double>,
			  WorldVector<double>,
971
			  std::vector<double>,
972 973
			  std::vector<WorldVector<double> > > *f_;

974
    AbstractFunction<WorldVector<double>,
975 976 977 978
		     WorldMatrix<double> > *divFct_;

    WorldVector<double> elementNormal_;

979 980
    mtl::dense_vector<WorldVector<double> > coordsAtQPs;

981
    std::vector<mtl::dense_vector<double> > vecsAtQPs;
982

983
    std::vector<mtl::dense_vector<WorldVector<double> > > gradsAtQPs_;
984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005

    bool symmetric_;
  };


  /**
   * \ingroup Assembler
   *
   * \brief
   * Laplace multiplied with a function which maps the gradient of a DOFVector
   * at each quadrature point to a double:
   * \f$ f(\nabla v(\vec{x})) \Delta u(\vec{x}) \f$
   */
  class VecAndGradAtQP_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    VecAndGradAtQP_SOT(DOFVectorBase<double> *dv,
		       BinaryAbstractFunction<double, double, WorldVector<double> > *af);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
1006
		     Quadrature *quad = NULL);
1007 1008 1009


    /// Implements SecondOrderTerm::getLALt().
1010
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
1011 1012 1013

    /// Implements SecondOrderTerm::eval().
    void eval(int nPoints,
1014 1015 1016 1017
	      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,
1018 1019 1020 1021 1022
	      double factor);

    /// Implements SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);
1023

1024 1025 1026 1027 1028
  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Vector v at quadrature points.
1029
    mtl::dense_vector<double> vecAtQPs;
1030 1031

    /// Gradient at quadrature points.
1032
    mtl::dense_vector<WorldVector<double> > gradAtQPs;
1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049

    /// Function for c.
    BinaryAbstractFunction<double, double, WorldVector<double> > *f;
  };


  /**
   * \ingroup Assembler
   *
   * \brief
   * Second order term: \f$ f(v(\vec{x}),\nabla w(\vec{x})) \Delta u(\vec{x})\f$
   */
  class VecGrad_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    VecGrad_SOT( DOFVectorBase<double> *dv1,
1050 1051
		 DOFVectorBase<double> *dv2,
		 BinaryAbstractFunction<double, double, WorldVector<double> > *f_)
1052
      : SecondOrderTerm(f_->getDegree()), vec1(dv1), vec2(dv2), f(f_)
1053 1054 1055 1056
    {}

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
1057
		     Quadrature *quad = NULL);
1058 1059

    /// Implements SecondOrderTerm::getLALt().
1060
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
1061 1062 1063

    /// Implements SecondOrderTerm::eval().
    void eval(int nPoints,
1064 1065 1066 1067
	      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,
1068 1069 1070 1071 1072
	      double factor);

    /// Implements SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);
1073

1074 1075 1076 1077 1078 1079
  protected:
    /// DOFVectors to be evaluated at quadrature points.
    DOFVectorBase<double>* vec1;
    DOFVectorBase<double>* vec2;

    /// Vector v at quadrature points.
1080
    mtl::dense_vector<double> vecAtQPs;
1081 1082

    /// Gradient at quadrature points.
1083
    mtl::dense_vector<WorldVector<double> > gradAtQPs;
1084 1085 1086 1087 1088 1089

    /// Function for c.
    BinaryAbstractFunction<double, double, WorldVector<double> > *f;
  };


1090
  /**
1091 1092 1093 1094 1095 1096 1097
   * \ingroup Assembler
   *
   * \brief
   * SecondOrderTerm where A(x) is a WorldMatrix having a in all positions
   * except possibly the position IJ
   * \f$ \nabla \cdot A(x) \nabla u(\vec{x}) \f$
   */
1098
  class CoordsAtQP_IJ_SOT : public SecondOrderTerm
1099 1100 1101 1102
  {
  public:
    /// Constructor.
    CoordsAtQP_IJ_SOT(AbstractFunction<double, WorldVector<double> > *af,
1103
		      int x_i, int x_j)
1104 1105 1106 1107 1108 1109 1110
      : SecondOrderTerm(af->getDegree()), g(af), xi(x_i), xj(x_j)
    {
      setSymmetric(xi == xj);
    }

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
1111
		     Quadrature *quad = NULL);
1112 1113

    /// Implements SecondOrderTerm::getLALt().
1114
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
1115 1116 1117

    /// Implements SecondOrderTerm::eval().
    void eval(int nPoints,
1118 1119 1120 1121
	      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,
1122 1123 1124 1125 1126
	      double factor);

    /// Implements SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);
1127

1128 1129
  private:
    /// Stores coordinates at quadrature points. Set in \ref initElement().
1130
    mtl::dense_vector<WorldVector<double> > coordsAtQPs;
1131 1132 1133 1134 1135 1136 1137 1138 1139 1140

    /// Function evaluated at quadrature points.
    AbstractFunction<double, WorldVector<double> > *g;

    /// Directions for the derivatives.
    int xi, xj;
  };

  /**
   * \ingroup Assembler
1141
   *
1142 1143 1144 1145 1146 1147 1148 1149 1150 1151
   * \brief
   * SecondOrderTerm where A is a WorldMatrix having a in all positions
   * except possibly the position IJ, multiplied with a function
   * evaluated at the quadrature points of a given DOFVector:
   * \f$ \nabla \cdot f(v(\vec{x})) A \nabla u(\vec{x}) \f$
   */
  class VecAtQP_IJ_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
1152 1153
    VecAtQP_IJ_SOT(DOFVectorBase<double> *dv,
		   AbstractFunction<double, double> *af,
1154 1155 1156 1157
		   int x_i, int x_j);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
1158
		     Quadrature *quad = NULL);
1159 1160

    /// Implements SecondOrderTerm::getLALt().
1161
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
1162 1163 1164

    /// Implements SecondOrderTerm::eval().
    void eval(int nPoints,
1165 1166 1167 1168
	      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,
1169 1170 1171 1172
	      double factor);

    /// Implements SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
1173 1174
		  std::vector<WorldVector<double> > &result);

1175 1176 1177 1178 1179
  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec;

    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
1180
    mtl::dense_vector<double> vecAtQPs;
1181 1182 1183 1184 1185 1186 1187 1188 1189

    /// Function evaluated at \ref vecAtQPs.
    AbstractFunction<double, double> *f;

  private:
    /// Directions for the derivatives.
    int xi, xj;
  };

1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209
  /**
   * \ingroup Assembler
   *
   * \brief
   * SecondOrderTerm where A is a WorldMatrix having a in all positions
   * except possibly the position IJ, multiplied with a function
   * evaluated at the quadrature points of two given DOFVectors:
   * \f$ \nabla \cdot f(v(\vec{x}),w(\vec{x})) A \nabla u(\vec{x}) \f$
   */
  class Vec2AtQP_IJ_SOT : public SecondOrderTerm
  {
  public:
    /// Constructor.
    Vec2AtQP_IJ_SOT(DOFVectorBase<double> *dv1,
		   DOFVectorBase<double> *dv2,
		   BinaryAbstractFunction<double, double, double> *af,
		   int x_i, int x_j);

    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
1210
		     Quadrature *quad = NULL);
1211 1212

    /// Implements SecondOrderTerm::getLALt().
1213
    void getLALt(const ElInfo *elInfo, std::vector<mtl::dense2D<double> > &LALt) const;
1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242

    /// Implements SecondOrderTerm::eval().
    void 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 factor);

    /// Implements SecondOrderTerm::weakEval().
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
		  std::vector<WorldVector<double> > &result);

  protected:
    /// DOFVector to be evaluated at quadrature points.
    DOFVectorBase<double>* vec1;
    DOFVectorBase<double>* vec2;

    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
    mtl::dense_vector<double> vecAtQPs1;
    mtl::dense_vector<double> vecAtQPs2;

    /// Function evaluated at \ref vecAtQPs.
    BinaryAbstractFunction<double, double, double> *f;

  private:
    /// Directions for the derivatives.
    int xi, xj;
  };
1243 1244 1245
}

#endif