Operator.h 106 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
573
574
575

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

    /** \brief
     * Implenetation of SecondOrderTerm::weakEval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
591
    void weakEval(int nPoints,
592
593
594
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const
    {
595
596
      if (grdUhAtQP)
	for (int iq = 0; iq < nPoints; iq++)
597
	  result[iq][xi] += (*factor) * grdUhAtQP[iq][xj];
598
    }
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621

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

    /** \brief
     * Pointer to the factor.
     */
    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
622
623
    /// Constructor.
    VecAtQP_SOT(DOFVectorBase<double> *dv, AbstractFunction<double, double> *af);
624

Thomas Witkowski's avatar
Thomas Witkowski committed
625
    /// Implementation of \ref OperatorTerm::initElement().
626
627
628
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

Thomas Witkowski's avatar
Thomas Witkowski committed
629
630
631
632
633
    /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
    void initElement(const ElInfo* smallElInfo,
		     const ElInfo* largeElInfo,
		     SubAssembler* subAssembler,
		     Quadrature *quad = NULL);
634

Thomas Witkowski's avatar
Thomas Witkowski committed
635
636
637
638
639
640
    /// 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,
641
642
643
644
645
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
646
647
    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(int nPoints,
648
649
650
651
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
655
    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
656
657
    double* vecAtQPs;

Thomas Witkowski's avatar
Thomas Witkowski committed
658
    /// Function evaluated at \ref vecAtQPs.
659
660
661
    AbstractFunction<double, double> *f;
  };

662
663
664
665
666
667
668
669
670
671
  /**
   * \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
672
    /// Constructor.
673
    Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, 
Thomas Witkowski's avatar
Thomas Witkowski committed
674
		 BinaryAbstractFunction<double, double, double> *af);
675

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

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

709
710
711
712
713
714
715
716
717
  /**
   * \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
718
719
720
    /// Constructor.
    CoordsAtQP_SOT(AbstractFunction<double, WorldVector<double> > *af)
      : SecondOrderTerm(af->getDegree()), g(af)
721
722
    {
      setSymmetric(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
723
    }
724
725
726
727
728
729
730
731
732
733

    /** \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
734
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;    
735
736
737
738

    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
739
740
    void eval(int nPoints,
	      const double *uhAtQP,
741
742
743
744
745
746
747
748
	      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
749
    void weakEval(int nPoints,
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;


  protected:
    /** \brief
     * Stores coordinates at quadrature points. Set in \ref initElement().
     */
    WorldVector<double>* coordsAtQPs;

    /** \brief
     * Function evaluated at quadrature points.
     */
    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
777
    /// Constructor.
778
    MatrixGradient_SOT(DOFVectorBase<double> *dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
779
780
781
		       AbstractFunction<WorldMatrix<double>, WorldVector<double> > *af,
		       AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
		       bool symm = false);
782
783
784
785
786
787
788
789
790
791

    /** \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
792
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
793
794
795
796

    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
797
798
    void eval(int nPoints,
	      const double *uhAtQP,
799
800
801
802
803
804
805
806
	      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
807
    void weakEval(int nPoints,
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
		  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
844
845
846
    /// Constructor.
    FctGradient_SOT(DOFVectorBase<double> *dv,  
		    AbstractFunction<double, WorldVector<double> > *af);
847
848
849
850
851
852
853
854
855
856

    /** \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
857
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
858
859
860
861

    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
862
863
    void eval(int nPoints,
	      const double *uhAtQP,
864
865
866
867
868
869
870
871
	      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
872
    void weakEval(int nPoints,
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
		  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
903
    /// Constructor.
904
    VecAndGradientAtQP_SOT(DOFVectorBase<double> *dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
905
			   BinaryAbstractFunction<double, double, WorldVector<double> > *af);
906
907
908
909
910
911
912
913
914
915

    /** \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
916
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
917
918
919
920

    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
921
922
    void eval(int nPoints,
	      const double *uhAtQP,
923
924
925
926
927
928
929
930
	      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
931
    void weakEval(int nPoints,
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
		  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;
  };

951
952
953
954
955
956
957
958
959
960
961
  /**
   * \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
962
    /// Constructor.
963
    VecMatrixGradientAtQP_SOT(DOFVectorBase<double> *dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
964
965
966
			      BinaryAbstractFunction<WorldMatrix<double>, double, WorldVector<double> > *af,
			      AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
			      bool symm = false);
967
968
969
970
971
972
973
974
975
976

    /** \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
977
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
978
979
980
981

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