oldroyd_theta_scheme1.h

Go to the documentation of this file.

template<class P>
oldroyd_theta_scheme<P>::oldroyd_theta_scheme()
: We(0), alpha(8./9), a(1), Re(1), delta_t(0.025), tol(1e-6), max_iter(500),
  Th(), Xh(), Qh(), b(), c(), d(), mt(), inv_mt(), mu(), mp(), th(), thb(),
  theta(), lambda(), eta(), nu(), c1(), c2(), c3(), c4(), c5(), stokes() {}
template<class P>
void oldroyd_theta_scheme<P>::reset(const geo& omega) {
  Th = space (omega, "P1d", "tensor");
  Xh = P::velocity_space (omega, "P2");
  Qh = space (omega, "P1d");
  theta = 1-1/sqrt(2.);
  lambda = Re/(theta*delta_t);
  eta = ((1 - alpha)*We + theta*delta_t)/(We + theta*delta_t);
  nu = 1/((1-2*theta)*delta_t);
  c1 = We/(We + theta*delta_t);
  c2 = - We*theta*delta_t/(We + theta*delta_t);
  c3 = alpha*theta*delta_t/(We + theta*delta_t);
  c4 = 1/((1-2*theta)*delta_t) - 1/We;
  c5 = alpha/We;
  trial u (Xh), tau(Th), p (Qh);
  test  v (Xh), xi (Th), q (Qh);
  mt = integrate (ddot(tau,xi));
  mu = integrate (dot(u,v));
  mp = integrate (p*q);
  integrate_option iopt;
  iopt.invert = true;
  inv_mt = integrate (ddot(tau,xi), iopt);
  b  = integrate (-ddot(tau,D(v)));
  c = integrate (lambda*dot(u,v) + 2*eta*ddot(D(u),D(v)));
  d  = integrate (-div(u)*q);
  stokes = problem_mixed (c, d);
  stokes.set_metric (mp);
}