~lukas-korous/hermes2d/master

« back to all changes in this revision

Viewing changes to hermes2d/test_examples/03-navier-stokes/definitions.cpp

  • Committer: Lukas Korous
  • Date: 2015-04-01 15:42:57 UTC
  • Revision ID: git-v1:212e78cb7fb92bfecaa4b337ebaf7d38b35e1982
Make stuff work and test drive Agros test-flow example.

Show diffs side-by-side

added added

removed removed

Lines of Context:
71
71
 
72
72
    double value(int n, double *wt, Func<double> *u_ext[], Func<double> *u, Func<double> *v, GeomVol<double> *e, Func<double> **ext) const{
73
73
      double result = int_grad_u_grad_v<double, double>(n, wt, u, v) / Reynolds;
74
 
      if(!Stokes)
75
 
        result += int_u_v<double, double>(n, wt, u, v) / time_step;
76
74
      return result;
77
75
    }
78
76
 
79
77
    Ord ord(int n, double *wt, Func<Ord> *u_ext[], Func<Ord> *u, Func<Ord> *v, GeomVol<Ord> *e, Func<Ord> **ext) const
80
78
    {
81
79
      Ord result = int_grad_u_grad_v<Ord, Ord>(n, wt, u, v) / Reynolds;
82
 
      if(!Stokes)
83
 
        result += int_u_v<Ord, Ord>(n, wt, u, v) / time_step;
84
80
      return result;
85
81
    }
86
82
 
185
181
 
186
182
    double value(int n, double *wt, Func<double> *u_ext[], Func<double> *v, GeomVol<double> *e, Func<double> **ext) const{
187
183
      double result = 0;
188
 
      if(!Stokes) {
189
 
        Func<double>* vel_prev_time = ext[2]; // this form is used with both velocity components
190
 
        result = int_u_v<double, double>(n, wt, vel_prev_time, v) / time_step;
191
 
      }
192
184
      return result;
193
185
    }
194
186
 
195
187
    Ord ord(int n, double *wt, Func<Ord> *u_ext[], Func<Ord> *v, GeomVol<Ord> *e, Func<Ord> **ext) const {
196
188
      Ord result = Ord(0);
197
 
      if(!Stokes) {
198
 
        Func<Ord>* vel_prev_time = ext[2]; // this form is used with both velocity components
199
 
        result = int_u_v<Ord, Ord>(n, wt, vel_prev_time, v) / time_step;
200
 
      }
201
189
      return result;
202
190
    }
203
191
    
268
256
 
269
257
    double value(int n, double *wt, Func<double> *u_ext[], Func<double> *u, Func<double> *v, GeomVol<double> *e, Func<double> **ext) const{
270
258
      double result = int_grad_u_grad_v<double, double>(n, wt, u, v) / Reynolds;
271
 
      if(!Stokes)
272
 
        result += int_u_v<double, double>(n, wt, u, v) / time_step;
273
259
      return result;
274
260
    }
275
261
 
510
496
        result += wt[i] * ((xvel_prev_newton->dx[i] * v->dx[i] + xvel_prev_newton->dy[i] * v->dy[i]) / Reynolds - (p_prev_newton->val[i] * v->dx[i]));
511
497
      if(!Stokes)
512
498
        for (int i = 0; i < n; i++)
513
 
          result += wt[i] * (((xvel_prev_newton->val[i] - xvel_prev_time->val[i]) * v->val[i] / time_step )
514
 
          + ((xvel_prev_newton->val[i] * xvel_prev_newton->dx[i] + yvel_prev_newton->val[i] * xvel_prev_newton->dy[i]) * v->val[i]));
 
499
          result += wt[i] * (xvel_prev_newton->val[i] * xvel_prev_newton->dx[i] + yvel_prev_newton->val[i] * xvel_prev_newton->dy[i]) * v->val[i];
515
500
      return result;
516
501
    }
517
502
 
526
511
        result += wt[i] * ((xvel_prev_newton->dx[i] * v->dx[i] + xvel_prev_newton->dy[i] * v->dy[i]) / Reynolds - (p_prev_newton->val[i] * v->dx[i]));
527
512
      if(!Stokes)
528
513
        for (int i = 0; i < n; i++)
529
 
          result += wt[i] * (((xvel_prev_newton->val[i] - xvel_prev_time->val[i]) * v->val[i] / time_step)
530
 
          + ((xvel_prev_newton->val[i] * xvel_prev_newton->dx[i] + yvel_prev_newton->val[i] * xvel_prev_newton->dy[i]) * v->val[i]));
 
514
          result += wt[i] * (xvel_prev_newton->val[i] * xvel_prev_newton->dx[i] + yvel_prev_newton->val[i] * xvel_prev_newton->dy[i]) * v->val[i];
531
515
      return result;
532
516
    }
533
517
    
559
543
        result += wt[i] * ((yvel_prev_newton->dx[i] * v->dx[i] + yvel_prev_newton->dy[i] * v->dy[i]) / Reynolds - (p_prev_newton->val[i] * v->dy[i]));
560
544
      if(!Stokes)
561
545
        for (int i = 0; i < n; i++)
562
 
          result += wt[i] * (((yvel_prev_newton->val[i] - yvel_prev_time->val[i]) * v->val[i] / time_step )
563
 
          + ((xvel_prev_newton->val[i] * yvel_prev_newton->dx[i] + yvel_prev_newton->val[i] * yvel_prev_newton->dy[i]) * v->val[i]));
 
546
          result += wt[i] * (xvel_prev_newton->val[i] * yvel_prev_newton->dx[i] + yvel_prev_newton->val[i] * yvel_prev_newton->dy[i]) * v->val[i];
564
547
      return result;
565
548
    }
566
549
 
648
631
      };
649
632
 
650
633
      virtual double value(double x, double y) const {
651
 
        double val_y = vel_inlet * y*(H-y) / (H/2.)/(H/2.);
652
 
        if (current_time <= startup_time) 
653
 
          return val_y * current_time/startup_time;
654
 
        else 
655
 
          return val_y;
 
634
        double val_y = -x*(1-x);
 
635
        return val_y / 4.;
656
636
      };
657
637
 
658
638
protected: