From a90b19a1f85ff432fa522a42fb22db6d67b617f4 Mon Sep 17 00:00:00 2001
From: Andrea Mola <mola.andrea@gmail.com>
Date: Mon, 28 Dec 2015 10:22:52 +0000
Subject: [PATCH] fixed a bug on the version that imposes water height at
 inflow

git-svn-id: svn://svn.sissa.it/openship@567 958231e9-8e66-0410-a4a0-f664109d2741
---
 include/restart_nonlinear_problem_diff.h |   1 +
 source/free_surface.cc                   | 116 ++++++++++++-----------
 waveBem.prm                              |   4 +-
 3 files changed, 65 insertions(+), 56 deletions(-)

diff --git a/include/restart_nonlinear_problem_diff.h b/include/restart_nonlinear_problem_diff.h
index e7d102d..5ae90f3 100644
--- a/include/restart_nonlinear_problem_diff.h
+++ b/include/restart_nonlinear_problem_diff.h
@@ -68,6 +68,7 @@ public NewtonArgument
       if ( (comp_dom.flags[i] & water) &&
            (comp_dom.flags[i] & near_boat)==0 &&
            (comp_dom.flags[i] & transom_on_water)==0 &&
+           (comp_dom.flags[i] & near_inflow)==0 &&
            (comp_dom.vector_constraints.is_constrained(3*i+2))==0 )
          {
          //cout<<"W "<<i<<" ("<<3*i+2<<")"<<endl;
diff --git a/source/free_surface.cc b/source/free_surface.cc
index 4ed121b..b462b08 100644
--- a/source/free_surface.cc
+++ b/source/free_surface.cc
@@ -4411,7 +4411,8 @@ std::cout<<"Preparing interpolated solution for restart"<<std::endl;
                      loc_y_smooth_res[i] += loc_stiffness_matrix[i][j]*(y_displs[j]-(1-blend_factor)*comp_dom.old_map_points(3*local_dof_indices[j]+1));
                      }
                  if ( !constraints.is_constrained(local_dof_indices[i]) &&
-                      !(comp_dom.flags[local_dof_indices[i]] & transom_on_water) )
+                      !(comp_dom.flags[local_dof_indices[i]] & transom_on_water) &&
+                      !(comp_dom.flags[local_dof_indices[i]] & near_inflow))
                     {
                     unsigned int ii = local_dof_indices[i];
                     eta_res[ii] += loc_eta_res[i].val();
@@ -4419,7 +4420,8 @@ std::cout<<"Preparing interpolated solution for restart"<<std::endl;
                                          
                     }
                  if ( !(constraints.is_constrained(local_dof_indices[i])) &&
-                      !(comp_dom.flags[local_dof_indices[i]] & edge) )
+                      !(comp_dom.flags[local_dof_indices[i]] & edge) &&
+                      !(comp_dom.flags[local_dof_indices[i]] & near_inflow))
                     {
                     unsigned int ii = local_dof_indices[i];
                     x_smoothing_res[ii] += loc_x_smooth_res[i].val();
@@ -4619,10 +4621,10 @@ std::cout<<"Preparing interpolated solution for restart"<<std::endl;
   RestartNonlinearProblemDiff rest_nonlin_prob_diff(*this,comp_dom,t,y,yp,jacobian_dot_matrix);
   std::map<unsigned int,unsigned int> &map_diff = rest_nonlin_prob_diff.indices_map;
 
-/*
+
   // these lines test the correctness of the jacobian for the
   // restart (reduced) nonlinear problem
-
+/*
   Vector<double> restart_prob_solution(rest_nonlin_prob_diff.n_dofs());
   Vector<double> restart_prob_residual(rest_nonlin_prob_diff.n_dofs());
   for (std::map<unsigned int, unsigned int>::iterator it = map_diff.begin(); it != map_diff.end(); ++it)
@@ -4646,7 +4648,7 @@ std::cout<<"Preparing interpolated solution for restart"<<std::endl;
      rest_nonlin_prob_diff.residual(restart_prob_residual,restart_prob_solution);
      cout<<"----------Test---------"<<endl;
      for (unsigned int i=0; i<rest_nonlin_prob_diff.n_dofs(); ++i)
-         if (fabs(restart_prob_residual(i)-delta_res(i)) > 1e-10)
+         if (fabs(restart_prob_residual(i)-delta_res(i)) > 1e-15)
          cout<<i<<"  "<<delta_res(i)<<" vs "<<restart_prob_residual(i)<<"      err "<<restart_prob_residual(i)-delta_res(i)<<"   "<<sys_comp(i)<<endl;
      delta_res*=-1;
      delta_res.add(restart_prob_residual);
@@ -4681,9 +4683,9 @@ std::cout<<"Preparing interpolated solution for restart"<<std::endl;
      //output_results(filename1, t, y, yp);
 
 
-   
+/*   
      // these lines test the jacobian of the DAE system
-/* 
+ 
      Vector<double> delta_y(this->n_dofs());
      Vector<double> delta_res(this->n_dofs());
 
@@ -4969,11 +4971,11 @@ int FreeSurface<dim>::jacobian(const double t,
   bem.solve_system(bem_phi, bem_dphi_dn, bem_bc);
   
 /*
-   cout<<"44!!!!!!"<<endl;
-    for (SparsityPattern::iterator col=jacobian_sparsity_pattern.begin(1484); col!=jacobian_sparsity_pattern.end(1484); ++col)
+   cout<<"JACOBIANO!!!!!!"<<endl;
+    for (SparsityPattern::iterator col=jacobian_sparsity_pattern.begin(80); col!=jacobian_sparsity_pattern.end(80); ++col)
         {
         unsigned int j = col->column(); 
-        cout<<j<<"("<<jacobian_matrix(1484,j)<<") ";
+        cout<<j<<"("<<jacobian_matrix(80,j)<<") ";
         }
     cout<<endl;
 //*/
@@ -5304,7 +5306,7 @@ int FreeSurface<dim>::residual_and_jacobian(const double t,
          if (comp_dom.flags[i] & near_inflow)
             {
             working_map_points(3*i+2) = comp_dom.old_map_points(3*i+2);
-            jacobian_matrix.add(3*i+2,3*i+2,1.0);
+            jacobian_matrix.set(3*i+2,3*i+2,1.0);
             }
          //cout<<"& "<<3*i<<" "<<3*i+1<<endl;
          }
@@ -7035,26 +7037,29 @@ int FreeSurface<dim>::residual_and_jacobian(const double t,
                      loc_y_smooth_res[i] += loc_stiffness_matrix[i][j]*(y_displs[j]-(1-blend_factor)*comp_dom.old_map_points(3*local_dof_indices[j]+1));
                      }
                  if ( !constraints.is_constrained(local_dof_indices[i]) &&
-                      !(comp_dom.flags[local_dof_indices[i]] & transom_on_water))
+                      !(comp_dom.flags[local_dof_indices[i]] & transom_on_water)  )
                     {
                     unsigned int ii = local_dof_indices[i];
-                    eta_res[ii] += loc_eta_res[i].val();
+                    if (!(comp_dom.flags[local_dof_indices[i]] & near_inflow))
+                       eta_res[ii] += loc_eta_res[i].val();
                     phi_res[ii] += loc_phi_res[i].val();
                     //cout<<"* "<<cell<<" "<<local_dof_indices[i]<<" "<<loc_phi_res[i]<<endl;
                     for (unsigned int j=0;j<dofs_per_cell;++j)
                         {
                         unsigned int jj = local_dof_indices[j];
-                        for (unsigned int k=0; k<dim; ++k)
-                            jacobian_matrix.add(3*ii+2,
-                                                3*jj+k,
-                                                loc_eta_res[i].fastAccessDx(3*j+k));
-                        jacobian_matrix.add(3*ii+2,
-                                            jj+comp_dom.vector_dh.n_dofs(),
-                                            loc_eta_res[i].fastAccessDx(3*dofs_per_cell+j));
-                        jacobian_matrix.add(3*ii+2,
-                                            jj+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs(),
-                                            loc_eta_res[i].fastAccessDx(4*dofs_per_cell+j));
-                       
+                        if (!(comp_dom.flags[local_dof_indices[i]] & near_inflow))
+                           {
+                           for (unsigned int k=0; k<dim; ++k)
+                               jacobian_matrix.add(3*ii+2,
+                                                   3*jj+k,
+                                                   loc_eta_res[i].fastAccessDx(3*j+k));
+                           jacobian_matrix.add(3*ii+2,
+                                               jj+comp_dom.vector_dh.n_dofs(),
+                                               loc_eta_res[i].fastAccessDx(3*dofs_per_cell+j));
+                           jacobian_matrix.add(3*ii+2,
+                                               jj+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs(),
+                                               loc_eta_res[i].fastAccessDx(4*dofs_per_cell+j));
+                           }
                         for (unsigned int k=0; k<dim; ++k)
                             jacobian_matrix.add(ii+comp_dom.vector_dh.n_dofs(),
                                                 3*jj+k,
@@ -7065,17 +7070,19 @@ int FreeSurface<dim>::residual_and_jacobian(const double t,
                         jacobian_matrix.add(ii+comp_dom.vector_dh.n_dofs(),
                                             jj+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs(),
                                             loc_phi_res[i].fastAccessDx(4*dofs_per_cell+j));
-                                                      
-                        for (unsigned int k=0; k<dim; ++k)
-                            jacobian_dot_matrix.add(3*ii+2,
-                                                    3*jj+k,
-                                                    loc_eta_res[i].fastAccessDx(5*dofs_per_cell+3*j+k));
-                        jacobian_dot_matrix.add(3*ii+2,
-                                                jj+comp_dom.vector_dh.n_dofs(),
-                                                loc_eta_res[i].fastAccessDx(8*dofs_per_cell+j));
-                        jacobian_dot_matrix.add(3*ii+2,
-                                                jj+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs(),
-                                                loc_eta_res[i].fastAccessDx(9*dofs_per_cell+j));
+                        if (!(comp_dom.flags[local_dof_indices[i]] & near_inflow))
+                           {                                                      
+                           for (unsigned int k=0; k<dim; ++k)
+                               jacobian_dot_matrix.add(3*ii+2,
+                                                          3*jj+k,
+                                                       loc_eta_res[i].fastAccessDx(5*dofs_per_cell+3*j+k));
+                           jacobian_dot_matrix.add(3*ii+2,
+                                                   jj+comp_dom.vector_dh.n_dofs(),
+                                                   loc_eta_res[i].fastAccessDx(8*dofs_per_cell+j));
+                           jacobian_dot_matrix.add(3*ii+2,
+                                                   jj+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs(),
+                                                   loc_eta_res[i].fastAccessDx(9*dofs_per_cell+j));
+                           }
                         for (unsigned int k=0; k<dim; ++k)
                             jacobian_dot_matrix.add(ii+comp_dom.vector_dh.n_dofs(),
                                                     3*jj+k,
@@ -7087,23 +7094,25 @@ int FreeSurface<dim>::residual_and_jacobian(const double t,
                                                 jj+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs(),
                                                 loc_phi_res[i].fastAccessDx(9*dofs_per_cell+j));
                         }
-
-                    for (unsigned int k=0; k<dim; ++k)
-                        jacobian_matrix.add(3*ii+2,
-                                            k+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
-                                            loc_eta_res[i].fastAccessDx(k+10*dofs_per_cell));
-                    for (unsigned int k=0; k<dim; ++k)
-                        jacobian_matrix.add(3*ii+2,
-                                            k+3+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
-                                            loc_eta_res[i].fastAccessDx(k+3+10*dofs_per_cell));
-                    for (unsigned int k=0; k<dim; ++k)
-                        jacobian_matrix.add(3*ii+2,
-                                            k+6+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
-                                            loc_eta_res[i].fastAccessDx(k+6+10*dofs_per_cell));
-                    for (unsigned int k=0; k<4; ++k)
-                        jacobian_matrix.add(3*ii+2,
-                                            k+9+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
-                                            loc_eta_res[i].fastAccessDx(k+9+10*dofs_per_cell));
+                    if (!(comp_dom.flags[local_dof_indices[i]] & near_inflow))
+                       {
+                       for (unsigned int k=0; k<dim; ++k)
+                           jacobian_matrix.add(3*ii+2,
+                                               k+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
+                                               loc_eta_res[i].fastAccessDx(k+10*dofs_per_cell));
+                       for (unsigned int k=0; k<dim; ++k)
+                           jacobian_matrix.add(3*ii+2,
+                                               k+3+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
+                                               loc_eta_res[i].fastAccessDx(k+3+10*dofs_per_cell));
+                       for (unsigned int k=0; k<dim; ++k)
+                           jacobian_matrix.add(3*ii+2,
+                                               k+6+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
+                                               loc_eta_res[i].fastAccessDx(k+6+10*dofs_per_cell));
+                       for (unsigned int k=0; k<4; ++k)
+                           jacobian_matrix.add(3*ii+2,
+                                               k+9+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
+                                               loc_eta_res[i].fastAccessDx(k+9+10*dofs_per_cell));
+                       }
                     for (unsigned int k=0; k<dim; ++k)
                         jacobian_matrix.add(ii+comp_dom.vector_dh.n_dofs(),
                                             k+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs(),
@@ -8112,8 +8121,7 @@ for (unsigned int i=0; i<this->n_dofs(); ++i)
         {
         unsigned int j = col->column();
         //cout<<" "<<i<<" "<<j<<" "<<jacobian_matrix(i,j)+alpha*jacobian_dot_matrix(i,j)<<endl;
-        //if ( (j > (int)i-preconditioner_band) &&
-        //     (j > (int)i+preconditioner_band)   ) 
+ 
            jacobian_preconditioner_matrix.set(i,j,jacobian_matrix(i,j)+alpha*jacobian_dot_matrix(i,j));
         }
 
diff --git a/waveBem.prm b/waveBem.prm
index 0a999e1..b5f9080 100644
--- a/waveBem.prm
+++ b/waveBem.prm
@@ -156,8 +156,8 @@ subsection Boundary Conditions ID Numbers
   set Wall Surface 1 ID = 3
   set Wall Surface 2 ID = 4
   set Wall Surface 3 ID = 0
-  set Inflow Surface 1 ID = 0
-  set Inflow Surface 2 ID = 0
+  set Inflow Surface 1 ID = 9
+  set Inflow Surface 2 ID = 10
   set Inflow Surface 3 ID = 0
   set Pressure Surface ID = 0
   set Free Surface Edge On Boat ID = 20