Operator.h 97.7 KB
Newer Older
3001
3002
3003
3004
3005
3006
3007
3008
3009
3010
3011
3012
3013
3014
3015
3016
3017
3018
3019
3020
3021
3022
3023
3024
3025
3026
3027
3028
3029
3030
3031
3032
3033
3034
3035
3036
3037
3038
3039
3040
3041
3042
3043
3044
3045
3046
3047
3048
3049
3050
3051
3052
3053
3054
3055
3056
3057
3058
3059
3060
3061
3062
3063
3064
3065
3066
3067
3068
3069
3070
3071
3072
3073
3074
3075
3076
3077
3078
3079
3080
3081
3082
3083
3084
3085
3086
3087
    /** \brief 
     * Function for c. 
     */ 
    BinaryAbstractFunction<double, double, WorldVector<double> > *f; 
  }; 

  // ============================================================================

  class VecAndGradVec2AtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
    /** \brief 
     * Constructor. 
     */ 
    VecAndGradVec2AtQP_ZOT(DOFVectorBase<double> *dv, 
			   DOFVectorBase<double> *dGrd1, DOFVectorBase<double> *dGrd2, 
			   TertiaryAbstractFunction<double, 
			   double, WorldVector<double>, WorldVector<double> > *f_) 
      : ZeroOrderTerm(f_->getDegree()), vec(dv), vecGrd1(dGrd1), vecGrd2(dGrd2), f(f_) 
    { 
    }; 

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
    void getC(const ElInfo *, int numPoints, double *C) const;

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
    void eval(int numPoints,
	      const double              *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

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

    /** \brief 
     * Vector v at quadrature points. 
     */ 
    double *vecAtQPs; 

    /** \brief 
     * First DOFVector whose gradient is evaluated at quadrature points. 
     */ 
    DOFVectorBase<double>* vecGrd1; 

    /** \brief 
     * Gradient of first vector at quadrature points. 
     */ 
    WorldVector<double> *grad1AtQPs; 

    /** \brief 
     * Second DOFVector whose gradient is evaluated at quadrature points. 
     */ 
    DOFVectorBase<double>* vecGrd2; 

    /** \brief 
     * Gradient of second vector at quadrature points. 
     */ 
    WorldVector<double> *grad2AtQPs; 

    /** \brief 
     * Function for c. 
     */ 
    TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *f; 
  }; 


  class VecOfDOFVecsAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
    /** \brief 
     * Constructor. 
     */ 
3088
3089
    VecOfDOFVecsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
			 AbstractFunction<double, std::vector<double> > *f_) 
3090
3091
3092
3093
3094
3095
3096
3097
3098
3099
3100
3101
3102
3103
3104
3105
3106
3107
3108
3109
3110
3111
3112
3113
3114
3115
3116
3117
3118
3119
      : ZeroOrderTerm(f_->getDegree()), vecs(dv), f(f_) 
    {
      vecsAtQPs.resize(vecs.size());
    }; 

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
    void getC(const ElInfo *, int numPoints, double *C) const;

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
    void eval(int numPoints,
	      const double              *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
    /** \brief 
     * Vector of DOFVectors to be evaluated at quadrature points. 
     */ 
3120
    std::vector<DOFVectorBase<double>*> vecs; 
3121
3122
3123
3124

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3125
    std::vector<double*> vecsAtQPs; 
3126
3127
3128
3129

    /** \brief 
     * Function for c. 
     */ 
3130
    AbstractFunction<double, std::vector<double> > *f; 
3131
3132
3133
3134
3135
3136
3137
3138
  }; 

  class VecOfGradientsAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
    /** \brief 
     * Constructor. 
     */ 
3139
3140
    VecOfGradientsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
			   AbstractFunction<double, std::vector<WorldVector<double>*> > *f_) 
3141
3142
3143
3144
3145
3146
3147
3148
3149
3150
3151
3152
3153
3154
3155
3156
3157
3158
3159
3160
3161
3162
3163
3164
3165
3166
3167
3168
3169
3170
      : ZeroOrderTerm(f_->getDegree()), vecs(dv), f(f_) 
    {
      gradsAtQPs.resize(vecs.size());
    }; 

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
    void getC(const ElInfo *, int numPoints, double *C) const;

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
    void eval(int numPoints,
	      const double              *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
    /** \brief 
     * Vector of DOFVectors to be evaluated at quadrature points. 
     */ 
3171
    std::vector<DOFVectorBase<double>*> vecs; 
3172
3173
3174
3175

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3176
    std::vector<WorldVector<double>*> gradsAtQPs; 
3177
3178
3179
3180

    /** \brief 
     * Function for c. 
     */ 
3181
    AbstractFunction<double, std::vector<WorldVector<double>*> > *f; 
3182
3183
3184
3185
3186
3187
3188
3189
3190
3191
3192
3193
3194
3195
3196
3197
3198
3199
3200
3201
3202
3203
3204
3205
3206
3207
3208
3209
3210
3211
3212
3213
3214
3215
3216
3217
3218
3219
3220
3221
3222
3223
3224
3225
3226
3227
3228
  };


  class VecDivergence_ZOT : public ZeroOrderTerm
  {
  public: 
    /** \brief 
     * Constructor. 
     */ 
    VecDivergence_ZOT(int numComponents,
		      DOFVectorBase<double> *vec0,
		      DOFVectorBase<double> *vec1 = NULL,
		      DOFVectorBase<double> *vec2 = NULL)
      : ZeroOrderTerm(0) 
    {
      vecs.resize(numComponents);
      gradsAtQPs.resize(numComponents);
      vecs[0] = vec0;
      vecs[1] = vec1;
      vecs[2] = vec2;
    }; 

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
    void getC(const ElInfo *, int numPoints, double *C) const;

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
    void eval(int numPoints,
	      const double              *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
    /** \brief 
     * Vector of DOFVectors to be evaluated at quadrature points. 
     */ 
3229
    std::vector<DOFVectorBase<double>*> vecs; 
3230
3231
3232
3233

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3234
    std::vector<WorldVector<double>*> gradsAtQPs; 
3235
3236
3237
3238
3239
3240
3241
3242
3243
3244
3245
3246
3247
3248
3249
3250
3251
3252
3253
3254
3255
3256
3257
3258
3259
3260
3261
3262
3263
3264
3265
3266
3267
3268
3269
3270
3271
3272
3273
3274
3275
3276
3277
3278
3279
3280
3281
3282
3283
3284
3285
3286
3287
3288
3289
3290
3291
3292
3293
3294
3295
3296
3297
3298
3299
3300
3301
3302
3303
3304
3305
3306
3307
3308
3309
3310
  };



  class VecAndVecOfGradientsAtQP_ZOT : public ZeroOrderTerm
  {
  public:
    /** \brief
     * Constructor.
     */
    VecAndVecOfGradientsAtQP_ZOT(DOFVector<double> *vec_,
				 const std::vector<DOFVector<double>*>& dv,
				 BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *f_)
      : ZeroOrderTerm(f_->getDegree()), vec(vec_), vecs(dv), f(f_)
    {
      gradsAtQPs.resize(vecs.size());
    };

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
    void getC(const ElInfo *, int numPoints, double *C) const;

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
    void eval(int numPoints,
	      const double              *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

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

    /** \brief
     * Vector v at quadrature points.
     */
    double *vecAtQPs;

    /** \brief
     * Vector of DOFVectors to be evaluated at quadrature points.
     */
    std::vector<DOFVector<double>*> vecs;

    /** \brief
     * Vectors at quadrature points.
     */
    std::vector<WorldVector<double>*> gradsAtQPs;

    /** \brief
     * Function for c.
     */
    BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *f;
  };



  // ============================================================================

  class General_ZOT : public ZeroOrderTerm
  {
  public:
    /** \brief
     * Constructor.
     */
3311
3312
    General_ZOT(std::vector<DOFVectorBase<double>*> vecs,
		std::vector<DOFVectorBase<double>*> grads,
3313
3314
		TertiaryAbstractFunction<double, 
		WorldVector<double>,
3315
3316
		std::vector<double>, 
		std::vector<WorldVector<double> > > *f)
3317
3318
3319
3320
3321
3322
3323
3324
3325
3326
3327
3328
3329
3330
3331
3332
3333
3334
3335
3336
3337
3338
3339
3340
3341
3342
3343
3344
3345
3346
3347
      : ZeroOrderTerm(f->getDegree()),
	vecs_(vecs),
	grads_(grads),
	f_(f)
    {
      vecsAtQPs_.resize(vecs_.size());
      gradsAtQPs_.resize(grads_.size());
    };

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
    void getC(const ElInfo *, int numPoints, double *C) const;

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
    void eval(int numPoints,
	      const double              *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
3348
    std::vector<DOFVectorBase<double>*> vecs_; 
3349

3350
    std::vector<DOFVectorBase<double>*> grads_;
3351
3352
3353

    TertiaryAbstractFunction<double, 
			     WorldVector<double>,
3354
3355
			     std::vector<double>, 
			     std::vector<WorldVector<double> > > *f_;
3356
3357
3358

    WorldVector<double> *coordsAtQPs_;

3359
    std::vector<double*> vecsAtQPs_;
3360

3361
    std::vector<WorldVector<double>*> gradsAtQPs_;
3362
3363
3364
3365
3366
3367
3368
3369
3370
3371
3372
  };


  // ============================================================================

  class GeneralParametric_ZOT : public ZeroOrderTerm
  {
  public:
    /** \brief
     * Constructor.
     */
3373
3374
    GeneralParametric_ZOT(std::vector<DOFVectorBase<double>*> vecs,
			  std::vector<DOFVectorBase<double>*> grads,
3375
3376
3377
			  QuartAbstractFunction<double, 
			  WorldVector<double>,
			  WorldVector<double>,
3378
3379
			  std::vector<double>, 
			  std::vector<WorldVector<double> > > *f)
3380
3381
3382
3383
3384
3385
3386
3387
3388
3389
3390
3391
3392
3393
3394
3395
3396
3397
3398
3399
3400
3401
3402
3403
3404
3405
3406
3407
3408
3409
3410
      : ZeroOrderTerm(f->getDegree()),
	vecs_(vecs),
	grads_(grads),
	f_(f)
    {
      vecsAtQPs_.resize(vecs_.size());
      gradsAtQPs_.resize(grads_.size());
    };

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
    void getC(const ElInfo *, int numPoints, double *C) const;

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
    void eval(int numPoints,
	      const double              *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
3411
    std::vector<DOFVectorBase<double>*> vecs_; 
3412

3413
    std::vector<DOFVectorBase<double>*> grads_;
3414
3415
3416
3417

    QuartAbstractFunction<double, 
			  WorldVector<double>,
			  WorldVector<double>,
3418
3419
			  std::vector<double>, 
			  std::vector<WorldVector<double> > > *f_;
3420
3421
3422
3423
3424

    WorldVector<double> *coordsAtQPs_;

    WorldVector<double> elementNormal_;

3425
    std::vector<double*> vecsAtQPs_;
3426

3427
    std::vector<WorldVector<double>*> gradsAtQPs_;
3428
3429
3430
3431
3432
3433
3434
3435
3436
3437
3438
3439
3440
3441
3442
3443
3444
3445
3446
3447
3448
3449
3450
3451
3452
3453
3454
3455
3456
3457
3458
3459
3460
3461
3462
3463
3464
3465
3466
3467
3468
3469
3470
  };



  /*****************************************************************************/
  /******      Operators for the least-square finite element method      *******/
  /*****************************************************************************/


  // ============================================================================
  // ===== class Operator =======================================================
  // ============================================================================

  /** \brief
   * An Operator holds all information needed to assemble the system matrix
   * and the right hand side. It consists of four OperatorTerm lists each storing
   * Terms of a specific order and type. You can define your own Operator by 
   * creating an empty Operator and than adding OperatorTerms to it.
   * An Operator can by of type MATRIX_OPERATOR, if it's used to assemble the
   * system matrix on the left hand side, or it can be of type VECTOR_OPERATOR,
   * if it's used to assemble the right hand side vector. If an Operator gives
   * contributions to both sides of the system it is a MATRIX_OPERATOR and a
   * VECTOR_OPERATOR in one instance. This allows to efficiently reuse element 
   * matrices once calculated.
   * By calling \ref getElementMatrix() or \ref getElementVector() one can 
   * initiate the assembling procedure. Therefor each Operator has its own
   * Assembler, especially created for this Operator, by the first call of
   * \ref getElementMatrix() or \ref getElementVector(). 
   */
  class Operator
  {
  public:
    MEMORY_MANAGED(Operator);

    /** \brief
     * Constructs an empty Operator of type operatorType for the given 
     * FiniteElemSpace. 
     * The type is given by a Flag that can contain the values MATRIX_OPERATOR,
     * VECTOR_OPERATOR, or MATRIX_OPERATOR | VECTOR_OPERATOR. This type specifies 
     * whether the Operator is used on the left hand side, the right hand side,
     * or on both sides of the system. 
     */
    Operator(Flag operatorType,
Thomas Witkowski's avatar
Thomas Witkowski committed
3471
3472
	     const FiniteElemSpace *rowFESpace,
	     const FiniteElemSpace *colFESpace = NULL);
3473
3474
3475
3476
3477
3478
3479
3480
3481
3482
3483
3484
3485
3486
3487
3488
3489
3490
3491
3492
3493
3494
3495

    /** \brief
     * Destructor.
     */
    virtual ~Operator() {};

    /** \brief
     * Sets \ref optimized.
     */
    inline void useOptimizedAssembler(bool opt) {
      optimized = opt;
    };

    /** \brief
     * Returns \ref optimized.
     */
    inline bool isOptimized() {
      return optimized;
    };

    /** \brief
     * Adds a SecondOrderTerm to the Operator
     */
3496
3497
    template <typename T>
    void addSecondOrderTerm(T *term);
3498
3499
3500
3501

    /** \brief
     * Adds a FirstOrderTerm to the Operator
     */
3502
3503
    template <typename T>
    void addFirstOrderTerm(T *term, 
3504
3505
3506
3507
			   FirstOrderType type = GRD_PHI);
    /** \brief
     * Adds a ZeroOrderTerm to the Operator
     */
3508
3509
    template <typename T>
    void addZeroOrderTerm(T *term);
3510
3511
3512
3513
3514
3515
3516
3517
3518
3519
3520
3521
3522
3523
3524
3525
3526
3527
3528
3529
3530
3531
3532
3533
3534
3535
3536
3537
3538
3539
3540
3541
3542
3543
3544
3545
3546
3547
3548
3549
3550
3551
3552
3553
3554
3555
3556
3557
3558


    /** \brief
     * Calculates the element matrix for this ElInfo and adds it multiplied by
     * factor to userMat.
     */
    virtual void getElementMatrix(const ElInfo *elInfo, 
				  ElementMatrix *userMat, 
				  double factor = 1.0);

    /** \brief
     * Calculates the element vector for this ElInfo and adds it multiplied by
     * factor to userVec.
     */
    virtual void getElementVector(const ElInfo *elInfo, 
				  ElementVector *userVec, 
				  double factor = 1.0);



    /** \brief
     * Returns \ref rowFESpace
     */
    inline const FiniteElemSpace *getRowFESpace() { 
      return rowFESpace; 
    };

    /** \brief
     * Returns \ref colFESpace
     */
    inline const FiniteElemSpace *getColFESpace() { 
      return colFESpace; 
    };

    /** \brief
     * Sets \ref uhOld.
     */
    void setUhOld(const DOFVectorBase<double> *uhOld);

    /** \brief
     * Returns \ref uhOld.
     */
    inline const DOFVectorBase<double> *getUhOld() {
      return uhOld;
    };

    /** \brief
     * Returns \ref assembler
     */
3559
    Assembler *getAssembler(int rank);
3560
3561
3562
3563

    /** \brief
     * Sets \ref assembler
     */
3564
    void setAssembler(int rank, Assembler *ass);
3565
3566
3567
3568
3569
3570
3571
3572
3573
3574
3575
3576
3577
3578
3579
3580
3581
3582
3583
3584
3585
3586
3587
3588
3589
3590
3591
3592
3593
3594
3595
3596

    /** \brief
     * Returns whether this is a matrix operator.
     */
    inline const bool isMatrixOperator() {
      return type.isSet(MATRIX_OPERATOR);
    };

    /** \brief
     * Returns whether this is a vector operator
     */
    inline const bool isVectorOperator() {
      return type.isSet(VECTOR_OPERATOR);
    };

    /** \brief
     * Sets \ref fillFlag, the flag used for this Operator while mesh traversal.
     */
    inline void setFillFlag(Flag f) { 
      fillFlag = f; 
    };

    /** \brief
     * Returns \ref fillFlag
     */
    inline Flag getFillFlag() { 
      return fillFlag; 
    };

    /** \brief
     * Initialization of the needed SubAssemblers using the given quadratures. 
     */
3597
3598
    void initAssembler(int rank,
		       Quadrature *quad2,
3599
3600
3601
3602
3603
3604
3605
3606
3607
3608
3609
3610
3611
3612
3613
3614
3615
3616
3617
3618
		       Quadrature *quad1GrdPsi,
		       Quadrature *quad1GrdPhi,
		       Quadrature *quad0);


    /** \brief
     * Calculates the needed quadrature degree for the given order. 
     */
    int getQuadratureDegree(int order, FirstOrderType firstOrderType = GRD_PHI);

    /** \brief
     * Evaluation of all terms in \ref zeroOrder. 
     */
    void evalZeroOrder(int numPoints,
		       const double              *uhAtQP,
		       const WorldVector<double> *grdUhAtQP,
		       const WorldMatrix<double> *D2UhAtQP,
		       double *result,
		       double factor) const
    {
3619
3620
      int myRank = omp_get_thread_num();

3621
      std::vector<OperatorTerm*>::const_iterator termIt;
3622
3623
3624
3625
3626
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
	   ++termIt) {
	(*termIt)->eval(numPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
      }
3627
3628
3629
3630
3631
3632
3633
3634
3635
3636
3637
3638
3639
    };


    /** \brief
     * Evaluation of all terms in \ref firstOrderGrdPsi. 
     */
    void evalFirstOrderGrdPsi(int numPoints,
			      const double              *uhAtQP,
			      const WorldVector<double> *grdUhAtQP,
			      const WorldMatrix<double> *D2UhAtQP,
			      double *result,
			      double factor) const
    {
3640
3641
      int myRank = omp_get_thread_num();

3642
      std::vector<OperatorTerm*>::const_iterator termIt;
3643
3644
3645
3646
3647
      for (termIt = firstOrderGrdPsi[myRank].begin(); 
	   termIt != firstOrderGrdPsi[myRank].end(); 
	   ++termIt) {
	(*termIt)->eval(numPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
      }
3648
3649
3650
3651
3652
3653
3654
3655
3656
3657
3658
3659
    };

    /** \brief
     * Evaluation of all terms in \ref firstOrderGrdPhi. 
     */
    void evalFirstOrderGrdPhi(int numPoints,
			      const double              *uhAtQP,
			      const WorldVector<double> *grdUhAtQP,
			      const WorldMatrix<double> *D2UhAtQP,
			      double *result,
			      double factor) const
    {
3660
3661
      int myRank = omp_get_thread_num();

3662
      std::vector<OperatorTerm*>::const_iterator termIt;
3663
3664
3665
3666
3667
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
	   ++termIt) {
	(*termIt)->eval(numPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
      }
3668
3669
3670
3671
3672
3673
3674
3675
3676
3677
3678
3679
3680
    };


    /** \brief
     * Evaluation of all terms in \ref secondOrder. 
     */
    void evalSecondOrder(int numPoints,
			 const double              *uhAtQP,
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
    {
3681
3682
      int myRank = omp_get_thread_num();

3683
      std::vector<OperatorTerm*>::const_iterator termIt;
3684
3685
3686
3687
3688
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
	   ++termIt) {
	(*termIt)->eval(numPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
      }
3689
3690
3691
3692
3693
3694
3695
3696
3697
    };

    /** \brief
     * Weak evaluation of all terms in \ref secondOrder. 
     */
    void weakEvalSecondOrder(int numPoints,
			     const WorldVector<double> *grdUhAtQP,
			     WorldVector<double> *result) const
    {
3698
3699
      int myRank = omp_get_thread_num();

3700
      std::vector<OperatorTerm*>::const_iterator termIt;
3701
3702
3703
3704
3705
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
	   ++termIt) {
	static_cast<SecondOrderTerm*>(*termIt)->weakEval(numPoints, grdUhAtQP, result);
      }
3706
3707
3708
3709
3710
3711
3712
3713
    };
  
    /** \brief
     * Calls getLALt() for each term in \ref secondOrder 
     * and adds the results to LALt.
     */
    void getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const
    {
3714
3715
      int myRank = omp_get_thread_num();

3716
      std::vector<OperatorTerm*>::const_iterator termIt;
3717
3718
3719
3720
3721
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
	   ++termIt) {
	static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, numPoints, LALt);
      }
3722
3723
3724
3725
3726
3727
3728
3729
    };
  
    /** \brief
     * Calls getLb() for each term in \ref firstOrderGrdPsi 
     * and adds the results to Lb.
     */
    void getLbGrdPsi(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const
    {
3730
3731
      int myRank = omp_get_thread_num();

3732
      std::vector<OperatorTerm*>::const_iterator termIt;
3733
3734
3735
3736
3737
      for (termIt = firstOrderGrdPsi[myRank].begin(); 
	   termIt != firstOrderGrdPsi[myRank].end(); 
	   ++termIt) {
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, numPoints, Lb);
      }
3738
3739
3740
3741
3742
3743
3744
3745
    };

    /** \brief
     * Calls getLb() for each term in \ref firstOrderGrdPhi 
     * and adds the results to Lb.
     */
    void getLbGrdPhi(const ElInfo *elInfo, int numPoints, VectorOfFixVecs<DimVec<double> >& Lb) const
    {
3746
3747
      int myRank = omp_get_thread_num();

3748
      std::vector<OperatorTerm*>::const_iterator termIt;
3749
3750
3751
3752
3753
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
	   ++termIt) {
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, numPoints, Lb);
      }
3754
3755
3756
3757
3758
3759
3760
3761
    };

    /** \brief
     * Calls getC() for each term in \ref zeroOrder
     * and adds the results to c.
     */
    void getC(const ElInfo *elInfo, int numPoints, double *c) const
    {
3762
3763
      int myRank = omp_get_thread_num();

3764
      std::vector<OperatorTerm*>::const_iterator termIt;
3765
3766
3767
3768
3769
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
	   ++termIt) {
	static_cast<ZeroOrderTerm*>(*termIt)->getC(elInfo, numPoints, c);
      }
3770
3771
3772
3773
3774
3775
    };

    /** \brief
     * Returns true, if there are second order terms. Returns false otherwise.
     */
    inline bool secondOrderTerms() {
3776
      return secondOrder[omp_get_thread_num()].size() != 0;
3777
3778
3779
3780
3781
3782
3783
    };

    /** \brief
     * Returns true, if there are first order terms (grdPsi). 
     * Returns false otherwise.
     */
    inline bool firstOrderTermsGrdPsi() {
3784
      return firstOrderGrdPsi[omp_get_thread_num()].size() != 0;
3785
3786
3787
3788
3789
3790
3791
    };

    /** \brief
     * Returns true, if there are first order terms (grdPhi). 
     * Returns false otherwise.
     */
    inline bool firstOrderTermsGrdPhi() {
3792
      return firstOrderGrdPhi[omp_get_thread_num()].size() != 0;
3793
3794
3795
3796
3797
3798
3799
    };

    /** \brief
     * Returns true, if there are zero order terms.
     * Returns false otherwise.
     */
    inline bool zeroOrderTerms() {
3800
      return zeroOrder[omp_get_thread_num()].size() != 0;
3801
3802
3803
3804
3805
3806
3807
3808
3809
3810
3811
3812
3813
3814
3815
3816
3817
3818
3819
3820
3821
3822
3823
3824
3825
3826
3827
3828
3829
3830
3831
3832
3833
3834
3835
3836
3837
3838
3839
3840
3841
3842
3843
3844
3845
3846
3847
3848
3849
    };

  public:
    /** \brief
     * Constant type flag for matrix operators
     */
    static const Flag MATRIX_OPERATOR;

    /** \brief
     * Constant type flag for vector operators
     */
    static const Flag VECTOR_OPERATOR;

  protected:
    /** \brief
     * FiniteElemSpace for matrix rows and element vector
     */
    const FiniteElemSpace *rowFESpace;

    /** \brief
     * FiniteElemSpace for matrix columns. Can be equal to \rowFESpace.
     */
    const FiniteElemSpace *colFESpace;

    /** \brief
     * Number of rows in the element matrix
     */
    int nRow;

    /** \brief
     * Number of columns in the element matrix
     */
    int nCol;

    /** \brief
     * Type of this Operator.
     */
    Flag type;

    /** \brief
     * Flag for mesh traversal
     */
    Flag fillFlag;

    /** \brief
     * Calculates the element matrix and/or the element vector. It is
     * created especially for this Operator, when \ref getElementMatrix()
     * or \ref getElementVector is called for the first time.
     */
3850
    std::vector<Assembler*> assembler;
3851
3852
3853
3854

    /** \brief
     * List of all second order terms
     */
3855
    std::vector< std::vector<OperatorTerm*> > secondOrder;
3856
3857
3858
3859

    /** \brief
     * List of all first order terms derived to psi
     */
3860
    std::vector< std::vector<OperatorTerm*> > firstOrderGrdPsi;
3861
3862
3863
3864

    /** \brief
     * List of all first order terms derived to phi
     */
3865
    std::vector< std::vector<OperatorTerm*> > firstOrderGrdPhi;
3866
3867
3868
3869

    /** \brief
     * List of all zero order terms
     */
3870
    std::vector< std::vector<OperatorTerm*> > zeroOrder;
3871
3872
3873
3874
3875
3876
3877
3878
3879
3880
3881
3882
3883
3884
3885
3886
3887
3888
3889
3890
3891
3892
3893

    /** \brief
     * Pointer to the solution of the last timestep. Can be used if the 
     * Operator is MATRIX_OPERATOR and VECTOR_OPERATOR for a more efficient
     * assemblage of the element vector when the element matrix was already
     * computed.
     */
    const DOFVectorBase<double> *uhOld;

    /** \brief
     * Spezifies whether optimized assemblers are used or not.
     */
    bool optimized;

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

}

3894
3895
#include "Operator.hh"

3896
#endif // AMDIS_OPERATOR_H
For faster browsing, not all history is shown. View entire blame