Operator.h 104 KB
Newer Older
3001
3002
3003
3004

    /** \brief
     * Implements SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3005
3006
    void eval(int nPoints,
	      const double *uhAtQP,
3007
3008
3009
3010
3011
3012
3013
3014
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor) const;

    /** \brief
     * Implements SecondOrderTerm::weakEval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3015
    void weakEval(int nPoints,
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
		  const WorldVector<double> *grdUhAtQP,
		  WorldVector<double> *result) const; 
  
  protected:
    /** \brief
     * DOFVector to be evaluated at quadrature points.
     */
    DOFVectorBase<double>* vec;

    /** \brief
     * Pointer to an array containing the DOFVector evaluated at quadrature
     * points.
     */
    double* vecAtQPs;

    /** \brief
     * Function evaluated at \ref vecAtQPs.
     */
    AbstractFunction<double, double> *f;

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

  class VecAndGradVecAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3046
    /// Constructor. 
3047
3048
    VecAndGradVecAtQP_ZOT(DOFVectorBase<double> *dv, 
			  DOFVectorBase<double> *dGrd,
Thomas Witkowski's avatar
Thomas Witkowski committed
3049
			  BinaryAbstractFunction<double, double, WorldVector<double> > *f);
3050
3051
3052
3053
3054
3055
3056
3057
3058
3059

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3060
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3061
3062
3063
3064

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3065
3066
    void eval(int nPoints,
	      const double *uhAtQP,
3067
3068
3069
3070
3071
3072
3073
3074
3075
3076
3077
3078
3079
3080
3081
3082
3083
3084
3085
3086
3087
3088
3089
3090
3091
3092
3093
3094
3095
3096
3097
3098
3099
3100
3101
	      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>* vecGrd; 

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

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

  class VecAndGradVec2AtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3102
    /// Constructor. 
3103
    VecAndGradVec2AtQP_ZOT(DOFVectorBase<double> *dv, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3104
3105
3106
			   DOFVectorBase<double> *dGrd1, 
			   DOFVectorBase<double> *dGrd2, 
			   TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *af);
3107
3108
3109
3110
3111
3112
3113
3114
3115
3116

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3117
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3118
3119
3120
3121

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3122
3123
    void eval(int nPoints,
	      const double *uhAtQP,
3124
3125
3126
3127
3128
3129
3130
3131
3132
3133
3134
3135
3136
3137
3138
3139
3140
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
	      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: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3170
    /// Constructor. 
3171
    VecOfDOFVecsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3172
			 AbstractFunction<double, std::vector<double> > *f);
3173
3174
3175
3176
3177
3178
3179
3180
3181
3182

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3183
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3184
3185
3186
3187

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3188
3189
    void eval(int nPoints,
	      const double *uhAtQP,
3190
3191
3192
3193
3194
3195
3196
3197
3198
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
    /** \brief 
     * Vector of DOFVectors to be evaluated at quadrature points. 
     */ 
3199
    std::vector<DOFVectorBase<double>*> vecs; 
3200
3201
3202
3203

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3204
    std::vector<double*> vecsAtQPs; 
3205
3206
3207
3208

    /** \brief 
     * Function for c. 
     */ 
3209
    AbstractFunction<double, std::vector<double> > *f; 
3210
3211
3212
3213
3214
  }; 

  class VecOfGradientsAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3215
    /// Constructor. 
3216
    VecOfGradientsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3217
			   AbstractFunction<double, std::vector<WorldVector<double>*> > *af);
3218
3219
3220
3221
3222
3223
3224
3225
3226
3227

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3228
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3229
3230
3231
3232

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3233
3234
    void eval(int nPoints,
	      const double *uhAtQP,
3235
3236
3237
3238
3239
3240
3241
3242
3243
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
    /** \brief 
     * Vector of DOFVectors to be evaluated at quadrature points. 
     */ 
3244
    std::vector<DOFVectorBase<double>*> vecs; 
3245
3246
3247
3248

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3249
    std::vector<WorldVector<double>*> gradsAtQPs; 
3250
3251
3252
3253

    /** \brief 
     * Function for c. 
     */ 
3254
    AbstractFunction<double, std::vector<WorldVector<double>*> > *f; 
3255
3256
3257
3258
3259
3260
  };


  class VecDivergence_ZOT : public ZeroOrderTerm
  {
  public: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3261
3262
    /// Constructor. 
    VecDivergence_ZOT(int nComponents,
3263
3264
		      DOFVectorBase<double> *vec0,
		      DOFVectorBase<double> *vec1 = NULL,
Thomas Witkowski's avatar
Thomas Witkowski committed
3265
		      DOFVectorBase<double> *vec2 = NULL);
3266
3267
3268
3269
3270
3271
3272
3273
3274
3275

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3276
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3277
3278
3279
3280

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3281
3282
    void eval(int nPoints,
	      const double *uhAtQP,
3283
3284
3285
3286
3287
3288
3289
3290
3291
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
    /** \brief 
     * Vector of DOFVectors to be evaluated at quadrature points. 
     */ 
3292
    std::vector<DOFVectorBase<double>*> vecs; 
3293
3294
3295
3296

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3297
    std::vector<WorldVector<double>*> gradsAtQPs; 
3298
3299
3300
3301
3302
3303
3304
  };



  class VecAndVecOfGradientsAtQP_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3305
3306
    /// Constructor.
    VecAndVecOfGradientsAtQP_ZOT(DOFVector<double> *v,
3307
				 const std::vector<DOFVector<double>*>& dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
3308
				 BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *af);
3309
3310
3311
3312
3313
3314
3315
3316
3317
3318

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3319
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3320
3321
3322
3323

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3324
3325
    void eval(int nPoints,
	      const double *uhAtQP,
3326
3327
3328
3329
3330
3331
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
3332
    /// DOFVector to be evaluated at quadrature points.
3333
3334
    DOFVector<double>* vec;

Thomas Witkowski's avatar
Thomas Witkowski committed
3335
    /// Vector v at quadrature points.
3336
3337
    double *vecAtQPs;

Thomas Witkowski's avatar
Thomas Witkowski committed
3338
    /// Vector of DOFVectors to be evaluated at quadrature points.
3339
3340
    std::vector<DOFVector<double>*> vecs;

Thomas Witkowski's avatar
Thomas Witkowski committed
3341
    /// Vectors at quadrature points.
3342
3343
    std::vector<WorldVector<double>*> gradsAtQPs;

Thomas Witkowski's avatar
Thomas Witkowski committed
3344
    /// Function for c.
3345
3346
3347
    BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *f;
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
3348
3349
3350
3351
3352
3353
3354
3355
3356
3357
3358
3359
3360
3361
3362
3363
3364
3365
3366
3367
3368
3369
3370
3371
3372
3373
3374
3375
3376
3377
3378
3379
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
3411
3412
3413
3414
3415

  class Vec2AndGrad2AtQP_ZOT : public ZeroOrderTerm
  {
  public:
    Vec2AndGrad2AtQP_ZOT(DOFVectorBase<double> *dv1, 
			 DOFVectorBase<double> *dv2,
			QuartAbstractFunction<double, double, double, WorldVector<double>, WorldVector<double> > *af);

    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;

    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
    DOFVectorBase<double>* vec1;
    DOFVectorBase<double>* vec2;

    double *vecAtQPs1;
    double *vecAtQPs2;

    WorldVector<double> *gradAtQPs1;
    WorldVector<double> *gradAtQPs2;

    QuartAbstractFunction<double, double, double, WorldVector<double>,WorldVector<double> > *f;
  };


  class Vec2AndGradVecAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
    Vec2AndGradVecAtQP_ZOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, 
			  DOFVectorBase<double> *dGrd,
			  TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f);

    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;

    void eval(int nPoints,
	      const double *uhAtQP,
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
    DOFVectorBase<double>* vec1; 
    DOFVectorBase<double>* vec2; 
    
    double *vec1AtQPs; 
    double *vec2AtQPs; 

    DOFVectorBase<double>* vecGrd; 
    WorldVector<double> *gradAtQPs; 

    TertiaryAbstractFunction<double, double, double, WorldVector<double> > *f; 
  }; 



3416
3417
3418
  class General_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3419
    /// Constructor.
3420
3421
    General_ZOT(std::vector<DOFVectorBase<double>*> vecs,
		std::vector<DOFVectorBase<double>*> grads,
Thomas Witkowski's avatar
Thomas Witkowski committed
3422
		TertiaryAbstractFunction<double, WorldVector<double>, std::vector<double>, std::vector<WorldVector<double> > > *f);
3423

Thomas Witkowski's avatar
Thomas Witkowski committed
3424
    /// Implementation of \ref OperatorTerm::initElement().
3425
3426
3427
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

Thomas Witkowski's avatar
Thomas Witkowski committed
3428
    /// Implements ZeroOrderTerm::getC().
3429
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3430

Thomas Witkowski's avatar
Thomas Witkowski committed
3431
    /// Implements ZeroOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
3432
3433
    void eval(int nPoints,
	      const double *uhAtQP,
3434
3435
3436
3437
3438
3439
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
3440
    std::vector<DOFVectorBase<double>*> vecs_; 
3441

3442
    std::vector<DOFVectorBase<double>*> grads_;
3443
3444
3445

    TertiaryAbstractFunction<double, 
			     WorldVector<double>,
3446
3447
			     std::vector<double>, 
			     std::vector<WorldVector<double> > > *f_;
3448
3449
3450

    WorldVector<double> *coordsAtQPs_;

3451
    std::vector<double*> vecsAtQPs_;
3452

3453
    std::vector<WorldVector<double>*> gradsAtQPs_;
3454
3455
3456
3457
3458
  };

  class GeneralParametric_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3459
    /// Constructor.
3460
3461
    GeneralParametric_ZOT(std::vector<DOFVectorBase<double>*> vecs,
			  std::vector<DOFVectorBase<double>*> grads,
3462
3463
3464
			  QuartAbstractFunction<double, 
			  WorldVector<double>,
			  WorldVector<double>,
3465
			  std::vector<double>, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3466
			  std::vector<WorldVector<double> > > *f);
3467
3468
3469
3470
3471
3472
3473
3474
3475
3476

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3477
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3478
3479
3480
3481

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3482
3483
    void eval(int nPoints,
	      const double *uhAtQP,
3484
3485
3486
3487
3488
3489
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
3490
    std::vector<DOFVectorBase<double>*> vecs_; 
3491

3492
    std::vector<DOFVectorBase<double>*> grads_;
3493
3494
3495
3496

    QuartAbstractFunction<double, 
			  WorldVector<double>,
			  WorldVector<double>,
3497
3498
			  std::vector<double>, 
			  std::vector<WorldVector<double> > > *f_;
3499
3500
3501
3502
3503

    WorldVector<double> *coordsAtQPs_;

    WorldVector<double> elementNormal_;

3504
    std::vector<double*> vecsAtQPs_;
3505

3506
    std::vector<WorldVector<double>*> gradsAtQPs_;
3507
3508
3509
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
  };


  /** \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:
    /** \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
3538
3539
	     const FiniteElemSpace *rowFESpace,
	     const FiniteElemSpace *colFESpace = NULL);
3540
3541
3542
3543

    /** \brief
     * Destructor.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3544
    virtual ~Operator() {}
3545
3546
3547
3548

    /** \brief
     * Sets \ref optimized.
     */
3549
3550
    inline void useOptimizedAssembler(bool opt) 
    {
3551
      optimized = opt;
Thomas Witkowski's avatar
Thomas Witkowski committed
3552
    }
3553
3554
3555
3556

    /** \brief
     * Returns \ref optimized.
     */
3557
3558
    inline bool isOptimized() 
    {
3559
      return optimized;
Thomas Witkowski's avatar
Thomas Witkowski committed
3560
    }
3561
3562
3563
3564

    /** \brief
     * Adds a SecondOrderTerm to the Operator
     */
3565
3566
    template <typename T>
    void addSecondOrderTerm(T *term);
3567
3568
3569
3570

    /** \brief
     * Adds a FirstOrderTerm to the Operator
     */
3571
3572
    template <typename T>
    void addFirstOrderTerm(T *term, 
3573
3574
3575
3576
			   FirstOrderType type = GRD_PHI);
    /** \brief
     * Adds a ZeroOrderTerm to the Operator
     */
3577
3578
    template <typename T>
    void addZeroOrderTerm(T *term);
3579
3580
3581
3582
3583
3584
3585


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

3589
3590
    virtual void getElementMatrix(const ElInfo *rowElInfo,
				  const ElInfo *colElInfo,
3591
3592
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
3593
				  ElementMatrix& userMat,
3594
3595
				  double factor = 1.0);

3596
3597
3598
3599
3600
    /** \brief
     * Calculates the element vector for this ElInfo and adds it multiplied by
     * factor to userVec.
     */
    virtual void getElementVector(const ElInfo *elInfo, 
3601
				  ElementVector& userVec, 
3602
3603
				  double factor = 1.0);

Thomas Witkowski's avatar
Thomas Witkowski committed
3604
3605
3606
3607
    virtual void getElementVector(const ElInfo *mainElInfo, 
				  const ElInfo *auxElInfo,
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
3608
				  ElementVector& userVec,
Thomas Witkowski's avatar
Thomas Witkowski committed
3609
				  double factor = 1.0);
3610

Thomas Witkowski's avatar
Thomas Witkowski committed
3611
3612
    /// That function must be called after one assembling cycle has been finished.
    void finishAssembling();
3613

Thomas Witkowski's avatar
Thomas Witkowski committed
3614
    /// Returns \ref rowFESpace.
3615
3616
    inline const FiniteElemSpace *getRowFESpace() 
    { 
3617
      return rowFESpace; 
3618
    }
3619

Thomas Witkowski's avatar
Thomas Witkowski committed
3620
    /// Returns \ref colFESpace.
3621
3622
    inline const FiniteElemSpace *getColFESpace() 
    { 
3623
      return colFESpace; 
3624
    }
3625

Thomas Witkowski's avatar
Thomas Witkowski committed
3626
    /// Returns \ref auxFESpaces.
3627
3628
    inline std::vector<const FiniteElemSpace*> getAuxFESpaces()
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
3629
3630
3631
      return auxFESpaces;
    }

3632
3633
3634
3635
3636
3637
3638
3639
    /** \brief
     * Sets \ref uhOld.
     */
    void setUhOld(const DOFVectorBase<double> *uhOld);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3645
    /// Returns \ref assembler
3646
    Assembler *getAssembler(int rank);
3647

Thomas Witkowski's avatar
Thomas Witkowski committed
3648
    /// Sets \ref assembler
3649
    void setAssembler(int rank, Assembler *ass);
3650

Thomas Witkowski's avatar
Thomas Witkowski committed
3651
    /// Returns whether this is a matrix operator.
3652
3653
    inline const bool isMatrixOperator() 
    {
3654
      return type.isSet(MATRIX_OPERATOR);
3655
    }
3656

Thomas Witkowski's avatar
Thomas Witkowski committed
3657
    /// Returns whether this is a vector operator
3658
3659
    inline const bool isVectorOperator() 
    {
3660
      return type.isSet(VECTOR_OPERATOR);
3661
    }
3662

Thomas Witkowski's avatar
Thomas Witkowski committed
3663
    /// Sets \ref fillFlag, the flag used for this Operator while mesh traversal.
3664
3665
    inline void setFillFlag(Flag f) 
    { 
3666
      fillFlag = f; 
3667
    }
3668

Thomas Witkowski's avatar
Thomas Witkowski committed
3669
    /// Sets \ref needDualTraverse.
3670
3671
    void setNeedDualTraverse(bool b) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
3672
3673
3674
3675
      needDualTraverse = b;
    }

    /// Returns \ref fillFlag
3676
3677
    inline Flag getFillFlag() 
    { 
3678
      return fillFlag; 
3679
    }
3680

Thomas Witkowski's avatar
Thomas Witkowski committed
3681
    /// Returns \ref needDualTraverse
3682
3683
    bool getNeedDualTraverse() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
3684
3685
3686
3687
      return needDualTraverse;
    }

    /// Initialization of the needed SubAssemblers using the given quadratures. 
3688
3689
    void initAssembler(int rank,
		       Quadrature *quad2,
3690
3691
3692
3693
3694
		       Quadrature *quad1GrdPsi,
		       Quadrature *quad1GrdPhi,
		       Quadrature *quad0);


Thomas Witkowski's avatar
Thomas Witkowski committed
3695
    /// Calculates the needed quadrature degree for the given order. 
3696
3697
    int getQuadratureDegree(int order, FirstOrderType firstOrderType = GRD_PHI);

Thomas Witkowski's avatar
Thomas Witkowski committed
3698
3699
3700
    /// Evaluation of all terms in \ref zeroOrder. 
    void evalZeroOrder(int nPoints,
		       const double *uhAtQP,
3701
3702
3703
3704
3705
		       const WorldVector<double> *grdUhAtQP,
		       const WorldMatrix<double> *D2UhAtQP,
		       double *result,
		       double factor) const
    {
3706
3707
      int myRank = omp_get_thread_num();

3708
      std::vector<OperatorTerm*>::const_iterator termIt;
3709
3710
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
3711
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3712
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
3713
    }
3714

3715
    /// Evaluation of all terms in \ref firstOrderGrdPsi. 
Thomas Witkowski's avatar
Thomas Witkowski committed
3716
3717
    void evalFirstOrderGrdPsi(int nPoints,
			      const double *uhAtQP,
3718
3719
3720
3721
3722
			      const WorldVector<double> *grdUhAtQP,
			      const WorldMatrix<double> *D2UhAtQP,
			      double *result,
			      double factor) const
    {
3723
3724
      int myRank = omp_get_thread_num();

3725
      std::vector<OperatorTerm*>::const_iterator termIt;
3726
3727
      for (termIt = firstOrderGrdPsi[myRank].begin(); 
	   termIt != firstOrderGrdPsi[myRank].end(); 
3728
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3729
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
3730
    }
3731

3732
    /// Evaluation of all terms in \ref firstOrderGrdPhi. 
Thomas Witkowski's avatar
Thomas Witkowski committed
3733
3734
    void evalFirstOrderGrdPhi(int nPoints,
			      const double *uhAtQP,
3735
3736
3737
3738
3739
			      const WorldVector<double> *grdUhAtQP,
			      const WorldMatrix<double> *D2UhAtQP,
			      double *result,
			      double factor) const
    {
3740
3741
      int myRank = omp_get_thread_num();

3742
      std::vector<OperatorTerm*>::const_iterator termIt;
3743
3744
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
3745
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3746
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
3747
    }
3748

3749
    /// Evaluation of all terms in \ref secondOrder. 
Thomas Witkowski's avatar
Thomas Witkowski committed
3750
3751
    void evalSecondOrder(int nPoints,
			 const double *uhAtQP,
3752
3753
3754
3755
3756
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
    {
3757
3758
      int myRank = omp_get_thread_num();

3759
      std::vector<OperatorTerm*>::const_iterator termIt;
3760
3761
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
3762
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3763
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
3764
    }
3765

3766
    /// Weak evaluation of all terms in \ref secondOrder. 
Thomas Witkowski's avatar
Thomas Witkowski committed
3767
    void weakEvalSecondOrder(int nPoints,
3768
3769
3770
			     const WorldVector<double> *grdUhAtQP,
			     WorldVector<double> *result) const
    {
3771
3772
      int myRank = omp_get_thread_num();

3773
      std::vector<OperatorTerm*>::const_iterator termIt;
3774
3775
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
3776
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3777
	static_cast<SecondOrderTerm*>(*termIt)->weakEval(nPoints, grdUhAtQP, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
3778
    }
3779
  
3780
    /// Calls getLALt() for each term in \ref secondOrder and adds the results to LALt.
Thomas Witkowski's avatar
Thomas Witkowski committed
3781
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
3782
    {
3783
3784
      int myRank = omp_get_thread_num();

3785
      std::vector<OperatorTerm*>::const_iterator termIt;
3786
3787
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
3788
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3789
	static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, nPoints, LALt);
Thomas Witkowski's avatar
Thomas Witkowski committed
3790
    }
3791
  
3792
3793
3794
3795
    /// Calls getLb() for each term in \ref firstOrderGrdPsi and adds the results to Lb.
    void getLbGrdPsi(const ElInfo *elInfo, 
		     int nPoints, 
		     VectorOfFixVecs<DimVec<double> >& Lb) const
3796
    {
3797
3798
      int myRank = omp_get_thread_num();

3799
      std::vector<OperatorTerm*>::const_iterator termIt;
3800
3801
      for (termIt = firstOrderGrdPsi[myRank].begin(); 
	   termIt != firstOrderGrdPsi[myRank].end(); 
3802
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3803
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
Thomas Witkowski's avatar
Thomas Witkowski committed
3804
    }
3805

3806
3807
3808
3809
    /// Calls getLb() for each term in \ref firstOrderGrdPhi and adds the results to Lb.
    void getLbGrdPhi(const ElInfo *elInfo, 
		     int nPoints, 
		     VectorOfFixVecs<DimVec<double> >& Lb) const
3810
    {
3811
3812
      int myRank = omp_get_thread_num();

3813
      std::vector<OperatorTerm*>::const_iterator termIt;
3814
3815
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
3816
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3817
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
Thomas Witkowski's avatar
Thomas Witkowski committed
3818
    }
3819

3820
    /// Calls getC() for each term in \ref zeroOrder and adds the results to c.
3821
    void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &c) const
3822
    {
3823
3824
      int myRank = omp_get_thread_num();

3825
      std::vector<OperatorTerm*>::const_iterator termIt;
3826
3827
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
3828
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3829
	static_cast<ZeroOrderTerm*>(*termIt)->getC(elInfo, nPoints, c);
Thomas Witkowski's avatar
Thomas Witkowski committed
3830
    }
3831

Thomas Witkowski's avatar
Thomas Witkowski committed
3832
    /// Returns true, if there are second order terms. Returns false otherwise.
3833
3834
    inline bool secondOrderTerms() 
    {
3835
      return secondOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3836
    }
3837

3838
    /// Returns true, if there are first order terms (grdPsi). Returns false otherwise.
3839
3840
    inline bool firstOrderTermsGrdPsi() 
    {
3841
      return firstOrderGrdPsi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3842
    }
3843

3844
    /// Returns true, if there are first order terms (grdPhi). Returns false otherwise.
3845
3846
    inline bool firstOrderTermsGrdPhi() 
    {
3847
      return firstOrderGrdPhi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3848
    }
3849

3850
    /// Returns true, if there are zero order terms. Returns false otherwise.
3851
3852
    inline bool zeroOrderTerms() 
    {
3853
      return zeroOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3854
    }
3855
3856

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3857
    /// Constant type flag for matrix operators
3858
3859
    static const Flag MATRIX_OPERATOR;

Thomas Witkowski's avatar
Thomas Witkowski committed
3860
    /// Constant type flag for vector operators
3861
3862
3863
    static const Flag VECTOR_OPERATOR;

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
3864
    /// FiniteElemSpace for matrix rows and element vector
3865
3866
    const FiniteElemSpace *rowFESpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
3867
    /// FiniteElemSpace for matrix columns. Can be equal to \rowFESpace.
3868
3869
    const FiniteElemSpace *colFESpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
3870
3871
3872
3873
    /// List of aux fe spaces, e.g., if a term is multiplied with a DOF vector
    std::vector<const FiniteElemSpace*> auxFESpaces;

    /// Number of rows in the element matrix
3874
3875
    int nRow;

Thomas Witkowski's avatar
Thomas Witkowski committed
3876
    /// Number of columns in the element matrix
3877
3878
    int nCol;

Thomas Witkowski's avatar
Thomas Witkowski committed
3879
    /// Type of this Operator.
3880
3881
    Flag type;

Thomas Witkowski's avatar
Thomas Witkowski committed
3882
    /// Flag for mesh traversal
3883
3884
    Flag fillFlag;

Thomas Witkowski's avatar
Thomas Witkowski committed
3885
3886
3887
    /// If true, the operator needs a dual traverse over two meshes for assembling.
    bool needDualTraverse;

3888
3889
3890
3891
3892
    /** \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.
     */
3893
    std::vector<Assembler*> assembler;
3894

Thomas Witkowski's avatar
Thomas Witkowski committed
3895
    /// List of all second order terms
3896
    std::vector< std::vector<OperatorTerm*> > secondOrder;
3897

Thomas Witkowski's avatar
Thomas Witkowski committed
3898
    /// List of all first order terms derived to psi
3899
    std::vector< std::vector<OperatorTerm*> > firstOrderGrdPsi;
3900

Thomas Witkowski's avatar
Thomas Witkowski committed
3901
    /// List of all first order terms derived to phi
3902
    std::vector< std::vector<OperatorTerm*> > firstOrderGrdPhi;
3903

Thomas Witkowski's avatar
Thomas Witkowski committed
3904
    /// List of all zero order terms
3905
    std::vector< std::vector<OperatorTerm*> > zeroOrder;
3906
3907
3908
3909
3910
3911
3912
3913
3914

    /** \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;

Thomas Witkowski's avatar
Thomas Witkowski committed
3915
    /// Spezifies whether optimized assemblers are used or not.
3916
3917
3918
3919
3920
3921
3922
3923
3924
3925
3926
    bool optimized;

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

}

3927
3928
#include "Operator.hh"

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