Operator.h 96.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file Operator.h */

#ifndef AMDIS_OPERATOR_H
#define AMDIS_OPERATOR_H

#include <vector>
#include "FixVec.h"
#include "Flag.h"
#include "MemoryManager.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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

namespace AMDiS {

  class Assembler;
  class ElInfo;
  class FiniteElemSpace;
  class Operator;
  class SubAssembler;
  class ElementMatrix;
  class ElementVector;
  class Quadrature;
  template<typename T> class DOFVectorBase;

  // ============================================================================
  // ===== class OperatorTerm ===================================================
  // ============================================================================

  /** 
   * \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:
    MEMORY_MANAGED(OperatorTerm);

    /** \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
70
    OperatorTerm(int deg) 
71
      : properties(0), 
Thomas Witkowski's avatar
Thomas Witkowski committed
72
73
	degree(deg),
	auxFESpaces(0)
74
    {}
75

Thomas Witkowski's avatar
Thomas Witkowski committed
76
77
    /// Destructor.
    virtual ~OperatorTerm() {}
78
79
80

    /** \brief
     * Virtual method. It's called by SubAssembler::initElement() for
Thomas Witkowski's avatar
Thomas Witkowski committed
81
82
     * each OperatorTerm belonging to this SubAssembler. E.g., vectors
     * and coordinates at quadrature points can be calculated here.
83
84
     */
    virtual void initElement(const ElInfo*, 
Thomas Witkowski's avatar
Thomas Witkowski committed
85
86
			     SubAssembler*,
			     Quadrature *quad = NULL) 
87
    {}
88

Thomas Witkowski's avatar
Thomas Witkowski committed
89
90
91
92
93
94
95
96
97
98
99
100
    virtual void initElement(const ElInfo* largeElInfo,
			     const ElInfo* smallElInfo,
			     SubAssembler*,
			     Quadrature *quad = NULL)
    {}

    /// Returs \auxFESpaces, the list of all aux fe spaces the operator makes use off.
    std::vector<const FiniteElemSpace*>& getAuxFESpaces() {
      return auxFESpaces;
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
103
    /// Returns true, if the term is piecewise constant, returns false otherwise
104
105
    inline bool isPWConst() { 
      return (degree == 0);
106
    }
107

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

Thomas Witkowski's avatar
Thomas Witkowski committed
111
    /// Returns \ref degree.
112
113
    inline int getDegree() { 
      return degree; 
114
    }
115

Thomas Witkowski's avatar
Thomas Witkowski committed
116
    /// Evaluation of the OperatorTerm at all quadrature points.
Thomas Witkowski's avatar
Thomas Witkowski committed
117
118
    virtual void eval(int nPoints,
		      const double *uhAtQP,
119
120
		      const WorldVector<double> *grdUhAtQP,
		      const WorldMatrix<double> *D2UhAtQP,
Thomas Witkowski's avatar
Thomas Witkowski committed
121
		      double *result,
122
123
124
		      double factor) const = 0;

    /** \brief
Thomas Witkowski's avatar
Thomas Witkowski committed
125
126
     * Determines the value of a dof vector at the quadrature points of a given 
     * element. It is used by all VecAtQP like operator terms.
127
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
    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);	

    WorldVector<double>* getGradientsAtQPs(DOFVectorBase<double>* vec,
					   const ElInfo* elInfo,
					   SubAssembler* subAssembler,
					   Quadrature *quad);

    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$.
158
    static void lalt(const DimVec<WorldVector<double> >& Lambda,
Thomas Witkowski's avatar
Thomas Witkowski committed
159
160
161
162
		     const WorldMatrix<double>& matrix,
		     DimMat<double>& LALt,
		     bool symm,
		     double factor);
163
164
165
166
167
168
169
170

    /** \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
171
172
			DimMat<double>& LALt,
			double factor);
173
174
175
176
177
178

    /** \brief
     * Evaluation of \f$ \Lambda \cdot A \cdot \Lambda^t\f$ for A equal to the 
     * identity.
     */
    static void l1lt(const DimVec<WorldVector<double> >& Lambda,
Thomas Witkowski's avatar
Thomas Witkowski committed
179
180
		     DimMat<double>& LALt,
		     double factor);
181

Thomas Witkowski's avatar
Thomas Witkowski committed
182
    /// Evaluation of \f$ \Lambda \cdot b\f$.
183
    static void lb(const DimVec<WorldVector<double> >& Lambda,
Thomas Witkowski's avatar
Thomas Witkowski committed
184
185
186
		   const WorldVector<double>& b,
		   DimVec<double>& Lb,
		   double factor);
187
188
189
190
191
192

    /** \brief
     * Evaluation of \f$ \Lambda \cdot b\f$ if b contains the value 1.0 in
     * each component.
     */
    static void l1(const DimVec<WorldVector<double> >& Lambda,
Thomas Witkowski's avatar
Thomas Witkowski committed
193
194
		   DimVec<double>& Lb,
		   double factor)
195
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
196
197
      int dim = Lb.getSize() - 1;
      static const int dimOfWorld = Global::getGeo(WORLD);
198
199
200
201
202
203
204
205
206
207

      for (int i = 0; i <= dim; i++) {
	double val = 0.0;
      
	for (int j = 0; j < dimOfWorld; j++) {
	  val += Lambda[i][j];
	}
	val *= factor;
	Lb[i] += val;
      }    
208
    }
209
210

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
211
    /// Stores the properties of this OperatorTerm
212
213
    Flag properties;

Thomas Witkowski's avatar
Thomas Witkowski committed
214
    /// Polynomial degree of the term. Used to detemine the degree of the quadrature.
215
216
    int degree;

Thomas Witkowski's avatar
Thomas Witkowski committed
217
218
219
220
    /// List off all fe spaces, the operator term makes use off.
    std::vector<const FiniteElemSpace*> auxFESpaces;

    /// Pointer to the Operator this OperatorTerm belongs to.
221
222
    Operator* operat;

Thomas Witkowski's avatar
Thomas Witkowski committed
223
    /// Flag for piecewise constant terms
224
225
    static const Flag PW_CONST;

Thomas Witkowski's avatar
Thomas Witkowski committed
226
    /// Flag for symmetric terms
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
    static const Flag SYMMETRIC;

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

  // ============================================================================
  // =====  class SecondOrderTerm ===============================================
  // ============================================================================

  /**
   * \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
249
    /// Constructor.
250
251
252
    SecondOrderTerm(int deg) 
      : OperatorTerm(deg) 
    {}
253
254
255
256

    /** \brief
     * Destructor.
     */
257
258
    virtual ~SecondOrderTerm() 
    {}
259
260
261
262
263

    /** \brief
     * Evaluation of \f$ \Lambda A \Lambda^t \f$ at all quadrature points.
     */
    virtual void getLALt(const ElInfo *elInfo, 
Thomas Witkowski's avatar
Thomas Witkowski committed
264
			 int nPoints, 
265
266
267
268
269
			 DimMat<double> **result) const = 0;

    /** \brief
     * Evaluation of \f$ A \nabla u(\vec{x}) \f$ at all quadrature points.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
270
    virtual void weakEval(int nPoints,
271
272
273
274
275
276
277
278
279
280
281
282
283
284
			  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
285
    /// Constructor.
286
287
288
289
    Laplace_SOT() 
      : SecondOrderTerm(0) 
    {
      setSymmetric(true);
290
    }
291
292
293
294

    /** \brief
     * Implenetation of SecondOrderTerm::getLALt().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
295
    inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
296
297
298
    {
      const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();

Thomas Witkowski's avatar
Thomas Witkowski committed
299
      for (int iq = 0; iq < nPoints; iq++) {
300
301
	l1lt(Lambda, *(LALt[iq]), 1.0);
      }
302
    }
303
304
305
306

    /** \brief
     * Implementation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
307
308
309
310
311
312
    inline void eval(int nPoints,
		     const double * ,    // uhAtQP
		     const WorldVector<double> * ,  // grdUhAtQP
		     const WorldMatrix<double> *D2UhAtQP,
		     double *result,
		     double factor) const
313
314
315
316
    {
      int dow = Global::getGeo(WORLD);
    
      if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
317
	for (int iq = 0; iq < nPoints; iq++) {
318
319
320
321
322
323
324
	  double resultQP = 0.0;
	  for (int i = 0; i < dow; i++) {
	    resultQP += D2UhAtQP[iq][i][i];
	  }
	  result[iq] += factor * resultQP;
	}
      }
325
    }
326
327
328
329

    /** \brief
     * Implenetation of SecondOrderTerm::weakEval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
330
    void weakEval(int nPoints,
331
332
333
334
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const
    {
      if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
335
	for (int iq = 0; iq < nPoints; iq++) { 
336
337
338
	  result[iq] += grdUhAtQP[iq];
	}
      }
339
340
    }
  };
341
342
343
344
345
346
347
348
349
350
351

  /**
   * \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
352
    /// Constructor.
353
    FactorLaplace_SOT(double f) 
354
355
      : SecondOrderTerm(0)  
    {
356
357
358
359
      factor = new double;
      *factor = f;

      setSymmetric(true);
360
    }
361

Thomas Witkowski's avatar
Thomas Witkowski committed
362
    /// Constructor.
363
    FactorLaplace_SOT(double *fptr) 
364
365
      : SecondOrderTerm(0), 
	factor(fptr)
366
367
    {
      setSymmetric(true);
368
    }
369
370
371
372

    /** \brief
     * Implements SecondOrderTerm::getLALt().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
373
    inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
374
    {
375
      const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
376
      for (int iq = 0; iq < nPoints; iq++) 
377
	l1lt(Lambda, *(LALt[iq]), (*factor));
378
    }
379
380
381
382
383


    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
384
    void eval(int nPoints,
385
	      const double *,
386
387
388
389
390
	      const WorldVector<double> *,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double f) const
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
391
      int dow = Global::getGeo(WORLD);
392

Thomas Witkowski's avatar
Thomas Witkowski committed
393
      if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
394
	for (int iq = 0; iq < nPoints; iq++) {
395
	  double resultQP = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
396
	  for (int i = 0; i < dow; i++) {
397
398
399
400
401
	    resultQP += D2UhAtQP[iq][i][i];
	  }
	  result[iq] += resultQP * f * (*factor);
	}
      }
402
    }
403
404
405
406

    /** \brief
     * Implenetation of SecondOrderTerm::weakEval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
407
    void weakEval(int nPoints,
408
409
410
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
411
      if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
412
	for (int iq = 0; iq < nPoints; iq++) {
413
414
415
	  axpy(*factor, grdUhAtQP[iq], result[iq]);
	}
      }
416
    }
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435

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

    /** \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
451
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
452
453
454
455
  
    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
456
457
    void eval(int nPoints,
	      const double *uhAtQP,
458
459
460
461
462
463
464
465
	      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
466
    void weakEval(int nPoints,
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;

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

    /** \brief
     * Pointer to the values of the DOFVector at quadrature points.
     */
    double* vecAtQPs;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
486
487
488
    /** \brief
     *
     */
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
    AbstractFunction<WorldVector<double>, WorldMatrix<double> >* divFct;

    /** \brief
     * True, if \ref matrixFct produces always symmetric matrices.
     */
    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
506
    /// Constructor
507
    Matrix_SOT(WorldMatrix<double> mat) 
Thomas Witkowski's avatar
Thomas Witkowski committed
508
      : SecondOrderTerm(0), matrix(mat)
509
510
511
    {
      symmetric = matrix.isSymmetric();
      setSymmetric(symmetric);
512
    }
513
514
515
516

    /** \brief
     * Implements SecondOrderTerm::getLALt().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
517
    inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const{
518
519
520
      const DimVec<WorldVector<double> >& Lambda     = elInfo->getGrdLambda();
      //double det = elInfo->getDet();
      int iq;
Thomas Witkowski's avatar
Thomas Witkowski committed
521
      for(iq = 0; iq < nPoints; iq++) 
522
	lalt(Lambda, matrix, *(LALt[iq]), symmetric, 1.0);
523
    }
524
525
526
527
  
    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
528
529
    void eval(int nPoints,
	      const double *uhAtQP,
530
531
532
533
534
535
536
537
	      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
538
    void weakEval(int nPoints,
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;

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

    /** \brief
     * True, if \ref matrix is symmetric.
     */
    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
565
    /// Constructor.
566
    FactorIJ_SOT(int x_i, int x_j, double f) 
Thomas Witkowski's avatar
Thomas Witkowski committed
567
      : SecondOrderTerm(0), xi(x_i), xj(x_j)
568
569
570
571
572
    {
      factor = new double;
      *factor = f;

      setSymmetric(xi == xj);
573
    }
574

Thomas Witkowski's avatar
Thomas Witkowski committed
575
    /// Constructor.
576
    FactorIJ_SOT(int x_i, int x_j, double *fptr) 
Thomas Witkowski's avatar
Thomas Witkowski committed
577
      : SecondOrderTerm(0), xi(x_i), xj(x_j), factor(fptr)
578
579
    {
      setSymmetric(xi == xj);
580
    }
581
582
583
584

    /** \brief
     * Implements SecondOrderTerm::getLALt().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
585
    inline void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
586
587
588
    {
      const DimVec<WorldVector<double> > &Lambda     = elInfo->getGrdLambda();
      int iq;
Thomas Witkowski's avatar
Thomas Witkowski committed
589
      for(iq = 0; iq < nPoints; iq++)
590
	lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*factor));
591
    }
592
593
594
595

    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
596
    void eval(int nPoints,
597
598
599
600
601
602
603
604
	      const double              *,
	      const WorldVector<double> *,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const
    {
      int iq;
      if(D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
605
	for(iq = 0; iq < nPoints; iq++) {
606
607
608
	  result[iq] += (*factor) * D2UhAtQP[iq][xi][xj] * fac;
	}
      }
609
    }
610
611
612
613

    /** \brief
     * Implenetation of SecondOrderTerm::weakEval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
614
    void weakEval(int nPoints,
615
616
617
618
619
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const
    {
      if(grdUhAtQP) {
	int iq;
Thomas Witkowski's avatar
Thomas Witkowski committed
620
	for(iq = 0; iq < nPoints; iq++) {
621
622
623
	  result[iq][xi] += (*factor) * grdUhAtQP[iq][xj];
	}
      }
624
    }
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647

  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
648
649
    /// Constructor.
    VecAtQP_SOT(DOFVectorBase<double> *dv, AbstractFunction<double, double> *af);
650

Thomas Witkowski's avatar
Thomas Witkowski committed
651
    /// Implementation of \ref OperatorTerm::initElement().
652
653
654
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

Thomas Witkowski's avatar
Thomas Witkowski committed
655
656
657
658
659
    /// Implementation of \ref OperatorTerm::initElement() for multilpe meshes.
    void initElement(const ElInfo* smallElInfo,
		     const ElInfo* largeElInfo,
		     SubAssembler* subAssembler,
		     Quadrature *quad = NULL);
660

Thomas Witkowski's avatar
Thomas Witkowski committed
661
662
663
664
665
666
    /// 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,
667
668
669
670
671
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor) const;

Thomas Witkowski's avatar
Thomas Witkowski committed
672
673
    /// Implenetation of SecondOrderTerm::weakEval().
    void weakEval(int nPoints,
674
675
676
677
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
681
    /// Pointer to an array containing the DOFVector evaluated at quadrature points.
682
683
    double* vecAtQPs;

Thomas Witkowski's avatar
Thomas Witkowski committed
684
    /// Function evaluated at \ref vecAtQPs.
685
686
687
    AbstractFunction<double, double> *f;
  };

688
689
690
691
692
693
694
695
696
697
  /**
   * \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
698
    /// Constructor.
699
    Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, 
Thomas Witkowski's avatar
Thomas Witkowski committed
700
		 BinaryAbstractFunction<double, double, double> *af);
701
702
703
704
705
706
707
708
709
710

    /** \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
711
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;    
712
713
714
715
    
    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
716
    void eval(int nPoints,
717
718
719
720
721
722
723
724
725
	      const double *uhAtQP,
	      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
726
    void weakEval(int nPoints,
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const;
    
  protected:
    /** \brief
     * DOFVector to be evaluated at quadrature points.
     */
    DOFVectorBase<double>* vec1;
    DOFVectorBase<double>* vec2;
    
    /** \brief
     * Pointer to an array containing the DOFVector evaluated at quadrature
     * points.
     */
    double* vecAtQPs1;
    double* vecAtQPs2;
    
    /** \brief
     * Function evaluated at \ref vecAtQPs.
     */
    BinaryAbstractFunction<double, double, double> *f;
  };

750
751
752
753
754
755
756
757
758
  /**
   * \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
759
760
761
    /// Constructor.
    CoordsAtQP_SOT(AbstractFunction<double, WorldVector<double> > *af)
      : SecondOrderTerm(af->getDegree()), g(af)
762
763
    {
      setSymmetric(true);
Thomas Witkowski's avatar
Thomas Witkowski committed
764
    }
765
766
767
768
769
770
771
772
773
774

    /** \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
775
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;    
776
777
778
779

    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
780
781
    void eval(int nPoints,
	      const double *uhAtQP,
782
783
784
785
786
787
788
789
	      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
790
    void weakEval(int nPoints,
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
		  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
818
    /// Constructor.
819
    MatrixGradient_SOT(DOFVectorBase<double> *dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
820
821
822
		       AbstractFunction<WorldMatrix<double>, WorldVector<double> > *af,
		       AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
		       bool symm = false);
823
824
825
826
827
828
829
830
831
832

    /** \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
833
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
834
835
836
837

    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
838
839
    void eval(int nPoints,
	      const double *uhAtQP,
840
841
842
843
844
845
846
847
	      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
848
    void weakEval(int nPoints,
849
850
851
852
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
883
884
		  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
885
886
887
    /// Constructor.
    FctGradient_SOT(DOFVectorBase<double> *dv,  
		    AbstractFunction<double, WorldVector<double> > *af);
888
889
890
891
892
893
894
895
896
897

    /** \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
898
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;
899
900
901
902

    /** \brief
     * Implenetation of SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
903
904
    void eval(int nPoints,
	      const double *uhAtQP,
905
906
907
908
909
910
911
912
	      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
913
    void weakEval(int nPoints,
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
		  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
944
    /// Constructor.
945
    VecAndGradientAtQP_SOT(DOFVectorBase<double> *dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
946
			   BinaryAbstractFunction<double, double, WorldVector<double> > *af);
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
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
		  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;
  };

992
993
994
995
996
997
998
999
1000
  /**
   * \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