Operator.h 100 KB
Newer Older
1
2
3
4
5
6
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
7
// ==  TU Dresden                                                            ==
8
// ==                                                                        ==
9
10
11
// ==  Institut fr Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
12
13
14
15
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
16
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
17
18
19
20
21
22
23
24
25
// ==                                                                        ==
// ============================================================================

/** \file Operator.h */

#ifndef AMDIS_OPERATOR_H
#define AMDIS_OPERATOR_H

#include <vector>
26
#include "AMDiS_fwd.h"
27
28
29
30
31
#include "FixVec.h"
#include "Flag.h"
#include "MatrixVector.h"
#include "ElInfo.h"
#include "AbstractFunction.h"
32
#include "OpenMP.h"
33
#include "SubAssembler.h"
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

namespace AMDiS {

  /** 
   * \ingroup Assembler
   * 
   * \brief
   * Base class for ZeroOrderTerm, FirstOrderTerm and SecondOrderTerm. 
   * OperatorTerms are the building blocks of an Operator. Each OperatorTerm
   * has its properties which are regarded, when constructing 
   * an Assembler for the corresponding Operator.
   */
  class OperatorTerm
  {
  public:
    /** \brief
     * Constructs an OperatorTerm with initially no properties.
     * degree_ is used to determine the degree of the needed quadrature
     * for the assemblage.  
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
54
    OperatorTerm(int deg) 
55
      : properties(0), 
Thomas Witkowski's avatar
Thomas Witkowski committed
56
	degree(deg),
57
58
	dimOfWorld(Global::getGeo(WORLD)),
	bOne(-1)
59
    {}
60

Thomas Witkowski's avatar
Thomas Witkowski committed
61
62
    /// Destructor.
    virtual ~OperatorTerm() {}
63
64
65

    /** \brief
     * Virtual method. It's called by SubAssembler::initElement() for
Thomas Witkowski's avatar
Thomas Witkowski committed
66
67
     * each OperatorTerm belonging to this SubAssembler. E.g., vectors
     * and coordinates at quadrature points can be calculated here.
68
     */
69
    virtual void initElement(const ElInfo*, SubAssembler*, Quadrature *quad = NULL) 
70
    {}
71

72
    virtual void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
73
			     SubAssembler*, Quadrature *quad = NULL)
74
75
    {}

Thomas Witkowski's avatar
Thomas Witkowski committed
76

77
78
    /// Returs \auxFeSpaces, the list of all aux fe spaces the operator makes use off.
    inline std::set<const FiniteElemSpace*>& getAuxFeSpaces() 
79
    {
80
      return auxFeSpaces;
Thomas Witkowski's avatar
Thomas Witkowski committed
81
82
83
    }

    /// Specifies whether the matrix of the term is symmetric
84
85
    void setSymmetric(bool symm);

Thomas Witkowski's avatar
Thomas Witkowski committed
86
    /// Returns true, if the term is piecewise constant, returns false otherwise
87
88
    inline bool isPWConst() 
    { 
89
      return (degree == 0);
90
    }
91

Thomas Witkowski's avatar
Thomas Witkowski committed
92
    /// Returns true, if the term has a symmetric matrix, returns false otherwise.
93
94
    bool isSymmetric();

Thomas Witkowski's avatar
Thomas Witkowski committed
95
    /// Returns \ref degree.
96
97
    inline int getDegree() 
    { 
98
      return degree; 
99
    }
100

101
102
103
104
105
106
    /// Sets one component of the b vector to be one. See \ref bOne.
    void setB(int b)
    {
      bOne = b;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
107
    /// Evaluation of the OperatorTerm at all quadrature points.
Thomas Witkowski's avatar
Thomas Witkowski committed
108
109
    virtual void eval(int nPoints,
		      const double *uhAtQP,
110
111
		      const WorldVector<double> *grdUhAtQP,
		      const WorldMatrix<double> *D2UhAtQP,
Thomas Witkowski's avatar
Thomas Witkowski committed
112
		      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
113
		      double factor) = 0;
114
115

    /** \brief
Thomas Witkowski's avatar
Thomas Witkowski committed
116
117
     * Determines the value of a dof vector at the quadrature points of a given 
     * element. It is used by all VecAtQP like operator terms.
118
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
    double *getVectorAtQPs(DOFVectorBase<double>* vec,
			   const ElInfo* elInfo, 
			   SubAssembler* subAssembler,
			   Quadrature *quad);		   

    /** \brief
     * Determines the value of a dof vector at the quadrature points of a given 
     * element. This function is used, if an operator is assembled on two different
     * meshes using the dual traverse. The element, i.e. the small or the large one,
     * is choosen, which corresponds to the mesh the dof vector is defined on.
     */
    double *getVectorAtQPs(DOFVectorBase<double>* vec,
			   const ElInfo* smallElInfo, 
			   const ElInfo* largeElInfo, 
			   SubAssembler* subAssembler,
			   Quadrature *quad);	

Thomas Witkowski's avatar
Thomas Witkowski committed
136
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
137
138
139
140
141
    WorldVector<double>* getGradientsAtQPs(DOFVectorBase<double>* vec,
					   const ElInfo* elInfo,
					   SubAssembler* subAssembler,
					   Quadrature *quad);

Thomas Witkowski's avatar
Thomas Witkowski committed
142
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
143
144
145
146
147
148
149
150
    WorldVector<double>* getGradientsAtQPs(DOFVectorBase<double>* vec,
					   const ElInfo* smallElInfo, 
					   const ElInfo* largeElInfo,
					   SubAssembler* subAssembler,
					   Quadrature *quad);

  protected:
    /// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$.
151
152
153
154
155
    void lalt(const DimVec<WorldVector<double> >& Lambda,
	      const WorldMatrix<double>& matrix,
	      DimMat<double>& LALt,
	      bool symm,
	      double factor) const;
156
157
158
159
160
161
162
163

    /** \brief
     * 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.
     */
    static void lalt_kl(const DimVec<WorldVector<double> >& Lambda,
			int k, int l,
Thomas Witkowski's avatar
Thomas Witkowski committed
164
165
			DimMat<double>& LALt,
			double factor);
166

167
    /// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to the identity.
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
    inline void l1lt(const DimVec<WorldVector<double> >& Lambda, 
		     DimMat<double>& LALt,
		     double factor) const
    {
      const int dim = LALt.getNumRows();
      
      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;
	}    
      }      
    }
190

Thomas Witkowski's avatar
Thomas Witkowski committed
191
    /// Evaluation of \f$ \Lambda \cdot b\f$.
192
193
194
195
    inline void lb(const DimVec<WorldVector<double> >& Lambda,
		   const WorldVector<double>& b,
		   DimVec<double>& Lb,
		   double factor) const
196
    {
197
      const int dim = Lambda.getSize();
198

199
      for (int i = 0; i < dim; i++) {
200
201
202
203
204
205
206
207
208
	double val = 0.0;
	
	for (int j = 0; j < dimOfWorld; j++)
	  val += Lambda[i][j] * b[j];
	
	Lb[i] += val * factor;
      }    
    }

209
210
211
212
213
    /// Evaluation of \f$ \Lambda \cdot b\f$.
    inline void lb_one(const DimVec<WorldVector<double> >& Lambda,
		       DimVec<double>& Lb,
		       double factor) const
    {
214
      const int dim = Lambda.getSize();
215
216
217
218
219

      for (int i = 0; i < dim; i++)
	Lb[i] += Lambda[i][bOne] * factor;
    }

220
221
222
223
    /** \brief
     * Evaluation of \f$ \Lambda \cdot b\f$ if b contains the value 1.0 in 
     * each component.
     */
224
225
226
    inline void l1(const DimVec<WorldVector<double> >& Lambda,
		   DimVec<double>& Lb,
		   double factor) const
227
    {
228
      const int dim = Lambda.getSize();
229

230
      for (int i = 0; i < dim; i++) {
231
232
	double val = 0.0;
      
233
	for (int j = 0; j < dimOfWorld; j++)
234
	  val += Lambda[i][j];
235

236
	Lb[i] += val * factor;
237
      }    
238
    }
239
240

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
241
    /// Stores the properties of this OperatorTerm
242
243
    Flag properties;

Thomas Witkowski's avatar
Thomas Witkowski committed
244
    /// Polynomial degree of the term. Used to detemine the degree of the quadrature.
245
246
    int degree;

247
248
249
    /// Stores the dimension of the world.
    int dimOfWorld;

Thomas Witkowski's avatar
Thomas Witkowski committed
250
    /// List off all fe spaces, the operator term makes use off.
251
    std::set<const FiniteElemSpace*> auxFeSpaces;
Thomas Witkowski's avatar
Thomas Witkowski committed
252
253

    /// Pointer to the Operator this OperatorTerm belongs to.
254
255
    Operator* operat;

256
257
258
259
260
261
262
263
264
    /** \brief
     * In many cases, the vector b in the evaluation \f$ \Lambda \cdot b\f$ has zeros
     * in all components expect one that is set to one. Using the function \ref lb is
     * then unnecessary time consuming. Instead, this variable defines the component
     * of the vector b to be one. The function \ref lb_one is used if this variable is
     * not -1.
     */
    int bOne;

Thomas Witkowski's avatar
Thomas Witkowski committed
265
    /// Flag for piecewise constant terms
266
267
    static const Flag PW_CONST;

Thomas Witkowski's avatar
Thomas Witkowski committed
268
    /// Flag for symmetric terms
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
    static const Flag SYMMETRIC;

    friend class SubAssembler;
    friend class ZeroOrderAssembler;
    friend class FirstOrderAssembler;
    friend class SecondOrderAssembler;
    friend class Operator;
  };

  /**
   * \ingroup Assembler
   * 
   * \brief
   * Describes the second order terms: \f$ \nabla \cdot A \nabla u(\vec{x}) \f$
   */
  class SecondOrderTerm : public OperatorTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
287
    /// Constructor.
288
289
290
    SecondOrderTerm(int deg) 
      : OperatorTerm(deg) 
    {}
291

292
    /// Destructor.
293
294
    virtual ~SecondOrderTerm() 
    {}
295

296
    /// Evaluation of \f$ \Lambda A \Lambda^t \f$ at all quadrature points.
297
    virtual void getLALt(const ElInfo *elInfo, 
Thomas Witkowski's avatar
Thomas Witkowski committed
298
			 int nPoints, 
299
300
			 DimMat<double> **result) const = 0;

301
    /// Evaluation of \f$ A \nabla u(\vec{x}) \f$ at all quadrature points.
302
    virtual void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
303
			  std::vector<WorldVector<double> > &result) = 0;
304
305
306
307
308
309
310
311
312
313
314
315

  };

  /**
   * \ingroup Assembler
   * 
   * \brief
   * Implements the laplace operator: \f$ \Delta u(\vec{x}) \f$
   */
  class Laplace_SOT : public SecondOrderTerm 
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
316
    /// Constructor.
317
318
319
320
    Laplace_SOT() 
      : SecondOrderTerm(0) 
    {
      setSymmetric(true);
321
    }
322

323
    /// Implenetation of SecondOrderTerm::getLALt().
Thomas Witkowski's avatar
Thomas Witkowski committed
324
    inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
325
326
327
    {
      const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();

328
      for (int iq = 0; iq < nPoints; iq++)
329
	l1lt(Lambda, *(LALt[iq]), 1.0);
330
    }
331

332
    /// Implementation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
333
    inline void eval(int nPoints,
334
335
		     const double *,
		     const WorldVector<double> *,
Thomas Witkowski's avatar
Thomas Witkowski committed
336
337
		     const WorldMatrix<double> *D2UhAtQP,
		     double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
338
		     double factor) 
339
340
341
342
    {
      int dow = Global::getGeo(WORLD);
    
      if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
343
	for (int iq = 0; iq < nPoints; iq++) {
344
	  double resultQP = 0.0;
345
	  for (int i = 0; i < dow; i++) 
346
347
348
349
	    resultQP += D2UhAtQP[iq][i][i];
	  result[iq] += factor * resultQP;
	}
      }
350
    }
351

352
    /// Implenetation of SecondOrderTerm::weakEval().
353
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
354
		  std::vector<WorldVector<double> > &result) 
355
    {
356
357
358
      int nPoints = grdUhAtQP.size();
      for (int iq = 0; iq < nPoints; iq++)
	result[iq] += grdUhAtQP[iq];
359
360
    }
  };
361
362
363
364
365
366
367
368
369
370
371

  /**
   * \ingroup Assembler
   * 
   * \brief
   * Implements the laplace operator multiplied with a scalar factor:
   * \f$ f \cdot \Delta u(\vec{x}) \f$
   */
  class FactorLaplace_SOT : public SecondOrderTerm 
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
372
    /// Constructor.
373
    FactorLaplace_SOT(double f) 
374
375
      : SecondOrderTerm(0)  
    {
376
377
378
379
      factor = new double;
      *factor = f;

      setSymmetric(true);
380
    }
381

Thomas Witkowski's avatar
Thomas Witkowski committed
382
    /// Constructor.
383
    FactorLaplace_SOT(double *fptr) 
384
385
      : SecondOrderTerm(0), 
	factor(fptr)
386
387
    {
      setSymmetric(true);
388
    }
389

390
    /// Implements SecondOrderTerm::getLALt().
Thomas Witkowski's avatar
Thomas Witkowski committed
391
    inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
392
    {
393
      const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
394
      for (int iq = 0; iq < nPoints; iq++) 
395
	l1lt(Lambda, *(LALt[iq]), (*factor));
396
    }
397

398
    /// Implenetation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
399
    void eval(int nPoints,
400
	      const double *,
401
402
403
	      const WorldVector<double> *,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
404
	      double f) 
405
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
406
      int dow = Global::getGeo(WORLD);
407

Thomas Witkowski's avatar
Thomas Witkowski committed
408
      if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
409
	for (int iq = 0; iq < nPoints; iq++) {
410
	  double resultQP = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
411
	  for (int i = 0; i < dow; i++) {
412
413
414
415
416
	    resultQP += D2UhAtQP[iq][i][i];
	  }
	  result[iq] += resultQP * f * (*factor);
	}
      }
417
    }
418

419
    /// Implenetation of SecondOrderTerm::weakEval().
420
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
421
		  std::vector<WorldVector<double> > &result) 
422
    {
423
424
425
      int nPoints = grdUhAtQP.size();
      for (int iq = 0; iq < nPoints; iq++)
	axpy(*factor, grdUhAtQP[iq], result[iq]);
426
    }
427
428

  private:
429
    /// Pointer to the factor.
430
431
432
433
434
435
436
437
438
439
440
441
442
443
    double *factor;
  };

  /**
   * \ingroup Assembler
   * 
   * \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$
   */
  class MatrixFct_SOT : public SecondOrderTerm 
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
444
    /// Constructor.
445
446
447
    MatrixFct_SOT(DOFVectorBase<double> *dv, 
		  AbstractFunction<WorldMatrix<double>, double> *fct,
		  AbstractFunction<WorldVector<double>, WorldMatrix<double> > *div,
Thomas Witkowski's avatar
Thomas Witkowski committed
448
		  bool sym = false);
449

450
    /// Implementation of \ref OperatorTerm::initElement().
451
452
453
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

454
    /// Implements SecondOrderTerm::getLALt().
Thomas Witkowski's avatar
Thomas Witkowski committed
455
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
456
  
457
    /// Implenetation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
458
459
    void eval(int nPoints,
	      const double *uhAtQP,
460
461
462
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
463
	      double factor);
464

465
    /// Implenetation of SecondOrderTerm::weakEval().
466
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
467
		  std::vector<WorldVector<double> > &result);
468
469

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

473
    /// Pointer to the values of the DOFVector at quadrature points.
474
475
    double* vecAtQPs;

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

479
    ///
480
481
    AbstractFunction<WorldVector<double>, WorldMatrix<double> >* divFct;

482
    /// True, if \ref matrixFct produces always symmetric matrices.
483
484
485
486
487
488
489
490
491
492
493
494
    bool symmetric;
  };

  /** 
   * \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:
Thomas Witkowski's avatar
Thomas Witkowski committed
495
    /// Constructor
496
    Matrix_SOT(WorldMatrix<double> mat) 
Thomas Witkowski's avatar
Thomas Witkowski committed
497
      : SecondOrderTerm(0), matrix(mat)
498
499
500
    {
      symmetric = matrix.isSymmetric();
      setSymmetric(symmetric);
501
    }
502

503
    /// Implements SecondOrderTerm::getLALt().
504
505
506
507
    inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
    {
      const DimVec<WorldVector<double> >& Lambda = elInfo->getGrdLambda();
      for (int iq = 0; iq < nPoints; iq++) 
508
	lalt(Lambda, matrix, *(LALt[iq]), symmetric, 1.0);
509
    }
510
  
511
    /// Implenetation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
512
513
    void eval(int nPoints,
	      const double *uhAtQP,
514
515
516
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
517
	      double factor);
518

519
    /// Implenetation of SecondOrderTerm::weakEval().
520
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
521
		  std::vector<WorldVector<double> > &result);
522
523

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

527
    /// True, if \ref matrix is symmetric.
528
529
530
531
532
533
534
535
536
537
538
539
540
541
    bool symmetric;
  };

  /** 
   * \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$
   */
  class FactorIJ_SOT : public SecondOrderTerm 
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
542
    /// Constructor.
543
    FactorIJ_SOT(int x_i, int x_j, double f) 
Thomas Witkowski's avatar
Thomas Witkowski committed
544
      : SecondOrderTerm(0), xi(x_i), xj(x_j)
545
546
547
548
549
    {
      factor = new double;
      *factor = f;

      setSymmetric(xi == xj);
550
    }
551

Thomas Witkowski's avatar
Thomas Witkowski committed
552
    /// Constructor.
553
    FactorIJ_SOT(int x_i, int x_j, double *fptr) 
Thomas Witkowski's avatar
Thomas Witkowski committed
554
      : SecondOrderTerm(0), xi(x_i), xj(x_j), factor(fptr)
555
556
    {
      setSymmetric(xi == xj);
557
    }
558

559
    /// Implements SecondOrderTerm::getLALt().
Thomas Witkowski's avatar
Thomas Witkowski committed
560
    inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
561
    {
562
      const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
563
      for (int iq = 0; iq < nPoints; iq++)
564
	lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*factor));
565
    }
566

Thomas Witkowski's avatar
Thomas Witkowski committed
567
    /// Implenetation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
568
    void eval(int nPoints,
Thomas Witkowski's avatar
Thomas Witkowski committed
569
	      const double *,
570
571
572
	      const WorldVector<double> *,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
573
	      double fac)
574
    {
575
576
      if (D2UhAtQP)
	for (int iq = 0; iq < nPoints; iq++)
577
	  result[iq] += (*factor) * D2UhAtQP[iq][xi][xj] * fac;
578
    }
579

Thomas Witkowski's avatar
Thomas Witkowski committed
580
    /// Implenetation of SecondOrderTerm::weakEval().
581
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
582
		  std::vector<WorldVector<double> > &result) 
583
    {
584
585
586
      int nPoints = grdUhAtQP.size();
      for (int iq = 0; iq < nPoints; iq++)
	result[iq][xi] += (*factor) * grdUhAtQP[iq][xj];
587
    }
588
589

  private:
Thomas Witkowski's avatar
Thomas Witkowski committed
590
    /// Directions for the derivatives.
591
592
    int xi, xj;

Thomas Witkowski's avatar
Thomas Witkowski committed
593
    /// Pointer to the factor.
594
595
596
597
598
599
600
601
602
603
604
605
606
    double *factor;
  };

  /**
   * \ingroup Assembler
   *  
   * \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 VecAtQP_SOT : public SecondOrderTerm {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
607
608
    /// Constructor.
    VecAtQP_SOT(DOFVectorBase<double> *dv, AbstractFunction<double, double> *af);
609

Thomas Witkowski's avatar
Thomas Witkowski committed
610
    /// Implementation of \ref OperatorTerm::initElement().
611
612
613
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

Thomas Witkowski's avatar
Thomas Witkowski committed
614
615
616
617
618
    /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
    void initElement(const ElInfo* smallElInfo,
		     const ElInfo* largeElInfo,
		     SubAssembler* subAssembler,
		     Quadrature *quad = NULL);
619

Thomas Witkowski's avatar
Thomas Witkowski committed
620
621
622
623
624
625
    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;    

    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
626
627
628
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
629
	      double factor);
630

Thomas Witkowski's avatar
Thomas Witkowski committed
631
    /// Implenetation of SecondOrderTerm::weakEval().
632
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
633
		  std::vector<WorldVector<double> > &result);
634
635

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
636
    /// DOFVector to be evaluated at quadrature points.
637
638
    DOFVectorBase<double>* vec;

Thomas Witkowski's avatar
Thomas Witkowski committed
639
    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
640
641
    double* vecAtQPs;

Thomas Witkowski's avatar
Thomas Witkowski committed
642
    /// Function evaluated at \ref vecAtQPs.
643
644
645
    AbstractFunction<double, double> *f;
  };

646
647
648
649
650
651
652
653
654
655
  /**
   * \ingroup Assembler
   *  
   * \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$
   */
  class Vec2AtQP_SOT : public SecondOrderTerm {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
656
    /// Constructor.
657
    Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, 
Thomas Witkowski's avatar
Thomas Witkowski committed
658
		 BinaryAbstractFunction<double, double, double> *af);
659

660
    /// Implementation of \ref OperatorTerm::initElement().
661
662
663
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

664
    /// Implements SecondOrderTerm::getLALt().
Thomas Witkowski's avatar
Thomas Witkowski committed
665
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;    
666
    
667
    /// Implenetation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
668
    void eval(int nPoints,
669
670
671
672
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
673
	      double factor);
674
    
675
    /// Implenetation of SecondOrderTerm::weakEval().
676
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
677
		  std::vector<WorldVector<double> > &result);
678
679
    
  protected:
680
    /// DOFVector to be evaluated at quadrature points.
681
682
683
    DOFVectorBase<double>* vec1;
    DOFVectorBase<double>* vec2;
    
684
    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
685
686
687
    double* vecAtQPs1;
    double* vecAtQPs2;
    
688
    /// Function evaluated at \ref vecAtQPs.
689
690
691
    BinaryAbstractFunction<double, double, double> *f;
  };

692
693
694
695
696
697
698
699
700
  /**
   * \ingroup Assembler
   * 
   * \brief
   * Laplace multiplied with a function evaluated at each quadrature point:
   * \f$ f(\vec{x}) \Delta u(\vec{x}) \f$
   */
  class CoordsAtQP_SOT : public SecondOrderTerm {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
701
702
703
    /// Constructor.
    CoordsAtQP_SOT(AbstractFunction<double, WorldVector<double> > *af)
      : SecondOrderTerm(af->getDegree()), g(af)
704
705
    {
      setSymmetric(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
706
    }
707

Thomas Witkowski's avatar
Thomas Witkowski committed
708
    /// Implementation of \ref OperatorTerm::initElement().
709
710
711
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

Thomas Witkowski's avatar
Thomas Witkowski committed
712
    /// Implements SecondOrderTerm::getLALt().
Thomas Witkowski's avatar
Thomas Witkowski committed
713
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;    
714

Thomas Witkowski's avatar
Thomas Witkowski committed
715
    /// Implenetation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
716
717
    void eval(int nPoints,
	      const double *uhAtQP,
718
719
720
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
721
	      double factor);
722

Thomas Witkowski's avatar
Thomas Witkowski committed
723
    /// Implenetation of SecondOrderTerm::weakEval().
724
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
725
		  std::vector<WorldVector<double> > &result);
726
727
728


  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
729
    /// Stores coordinates at quadrature points. Set in \ref initElement().
730
731
    WorldVector<double>* coordsAtQPs;

Thomas Witkowski's avatar
Thomas Witkowski committed
732
    /// Function evaluated at quadrature points.
733
734
735
736
737
738
739
740
741
742
743
744
745
746
    AbstractFunction<double, WorldVector<double> > *g;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * SecondOrderTerm where A is a function which maps the gradient of a 
   * 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:
Thomas Witkowski's avatar
Thomas Witkowski committed
747
    /// Constructor.
748
    MatrixGradient_SOT(DOFVectorBase<double> *dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
749
750
751
		       AbstractFunction<WorldMatrix<double>, WorldVector<double> > *af,
		       AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
		       bool symm = false);
752

753
    /// Implementation of \ref OperatorTerm::initElement().
754
755
756
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

757
758
759
760
761
762
    /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
    void initElement(const ElInfo* smallElInfo, 
		     const ElInfo* largeElInfo,
		     SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

763
    /// Implements SecondOrderTerm::getLALt().
Thomas Witkowski's avatar
Thomas Witkowski committed
764
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
765

766
    /// Implenetation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
767
768
    void eval(int nPoints,
	      const double *uhAtQP,
769
770
771
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
772
	      double factor);
773

774
    /// Implenetation of SecondOrderTerm::weakEval().
775
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
776
		  std::vector<WorldVector<double> > &result);
777
778
779
780

  protected:
    DOFVectorBase<double>* vec;

781
    /// Function which maps each entry in \ref gradAtQPs to a WorldMatrix<double>.
782
783
784
785
786
787
788
789
790
791
    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.
     */
    WorldVector<double>* gradAtQPs;

792
    /// True, if \ref f provides always a symmetric WorldMatrix<double>.
793
    bool symmetric;
794
795
796
797

  private:
    /// Temporary matrix used in \ref MatrixGradient_SOT::weakEval.
    WorldMatrix<double> tmpMat;
798
799
800
801
802
803
804
805
806
807
808
809
810
  };

  /**
   * \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:
Thomas Witkowski's avatar
Thomas Witkowski committed
811
812
813
    /// Constructor.
    FctGradient_SOT(DOFVectorBase<double> *dv,  
		    AbstractFunction<double, WorldVector<double> > *af);
814

815
    /// Implementation of \ref OperatorTerm::initElement().
816
817
818
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

819
    /// Implements SecondOrderTerm::getLALt().
Thomas Witkowski's avatar
Thomas Witkowski committed
820
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
821

822
    /// Implenetation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
823
824
    void eval(int nPoints,
	      const double *uhAtQP,
825
826
827
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
828
	      double factor);
829

830
    /// Implenetation of SecondOrderTerm::weakEval().
831
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
832
		  std::vector<WorldVector<double> > &result);
833
834
835
836

  protected:
    DOFVectorBase<double>* vec;

837
    /// Function wich maps \ref gradAtQPs to a double.
838
839
840
841
842
843
844
845
846
847
    AbstractFunction<double, WorldVector<double> > *f;

    /** \brief
     * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
     * at each quadrature point.
     */
    WorldVector<double>* gradAtQPs;
  };


848
849
850
851
852
853
854
855
856
857
858
  /**
   * \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:
Thomas Witkowski's avatar
Thomas Witkowski committed
859
    /// Constructor.
860
    VecMatrixGradientAtQP_SOT(DOFVectorBase<double> *dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
861
862
863
			      BinaryAbstractFunction<WorldMatrix<double>, double, WorldVector<double> > *af,
			      AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
			      bool symm = false);
864

865
    /// Implementation of \ref OperatorTerm::initElement().
866
867
868
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

869
    /// Implements SecondOrderTerm::getLALt().
Thomas Witkowski's avatar
Thomas Witkowski committed
870
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
871

872
    /// Implenetation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
873
874
    void eval(int nPoints,
	      const double *uhAtQP,
875
876
877
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
878
	      double factor);
879

880
    /// Implenetation of SecondOrderTerm::weakEval().
881
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
882
		  std::vector<WorldVector<double> > &result);
Thomas Witkowski's avatar
Thomas Witkowski committed
883
    
884
885
  protected:
    DOFVectorBase<double>* vec;
Thomas Witkowski's avatar
Thomas Witkowski committed
886
    
887
    /// Function wich maps \ref gradAtQPs to a double.
888
889
    BinaryAbstractFunction<WorldMatrix<double>, 
			   double, WorldVector<double> > *f;
Thomas Witkowski's avatar
Thomas Witkowski committed
890
891
892
893
894
895
896
897
898
899
900
901
    
    AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divFct;
    
    /** \brief
     * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
     * at each quadrature point.
     */
    double* vecAtQPs;
    WorldVector<double>* gradAtQPs;
    
    bool symmetric;
  };
902
903
904
905
906
907
908
909
910
911
912
913

  /**
   * \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:
Thomas Witkowski's avatar
Thomas Witkowski committed
914
    /// Constructor.
915
916
    VecGradCoordsAtQP_SOT(DOFVectorBase<double> *dv,
			  TertiaryAbstractFunction<double, double,
Thomas Witkowski's avatar
Thomas Witkowski committed
917
			  WorldVector<double>, WorldVector<double> > *af);
918

919
    /// Implementation of \ref OperatorTerm::initElement().
920
921
922
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

923
    /// Implements SecondOrderTerm::getLALt().
Thomas Witkowski's avatar
Thomas Witkowski committed
924
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
925

926
    /// Implenetation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
927
928
    void eval(int nPoints,
	      const double *uhAtQP,
929
930
931
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
932
	      double factor);
933

934
    /// Implenetation of SecondOrderTerm::weakEval().
935
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
936
		  std::vector<WorldVector<double> > &result);
937
938
939
940

  protected:
    DOFVectorBase<double>* vec;

941
    /// Function wich maps \ref gradAtQPs to a double.
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
    TertiaryAbstractFunction<double, double, WorldVector<double>,
			     WorldVector<double> > *f;

    /** \brief
     * Pointer to a WorldVector<double> array containing the gradients of the DOFVector
     * at each quadrature point.
     */
    double* vecAtQPs;
    WorldVector<double>* gradAtQPs;
    WorldVector<double>* coordsAtQPs;
  
  };


  /**
   * \ingroup Assembler
   *  
   * \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:
Thomas Witkowski's avatar
Thomas Witkowski committed
966
    /// Constructor.
967
    VecAndCoordsAtQP_SOT(DOFVectorBase<double> *dv, 
Thomas Witkowski's avatar
Thomas Witkowski committed
968
			 BinaryAbstractFunction<double, double, WorldVector<double> > *af);
969

970
    /// Implementation of \ref OperatorTerm::initElement().
971
972
973
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

974
    /// Implements SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
975
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;    
976

Thomas Witkowski's avatar
Thomas Witkowski committed
977
978
    void eval(int nPoints,
	      const double *uhAtQP,
979
980
981
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
982
	      double factor);
983

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

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

991
    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
992
993
994
995
    double* vecAtQPs;
  
    WorldVector<double>* coordsAtQPs;

996
    /// Function evaluated at \ref vecAtQPs.
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
    BinaryAbstractFunction<double, double,  WorldVector<double> > *f;
  };

  /**
   * \ingroup Assembler
   *
   * \brief
   * SecondOrderTerm where A is a function which maps the gradient of a 
   * 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:
Thomas Witkowski's avatar
Thomas Witkowski committed
1011
    /// Constructor.
1012
1013
1014
    MatrixGradientAndCoords_SOT(DOFVectorBase<double> *dv,
				BinaryAbstractFunction<WorldMatrix<double>,
				WorldVector<double> , 
Thomas Witkowski's avatar
Thomas Witkowski committed
1015
				WorldVector<double> > *af,
1016
				AbstractFunction<WorldVector<double>,
Thomas Witkowski's avatar
Thomas Witkowski committed
1017
1018
				WorldMatrix<double> > *divAf,
				bool symm = false);
1019

1020
    /// Implementation of \ref OperatorTerm::initElement().
1021
1022
1023
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

1024
    /// Implements SecondOrderTerm::getLALt().
Thomas Witkowski's avatar
Thomas Witkowski committed
1025
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
1026

1027
    /// Implenetation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
1028
1029
    void eval(int nPoints,
	      const double *uhAtQP,
1030
1031
1032
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
1033
	      double factor);
1034

1035
    /// Implenetation of SecondOrderTerm::weakEval().
1036
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
1037
		  std::vector<WorldVector<double> > &result);
1038
1039
1040
1041

  protected:
    DOFVectorBase<double>* vec;

1042
    /// Function which maps each entry in \ref gradAtQPs to a WorldMatrix<double>.
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
    BinaryAbstractFunction<WorldMatrix<double>,
			   WorldVector<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.
     */
    WorldVector<double>* gradAtQPs;

    WorldVector<double>* coordsAtQPs;

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

1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
  class MatrixVec2_SOT : public SecondOrderTerm 
  {
  public:
    /// Constructor.
    MatrixVec2_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,  
		  BinaryAbstractFunction<double, double, double> *f,
		  WorldMatrix<double> Af,
		  bool sym = false);

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

    /// Implements SecondOrderTerm::getLALt().
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
  
    /// Implenetation of SecondOrderTerm::eval().
    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
Thomas Witkowski's avatar
Thomas Witkowski committed
1082
	      double factor);
1083
1084

    /// Implenetation of SecondOrderTerm::weakEval().
1085
    void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
1086
		  std::vector<WorldVector<double> > &result);
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107

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

    /// Pointer to the values of the DOFVector at quadrature points.
    double* vec1AtQPs;
    double* vec2AtQPs;

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

    ///
    WorldMatrix<double> A;

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


1108
1109
1110
  class General_SOT : public SecondOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
1111
    /// Constructor.
1112
1113
    General_SOT(std::vector<DOFVectorBase<double>*> vecs,
		std::vector<DOFVectorBase<double>*> grads,
1114
1115
		TertiaryAbstractFunction<WorldMatrix<double>, 
		  WorldVector<double>,
1116
1117
		  std::vector<double>, 
		  std::vector<WorldVector<double> > > *f,
1118
1119
		AbstractFunction<WorldVector<double>, 
		  WorldMatrix<double> > *divFct,
Thomas Witkowski's avatar
Thomas Witkowski committed
1120
		bool symmetric);
1121

1122
1123
    /// Implementation of \ref OperatorTerm::initElement().
    void initElement(const ElInfo*, SubAssembler*, Quadrature *quad= NULL);
1124

1125
    /// Implements SecondOrderTerm::getLALt().
Thomas Witkowski's avatar
Thomas Witkowski committed
1126
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
1127