Operator.h 103 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
3001
3002
    void eval(int nPoints,
	      const double *uhAtQP,
3003
3004
3005
3006
3007
3008
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3012
    /// Vector v at quadrature points. 
3013
3014
    double *vecAtQPs; 

Thomas Witkowski's avatar
Thomas Witkowski committed
3015
    /// First DOFVector whose gradient is evaluated at quadrature points. 
3016
3017
    DOFVectorBase<double>* vecGrd1; 

Thomas Witkowski's avatar
Thomas Witkowski committed
3018
    /// Gradient of first vector at quadrature points. 
3019
3020
    WorldVector<double> *grad1AtQPs; 

Thomas Witkowski's avatar
Thomas Witkowski committed
3021
    /// Second DOFVector whose gradient is evaluated at quadrature points. 
3022
3023
    DOFVectorBase<double>* vecGrd2; 

Thomas Witkowski's avatar
Thomas Witkowski committed
3024
    /// Gradient of second vector at quadrature points. 
3025
3026
    WorldVector<double> *grad2AtQPs; 

Thomas Witkowski's avatar
Thomas Witkowski committed
3027
    /// Function for c. 
3028
3029
3030
3031
3032
3033
3034
    TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *f; 
  }; 


  class VecOfDOFVecsAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3035
    /// Constructor. 
3036
    VecOfDOFVecsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3037
			 AbstractFunction<double, std::vector<double> > *f);
3038

Thomas Witkowski's avatar
Thomas Witkowski committed
3039
    /// Implementation of \ref OperatorTerm::initElement().
3040
3041
3042
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

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

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

  protected: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3055
    /// Vector of DOFVectors to be evaluated at quadrature points. 
3056
    std::vector<DOFVectorBase<double>*> vecs; 
3057

Thomas Witkowski's avatar
Thomas Witkowski committed
3058
    /// Vectors at quadrature points. 
3059
    std::vector<double*> vecsAtQPs; 
3060

Thomas Witkowski's avatar
Thomas Witkowski committed
3061
    /// Function for c. 
3062
    AbstractFunction<double, std::vector<double> > *f; 
3063
3064
3065
3066
3067
  }; 

  class VecOfGradientsAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3068
    /// Constructor. 
3069
    VecOfGradientsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3070
			   AbstractFunction<double, std::vector<WorldVector<double>*> > *af);
3071

Thomas Witkowski's avatar
Thomas Witkowski committed
3072
    /// Implementation of \ref OperatorTerm::initElement().
3073
3074
3075
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3079
    /// Implements ZeroOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
3080
3081
    void eval(int nPoints,
	      const double *uhAtQP,
3082
3083
3084
3085
3086
3087
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3088
    /// Vector of DOFVectors to be evaluated at quadrature points. 
3089
    std::vector<DOFVectorBase<double>*> vecs; 
3090

Thomas Witkowski's avatar
Thomas Witkowski committed
3091
    /// Vectors at quadrature points. 
3092
    std::vector<WorldVector<double>*> gradsAtQPs; 
3093

Thomas Witkowski's avatar
Thomas Witkowski committed
3094
    /// Function for c. 
3095
    AbstractFunction<double, std::vector<WorldVector<double>*> > *f; 
3096
3097
3098
3099
3100
3101
  };


  class VecDivergence_ZOT : public ZeroOrderTerm
  {
  public: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3102
3103
    /// Constructor. 
    VecDivergence_ZOT(int nComponents,
3104
3105
		      DOFVectorBase<double> *vec0,
		      DOFVectorBase<double> *vec1 = NULL,
Thomas Witkowski's avatar
Thomas Witkowski committed
3106
		      DOFVectorBase<double> *vec2 = NULL);
3107

Thomas Witkowski's avatar
Thomas Witkowski committed
3108
    /// Implementation of \ref OperatorTerm::initElement().
3109
3110
3111
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3115
    /// Implements ZeroOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
3116
3117
    void eval(int nPoints,
	      const double *uhAtQP,
3118
3119
3120
3121
3122
3123
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3124
    /// Vector of DOFVectors to be evaluated at quadrature points. 
3125
    std::vector<DOFVectorBase<double>*> vecs; 
3126

Thomas Witkowski's avatar
Thomas Witkowski committed
3127
    /// Vectors at quadrature points. 
3128
    std::vector<WorldVector<double>*> gradsAtQPs; 
3129
3130
3131
3132
3133
3134
3135
  };



  class VecAndVecOfGradientsAtQP_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3136
3137
    /// Constructor.
    VecAndVecOfGradientsAtQP_ZOT(DOFVector<double> *v,
3138
				 const std::vector<DOFVector<double>*>& dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
3139
				 BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *af);
3140

Thomas Witkowski's avatar
Thomas Witkowski committed
3141
    /// Implementation of \ref OperatorTerm::initElement().
3142
3143
3144
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3148
    /// Implements ZeroOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
3149
3150
    void eval(int nPoints,
	      const double *uhAtQP,
3151
3152
3153
3154
3155
3156
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
3157
    /// DOFVector to be evaluated at quadrature points.
3158
3159
    DOFVector<double>* vec;

Thomas Witkowski's avatar
Thomas Witkowski committed
3160
    /// Vector v at quadrature points.
3161
3162
    double *vecAtQPs;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3166
    /// Vectors at quadrature points.
3167
3168
    std::vector<WorldVector<double>*> gradsAtQPs;

Thomas Witkowski's avatar
Thomas Witkowski committed
3169
    /// Function for c.
3170
3171
3172
    BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *f;
  };

Thomas Witkowski's avatar
Thomas Witkowski committed
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
3212
3213
3214
3215
3216
3217
3218
3219
3220
3221
3222
3223
3224
3225
3226
3227
3228
3229
3230
3231
3232
3233
3234
3235
3236
3237
3238
3239
3240

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



3241
3242
3243
  class General_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3244
    /// Constructor.
3245
3246
    General_ZOT(std::vector<DOFVectorBase<double>*> vecs,
		std::vector<DOFVectorBase<double>*> grads,
Thomas Witkowski's avatar
Thomas Witkowski committed
3247
		TertiaryAbstractFunction<double, WorldVector<double>, std::vector<double>, std::vector<WorldVector<double> > > *f);
3248

Thomas Witkowski's avatar
Thomas Witkowski committed
3249
    /// Implementation of \ref OperatorTerm::initElement().
3250
3251
3252
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3256
    /// Implements ZeroOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
3257
3258
    void eval(int nPoints,
	      const double *uhAtQP,
3259
3260
3261
3262
3263
3264
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
3265
    std::vector<DOFVectorBase<double>*> vecs_; 
3266

3267
    std::vector<DOFVectorBase<double>*> grads_;
3268
3269
3270

    TertiaryAbstractFunction<double, 
			     WorldVector<double>,
3271
3272
			     std::vector<double>, 
			     std::vector<WorldVector<double> > > *f_;
3273
3274
3275

    WorldVector<double> *coordsAtQPs_;

3276
    std::vector<double*> vecsAtQPs_;
3277

3278
    std::vector<WorldVector<double>*> gradsAtQPs_;
3279
3280
3281
3282
3283
  };

  class GeneralParametric_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3284
    /// Constructor.
3285
3286
    GeneralParametric_ZOT(std::vector<DOFVectorBase<double>*> vecs,
			  std::vector<DOFVectorBase<double>*> grads,
3287
3288
3289
			  QuartAbstractFunction<double, 
			  WorldVector<double>,
			  WorldVector<double>,
3290
			  std::vector<double>, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3291
			  std::vector<WorldVector<double> > > *f);
3292

Thomas Witkowski's avatar
Thomas Witkowski committed
3293
    /// Implementation of \ref OperatorTerm::initElement().
3294
3295
3296
    void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
		     Quadrature *quad = NULL);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3300
    /// Implements ZeroOrderTerm::eval().
Thomas Witkowski's avatar
Thomas Witkowski committed
3301
3302
    void eval(int nPoints,
	      const double *uhAtQP,
3303
3304
3305
3306
3307
3308
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
3309
    std::vector<DOFVectorBase<double>*> vecs_; 
3310

3311
    std::vector<DOFVectorBase<double>*> grads_;
3312
3313
3314
3315

    QuartAbstractFunction<double, 
			  WorldVector<double>,
			  WorldVector<double>,
3316
3317
			  std::vector<double>, 
			  std::vector<WorldVector<double> > > *f_;
3318
3319
3320
3321
3322

    WorldVector<double> *coordsAtQPs_;

    WorldVector<double> elementNormal_;

3323
    std::vector<double*> vecsAtQPs_;
3324

3325
    std::vector<WorldVector<double>*> gradsAtQPs_;
3326
3327
3328
3329
3330
3331
3332
3333
3334
3335
3336
3337
3338
3339
3340
3341
3342
3343
3344
3345
3346
3347
3348
3349
3350
3351
3352
3353
3354
3355
3356
  };


  /** \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
3357
3358
	     const FiniteElemSpace *rowFESpace,
	     const FiniteElemSpace *colFESpace = NULL);
3359

3360
    /// Destructor.
Thomas Witkowski's avatar
Thomas Witkowski committed
3361
    virtual ~Operator() {}
3362

3363
    /// Sets \ref optimized.
3364
3365
    inline void useOptimizedAssembler(bool opt) 
    {
3366
      optimized = opt;
Thomas Witkowski's avatar
Thomas Witkowski committed
3367
    }
3368

3369
    /// Returns \ref optimized.
3370
3371
    inline bool isOptimized() 
    {
3372
      return optimized;
Thomas Witkowski's avatar
Thomas Witkowski committed
3373
    }
3374

3375
    /// Adds a SecondOrderTerm to the Operator
3376
3377
    template <typename T>
    void addSecondOrderTerm(T *term);
3378

3379
    /// Adds a FirstOrderTerm to the Operator
3380
    template <typename T>
3381
3382
3383
    void addFirstOrderTerm(T *term, FirstOrderType type = GRD_PHI);

    /// Adds a ZeroOrderTerm to the Operator
3384
3385
    template <typename T>
    void addZeroOrderTerm(T *term);
3386
3387
3388
3389
3390
3391
3392


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

3396
3397
    virtual void getElementMatrix(const ElInfo *rowElInfo,
				  const ElInfo *colElInfo,
3398
3399
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
3400
				  ElementMatrix& userMat,
3401
3402
				  double factor = 1.0);

3403
3404
3405
3406
3407
    /** \brief
     * Calculates the element vector for this ElInfo and adds it multiplied by
     * factor to userVec.
     */
    virtual void getElementVector(const ElInfo *elInfo, 
3408
				  ElementVector& userVec, 
3409
3410
				  double factor = 1.0);

Thomas Witkowski's avatar
Thomas Witkowski committed
3411
3412
3413
3414
    virtual void getElementVector(const ElInfo *mainElInfo, 
				  const ElInfo *auxElInfo,
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
3415
				  ElementVector& userVec,
Thomas Witkowski's avatar
Thomas Witkowski committed
3416
				  double factor = 1.0);
3417

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3421
    /// Returns \ref rowFESpace.
3422
3423
    inline const FiniteElemSpace *getRowFESpace() 
    { 
3424
      return rowFESpace; 
3425
    }
3426

Thomas Witkowski's avatar
Thomas Witkowski committed
3427
    /// Returns \ref colFESpace.
3428
3429
    inline const FiniteElemSpace *getColFESpace() 
    { 
3430
      return colFESpace; 
3431
    }
3432

Thomas Witkowski's avatar
Thomas Witkowski committed
3433
    /// Returns \ref auxFESpaces.
3434
3435
    inline std::vector<const FiniteElemSpace*> getAuxFESpaces()
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
3436
3437
3438
      return auxFESpaces;
    }

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

3442
    /// Returns \ref uhOld.
3443
3444
    inline const DOFVectorBase<double> *getUhOld() 
    {
3445
      return uhOld;
3446
    }    
3447

Thomas Witkowski's avatar
Thomas Witkowski committed
3448
    /// Returns \ref assembler
3449
    Assembler *getAssembler(int rank);
3450

Thomas Witkowski's avatar
Thomas Witkowski committed
3451
    /// Sets \ref assembler
3452
    void setAssembler(int rank, Assembler *ass);
3453

Thomas Witkowski's avatar
Thomas Witkowski committed
3454
    /// Returns whether this is a matrix operator.
3455
3456
    inline const bool isMatrixOperator() 
    {
3457
      return type.isSet(MATRIX_OPERATOR);
3458
    }
3459

Thomas Witkowski's avatar
Thomas Witkowski committed
3460
    /// Returns whether this is a vector operator
3461
3462
    inline const bool isVectorOperator() 
    {
3463
      return type.isSet(VECTOR_OPERATOR);
3464
    }
3465

Thomas Witkowski's avatar
Thomas Witkowski committed
3466
    /// Sets \ref fillFlag, the flag used for this Operator while mesh traversal.
3467
3468
    inline void setFillFlag(Flag f) 
    { 
3469
      fillFlag = f; 
3470
    }
3471

Thomas Witkowski's avatar
Thomas Witkowski committed
3472
    /// Sets \ref needDualTraverse.
3473
3474
    void setNeedDualTraverse(bool b) 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
3475
3476
3477
3478
      needDualTraverse = b;
    }

    /// Returns \ref fillFlag
3479
3480
    inline Flag getFillFlag() 
    { 
3481
      return fillFlag; 
3482
    }
3483

Thomas Witkowski's avatar
Thomas Witkowski committed
3484
    /// Returns \ref needDualTraverse
3485
3486
    bool getNeedDualTraverse() 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
3487
3488
3489
3490
      return needDualTraverse;
    }

    /// Initialization of the needed SubAssemblers using the given quadratures. 
3491
3492
    void initAssembler(int rank,
		       Quadrature *quad2,
3493
3494
3495
3496
3497
		       Quadrature *quad1GrdPsi,
		       Quadrature *quad1GrdPhi,
		       Quadrature *quad0);


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

Thomas Witkowski's avatar
Thomas Witkowski committed
3501
3502
3503
    /// Evaluation of all terms in \ref zeroOrder. 
    void evalZeroOrder(int nPoints,
		       const double *uhAtQP,
3504
3505
3506
3507
3508
		       const WorldVector<double> *grdUhAtQP,
		       const WorldMatrix<double> *D2UhAtQP,
		       double *result,
		       double factor) const
    {
3509
3510
      int myRank = omp_get_thread_num();

3511
      std::vector<OperatorTerm*>::const_iterator termIt;
3512
3513
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
3514
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3515
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
3516
    }
3517

3518
    /// Evaluation of all terms in \ref firstOrderGrdPsi. 
Thomas Witkowski's avatar
Thomas Witkowski committed
3519
3520
    void evalFirstOrderGrdPsi(int nPoints,
			      const double *uhAtQP,
3521
3522
3523
3524
3525
			      const WorldVector<double> *grdUhAtQP,
			      const WorldMatrix<double> *D2UhAtQP,
			      double *result,
			      double factor) const
    {
3526
3527
      int myRank = omp_get_thread_num();

3528
      std::vector<OperatorTerm*>::const_iterator termIt;
3529
3530
      for (termIt = firstOrderGrdPsi[myRank].begin(); 
	   termIt != firstOrderGrdPsi[myRank].end(); 
3531
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3532
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
3533
    }
3534

3535
    /// Evaluation of all terms in \ref firstOrderGrdPhi. 
Thomas Witkowski's avatar
Thomas Witkowski committed
3536
3537
    void evalFirstOrderGrdPhi(int nPoints,
			      const double *uhAtQP,
3538
3539
3540
3541
3542
			      const WorldVector<double> *grdUhAtQP,
			      const WorldMatrix<double> *D2UhAtQP,
			      double *result,
			      double factor) const
    {
3543
3544
      int myRank = omp_get_thread_num();

3545
      std::vector<OperatorTerm*>::const_iterator termIt;
3546
3547
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
3548
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3549
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
3550
    }
3551

3552
    /// Evaluation of all terms in \ref secondOrder. 
Thomas Witkowski's avatar
Thomas Witkowski committed
3553
3554
    void evalSecondOrder(int nPoints,
			 const double *uhAtQP,
3555
3556
3557
3558
3559
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
    {
3560
3561
      int myRank = omp_get_thread_num();

3562
      std::vector<OperatorTerm*>::const_iterator termIt;
3563
3564
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
3565
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3566
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
Thomas Witkowski's avatar
Thomas Witkowski committed
3567
    }
3568

3569
    /// Weak evaluation of all terms in \ref secondOrder. 
Thomas Witkowski's avatar
Thomas Witkowski committed
3570
    void weakEvalSecondOrder(int nPoints,
3571
3572
3573
			     const WorldVector<double> *grdUhAtQP,
			     WorldVector<double> *result) const
    {
3574
3575
      int myRank = omp_get_thread_num();

3576
      std::vector<OperatorTerm*>::const_iterator termIt;
3577
3578
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
3579
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3580
	static_cast<SecondOrderTerm*>(*termIt)->weakEval(nPoints, grdUhAtQP, result);
Thomas Witkowski's avatar
Thomas Witkowski committed
3581
    }
3582
  
3583
    /// Calls getLALt() for each term in \ref secondOrder and adds the results to LALt.
Thomas Witkowski's avatar
Thomas Witkowski committed
3584
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
3585
    {
3586
3587
      int myRank = omp_get_thread_num();

3588
      std::vector<OperatorTerm*>::const_iterator termIt;
3589
3590
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
3591
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3592
	static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, nPoints, LALt);
Thomas Witkowski's avatar
Thomas Witkowski committed
3593
    }
3594
  
3595
3596
3597
3598
    /// 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
3599
    {
3600
3601
      int myRank = omp_get_thread_num();

3602
      std::vector<OperatorTerm*>::const_iterator termIt;
3603
3604
      for (termIt = firstOrderGrdPsi[myRank].begin(); 
	   termIt != firstOrderGrdPsi[myRank].end(); 
3605
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3606
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
Thomas Witkowski's avatar
Thomas Witkowski committed
3607
    }
3608

3609
3610
3611
3612
    /// 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
3613
    {
3614
3615
      int myRank = omp_get_thread_num();

3616
      std::vector<OperatorTerm*>::const_iterator termIt;
3617
3618
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
3619
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3620
	static_cast<FirstOrderTerm*>(*termIt)->getLb(elInfo, nPoints, Lb);
Thomas Witkowski's avatar
Thomas Witkowski committed
3621
    }
3622

3623
    /// Calls getC() for each term in \ref zeroOrder and adds the results to c.
3624
    void getC(const ElInfo *elInfo, int nPoints, std::vector<double> &c) const
3625
    {
3626
3627
      int myRank = omp_get_thread_num();

3628
      std::vector<OperatorTerm*>::const_iterator termIt;
3629
3630
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
3631
	   ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
3632
	static_cast<ZeroOrderTerm*>(*termIt)->getC(elInfo, nPoints, c);
Thomas Witkowski's avatar
Thomas Witkowski committed
3633
    }
3634

Thomas Witkowski's avatar
Thomas Witkowski committed
3635
    /// Returns true, if there are second order terms. Returns false otherwise.
3636
3637
    inline bool secondOrderTerms() 
    {
3638
      return secondOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3639
    }
3640

3641
    /// Returns true, if there are first order terms (grdPsi). Returns false otherwise.
3642
3643
    inline bool firstOrderTermsGrdPsi() 
    {
3644
      return firstOrderGrdPsi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3645
    }
3646

3647
    /// Returns true, if there are first order terms (grdPhi). Returns false otherwise.
3648
3649
    inline bool firstOrderTermsGrdPhi() 
    {
3650
      return firstOrderGrdPhi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3651
    }
3652

3653
    /// Returns true, if there are zero order terms. Returns false otherwise.
3654
3655
    inline bool zeroOrderTerms() 
    {
3656
      return zeroOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3657
    }
3658
3659

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3660
    /// Constant type flag for matrix operators
3661
3662
    static const Flag MATRIX_OPERATOR;

Thomas Witkowski's avatar
Thomas Witkowski committed
3663
    /// Constant type flag for vector operators
3664
3665
3666
    static const Flag VECTOR_OPERATOR;

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
3667
    /// FiniteElemSpace for matrix rows and element vector
3668
3669
    const FiniteElemSpace *rowFESpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
3670
    /// FiniteElemSpace for matrix columns. Can be equal to \rowFESpace.
3671
3672
    const FiniteElemSpace *colFESpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
3673
3674
3675
3676
    /// 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
3677
3678
    int nRow;

Thomas Witkowski's avatar
Thomas Witkowski committed
3679
    /// Number of columns in the element matrix
3680
3681
    int nCol;

Thomas Witkowski's avatar
Thomas Witkowski committed
3682
    /// Type of this Operator.
3683
3684
    Flag type;

Thomas Witkowski's avatar
Thomas Witkowski committed
3685
    /// Flag for mesh traversal
3686
3687
    Flag fillFlag;

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

3691
3692
3693
3694
3695
    /** \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.
     */
3696
    std::vector<Assembler*> assembler;
3697

Thomas Witkowski's avatar
Thomas Witkowski committed
3698
    /// List of all second order terms
3699
    std::vector< std::vector<OperatorTerm*> > secondOrder;
3700

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3707
    /// List of all zero order terms
3708
    std::vector< std::vector<OperatorTerm*> > zeroOrder;
3709
3710
3711
3712
3713
3714
3715
3716
3717

    /** \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
3718
    /// Spezifies whether optimized assemblers are used or not.
3719
3720
3721
3722
3723
3724
3725
3726
3727
3728
3729
    bool optimized;

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

}

3730
3731
#include "Operator.hh"

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