Operator.h 96.8 KB
Newer Older
3001
3002
3003
3004
3005
3006
3007
3008
     */ 
    TertiaryAbstractFunction<double, double, WorldVector<double>, WorldVector<double> > *f; 
  }; 


  class VecOfDOFVecsAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3009
    /// Constructor. 
3010
    VecOfDOFVecsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3011
			 AbstractFunction<double, std::vector<double> > *f);
3012
3013
3014
3015
3016
3017
3018
3019
3020
3021

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3022
    void getC(const ElInfo *, int nPoints, double *C) const;
3023
3024
3025
3026

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3027
3028
    void eval(int nPoints,
	      const double *uhAtQP,
3029
3030
3031
3032
3033
3034
3035
3036
3037
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
    /** \brief 
     * Vector of DOFVectors to be evaluated at quadrature points. 
     */ 
3038
    std::vector<DOFVectorBase<double>*> vecs; 
3039
3040
3041
3042

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3043
    std::vector<double*> vecsAtQPs; 
3044
3045
3046
3047

    /** \brief 
     * Function for c. 
     */ 
3048
    AbstractFunction<double, std::vector<double> > *f; 
3049
3050
3051
3052
3053
  }; 

  class VecOfGradientsAtQP_ZOT : public ZeroOrderTerm 
  { 
  public: 
Thomas Witkowski's avatar
Thomas Witkowski committed
3054
    /// Constructor. 
3055
    VecOfGradientsAtQP_ZOT(const std::vector<DOFVectorBase<double>*>& dv, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3056
			   AbstractFunction<double, std::vector<WorldVector<double>*> > *af);
3057
3058
3059
3060
3061
3062
3063
3064
3065
3066

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3067
    void getC(const ElInfo *, int nPoints, double *C) const;
3068
3069
3070
3071

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3072
3073
    void eval(int nPoints,
	      const double *uhAtQP,
3074
3075
3076
3077
3078
3079
3080
3081
3082
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
    /** \brief 
     * Vector of DOFVectors to be evaluated at quadrature points. 
     */ 
3083
    std::vector<DOFVectorBase<double>*> vecs; 
3084
3085
3086
3087

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3088
    std::vector<WorldVector<double>*> gradsAtQPs; 
3089
3090
3091
3092

    /** \brief 
     * Function for c. 
     */ 
3093
    AbstractFunction<double, std::vector<WorldVector<double>*> > *f; 
3094
3095
3096
3097
3098
3099
  };


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

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3115
    void getC(const ElInfo *, int nPoints, double *C) const;
3116
3117
3118
3119

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3120
3121
    void eval(int nPoints,
	      const double *uhAtQP,
3122
3123
3124
3125
3126
3127
3128
3129
3130
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected: 
    /** \brief 
     * Vector of DOFVectors to be evaluated at quadrature points. 
     */ 
3131
    std::vector<DOFVectorBase<double>*> vecs; 
3132
3133
3134
3135

    /** \brief 
     * Vectors at quadrature points. 
     */ 
3136
    std::vector<WorldVector<double>*> gradsAtQPs; 
3137
3138
3139
3140
3141
3142
3143
  };



  class VecAndVecOfGradientsAtQP_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3144
3145
    /// Constructor.
    VecAndVecOfGradientsAtQP_ZOT(DOFVector<double> *v,
3146
				 const std::vector<DOFVector<double>*>& dv,
Thomas Witkowski's avatar
Thomas Witkowski committed
3147
				 BinaryAbstractFunction<double, double, std::vector<WorldVector<double>*> > *af);
3148
3149
3150
3151
3152
3153
3154
3155
3156
3157

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3158
    void getC(const ElInfo *, int nPoints, double *C) const;
3159
3160
3161
3162

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3163
3164
    void eval(int nPoints,
	      const double *uhAtQP,
3165
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
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

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

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

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

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

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

  class General_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3200
    /// Constructor.
3201
3202
    General_ZOT(std::vector<DOFVectorBase<double>*> vecs,
		std::vector<DOFVectorBase<double>*> grads,
Thomas Witkowski's avatar
Thomas Witkowski committed
3203
		TertiaryAbstractFunction<double, WorldVector<double>, std::vector<double>, std::vector<WorldVector<double> > > *f);
3204
3205
3206
3207
3208
3209
3210
3211
3212
3213

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3214
    void getC(const ElInfo *, int nPoints, double *C) const;
3215
3216
3217
3218

    /** \brief
     * Implements ZeroOrderTerm::eval().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3219
3220
    void eval(int nPoints,
	      const double *uhAtQP,
3221
3222
3223
3224
3225
3226
	      const WorldVector<double> *grdUhAtQP,
	      const WorldMatrix<double> *D2UhAtQP,
	      double *result,
	      double fac) const;

  protected:
3227
    std::vector<DOFVectorBase<double>*> vecs_; 
3228

3229
    std::vector<DOFVectorBase<double>*> grads_;
3230
3231
3232

    TertiaryAbstractFunction<double, 
			     WorldVector<double>,
3233
3234
			     std::vector<double>, 
			     std::vector<WorldVector<double> > > *f_;
3235
3236
3237

    WorldVector<double> *coordsAtQPs_;

3238
    std::vector<double*> vecsAtQPs_;
3239

3240
    std::vector<WorldVector<double>*> gradsAtQPs_;
3241
3242
3243
3244
3245
  };

  class GeneralParametric_ZOT : public ZeroOrderTerm
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3246
    /// Constructor.
3247
3248
    GeneralParametric_ZOT(std::vector<DOFVectorBase<double>*> vecs,
			  std::vector<DOFVectorBase<double>*> grads,
3249
3250
3251
			  QuartAbstractFunction<double, 
			  WorldVector<double>,
			  WorldVector<double>,
3252
			  std::vector<double>, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3253
			  std::vector<WorldVector<double> > > *f);
3254
3255
3256
3257
3258
3259
3260
3261
3262
3263

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

    /** \brief
     * Implements ZeroOrderTerm::getC().
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3264
    void getC(const ElInfo *, int nPoints, double *C) const;
3265
3266
3267
3268

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

  protected:
3277
    std::vector<DOFVectorBase<double>*> vecs_; 
3278

3279
    std::vector<DOFVectorBase<double>*> grads_;
3280
3281
3282
3283

    QuartAbstractFunction<double, 
			  WorldVector<double>,
			  WorldVector<double>,
3284
3285
			  std::vector<double>, 
			  std::vector<WorldVector<double> > > *f_;
3286
3287
3288
3289
3290

    WorldVector<double> *coordsAtQPs_;

    WorldVector<double> elementNormal_;

3291
    std::vector<double*> vecsAtQPs_;
3292

3293
    std::vector<WorldVector<double>*> gradsAtQPs_;
3294
3295
3296
3297
3298
3299
3300
3301
3302
3303
3304
3305
3306
3307
3308
3309
3310
3311
3312
3313
3314
3315
3316
3317
3318
3319
3320
3321
3322
3323
3324
3325
3326
3327
3328
3329
3330
  };


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

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

    /** \brief
     * Constructs an empty Operator of type operatorType for the given 
     * FiniteElemSpace. 
     * The type is given by a Flag that can contain the values MATRIX_OPERATOR,
     * VECTOR_OPERATOR, or MATRIX_OPERATOR | VECTOR_OPERATOR. This type specifies 
     * whether the Operator is used on the left hand side, the right hand side,
     * or on both sides of the system. 
     */
    Operator(Flag operatorType,
Thomas Witkowski's avatar
Thomas Witkowski committed
3331
3332
	     const FiniteElemSpace *rowFESpace,
	     const FiniteElemSpace *colFESpace = NULL);
3333
3334
3335
3336

    /** \brief
     * Destructor.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3337
    virtual ~Operator() {}
3338
3339
3340
3341
3342
3343

    /** \brief
     * Sets \ref optimized.
     */
    inline void useOptimizedAssembler(bool opt) {
      optimized = opt;
Thomas Witkowski's avatar
Thomas Witkowski committed
3344
    }
3345
3346
3347
3348
3349
3350

    /** \brief
     * Returns \ref optimized.
     */
    inline bool isOptimized() {
      return optimized;
Thomas Witkowski's avatar
Thomas Witkowski committed
3351
    }
3352
3353
3354
3355

    /** \brief
     * Adds a SecondOrderTerm to the Operator
     */
3356
3357
    template <typename T>
    void addSecondOrderTerm(T *term);
3358
3359
3360
3361

    /** \brief
     * Adds a FirstOrderTerm to the Operator
     */
3362
3363
    template <typename T>
    void addFirstOrderTerm(T *term, 
3364
3365
3366
3367
			   FirstOrderType type = GRD_PHI);
    /** \brief
     * Adds a ZeroOrderTerm to the Operator
     */
3368
3369
    template <typename T>
    void addZeroOrderTerm(T *term);
3370
3371
3372
3373
3374
3375
3376
3377
3378
3379


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

3380
3381
    virtual void getElementMatrix(const ElInfo *rowElInfo,
				  const ElInfo *colElInfo,
3382
3383
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
3384
3385
3386
				  ElementMatrix *userMat,
				  double factor = 1.0);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3395
3396
3397
3398
3399
3400
    virtual void getElementVector(const ElInfo *mainElInfo, 
				  const ElInfo *auxElInfo,
				  const ElInfo *smallElInfo,
				  const ElInfo *largeElInfo,
				  ElementVector *userVec,
				  double factor = 1.0);
3401
3402


Thomas Witkowski's avatar
Thomas Witkowski committed
3403
    /// Returns \ref rowFESpace.
3404
3405
    inline const FiniteElemSpace *getRowFESpace() { 
      return rowFESpace; 
3406
    }
3407

Thomas Witkowski's avatar
Thomas Witkowski committed
3408
    /// Returns \ref colFESpace.
3409
3410
    inline const FiniteElemSpace *getColFESpace() { 
      return colFESpace; 
3411
    }
3412

Thomas Witkowski's avatar
Thomas Witkowski committed
3413
3414
3415
3416
3417
    /// Returns \ref auxFESpaces.
    inline std::vector<const FiniteElemSpace*> getAuxFESpaces() {
      return auxFESpaces;
    }

3418
3419
3420
3421
3422
3423
3424
3425
3426
3427
    /** \brief
     * Sets \ref uhOld.
     */
    void setUhOld(const DOFVectorBase<double> *uhOld);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3430
    /// Returns \ref assembler
3431
    Assembler *getAssembler(int rank);
3432

Thomas Witkowski's avatar
Thomas Witkowski committed
3433
    /// Sets \ref assembler
3434
    void setAssembler(int rank, Assembler *ass);
3435

Thomas Witkowski's avatar
Thomas Witkowski committed
3436
    /// Returns whether this is a matrix operator.
3437
3438
    inline const bool isMatrixOperator() {
      return type.isSet(MATRIX_OPERATOR);
3439
    }
3440

Thomas Witkowski's avatar
Thomas Witkowski committed
3441
    /// Returns whether this is a vector operator
3442
3443
    inline const bool isVectorOperator() {
      return type.isSet(VECTOR_OPERATOR);
3444
    }
3445

Thomas Witkowski's avatar
Thomas Witkowski committed
3446
    /// Sets \ref fillFlag, the flag used for this Operator while mesh traversal.
3447
3448
    inline void setFillFlag(Flag f) { 
      fillFlag = f; 
3449
    }
3450

Thomas Witkowski's avatar
Thomas Witkowski committed
3451
3452
3453
3454
3455
3456
    /// Sets \ref needDualTraverse.
    void setNeedDualTraverse(bool b) {
      needDualTraverse = b;
    }

    /// Returns \ref fillFlag
3457
3458
    inline Flag getFillFlag() { 
      return fillFlag; 
3459
    }
3460

Thomas Witkowski's avatar
Thomas Witkowski committed
3461
3462
3463
3464
3465
3466
    /// Returns \ref needDualTraverse
    bool getNeedDualTraverse() {
      return needDualTraverse;
    }

    /// Initialization of the needed SubAssemblers using the given quadratures. 
3467
3468
    void initAssembler(int rank,
		       Quadrature *quad2,
3469
3470
3471
3472
3473
		       Quadrature *quad1GrdPsi,
		       Quadrature *quad1GrdPhi,
		       Quadrature *quad0);


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

Thomas Witkowski's avatar
Thomas Witkowski committed
3477
3478
3479
    /// Evaluation of all terms in \ref zeroOrder. 
    void evalZeroOrder(int nPoints,
		       const double *uhAtQP,
3480
3481
3482
3483
3484
		       const WorldVector<double> *grdUhAtQP,
		       const WorldMatrix<double> *D2UhAtQP,
		       double *result,
		       double factor) const
    {
3485
3486
      int myRank = omp_get_thread_num();

3487
      std::vector<OperatorTerm*>::const_iterator termIt;
3488
3489
3490
      for (termIt = zeroOrder[myRank].begin(); 
	   termIt != zeroOrder[myRank].end(); 
	   ++termIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3491
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
3492
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
3493
    }
3494
3495
3496
3497
3498


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

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

    /** \brief
     * Evaluation of all terms in \ref firstOrderGrdPhi. 
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3519
3520
    void evalFirstOrderGrdPhi(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
3531
      for (termIt = firstOrderGrdPhi[myRank].begin(); 
	   termIt != firstOrderGrdPhi[myRank].end(); 
	   ++termIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3532
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
3533
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
3534
    }
3535
3536
3537
3538
3539


    /** \brief
     * Evaluation of all terms in \ref secondOrder. 
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3540
3541
    void evalSecondOrder(int nPoints,
			 const double *uhAtQP,
3542
3543
3544
3545
3546
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
    {
3547
3548
      int myRank = omp_get_thread_num();

3549
      std::vector<OperatorTerm*>::const_iterator termIt;
3550
3551
3552
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
	   ++termIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3553
	(*termIt)->eval(nPoints, uhAtQP, grdUhAtQP, D2UhAtQP, result, factor);
3554
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
3555
    }
3556
3557
3558
3559

    /** \brief
     * Weak evaluation of all terms in \ref secondOrder. 
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3560
    void weakEvalSecondOrder(int nPoints,
3561
3562
3563
			     const WorldVector<double> *grdUhAtQP,
			     WorldVector<double> *result) const
    {
3564
3565
      int myRank = omp_get_thread_num();

3566
      std::vector<OperatorTerm*>::const_iterator termIt;
3567
3568
3569
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
	   ++termIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3570
	static_cast<SecondOrderTerm*>(*termIt)->weakEval(nPoints, grdUhAtQP, result);
3571
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
3572
    }
3573
3574
3575
3576
3577
  
    /** \brief
     * Calls getLALt() for each term in \ref secondOrder 
     * and adds the results to LALt.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3578
    void getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const
3579
    {
3580
3581
      int myRank = omp_get_thread_num();

3582
      std::vector<OperatorTerm*>::const_iterator termIt;
3583
3584
3585
      for (termIt = secondOrder[myRank].begin(); 
	   termIt != secondOrder[myRank].end(); 
	   ++termIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3586
	static_cast<SecondOrderTerm*>(*termIt)->getLALt(elInfo, nPoints, LALt);
3587
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
3588
    }
3589
3590
3591
3592
3593
  
    /** \brief
     * Calls getLb() for each term in \ref firstOrderGrdPsi 
     * and adds the results to Lb.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3594
    void getLbGrdPsi(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const
3595
    {
3596
3597
      int myRank = omp_get_thread_num();

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

    /** \brief
     * Calls getLb() for each term in \ref firstOrderGrdPhi 
     * and adds the results to Lb.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3610
    void getLbGrdPhi(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const
3611
    {
3612
3613
      int myRank = omp_get_thread_num();

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

    /** \brief
     * Calls getC() for each term in \ref zeroOrder
     * and adds the results to c.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
3626
    void getC(const ElInfo *elInfo, int nPoints, double *c) const
3627
    {
3628
3629
      int myRank = omp_get_thread_num();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3638
    /// Returns true, if there are second order terms. Returns false otherwise.
3639
    inline bool secondOrderTerms() {
3640
      return secondOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3641
    }
3642
3643
3644
3645
3646
3647

    /** \brief
     * Returns true, if there are first order terms (grdPsi). 
     * Returns false otherwise.
     */
    inline bool firstOrderTermsGrdPsi() {
3648
      return firstOrderGrdPsi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3649
    }
3650
3651
3652
3653
3654
3655

    /** \brief
     * Returns true, if there are first order terms (grdPhi). 
     * Returns false otherwise.
     */
    inline bool firstOrderTermsGrdPhi() {
3656
      return firstOrderGrdPhi[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3657
    }
3658
3659
3660
3661
3662
3663

    /** \brief
     * Returns true, if there are zero order terms.
     * Returns false otherwise.
     */
    inline bool zeroOrderTerms() {
3664
      return zeroOrder[omp_get_thread_num()].size() != 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
3665
    }
3666
3667

  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
3668
    /// Constant type flag for matrix operators
3669
3670
    static const Flag MATRIX_OPERATOR;

Thomas Witkowski's avatar
Thomas Witkowski committed
3671
    /// Constant type flag for vector operators
3672
3673
3674
    static const Flag VECTOR_OPERATOR;

  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
3675
    /// FiniteElemSpace for matrix rows and element vector
3676
3677
    const FiniteElemSpace *rowFESpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
3678
    /// FiniteElemSpace for matrix columns. Can be equal to \rowFESpace.
3679
3680
    const FiniteElemSpace *colFESpace;

Thomas Witkowski's avatar
Thomas Witkowski committed
3681
3682
3683
3684
    /// 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
3685
3686
    int nRow;

Thomas Witkowski's avatar
Thomas Witkowski committed
3687
    /// Number of columns in the element matrix
3688
3689
    int nCol;

Thomas Witkowski's avatar
Thomas Witkowski committed
3690
    /// Type of this Operator.
3691
3692
    Flag type;

Thomas Witkowski's avatar
Thomas Witkowski committed
3693
    /// Flag for mesh traversal
3694
3695
    Flag fillFlag;

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

3699
3700
3701
3702
3703
    /** \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.
     */
3704
    std::vector<Assembler*> assembler;
3705

Thomas Witkowski's avatar
Thomas Witkowski committed
3706
    /// List of all second order terms
3707
    std::vector< std::vector<OperatorTerm*> > secondOrder;
3708

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3715
    /// List of all zero order terms
3716
    std::vector< std::vector<OperatorTerm*> > zeroOrder;
3717
3718
3719
3720
3721
3722
3723
3724
3725

    /** \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
3726
    /// Spezifies whether optimized assemblers are used or not.
3727
3728
3729
3730
3731
3732
3733
3734
3735
3736
3737
    bool optimized;

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

}

3738
3739
#include "Operator.hh"

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