Operator.h 103 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
59
	dimOfWorld(Global::getGeo(WORLD)),
	auxFESpaces(0),
	bOne(-1)
60
    {}
61

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

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

73
74
    virtual void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
			     SubAssembler*, 
Thomas Witkowski's avatar
Thomas Witkowski committed
75
76
77
78
			     Quadrature *quad = NULL)
    {}

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

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

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

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

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

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

    /** \brief
Thomas Witkowski's avatar
Thomas Witkowski committed
117
118
     * Determines the value of a dof vector at the quadrature points of a given 
     * element. It is used by all VecAtQP like operator terms.
119
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
    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
137
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
138
139
140
141
142
    WorldVector<double>* getGradientsAtQPs(DOFVectorBase<double>* vec,
					   const ElInfo* elInfo,
					   SubAssembler* subAssembler,
					   Quadrature *quad);

Thomas Witkowski's avatar
Thomas Witkowski committed
143
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
144
145
146
147
148
149
150
151
    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$.
152
153
154
155
156
    void lalt(const DimVec<WorldVector<double> >& Lambda,
	      const WorldMatrix<double>& matrix,
	      DimMat<double>& LALt,
	      bool symm,
	      double factor) const;
157
158
159
160
161
162
163
164

    /** \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
165
166
			DimMat<double>& LALt,
			double factor);
167

168
    /// Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to the identity.
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
    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;
	}    
      }      
    }
191

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

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

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

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

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

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

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

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

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

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

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

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

257
258
259
260
261
262
263
264
265
    /** \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
266
    /// Flag for piecewise constant terms
267
268
    static const Flag PW_CONST;

Thomas Witkowski's avatar
Thomas Witkowski committed
269
    /// Flag for symmetric terms
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
    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
288
    /// Constructor.
289
290
291
    SecondOrderTerm(int deg) 
      : OperatorTerm(deg) 
    {}
292

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

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

302
    /// Evaluation of \f$ A \nabla u(\vec{x}) \f$ at all quadrature points.
Thomas Witkowski's avatar
Thomas Witkowski committed
303
    virtual void weakEval(int nPoints,
304
305
306
307
308
309
310
311
312
313
314
315
316
317
			  const WorldVector<double> *grdUhAtQP,
			  WorldVector<double> *result) const = 0;

  };

  /**
   * \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
318
    /// Constructor.
319
320
321
322
    Laplace_SOT() 
      : SecondOrderTerm(0) 
    {
      setSymmetric(true);
323
    }
324

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

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

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

354
    /// Implenetation of SecondOrderTerm::weakEval().
Thomas Witkowski's avatar
Thomas Witkowski committed
355
    void weakEval(int nPoints,
356
357
358
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const
    {
359
      if (grdUhAtQP)
360
	for (int iq = 0; iq < nPoints; iq++)
361
	  result[iq] += grdUhAtQP[iq];
362
363
    }
  };
364
365
366
367
368
369
370
371
372
373
374

  /**
   * \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
375
    /// Constructor.
376
    FactorLaplace_SOT(double f) 
377
378
      : SecondOrderTerm(0)  
    {
379
380
381
382
      factor = new double;
      *factor = f;

      setSymmetric(true);
383
    }
384

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

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

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

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

422
    /// Implenetation of SecondOrderTerm::weakEval().
Thomas Witkowski's avatar
Thomas Witkowski committed
423
    void weakEval(int nPoints,
424
425
426
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const
    {
427
428
      if (grdUhAtQP)
	for (int iq = 0; iq < nPoints; iq++)
429
	  axpy(*factor, grdUhAtQP[iq], result[iq]);
430
    }
431
432

  private:
433
    /// Pointer to the factor.
434
435
436
437
438
439
440
441
442
443
444
445
446
447
    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
448
    /// Constructor.
449
450
451
    MatrixFct_SOT(DOFVectorBase<double> *dv, 
		  AbstractFunction<WorldMatrix<double>, double> *fct,
		  AbstractFunction<WorldVector<double>, WorldMatrix<double> > *div,
Thomas Witkowski's avatar
Thomas Witkowski committed
452
		  bool sym = false);
453

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

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

469
    /// Implenetation of SecondOrderTerm::weakEval().
Thomas Witkowski's avatar
Thomas Witkowski committed
470
    void weakEval(int nPoints,
471
472
473
474
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;

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

478
    /// Pointer to the values of the DOFVector at quadrature points.
479
480
    double* vecAtQPs;

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

484
    ///
485
486
    AbstractFunction<WorldVector<double>, WorldMatrix<double> >* divFct;

487
    /// True, if \ref matrixFct produces always symmetric matrices.
488
489
490
491
492
493
494
495
496
497
498
499
    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
500
    /// Constructor
501
    Matrix_SOT(WorldMatrix<double> mat) 
Thomas Witkowski's avatar
Thomas Witkowski committed
502
      : SecondOrderTerm(0), matrix(mat)
503
504
505
    {
      symmetric = matrix.isSymmetric();
      setSymmetric(symmetric);
506
    }
507

508
    /// Implements SecondOrderTerm::getLALt().
509
510
511
512
    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++) 
513
	lalt(Lambda, matrix, *(LALt[iq]), symmetric, 1.0);
514
    }
515
  
516
    /// Implenetation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
517
518
    void eval(int nPoints,
	      const double *uhAtQP,
519
520
521
522
523
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor) const;

524
    /// Implenetation of SecondOrderTerm::weakEval().
Thomas Witkowski's avatar
Thomas Witkowski committed
525
    void weakEval(int nPoints,
526
527
528
529
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;

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

533
    /// True, if \ref matrix is symmetric.
534
535
536
537
538
539
540
541
542
543
544
545
546
547
    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
548
    /// Constructor.
549
    FactorIJ_SOT(int x_i, int x_j, double f) 
Thomas Witkowski's avatar
Thomas Witkowski committed
550
      : SecondOrderTerm(0), xi(x_i), xj(x_j)
551
552
553
554
555
    {
      factor = new double;
      *factor = f;

      setSymmetric(xi == xj);
556
    }
557

Thomas Witkowski's avatar
Thomas Witkowski committed
558
    /// Constructor.
559
    FactorIJ_SOT(int x_i, int x_j, double *fptr) 
Thomas Witkowski's avatar
Thomas Witkowski committed
560
      : SecondOrderTerm(0), xi(x_i), xj(x_j), factor(fptr)
561
562
    {
      setSymmetric(xi == xj);
563
    }
564

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
586
    /// Implenetation of SecondOrderTerm::weakEval().
Thomas Witkowski's avatar
Thomas Witkowski committed
587
    void weakEval(int nPoints,
588
589
590
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const
    {
591
592
      if (grdUhAtQP)
	for (int iq = 0; iq < nPoints; iq++)
593
	  result[iq][xi] += (*factor) * grdUhAtQP[iq][xj];
594
    }
595
596

  private:
Thomas Witkowski's avatar
Thomas Witkowski committed
597
    /// Directions for the derivatives.
598
599
    int xi, xj;

Thomas Witkowski's avatar
Thomas Witkowski committed
600
    /// Pointer to the factor.
601
602
603
604
605
606
607
608
609
610
611
612
613
    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
614
615
    /// Constructor.
    VecAtQP_SOT(DOFVectorBase<double> *dv, AbstractFunction<double, double> *af);
616

Thomas Witkowski's avatar
Thomas Witkowski committed
617
    /// Implementation of \ref OperatorTerm::initElement().
618
619
620
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

Thomas Witkowski's avatar
Thomas Witkowski committed
621
622
623
624
625
    /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
    void initElement(const ElInfo* smallElInfo,
		     const ElInfo* largeElInfo,
		     SubAssembler* subAssembler,
		     Quadrature *quad = NULL);
626

Thomas Witkowski's avatar
Thomas Witkowski committed
627
628
629
630
631
632
    /// 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,
633
634
635
636
637
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
638
639
    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(int nPoints,
640
641
642
643
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
647
    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
648
649
    double* vecAtQPs;

Thomas Witkowski's avatar
Thomas Witkowski committed
650
    /// Function evaluated at \ref vecAtQPs.
651
652
653
    AbstractFunction<double, double> *f;
  };

654
655
656
657
658
659
660
661
662
663
  /**
   * \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
664
    /// Constructor.
665
    Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, 
Thomas Witkowski's avatar
Thomas Witkowski committed
666
		 BinaryAbstractFunction<double, double, double> *af);
667

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

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

701
702
703
704
705
706
707
708
709
  /**
   * \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
710
711
712
    /// Constructor.
    CoordsAtQP_SOT(AbstractFunction<double, WorldVector<double> > *af)
      : SecondOrderTerm(af->getDegree()), g(af)
713
714
    {
      setSymmetric(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
715
    }
716

Thomas Witkowski's avatar
Thomas Witkowski committed
717
    /// Implementation of \ref OperatorTerm::initElement().
718
719
720
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
724
    /// Implenetation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
725
726
    void eval(int nPoints,
	      const double *uhAtQP,
727
728
729
730
731
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
732
    /// Implenetation of SecondOrderTerm::weakEval().
Thomas Witkowski's avatar
Thomas Witkowski committed
733
    void weakEval(int nPoints,
734
735
736
737
738
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;


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

Thomas Witkowski's avatar
Thomas Witkowski committed
742
    /// Function evaluated at quadrature points.
743
744
745
746
747
748
749
750
751
752
753
754
755
756
    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
757
    /// Constructor.
758
    MatrixGradient_SOT(DOFVectorBase<double> *dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
759
760
761
		       AbstractFunction<WorldMatrix<double>, WorldVector<double> > *af,
		       AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
		       bool symm = false);
762
763
764
765
766
767
768
769
770
771

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

    /** \brief
     * Implements SecondOrderTerm::getLALt().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
772
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
773
774
775
776

    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
777
778
    void eval(int nPoints,
	      const double *uhAtQP,
779
780
781
782
783
784
785
786
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor) const;

    /** \brief
     * Implenetation of SecondOrderTerm::weakEval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
787
    void weakEval(int nPoints,
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;

  protected:
    DOFVectorBase<double>* vec;

    /** \brief
     * 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.
     */
    WorldVector<double>* gradAtQPs;

    /** \brief
     * True, if \ref f provides always a symmetric WorldMatrix<double>.
     */
    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 FctGradient_SOT : public SecondOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
824
825
826
    /// Constructor.
    FctGradient_SOT(DOFVectorBase<double> *dv,  
		    AbstractFunction<double, WorldVector<double> > *af);
827
828
829
830
831
832
833
834
835
836

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

    /** \brief
     * Implements SecondOrderTerm::getLALt().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
837
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
838
839
840
841

    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
842
843
    void eval(int nPoints,
	      const double *uhAtQP,
844
845
846
847
848
849
850
851
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor) const;

    /** \brief
     * Implenetation of SecondOrderTerm::weakEval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
852
    void weakEval(int nPoints,
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;

  protected:
    DOFVectorBase<double>* vec;

    /** \brief
     * 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.
     */
    WorldVector<double>* gradAtQPs;
  };


  /**
   * \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 VecAndGradientAtQP_SOT : public SecondOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
883
    /// Constructor.
884
    VecAndGradientAtQP_SOT(DOFVectorBase<double> *dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
885
			   BinaryAbstractFunction<double, double, WorldVector<double> > *af);
886
887
888
889
890
891
892
893
894
895

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

    /** \brief
     * Implements SecondOrderTerm::getLALt().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
896
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
897
898
899
900

    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
901
902
    void eval(int nPoints,
	      const double *uhAtQP,
903
904
905
906
907
908
909
910
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor) const;

    /** \brief
     * Implenetation of SecondOrderTerm::weakEval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
911
    void weakEval(int nPoints,
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;

  protected:
    DOFVectorBase<double>* vec;

    /** \brief
     * Function wich maps \ref gradAtQPs to a double.
     */
    BinaryAbstractFunction<double, 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;
  };

931
932
933
934
935
936
937
938
939
940
941
  /**
   * \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
942
    /// Constructor.
943
    VecMatrixGradientAtQP_SOT(DOFVectorBase<double> *dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
944
945
946
			      BinaryAbstractFunction<WorldMatrix<double>, double, WorldVector<double> > *af,
			      AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
			      bool symm = false);
947
948
949
950
951
952
953
954
955
956

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

    /** \brief
     * Implements SecondOrderTerm::getLALt().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
957
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
958
959
960
961

    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
962
963
    void eval(int nPoints,
	      const double *uhAtQP,
964
965
966
967
968
969
970
971
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor) const;

    /** \brief
     * Implenetation of SecondOrderTerm::weakEval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
972
    void weakEval(int nPoints,
973
974
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
975
    
976
977
  protected:
    DOFVectorBase<double>* vec;
Thomas Witkowski's avatar
Thomas Witkowski committed
978
    
979
980
981
982
983
    /** \brief
     * Function wich maps \ref gradAtQPs to a double.
     */
    BinaryAbstractFunction<WorldMatrix<double>, 
			   double, WorldVector<double> > *f;
Thomas Witkowski's avatar
Thomas Witkowski committed
984
985
986
987
988
989
990
991
992
993
994
995
    
    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;
  };
996
997
998
999
1000

  /**
   * \ingroup Assembler
   *
   * \brief