Operator.h 106 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
    /** \brief
     * Stores coordinates at quadrature points. Set in \ref initElement().
     */
    WorldVector<double>* coordsAtQPs;

    /** \brief
     * Function evaluated at quadrature points.
     */
    AbstractFunction<double, WorldVector<double> > *g;

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

  /**
   * \ingroup Assembler
   *  
   * \brief
   * SecondOrderTerm where A is a WorldMatrix having a in all positions
   * except possibly the position IJ, multiplied with a function
   * evaluated at the quadrature points of a given DOFVector:
   * \f$ \nabla \cdot f(v(\vec{x})) A \nabla u(\vec{x}) \f$
   */
  class VecAtQP_IJ_SOT : public SecondOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3029
3030
3031
    /// Constructor.
    VecAtQP_IJ_SOT(DOFVectorBase<double> *dv, AbstractFunction<double, double> *af,
		   int x_i, int x_j);
3032
3033
3034
3035
3036
3037
3038
3039
3040
3041

    /** \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
3042
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const;    
3043
3044
3045
3046

    /** \brief
     * Implements SecondOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3047
3048
    void eval(int nPoints,
	      const double *uhAtQP,
3049
3050
3051
3052
3053
3054
3055
3056
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double factor) const;

    /** \brief
     * Implements SecondOrderTerm::weakEval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3057
    void weakEval(int nPoints,
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
		  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
3088
    /// Constructor. 
3089
3090
    VecAndGradVecAtQP_ZOT(DOFVectorBase<double> *dv, 
			  DOFVectorBase<double> *dGrd,
Thomas Witkowski's avatar
Thomas Witkowski committed
3091
			  BinaryAbstractFunction<double, double, WorldVector<double> > *f);
3092
3093
3094
3095
3096
3097
3098
3099
3100
3101

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3102
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3103
3104
3105
3106

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3107
3108
    void eval(int nPoints,
	      const double *uhAtQP,
3109
3110
3111
3112
3113
3114
3115
3116
3117
3118
3119
3120
3121
3122
3123
3124
3125
3126
3127
3128
3129
3130
3131
3132
3133
3134
3135
3136
3137
3138
3139
3140
3141
3142
3143
	      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
3144
    /// Constructor. 
3145
    VecAndGradVec2AtQP_ZOT(DOFVectorBase<double> *dv, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3146
3147
3148
			   DOFVectorBase<double> *dGrd1, 
			   DOFVectorBase<double> *dGrd2, 
			   TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *af);
3149
3150
3151
3152
3153
3154
3155
3156
3157
3158

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3159
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3160
3161
3162
3163

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3164
3165
    void eval(int nPoints,
	      const double *uhAtQP,
3166
3167
3168
3169
3170
3171
3172
3173
3174
3175
3176
3177
3178
3179
3180
3181
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
	      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
3212
    /// Constructor. 
3213
    VecOfDOFVecsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3214
			 AbstractFunction<double, std::vector<double> > *f);
3215
3216
3217
3218
3219
3220
3221
3222
3223
3224

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

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

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

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

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3246
    std::vector<double*> vecsAtQPs; 
3247
3248
3249
3250

    /** \brief 
     * Function for c. 
     */ 
3251
    AbstractFunction<double, std::vector<double> > *f; 
3252
3253
3254
3255
3256
  }; 

  class VecOfGradientsAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3257
    /// Constructor. 
3258
    VecOfGradientsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3259
			   AbstractFunction<double, std::vector<WorldVector<double>*> > *af);
3260
3261
3262
3263
3264
3265
3266
3267
3268
3269

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3270
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3271
3272
3273
3274

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3275
3276
    void eval(int nPoints,
	      const double *uhAtQP,
3277
3278
3279
3280
3281
3282
3283
3284
3285
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
    /** \brief 
     * Vector of DOFVectors to be evaluated at quadrature points. 
     */ 
3286
    std::vector<DOFVectorBase<double>*> vecs; 
3287
3288
3289
3290

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3291
    std::vector<WorldVector<double>*> gradsAtQPs; 
3292
3293
3294
3295

    /** \brief 
     * Function for c. 
     */ 
3296
    AbstractFunction<double, std::vector<WorldVector<double>*> > *f; 
3297
3298
3299
3300
3301
3302
  };


  class VecDivergence_ZOT : public ZeroOrderTerm
  {
  public: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3303
3304
    /// Constructor. 
    VecDivergence_ZOT(int nComponents,
3305
3306
		      DOFVectorBase<double> *vec0,
		      DOFVectorBase<double> *vec1 = NULL,
Thomas Witkowski's avatar
Thomas Witkowski committed
3307
		      DOFVectorBase<double> *vec2 = NULL);
3308
3309
3310
3311
3312
3313
3314
3315
3316
3317

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

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

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

  protected: 
    /** \brief 
     * Vector of DOFVectors to be evaluated at quadrature points. 
     */ 
3334
    std::vector<DOFVectorBase<double>*> vecs; 
3335
3336
3337
3338

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3339
    std::vector<WorldVector<double>*> gradsAtQPs; 
3340
3341
3342
3343
3344
3345
3346
  };



  class VecAndVecOfGradientsAtQP_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3347
3348
    /// Constructor.
    VecAndVecOfGradientsAtQP_ZOT(DOFVector<double> *v,
3349
				 const std::vector<DOFVector<double>*>& dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
3350
				 BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *af);
3351
3352
3353
3354
3355
3356
3357
3358
3359
3360

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3361
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3362
3363
3364
3365

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3366
3367
    void eval(int nPoints,
	      const double *uhAtQP,
3368
3369
3370
3371
3372
3373
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
3374
    /// DOFVector to be evaluated at quadrature points.
3375
3376
    DOFVector<double>* vec;

Thomas Witkowski's avatar
Thomas Witkowski committed
3377
    /// Vector v at quadrature points.
3378
3379
    double *vecAtQPs;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3383
    /// Vectors at quadrature points.
3384
3385
    std::vector<WorldVector<double>*> gradsAtQPs;

Thomas Witkowski's avatar
Thomas Witkowski committed
3386
    /// Function for c.
3387
3388
3389
    BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *f;
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
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
3416
3417
3418
3419
3420
3421
3422
3423
3424
3425
3426
3427
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

  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; 
  }; 



3458
3459
3460
  class General_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3461
    /// Constructor.
3462
3463
    General_ZOT(std::vector<DOFVectorBase<double>*> vecs,
		std::vector<DOFVectorBase<double>*> grads,
Thomas Witkowski's avatar
Thomas Witkowski committed
3464
		TertiaryAbstractFunction<double, WorldVector<double>, std::vector<double>, std::vector<WorldVector<double> > > *f);
3465

Thomas Witkowski's avatar
Thomas Witkowski committed
3466
    /// Implementation of \ref OperatorTerm::initElement().
3467
3468
3469
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3473
    /// Implements ZeroOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
3474
3475
    void eval(int nPoints,
	      const double *uhAtQP,
3476
3477
3478
3479
3480
3481
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
3482
    std::vector<DOFVectorBase<double>*> vecs_; 
3483

3484
    std::vector<DOFVectorBase<double>*> grads_;
3485
3486
3487

    TertiaryAbstractFunction<double, 
			     WorldVector<double>,
3488
3489
			     std::vector<double>, 
			     std::vector<WorldVector<double> > > *f_;
3490
3491
3492

    WorldVector<double> *coordsAtQPs_;

3493
    std::vector<double*> vecsAtQPs_;
3494

3495
    std::vector<WorldVector<double>*> gradsAtQPs_;
3496
3497
3498
3499
3500
  };

  class GeneralParametric_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3501
    /// Constructor.
3502
3503
    GeneralParametric_ZOT(std::vector<DOFVectorBase<double>*> vecs,
			  std::vector<DOFVectorBase<double>*> grads,
3504
3505
3506
			  QuartAbstractFunction<double, 
			  WorldVector<double>,
			  WorldVector<double>,
3507
			  std::vector<double>, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3508
			  std::vector<WorldVector<double> > > *f);
3509
3510
3511
3512
3513
3514
3515
3516
3517
3518

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3519
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3520
3521
3522
3523

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3524
3525
    void eval(int nPoints,
	      const double *uhAtQP,
3526
3527
3528
3529
3530
3531
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
3532
    std::vector<DOFVectorBase<double>*> vecs_; 
3533

3534
    std::vector<DOFVectorBase<double>*> grads_;
3535
3536
3537
3538

    QuartAbstractFunction<double, 
			  WorldVector<double>,
			  WorldVector<double>,
3539
3540
			  std::vector<double>, 
			  std::vector<WorldVector<double> > > *f_;
3541
3542
3543
3544
3545

    WorldVector<double> *coordsAtQPs_;

    WorldVector<double> elementNormal_;

3546
    std::vector<double*> vecsAtQPs_;
3547

3548
    std::vector<WorldVector<double>*> gradsAtQPs_;
3549
3550
3551
3552
3553
3554
3555
3556
3557
3558
3559
3560
3561
3562
3563
3564
3565
3566
3567
3568
3569
3570
3571
3572
3573
3574
3575
3576
3577
3578
3579
  };


  /** \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
3580
3581
	     const FiniteElemSpace *rowFESpace,
	     const FiniteElemSpace *colFESpace = NULL);
3582

3583
    /// Destructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
3584
    virtual ~Operator() {}
3585

3586
    /// Sets \ref optimized.
3587
3588
    inline void useOptimizedAssembler(bool opt) 
    {
3589
      optimized = opt;
Thomas Witkowski's avatar
Thomas Witkowski committed
3590
    }
3591

3592
    /// Returns \ref optimized.
3593
3594
    inline bool isOptimized() 
    {
3595
      return optimized;
Thomas Witkowski's avatar
Thomas Witkowski committed
3596
    }
3597

3598
    /// Adds a SecondOrderTerm to the Operator
3599
3600
    template <typename T>
    void addSecondOrderTerm(T *term);
3601

3602
    /// Adds a FirstOrderTerm to the Operator
3603
    template <typename T>
3604
3605
3606
    void addFirstOrderTerm(T *term, FirstOrderType type = GRD_PHI);

    /// Adds a ZeroOrderTerm to the Operator
3607
3608
    template <typename T>
    void addZeroOrderTerm(T *term);
3609
3610
3611
3612
3613
3614
3615


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

3619
3620
    virtual void getElementMatrix(const ElInfo *rowElInfo,
				  const ElInfo *colElInfo,
3621
3622
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
3623
				  ElementMatrix& userMat,
3624
3625
				  double factor = 1.0);

3626
3627
3628
3629
3630
    /** \brief
     * Calculates the element vector for this ElInfo and adds it multiplied by
     * factor to userVec.
     */
    virtual void getElementVector(const ElInfo *elInfo, 
3631
				  ElementVector& userVec, 
3632
3633
				  double factor = 1.0);

Thomas Witkowski's avatar
Thomas Witkowski committed
3634
3635
3636
3637
    virtual void getElementVector(const ElInfo *mainElInfo, 
				  const ElInfo *auxElInfo,
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
3638
				  ElementVector& userVec,
Thomas Witkowski's avatar
Thomas Witkowski committed
3639
				  double factor = 1.0);
3640

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3644
    /// Returns \ref rowFESpace.
3645
3646
    inline const FiniteElemSpace *getRowFESpace() 
    { 
3647
      return rowFESpace; 
3648
    }
3649

Thomas Witkowski's avatar
Thomas Witkowski committed
3650
    /// Returns \ref colFESpace.
3651
3652
    inline const FiniteElemSpace *getColFESpace() 
    { 
3653
      return colFESpace; 
3654
    }
3655

Thomas Witkowski's avatar
Thomas Witkowski committed
3656
    /// Returns \ref auxFESpaces.
3657
3658
    inline std::vector<const FiniteElemSpace*> getAuxFESpaces()
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
3659
3660
3661
      return auxFESpaces;
    }

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

3665
    /// Returns \ref uhOld.
3666
3667
    inline const DOFVectorBase<double> *getUhOld() 
    {
3668
      return uhOld;
3669
    }    
3670

Thomas Witkowski's avatar
Thomas Witkowski committed
3671
    /// Returns \ref assembler
3672
    Assembler *getAssembler(int rank);
3673

Thomas Witkowski's avatar
Thomas Witkowski committed
3674
    /// Sets \ref assembler
3675
    void setAssembler(int rank, Assembler *ass);
3676

Thomas Witkowski's avatar
Thomas Witkowski committed
3677
    /// Returns whether this is a matrix operator.
3678
3679
    inline const bool isMatrixOperator() 
    {
3680
      return type.isSet(MATRIX_OPERATOR);
3681
    }
3682

Thomas Witkowski's avatar
Thomas Witkowski committed
3683
    /// Returns whether this is a vector operator
3684
3685
    inline const bool isVectorOperator() 
    {
3686
      return type.isSet(VECTOR_OPERATOR);
3687
    }
3688

Thomas Witkowski's avatar
Thomas Witkowski committed
3689
    /// Sets \ref fillFlag, the flag used for this Operator while mesh traversal.
3690
3691
    inline void setFillFlag(Flag f) 
    { 
3692
      fillFlag = f; 
3693
    }
3694

Thomas Witkowski's avatar
Thomas Witkowski committed
3695
    /// Sets \ref needDualTraverse.
3696
3697
    void setNeedDualTraverse(bool b) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
3698
3699
3700
3701
      needDualTraverse = b;
    }

    /// Returns \ref fillFlag
3702
3703
    inline Flag getFillFlag() 
    { 
3704
      return fillFlag; 
3705
    }
3706

Thomas Witkowski's avatar
Thomas Witkowski committed
3707
    /// Returns \ref needDualTraverse
3708
3709
    bool getNeedDualTraverse() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
3710
3711
3712
3713
      return needDualTraverse;
    }

    /// Initialization of the needed SubAssemblers using the given quadratures. 
3714
3715
    void initAssembler(int rank,
		       Quadrature *quad2,
3716
3717
3718
3719
3720
		       Quadrature *quad1GrdPsi,
		       Quadrature *quad1GrdPhi,
		       Quadrature *quad0);


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

Thomas Witkowski's avatar
Thomas Witkowski committed
3724
3725
3726
    /// Evaluation of all terms in \ref zeroOrder. 
    void evalZeroOrder(int nPoints,
		       const double *uhAtQP,
3727
3728
3729
3730
3731
		       const WorldVector<double> *grdUhAtQP,
		       const WorldMatrix<double> *D2UhAtQP,
		       double *result,
		       double factor) const
    {
3732
3733
      int myRank = omp_get_thread_num();

3734
      std::vector<OperatorTerm*>::const_iterator termIt;
3735
3736
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
3737
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3738
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
3739
    }
3740

3741
    /// Evaluation of all terms in \ref firstOrderGrdPsi. 
Thomas Witkowski's avatar
Thomas Witkowski committed
3742
3743
    void evalFirstOrderGrdPsi(int nPoints,
			      const double *uhAtQP,
3744
3745
3746
3747
3748
			      const WorldVector<double> *grdUhAtQP,
			      const WorldMatrix<double> *D2UhAtQP,
			      double *result,
			      double factor) const
    {
3749
3750
      int myRank = omp_get_thread_num();

3751
      std::vector<OperatorTerm*>::const_iterator termIt;
3752
3753
      for (termIt = firstOrderGrdPsi[myRank].begin(); 
	   termIt != firstOrderGrdPsi[myRank].end(); 
3754
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3755
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
3756
    }
3757

3758
    /// Evaluation of all terms in \ref firstOrderGrdPhi. 
Thomas Witkowski's avatar
Thomas Witkowski committed
3759
3760
    void evalFirstOrderGrdPhi(int nPoints,
			      const double *uhAtQP,
3761
3762
3763
3764
3765
			      const WorldVector<double> *grdUhAtQP,
			      const WorldMatrix<double> *D2UhAtQP,
			      double *result,
			      double factor) const
    {
3766
3767
      int myRank = omp_get_thread_num();

3768
      std::vector<OperatorTerm*>::const_iterator termIt;
3769
3770
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
3771
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3772
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
3773
    }
3774

3775
    /// Evaluation of all terms in \ref secondOrder. 
Thomas Witkowski's avatar
Thomas Witkowski committed
3776
3777
    void evalSecondOrder(int nPoints,
			 const double *uhAtQP,
3778
3779
3780
3781
3782
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
    {
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
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
3790
    }
3791

3792
    /// Weak evaluation of all terms in \ref secondOrder. 
Thomas Witkowski's avatar
Thomas Witkowski committed
3793
    void weakEvalSecondOrder(int nPoints,
3794
3795
3796
			     const WorldVector<double> *grdUhAtQP,
			     WorldVector<double> *result) const
    {
3797
3798
      int myRank = omp_get_thread_num();

3799
      std::vector<OperatorTerm*>::const_iterator termIt;
3800
3801
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
3802
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3803
	static_cast<SecondOrderTerm*>(*termIt)->weakEval(nPoints, grdUhAtQP, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
3804
    }
3805
  
3806
    /// Calls getLALt() for each term in \ref secondOrder and adds the results to LALt.
Thomas Witkowski's avatar
Thomas Witkowski committed
3807
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
3808
    {
3809
3810
      int myRank = omp_get_thread_num();

3811
      std::vector<OperatorTerm*>::const_iterator termIt;
3812
3813
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
3814
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3815
	static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, nPoints, LALt);
Thomas Witkowski's avatar
Thomas Witkowski committed
3816
    }
3817
  
3818
3819
3820
3821
    /// 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
3822
    {
3823
3824
      int myRank = omp_get_thread_num();

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

3832
3833
3834
3835
    /// 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
3836
    {
3837
3838
      int myRank = omp_get_thread_num();

3839
      std::vector<OperatorTerm*>::const_iterator termIt;
3840
3841
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
3842
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3843
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
Thomas Witkowski's avatar
Thomas Witkowski committed
3844
    }
3845

3846
    /// Calls getC() for each term in \ref zeroOrder and adds the results to c.
3847
    void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &c) const
3848
    {
3849
3850
      int myRank = omp_get_thread_num();

3851
      std::vector<OperatorTerm*>::const_iterator termIt;
3852
3853
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
3854
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3855
	static_cast<ZeroOrderTerm*>(*termIt)->getC(elInfo, nPoints, c);
Thomas Witkowski's avatar
Thomas Witkowski committed
3856
    }
3857

Thomas Witkowski's avatar
Thomas Witkowski committed
3858
    /// Returns true, if there are second order terms. Returns false otherwise.
3859
3860
    inline bool secondOrderTerms() 
    {
3861
      return secondOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3862
    }
3863

3864
    /// Returns true, if there are first order terms (grdPsi). Returns false otherwise.
3865
3866
    inline bool firstOrderTermsGrdPsi() 
    {
3867
      return firstOrderGrdPsi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3868
    }
3869

3870
    /// Returns true, if there are first order terms (grdPhi). Returns false otherwise.
3871
3872
    inline bool firstOrderTermsGrdPhi() 
    {
3873
      return firstOrderGrdPhi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3874
    }
3875

3876
    /// Returns true, if there are zero order terms. Returns false otherwise.
3877
3878
    inline bool zeroOrderTerms() 
    {
3879
      return zeroOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3880
    }
3881
3882

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3883
    /// Constant type flag for matrix operators
3884
3885
    static const Flag MATRIX_OPERATOR;

Thomas Witkowski's avatar
Thomas Witkowski committed
3886
    /// Constant type flag for vector operators
3887
3888
3889
    static const Flag VECTOR_OPERATOR;

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
3890
    /// FiniteElemSpace for matrix rows and element vector
3891
3892
    const FiniteElemSpace *rowFESpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
3893
    /// FiniteElemSpace for matrix columns. Can be equal to \rowFESpace.
3894
3895
    const FiniteElemSpace *colFESpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
3896
3897
3898
3899
    /// 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
3900
3901
    int nRow;

Thomas Witkowski's avatar
Thomas Witkowski committed
3902
    /// Number of columns in the element matrix
3903
3904
    int nCol;

Thomas Witkowski's avatar
Thomas Witkowski committed
3905
    /// Type of this Operator.
3906
3907
    Flag type;

Thomas Witkowski's avatar
Thomas Witkowski committed
3908
    /// Flag for mesh traversal
3909
3910
    Flag fillFlag;

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

3914
3915
3916
3917
3918
    /** \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.
     */
3919
    std::vector<Assembler*> assembler;
3920

Thomas Witkowski's avatar
Thomas Witkowski committed
3921
    /// List of all second order terms
3922
    std::vector< std::vector<OperatorTerm*> > secondOrder;
3923

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3930
    /// List of all zero order terms
3931
    std::vector< std::vector<OperatorTerm*> > zeroOrder;
3932
3933
3934
3935
3936
3937
3938
3939
3940

    /** \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
3941
    /// Spezifies whether optimized assemblers are used or not.
3942
3943
3944
3945
3946
3947
3948
3949
3950
3951
3952
    bool optimized;

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

}

3953
3954
#include "Operator.hh"

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