Improvements.

Author: l-korous

hermes2d/include/discrete_problem/discrete_problem.h

@@ -52,21 +52,21 @@ namespace Hermes
       /// 2 - allows for assembling Dirichlet boundary conditions using a Dirichlet lift.
       /// \param[in] dirichlet_lift_accordingly If true, the appropriate settings for (linear / nonlinear)
       /// problem will be used (use Dirichlet lift iff the problem is linear). If false, the other way round.
-      DiscreteProblem(WeakFormSharedPtr<Scalar> wf, std::vector<SpaceSharedPtr<Scalar> > spaces, bool linear = false, bool dirichlet_lift_accordingly = true);
+      DiscreteProblem(WeakFormSharedPtr<Scalar> wf, std::vector<SpaceSharedPtr<Scalar> > spaces, bool linear = false, bool dirichlet_lift_accordingly = true, bool use_direct_for_Dirichlet_lift = false);
       /// Constructor for one equation.
       /// Making this DiscreteProblem linear does 2 things
       /// 1 - turns off initialization of previous iterations for nonlinear solvers.
       /// 2 - allows for assembling Dirichlet boundary conditions using a Dirichlet lift.
       /// \param[in] dirichlet_lift_accordingly If true, the appropriate settings for (linear / nonlinear)
       /// problem will be used (use Dirichlet lift iff the problem is linear). If false, the other way round.
-      DiscreteProblem(WeakFormSharedPtr<Scalar> wf, SpaceSharedPtr<Scalar> space, bool linear = false, bool dirichlet_lift_accordingly = true);
+      DiscreteProblem(WeakFormSharedPtr<Scalar> wf, SpaceSharedPtr<Scalar> space, bool linear = false, bool dirichlet_lift_accordingly = true, bool use_direct_for_Dirichlet_lift = false);
       /// Non-parameterized constructor.
       /// Making this DiscreteProblem linear does 2 things
       /// 1 - turns off initialization of previous iterations for nonlinear solvers.
       /// 2 - allows for assembling Dirichlet boundary conditions using a Dirichlet lift.
       /// \param[in] dirichlet_lift_accordingly If true, the appropriate settings for (linear / nonlinear)
       /// problem will be used (use Dirichlet lift iff the problem is linear). If false, the other way round.
-      DiscreteProblem(bool linear = false, bool dirichlet_lift_accordingly = true);
+      DiscreteProblem(bool linear = false, bool dirichlet_lift_accordingly = true, bool use_direct_for_Dirichlet_lift = false);
       /// Destuctor.
       virtual ~DiscreteProblem();
 
@@ -128,7 +128,7 @@ namespace Hermes
       inline std::string getClassName() const { return "DiscreteProblem"; }
 
       /// Init function. Common code for the constructors.
-      void init(bool linear, bool dirichlet_lift_accordingly);
+      void init(bool linear, bool dirichlet_lift_accordingly, bool use_direct_for_Dirichlet_lift);
 
       /// Space instances for all equations in the system.
       std::vector<SpaceSharedPtr<Scalar> > spaces;
@@ -138,7 +138,7 @@ namespace Hermes
       Vector<Scalar>* dirichlet_lift_rhs;
 
       /// Internal.
-      bool nonlinear, add_dirichlet_lift;
+      bool nonlinear, add_dirichlet_lift, use_direct_for_Dirichlet_lift;
 
       /// DiscreteProblemMatrixVector methods.
       bool set_matrix(SparseMatrix<Scalar>* mat);

hermes2d/src/discrete_problem/discrete_problem.cpp

@@ -29,29 +29,29 @@ namespace Hermes
   namespace Hermes2D
   {
     template<typename Scalar>
-    DiscreteProblem<Scalar>::DiscreteProblem(WeakFormSharedPtr<Scalar> wf_, std::vector<SpaceSharedPtr<Scalar> > spaces, bool to_set, bool dirichlet_lift_accordingly)
+    DiscreteProblem<Scalar>::DiscreteProblem(WeakFormSharedPtr<Scalar> wf_, std::vector<SpaceSharedPtr<Scalar> > spaces, bool to_set, bool dirichlet_lift_accordingly, bool use_direct_for_Dirichlet_lift)
     {
-      this->init(to_set, dirichlet_lift_accordingly);
+      this->init(to_set, dirichlet_lift_accordingly, use_direct_for_Dirichlet_lift);
       this->set_spaces(spaces);
       this->set_weak_formulation(wf_);
     }
 
     template<typename Scalar>
-    DiscreteProblem<Scalar>::DiscreteProblem(WeakFormSharedPtr<Scalar> wf_, SpaceSharedPtr<Scalar> space, bool to_set, bool dirichlet_lift_accordingly)
+    DiscreteProblem<Scalar>::DiscreteProblem(WeakFormSharedPtr<Scalar> wf_, SpaceSharedPtr<Scalar> space, bool to_set, bool dirichlet_lift_accordingly, bool use_direct_for_Dirichlet_lift)
     {
-      this->init(to_set, dirichlet_lift_accordingly);
+      this->init(to_set, dirichlet_lift_accordingly, use_direct_for_Dirichlet_lift);
       this->set_space(space);
       this->set_weak_formulation(wf_);
     }
 
     template<typename Scalar>
-    DiscreteProblem<Scalar>::DiscreteProblem(bool to_set, bool dirichlet_lift_accordingly)
+    DiscreteProblem<Scalar>::DiscreteProblem(bool to_set, bool dirichlet_lift_accordingly, bool use_direct_for_Dirichlet_lift)
     {
-      init(to_set, dirichlet_lift_accordingly);
+      init(to_set, dirichlet_lift_accordingly, use_direct_for_Dirichlet_lift);
     }
 
     template<typename Scalar>
-    void DiscreteProblem<Scalar>::init(bool to_set, bool dirichlet_lift_accordingly)
+    void DiscreteProblem<Scalar>::init(bool to_set, bool dirichlet_lift_accordingly, bool use_direct_for_Dirichlet_lift)
     {
       this->reassembled_states_reuse_linear_system = nullptr;
 
@@ -64,7 +64,7 @@ namespace Hermes
         this->add_dirichlet_lift = this->nonlinear;
 
       if (this->add_dirichlet_lift)
-        this->dirichlet_lift_rhs = create_vector<Scalar>(false);
+        this->dirichlet_lift_rhs = create_vector<Scalar>(use_direct_for_Dirichlet_lift);
       else
         this->dirichlet_lift_rhs = nullptr;
 

hermes2d/src/projections/ogprojection.cpp

@@ -31,7 +31,7 @@ namespace Hermes
         throw Exceptions::NullException(2);
 
       // Initialize linear solver.
-      Hermes::Hermes2D::LinearSolver<Scalar> linear_solver(wf, space);
+      Hermes::Hermes2D::LinearSolver<Scalar> linear_solver(wf, space, true);
       linear_solver.set_verbose_output(false);
 
       // Perform Newton iteration.

hermes2d/src/solver/linear_solver.cpp

@@ -43,7 +43,7 @@ namespace Hermes
     template<typename Scalar>
     LinearSolver<Scalar>::LinearSolver(WeakFormSharedPtr<Scalar> wf, SpaceSharedPtr<Scalar> space, bool force_use_direct_solver) : Solver<Scalar>(false), Hermes::Solvers::MatrixSolver<Scalar>(force_use_direct_solver)
     {
-      this->dp = new DiscreteProblem<Scalar>(wf, space, true);
+      this->dp = new DiscreteProblem<Scalar>(wf, space, true, true, true);
       this->own_dp = true;
     }
 

hermes2d/test_examples/03-navier-stokes/definitions.cpp

@@ -71,12 +71,16 @@ class WeakFormNSSimpleLinearization : public WeakForm<double>
 
     double value(int n, double *wt, Func<double> *u_ext[], Func<double> *u, Func<double> *v, GeomVol<double> *e, Func<double> **ext) const{
       double result = int_grad_u_grad_v<double, double>(n, wt, u, v) / Reynolds;
+      if(!Stokes)
+        result += int_u_v<double, double>(n, wt, u, v) / time_step;
       return result;
     }
 
     Ord ord(int n, double *wt, Func<Ord> *u_ext[], Func<Ord> *u, Func<Ord> *v, GeomVol<Ord> *e, Func<Ord> **ext) const
     {
       Ord result = int_grad_u_grad_v<Ord, Ord>(n, wt, u, v) / Reynolds;
+      if(!Stokes)
+        result += int_u_v<Ord, Ord>(n, wt, u, v) / time_step;
       return result;
     }
 
@@ -179,13 +183,21 @@ class WeakFormNSSimpleLinearization : public WeakForm<double>
     VectorFormVolVel(int i, bool Stokes, double time_step) : VectorFormVol<double>(i), Stokes(Stokes), time_step(time_step) {
     }
 
-    double value(int n, double *wt, Func<double> *u_ext[], Func<double> *v, GeomVol<double> *e, Func<double> **ext) const{
+    double value(int n, double *wt, Func<double> *u_ext[], Func<double> *v, Geom<double> *e, Func<double> **ext) const{
       double result = 0;
+      if(!Stokes) {
+        Func<double>* vel_prev_time = ext[2]; // this form is used with both velocity components
+        result = int_u_v<double, double>(n, wt, vel_prev_time, v) / time_step;
+      }
       return result;
     }
 
-    Ord ord(int n, double *wt, Func<Ord> *u_ext[], Func<Ord> *v, GeomVol<Ord> *e, Func<Ord> **ext) const {
+    Ord ord(int n, double *wt, Func<Ord> *u_ext[], Func<Ord> *v, Geom<Ord> *e, Func<Ord> **ext) const {
       Ord result = Ord(0);
+      if(!Stokes) {
+        Func<Ord>* vel_prev_time = ext[2]; // this form is used with both velocity components
+        result = int_u_v<Ord, Ord>(n, wt, vel_prev_time, v) / time_step;
+      }
       return result;
     }
 
@@ -256,6 +268,8 @@ class WeakFormNSNewton : public WeakForm<double>
 
     double value(int n, double *wt, Func<double> *u_ext[], Func<double> *u, Func<double> *v, GeomVol<double> *e, Func<double> **ext) const{
       double result = int_grad_u_grad_v<double, double>(n, wt, u, v) / Reynolds;
+      if(!Stokes)
+        result += int_u_v<double, double>(n, wt, u, v) / time_step;
       return result;
     }
 
@@ -496,7 +510,8 @@ class WeakFormNSNewton : public WeakForm<double>
         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]));
       if(!Stokes)
         for (int i = 0; i < n; i++)
-          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];
+          result += wt[i] * (((xvel_prev_newton->val[i] - xvel_prev_time->val[i]) * v->val[i] / time_step )
+          + ((xvel_prev_newton->val[i] * xvel_prev_newton->dx[i] + yvel_prev_newton->val[i] * xvel_prev_newton->dy[i]) * v->val[i]));
       return result;
     }
 
@@ -511,7 +526,8 @@ class WeakFormNSNewton : public WeakForm<double>
         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]));
       if(!Stokes)
         for (int i = 0; i < n; i++)
-          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];
+          result += wt[i] * (((xvel_prev_newton->val[i] - xvel_prev_time->val[i]) * v->val[i] / time_step)
+          + ((xvel_prev_newton->val[i] * xvel_prev_newton->dx[i] + yvel_prev_newton->val[i] * xvel_prev_newton->dy[i]) * v->val[i]));
       return result;
     }
 
@@ -543,7 +559,8 @@ class WeakFormNSNewton : public WeakForm<double>
         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]));
       if(!Stokes)
         for (int i = 0; i < n; i++)
-          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];
+          result += wt[i] * (((yvel_prev_newton->val[i] - yvel_prev_time->val[i]) * v->val[i] / time_step )
+          + ((xvel_prev_newton->val[i] * yvel_prev_newton->dx[i] + yvel_prev_newton->val[i] * yvel_prev_newton->dy[i]) * v->val[i]));
       return result;
     }
 
@@ -631,8 +648,11 @@ class EssentialBCNonConst : public EssentialBoundaryCondition<double>
       };
 
       virtual double value(double x, double y) const {
-        double val_y = -x*(1-x);
-        return val_y / 4.;
+        double val_y = vel_inlet * y*(H-y) / (H/2.)/(H/2.);
+        if (current_time <= startup_time) 
+          return val_y * current_time/startup_time;
+        else 
+          return val_y;
       };
 
 protected: