Operator.h 101 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
    std::vector<const FiniteElemSpace*>& getAuxFeSpaces() 
80
    {
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
    /// Implementation of \ref OperatorTerm::initElement().
764
765
766
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

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

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

778
    /// Implenetation of SecondOrderTerm::weakEval().
Thomas Witkowski's avatar
Thomas Witkowski committed
779
    void weakEval(int nPoints,
780
781
782
783
784
785
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;

  protected:
    DOFVectorBase<double>* vec;

786
    /// Function which maps each entry in \ref gradAtQPs to a WorldMatrix<double>.
787
788
789
790
791
792
793
794
795
796
    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;

797
    /// True, if \ref f provides always a symmetric WorldMatrix<double>.
798
799
800
801
802
803
804
805
806
807
808
809
810
811
    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
812
813
814
    /// Constructor.
    FctGradient_SOT(DOFVectorBase<double> *dv,  
		    AbstractFunction<double, WorldVector<double> > *af);
815

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

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

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

831
    /// Implenetation of SecondOrderTerm::weakEval().
Thomas Witkowski's avatar
Thomas Witkowski committed
832
    void weakEval(int nPoints,
833
834
835
836
837
838
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;

  protected:
    DOFVectorBase<double>* vec;

839
    /// Function wich maps \ref gradAtQPs to a double.
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
    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
861
    /// Constructor.
862
    VecAndGradientAtQP_SOT(DOFVectorBase<double> *dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
863
			   BinaryAbstractFunction<double, double, WorldVector<double> > *af);
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
878
879
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor) const;

880
    /// Implenetation of SecondOrderTerm::weakEval().
Thomas Witkowski's avatar
Thomas Witkowski committed
881
    void weakEval(int nPoints,
882
883
884
885
886
887
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;

  protected:
    DOFVectorBase<double>* vec;

888
    /// Function wich maps \ref gradAtQPs to a double.
889
890
891
892
893
894
895
896
897
898
    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;
  };

899
900
901
902
903
904
905
906
907
908
909
  /**
   * \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
910
    /// Constructor.
911
    VecMatrixGradientAtQP_SOT(DOFVectorBase<double> *dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
912
913
914
			      BinaryAbstractFunction<WorldMatrix<double>, double, WorldVector<double> > *af,
			      AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
			      bool symm = false);
915

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

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

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

931
    /// Implenetation of SecondOrderTerm::weakEval().
Thomas Witkowski's avatar
Thomas Witkowski committed
932
    void weakEval(int nPoints,
933
934
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;
Thomas Witkowski's avatar
Thomas Witkowski committed
935
    
936
937
  protected:
    DOFVectorBase<double>* vec;
Thomas Witkowski's avatar
Thomas Witkowski committed
938
    
939
    /// Function wich maps \ref gradAtQPs to a double.
940
941
    BinaryAbstractFunction<WorldMatrix<double>, 
			   double, WorldVector<double> > *f;
Thomas Witkowski's avatar
Thomas Witkowski committed
942
943
944
945
946
947
948
949
950
951
952
953
    
    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;
  };
954
955
956
957
958
959
960
961
962
963
964
965

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

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

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

978
    /// Implenetation of SecondOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
979
980
    void eval(int nPoints,
	      const double *uhAtQP,
981
982
983
984
985
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor) const;

986
    /// Implenetation of SecondOrderTerm::weakEval().
Thomas Witkowski's avatar
Thomas Witkowski committed
987
    void weakEval(int nPoints,
988
989
990
991
992
993
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;

  protected:
    DOFVectorBase<double>* vec;

994
    /// Function wich maps \ref gradAtQPs to a double.
995
996
997
998
999
1000
    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.