Operator.h 104 KB
Newer Older
3001
3002
3003
3004

    /** \brief
     * Implements SecondOrderTerm::weakEval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3005
    void weakEval(int nPoints,
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
		  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
3036
    /// Constructor. 
3037
3038
    VecAndGradVecAtQP_ZOT(DOFVectorBase<double> *dv, 
			  DOFVectorBase<double> *dGrd,
Thomas Witkowski's avatar
Thomas Witkowski committed
3039
			  BinaryAbstractFunction<double, double, WorldVector<double> > *f);
3040
3041
3042
3043
3044
3045
3046
3047
3048
3049

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3050
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3051
3052
3053
3054

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3055
3056
    void eval(int nPoints,
	      const double *uhAtQP,
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
3088
3089
3090
3091
	      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
3092
    /// Constructor. 
3093
    VecAndGradVec2AtQP_ZOT(DOFVectorBase<double> *dv, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3094
3095
3096
			   DOFVectorBase<double> *dGrd1, 
			   DOFVectorBase<double> *dGrd2, 
			   TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *af);
3097
3098
3099
3100
3101
3102
3103
3104
3105
3106

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3107
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3108
3109
3110
3111

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3112
3113
    void eval(int nPoints,
	      const double *uhAtQP,
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
3144
3145
3146
3147
3148
3149
3150
3151
3152
3153
3154
3155
3156
3157
3158
3159
	      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
3160
    /// Constructor. 
3161
    VecOfDOFVecsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3162
			 AbstractFunction<double, std::vector<double> > *f);
3163
3164
3165
3166
3167
3168
3169
3170
3171
3172

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3173
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3174
3175
3176
3177

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3178
3179
    void eval(int nPoints,
	      const double *uhAtQP,
3180
3181
3182
3183
3184
3185
3186
3187
3188
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
    /** \brief 
     * Vector of DOFVectors to be evaluated at quadrature points. 
     */ 
3189
    std::vector<DOFVectorBase<double>*> vecs; 
3190
3191
3192
3193

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3194
    std::vector<double*> vecsAtQPs; 
3195
3196
3197
3198

    /** \brief 
     * Function for c. 
     */ 
3199
    AbstractFunction<double, std::vector<double> > *f; 
3200
3201
3202
3203
3204
  }; 

  class VecOfGradientsAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3205
    /// Constructor. 
3206
    VecOfGradientsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3207
			   AbstractFunction<double, std::vector<WorldVector<double>*> > *af);
3208
3209
3210
3211
3212
3213
3214
3215
3216
3217

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3218
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3219
3220
3221
3222

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3223
3224
    void eval(int nPoints,
	      const double *uhAtQP,
3225
3226
3227
3228
3229
3230
3231
3232
3233
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
    /** \brief 
     * Vector of DOFVectors to be evaluated at quadrature points. 
     */ 
3234
    std::vector<DOFVectorBase<double>*> vecs; 
3235
3236
3237
3238

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3239
    std::vector<WorldVector<double>*> gradsAtQPs; 
3240
3241
3242
3243

    /** \brief 
     * Function for c. 
     */ 
3244
    AbstractFunction<double, std::vector<WorldVector<double>*> > *f; 
3245
3246
3247
3248
3249
3250
  };


  class VecDivergence_ZOT : public ZeroOrderTerm
  {
  public: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3251
3252
    /// Constructor. 
    VecDivergence_ZOT(int nComponents,
3253
3254
		      DOFVectorBase<double> *vec0,
		      DOFVectorBase<double> *vec1 = NULL,
Thomas Witkowski's avatar
Thomas Witkowski committed
3255
		      DOFVectorBase<double> *vec2 = NULL);
3256
3257
3258
3259
3260
3261
3262
3263
3264
3265

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3266
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3267
3268
3269
3270

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

  protected: 
    /** \brief 
     * Vector of DOFVectors to be evaluated at quadrature points. 
     */ 
3282
    std::vector<DOFVectorBase<double>*> vecs; 
3283
3284
3285
3286

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3287
    std::vector<WorldVector<double>*> gradsAtQPs; 
3288
3289
3290
3291
3292
3293
3294
  };



  class VecAndVecOfGradientsAtQP_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3295
3296
    /// Constructor.
    VecAndVecOfGradientsAtQP_ZOT(DOFVector<double> *v,
3297
				 const std::vector<DOFVector<double>*>& dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
3298
				 BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *af);
3299
3300
3301
3302
3303
3304
3305
3306
3307
3308

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3309
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3310
3311
3312
3313

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3314
3315
    void eval(int nPoints,
	      const double *uhAtQP,
3316
3317
3318
3319
3320
3321
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
3322
    /// DOFVector to be evaluated at quadrature points.
3323
3324
    DOFVector<double>* vec;

Thomas Witkowski's avatar
Thomas Witkowski committed
3325
    /// Vector v at quadrature points.
3326
3327
    double *vecAtQPs;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3331
    /// Vectors at quadrature points.
3332
3333
    std::vector<WorldVector<double>*> gradsAtQPs;

Thomas Witkowski's avatar
Thomas Witkowski committed
3334
    /// Function for c.
3335
3336
3337
    BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *f;
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
3338
3339
3340
3341
3342
3343
3344
3345
3346
3347
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

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



3406
3407
3408
  class General_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3409
    /// Constructor.
3410
3411
    General_ZOT(std::vector<DOFVectorBase<double>*> vecs,
		std::vector<DOFVectorBase<double>*> grads,
Thomas Witkowski's avatar
Thomas Witkowski committed
3412
		TertiaryAbstractFunction<double, WorldVector<double>, std::vector<double>, std::vector<WorldVector<double> > > *f);
3413

Thomas Witkowski's avatar
Thomas Witkowski committed
3414
    /// Implementation of \ref OperatorTerm::initElement().
3415
3416
3417
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3421
    /// Implements ZeroOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
3422
3423
    void eval(int nPoints,
	      const double *uhAtQP,
3424
3425
3426
3427
3428
3429
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
3430
    std::vector<DOFVectorBase<double>*> vecs_; 
3431

3432
    std::vector<DOFVectorBase<double>*> grads_;
3433
3434
3435

    TertiaryAbstractFunction<double, 
			     WorldVector<double>,
3436
3437
			     std::vector<double>, 
			     std::vector<WorldVector<double> > > *f_;
3438
3439
3440

    WorldVector<double> *coordsAtQPs_;

3441
    std::vector<double*> vecsAtQPs_;
3442

3443
    std::vector<WorldVector<double>*> gradsAtQPs_;
3444
3445
3446
3447
3448
  };

  class GeneralParametric_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3449
    /// Constructor.
3450
3451
    GeneralParametric_ZOT(std::vector<DOFVectorBase<double>*> vecs,
			  std::vector<DOFVectorBase<double>*> grads,
3452
3453
3454
			  QuartAbstractFunction<double, 
			  WorldVector<double>,
			  WorldVector<double>,
3455
			  std::vector<double>, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3456
			  std::vector<WorldVector<double> > > *f);
3457
3458
3459
3460
3461
3462
3463
3464
3465
3466

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
3467
    void getC(const ElInfo *, int nPoints, std::vector<double> &C) const;
3468
3469
3470
3471

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

  protected:
3480
    std::vector<DOFVectorBase<double>*> vecs_; 
3481

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

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

    WorldVector<double> *coordsAtQPs_;

    WorldVector<double> elementNormal_;

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

3496
    std::vector<WorldVector<double>*> gradsAtQPs_;
3497
3498
3499
3500
3501
3502
3503
3504
3505
3506
3507
3508
3509
3510
3511
3512
3513
3514
3515
3516
3517
3518
3519
3520
3521
3522
3523
3524
3525
3526
3527
  };


  /** \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
3528
3529
	     const FiniteElemSpace *rowFESpace,
	     const FiniteElemSpace *colFESpace = NULL);
3530
3531
3532
3533

    /** \brief
     * Destructor.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3534
    virtual ~Operator() {}
3535
3536
3537
3538

    /** \brief
     * Sets \ref optimized.
     */
3539
3540
    inline void useOptimizedAssembler(bool opt) 
    {
3541
      optimized = opt;
Thomas Witkowski's avatar
Thomas Witkowski committed
3542
    }
3543
3544
3545
3546

    /** \brief
     * Returns \ref optimized.
     */
3547
3548
    inline bool isOptimized() 
    {
3549
      return optimized;
Thomas Witkowski's avatar
Thomas Witkowski committed
3550
    }
3551
3552
3553
3554

    /** \brief
     * Adds a SecondOrderTerm to the Operator
     */
3555
3556
    template <typename T>
    void addSecondOrderTerm(T *term);
3557
3558
3559
3560

    /** \brief
     * Adds a FirstOrderTerm to the Operator
     */
3561
3562
    template <typename T>
    void addFirstOrderTerm(T *term, 
3563
3564
3565
3566
			   FirstOrderType type = GRD_PHI);
    /** \brief
     * Adds a ZeroOrderTerm to the Operator
     */
3567
3568
    template <typename T>
    void addZeroOrderTerm(T *term);
3569
3570
3571
3572
3573
3574
3575


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

3579
3580
    virtual void getElementMatrix(const ElInfo *rowElInfo,
				  const ElInfo *colElInfo,
3581
3582
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
3583
				  ElementMatrix& userMat,
3584
3585
				  double factor = 1.0);

3586
3587
3588
3589
3590
    /** \brief
     * Calculates the element vector for this ElInfo and adds it multiplied by
     * factor to userVec.
     */
    virtual void getElementVector(const ElInfo *elInfo, 
3591
				  ElementVector& userVec, 
3592
3593
				  double factor = 1.0);

Thomas Witkowski's avatar
Thomas Witkowski committed
3594
3595
3596
3597
    virtual void getElementVector(const ElInfo *mainElInfo, 
				  const ElInfo *auxElInfo,
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
3598
				  ElementVector& userVec,
Thomas Witkowski's avatar
Thomas Witkowski committed
3599
				  double factor = 1.0);
3600

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3604
    /// Returns \ref rowFESpace.
3605
3606
    inline const FiniteElemSpace *getRowFESpace() 
    { 
3607
      return rowFESpace; 
3608
    }
3609

Thomas Witkowski's avatar
Thomas Witkowski committed
3610
    /// Returns \ref colFESpace.
3611
3612
    inline const FiniteElemSpace *getColFESpace() 
    { 
3613
      return colFESpace; 
3614
    }
3615

Thomas Witkowski's avatar
Thomas Witkowski committed
3616
    /// Returns \ref auxFESpaces.
3617
3618
    inline std::vector<const FiniteElemSpace*> getAuxFESpaces()
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
3619
3620
3621
      return auxFESpaces;
    }

3622
3623
3624
3625
3626
3627
3628
3629
    /** \brief
     * Sets \ref uhOld.
     */
    void setUhOld(const DOFVectorBase<double> *uhOld);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3635
    /// Returns \ref assembler
3636
    Assembler *getAssembler(int rank);
3637

Thomas Witkowski's avatar
Thomas Witkowski committed
3638
    /// Sets \ref assembler
3639
    void setAssembler(int rank, Assembler *ass);
3640

Thomas Witkowski's avatar
Thomas Witkowski committed
3641
    /// Returns whether this is a matrix operator.
3642
3643
    inline const bool isMatrixOperator() 
    {
3644
      return type.isSet(MATRIX_OPERATOR);
3645
    }
3646

Thomas Witkowski's avatar
Thomas Witkowski committed
3647
    /// Returns whether this is a vector operator
3648
3649
    inline const bool isVectorOperator() 
    {
3650
      return type.isSet(VECTOR_OPERATOR);
3651
    }
3652

Thomas Witkowski's avatar
Thomas Witkowski committed
3653
    /// Sets \ref fillFlag, the flag used for this Operator while mesh traversal.
3654
3655
    inline void setFillFlag(Flag f) 
    { 
3656
      fillFlag = f; 
3657
    }
3658

Thomas Witkowski's avatar
Thomas Witkowski committed
3659
    /// Sets \ref needDualTraverse.
3660
3661
    void setNeedDualTraverse(bool b) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
3662
3663
3664
3665
      needDualTraverse = b;
    }

    /// Returns \ref fillFlag
3666
3667
    inline Flag getFillFlag() 
    { 
3668
      return fillFlag; 
3669
    }
3670

Thomas Witkowski's avatar
Thomas Witkowski committed
3671
    /// Returns \ref needDualTraverse
3672
3673
    bool getNeedDualTraverse() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
3674
3675
3676
3677
      return needDualTraverse;
    }

    /// Initialization of the needed SubAssemblers using the given quadratures. 
3678
3679
    void initAssembler(int rank,
		       Quadrature *quad2,
3680
3681
3682
3683
3684
		       Quadrature *quad1GrdPsi,
		       Quadrature *quad1GrdPhi,
		       Quadrature *quad0);


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

Thomas Witkowski's avatar
Thomas Witkowski committed
3688
3689
3690
    /// Evaluation of all terms in \ref zeroOrder. 
    void evalZeroOrder(int nPoints,
		       const double *uhAtQP,
3691
3692
3693
3694
3695
		       const WorldVector<double> *grdUhAtQP,
		       const WorldMatrix<double> *D2UhAtQP,
		       double *result,
		       double factor) const
    {
3696
3697
      int myRank = omp_get_thread_num();

3698
      std::vector<OperatorTerm*>::const_iterator termIt;
3699
3700
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
3701
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3702
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
3703
    }
3704

3705
    /// Evaluation of all terms in \ref firstOrderGrdPsi. 
Thomas Witkowski's avatar
Thomas Witkowski committed
3706
3707
    void evalFirstOrderGrdPsi(int nPoints,
			      const double *uhAtQP,
3708
3709
3710
3711
3712
			      const WorldVector<double> *grdUhAtQP,
			      const WorldMatrix<double> *D2UhAtQP,
			      double *result,
			      double factor) const
    {
3713
3714
      int myRank = omp_get_thread_num();

3715
      std::vector<OperatorTerm*>::const_iterator termIt;
3716
3717
      for (termIt = firstOrderGrdPsi[myRank].begin(); 
	   termIt != firstOrderGrdPsi[myRank].end(); 
3718
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3719
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
3720
    }
3721

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

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

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

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

3756
    /// Weak evaluation of all terms in \ref secondOrder. 
Thomas Witkowski's avatar
Thomas Witkowski committed
3757
    void weakEvalSecondOrder(int nPoints,
3758
3759
3760
			     const WorldVector<double> *grdUhAtQP,
			     WorldVector<double> *result) const
    {
3761
3762
      int myRank = omp_get_thread_num();

3763
      std::vector<OperatorTerm*>::const_iterator termIt;
3764
3765
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
3766
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3767
	static_cast<SecondOrderTerm*>(*termIt)->weakEval(nPoints, grdUhAtQP, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
3768
    }
3769
  
3770
    /// Calls getLALt() for each term in \ref secondOrder and adds the results to LALt.
Thomas Witkowski's avatar
Thomas Witkowski committed
3771
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
3772
    {
3773
3774
      int myRank = omp_get_thread_num();

3775
      std::vector<OperatorTerm*>::const_iterator termIt;
3776
3777
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
3778
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3779
	static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, nPoints, LALt);
Thomas Witkowski's avatar
Thomas Witkowski committed
3780
    }
3781
  
3782
3783
3784
3785
    /// 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
3786
    {
3787
3788
      int myRank = omp_get_thread_num();

3789
      std::vector<OperatorTerm*>::const_iterator termIt;
3790
3791
      for (termIt = firstOrderGrdPsi[myRank].begin(); 
	   termIt != firstOrderGrdPsi[myRank].end(); 
3792
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3793
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
Thomas Witkowski's avatar
Thomas Witkowski committed
3794
    }
3795

3796
3797
3798
3799
    /// 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
3800
    {
3801
3802
      int myRank = omp_get_thread_num();

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

3810
    /// Calls getC() for each term in \ref zeroOrder and adds the results to c.
3811
    void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &c) const
3812
    {
3813
3814
      int myRank = omp_get_thread_num();

3815
      std::vector<OperatorTerm*>::const_iterator termIt;
3816
3817
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
3818
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3819
	static_cast<ZeroOrderTerm*>(*termIt)->getC(elInfo, nPoints, c);
Thomas Witkowski's avatar
Thomas Witkowski committed
3820
    }
3821

Thomas Witkowski's avatar
Thomas Witkowski committed
3822
    /// Returns true, if there are second order terms. Returns false otherwise.
3823
3824
    inline bool secondOrderTerms() 
    {
3825
      return secondOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3826
    }
3827

3828
    /// Returns true, if there are first order terms (grdPsi). Returns false otherwise.
3829
3830
    inline bool firstOrderTermsGrdPsi() 
    {
3831
      return firstOrderGrdPsi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3832
    }
3833

3834
    /// Returns true, if there are first order terms (grdPhi). Returns false otherwise.
3835
3836
    inline bool firstOrderTermsGrdPhi() 
    {
3837
      return firstOrderGrdPhi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3838
    }
3839

3840
    /// Returns true, if there are zero order terms. Returns false otherwise.
3841
3842
    inline bool zeroOrderTerms() 
    {
3843
      return zeroOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3844
    }
3845
3846

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3847
    /// Constant type flag for matrix operators
3848
3849
    static const Flag MATRIX_OPERATOR;

Thomas Witkowski's avatar
Thomas Witkowski committed
3850
    /// Constant type flag for vector operators
3851
3852
3853
    static const Flag VECTOR_OPERATOR;

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
3854
    /// FiniteElemSpace for matrix rows and element vector
3855
3856
    const FiniteElemSpace *rowFESpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
3857
    /// FiniteElemSpace for matrix columns. Can be equal to \rowFESpace.
3858
3859
    const FiniteElemSpace *colFESpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
3860
3861
3862
3863
    /// 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
3864
3865
    int nRow;

Thomas Witkowski's avatar
Thomas Witkowski committed
3866
    /// Number of columns in the element matrix
3867
3868
    int nCol;

Thomas Witkowski's avatar
Thomas Witkowski committed
3869
    /// Type of this Operator.
3870
3871
    Flag type;

Thomas Witkowski's avatar
Thomas Witkowski committed
3872
    /// Flag for mesh traversal
3873
3874
    Flag fillFlag;

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

3878
3879
3880
3881
3882
    /** \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.
     */
3883
    std::vector<Assembler*> assembler;
3884

Thomas Witkowski's avatar
Thomas Witkowski committed
3885
    /// List of all second order terms
3886
    std::vector< std::vector<OperatorTerm*> > secondOrder;
3887

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3894
    /// List of all zero order terms
3895
    std::vector< std::vector<OperatorTerm*> > zeroOrder;
3896
3897
3898
3899
3900
3901
3902
3903
3904

    /** \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
3905
    /// Spezifies whether optimized assemblers are used or not.
3906
3907
3908
3909
3910
3911
3912
3913
3914
3915
3916
    bool optimized;

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

}

3917
3918
#include "Operator.hh"

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