first WaveBEM version with FSI up and running. Rotations are yet to be added, plus we need to add rigid motions in output files and correct hull nodes velocities in .vtu outputs. But apparently it is working

git-svn-id: svn://svn.sissa.it/openship@560 958231e9-8e66-0410-a4a0-f664109d2741

Author: andreamola

include/boat_model.h

@@ -43,7 +43,7 @@ class BoatModel{
 
     Point<3> compute_hydrostatic_moment(const double &sink);
 
-    gp_Trsf set_current_position(const double &sink);
+    gp_Trsf set_current_position(const Point<3> &translation_vect, const double &trim);
 
 //private: 
 				     // keel intersection with undisturbed free surface at bow
@@ -143,6 +143,11 @@ class BoatModel{
                                      // flag to determine if we have a transom stern or not
     bool is_transom;
 
+    double boat_mass;
+
+    Point<3> current_left_transom_tangent;
+
+    Point<3> current_right_transom_tangent;
 
 
 };

include/free_surface.h

@@ -309,8 +309,14 @@ public NewtonArgument
   Functions::ParsedFunction<dim> initial_wave_potential;
   Functions::ParsedFunction<dim> initial_norm_potential_grad;
   Functions::ParsedFunction<dim> inflow_norm_potential_grad;
+  Functions::ParsedFunction<1> hull_x_axis_translation;
+  Functions::ParsedFunction<1> hull_y_axis_translation;
+  Functions::ParsedFunction<1> hull_z_axis_translation;
+
+  bool is_hull_x_translation_imposed;
+  bool is_hull_y_translation_imposed;
+  bool is_hull_z_translation_imposed;
 
-  
   std::string node_displacement_type;
 
   SolverControl solver_control;
@@ -454,7 +460,12 @@ public NewtonArgument
   double adaptive_ref_limit;
 
   TopLoc_Location restart_hull_location;
-
+  Point<3> restart_hull_displacement;
+  Point<3> restart_transom_center_point;
+  Point<3> restart_transom_left_point;
+  Point<3> restart_transom_right_point;
+  Point<3> restart_transom_left_tangent;
+  Point<3> restart_transom_right_tangent;
 
 };
 

include/restart_nonlinear_problem_diff.h

@@ -36,6 +36,8 @@ public NewtonArgument
   bow_stern_indices.clear();
   water_indices.clear();
   phi_indices.clear();
+  rigid_modes_indices.clear();
+
   if (!comp_dom.no_boat)
      for (unsigned int i=0; i<comp_dom.dh.n_dofs(); ++i)
          {
@@ -89,7 +91,12 @@ public NewtonArgument
                 }
              }
           }
-
+  
+        // we now add the rigid modes dofs
+  for (unsigned int d=0; d<6; ++d)
+      {
+      rigid_modes_indices.insert(d);
+      }
 
   free_surf_jac_x_delta.reinit(free_surf_y.size());
   free_surf_res.reinit(free_surf_y.size());
@@ -130,6 +137,13 @@ public NewtonArgument
       count++;
       }  
 
+  for (std::set <unsigned int>::iterator pos = rigid_modes_indices.begin(); pos != rigid_modes_indices.end(); ++pos)
+      {
+      unsigned int i=*pos;
+      indices_map[i+comp_dom.vector_dh.n_dofs()+comp_dom.dh.n_dofs()+comp_dom.dh.n_dofs()] = count;
+      count++;
+      }
+
 
   //cout<<" MAP "<<endl;
   //for (std::map <unsigned int,unsigned int>::iterator pos = indices_map.begin(); pos != indices_map.end(); ++pos)
@@ -143,9 +157,20 @@ public NewtonArgument
   //cout<<water_indices.size()<<endl;
   //cout<<phi_indices.size()<<endl;
   //cout<<2*water_line_indices.size()+3*bow_stern_indices.size()+water_indices.size()+phi_indices.size()+dphi_dn_indices.size()<<endl;
-  jacobian_sparsity_pattern.reinit(2*water_line_indices.size()+3*bow_stern_indices.size()+water_indices.size()+phi_indices.size(),
-                                   2*water_line_indices.size()+3*bow_stern_indices.size()+water_indices.size()+phi_indices.size(),
-                                   30);
+  std::vector<unsigned int> line_lengths(2*water_line_indices.size()+
+                                         3*bow_stern_indices.size()+
+                                         water_indices.size()+
+                                         phi_indices.size()+
+                                         rigid_modes_indices.size());
+  for (unsigned int i=0; i<2*water_line_indices.size()+3*bow_stern_indices.size()+water_indices.size()+phi_indices.size(); ++i)
+      line_lengths[i] = 30;
+  for (unsigned int i=0; i<rigid_modes_indices.size(); ++i)
+      line_lengths[i+2*water_line_indices.size()+3*bow_stern_indices.size()+water_indices.size()+phi_indices.size()] = 
+                   2*water_line_indices.size()+3*bow_stern_indices.size()+water_indices.size()+phi_indices.size()+rigid_modes_indices.size();
+
+  jacobian_sparsity_pattern.reinit(2*water_line_indices.size()+3*bow_stern_indices.size()+water_indices.size()+phi_indices.size()+rigid_modes_indices.size(),
+                                   2*water_line_indices.size()+3*bow_stern_indices.size()+water_indices.size()+phi_indices.size()+rigid_modes_indices.size(),
+                 line_lengths );
 
   count = 0;
   for (std::set <unsigned int>::iterator pos = water_line_indices.begin(); pos != water_line_indices.end(); ++pos)
@@ -196,7 +221,15 @@ public NewtonArgument
       count++;
       }
 
-
+  for (std::set <unsigned int>::iterator pos = rigid_modes_indices.begin(); pos != rigid_modes_indices.end(); ++pos)
+      {
+      unsigned int i=*pos;
+      for (SparsityPattern::iterator col=free_surf_jacobian_dot.get_sparsity_pattern().begin(i+comp_dom.vector_dh.n_dofs());
+           col!=free_surf_jacobian_dot.get_sparsity_pattern().end(i+comp_dom.vector_dh.n_dofs()); ++col)
+          if ( indices_map.count(col->column()) )
+             jacobian_sparsity_pattern.add(count,indices_map.find(col->column())->second);
+      count++;
+      }
   jacobian_sparsity_pattern.compress();
   jacobian_matrix.reinit(jacobian_sparsity_pattern);
 
@@ -250,6 +283,7 @@ public NewtonArgument
   std::set<unsigned int> bow_stern_indices;
   std::set<unsigned int> water_indices;
   std::set<unsigned int> phi_indices;
+  std::set<unsigned int> rigid_modes_indices;
   SparsityPattern      jacobian_sparsity_pattern;
   SparseMatrix<double> jacobian_matrix;
 

source/boat_model.cc

@@ -288,6 +288,7 @@ void BoatModel::start_iges_model(std::string igesFileName,
      }
 
      Point<3> hs_force = compute_hydrostatic_force(0.0);
+     boat_mass = hs_force(2)/9.81;
      Point<3> hs_moment = compute_hydrostatic_moment(0.0);
 
      double rot_axis_x_coor = -hs_moment(1)/hs_force(2);
@@ -804,26 +805,26 @@ cout<<"UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU"<<endl;
       return hydrostatic_force;
 }
 
-gp_Trsf BoatModel::set_current_position(const double &sink)
+gp_Trsf BoatModel::set_current_position(const Point<3> &translation_vect, const double &trim)
 {
 
                                   //here we prepare the rotation of the boat of the requested trim angle     
-     //gp_Pnt rot_center(0.0,0.0,0.0);
-     //gp_Dir rot_dir(0.0,1.0,0.0);
-     //gp_Ax1 rot_axis(rot_center, rot_dir);
-     //gp_Trsf rotation;
-     //rotation.SetRotation(rot_axis,assigned_trim);
+     gp_Pnt rot_center(0.0,0.0,0.0);
+     gp_Dir rot_dir(0.0,1.0,0.0);
+     gp_Ax1 rot_axis(rot_center, rot_dir);
+     gp_Trsf rotation;
+     rotation.SetRotation(rot_axis,trim);
                                   //we first get the full transformation currently applied to the shape
      TopLoc_Location prev_L = reference_loc; 
      gp_Trsf prev_Transf = prev_L.Transformation();
 
                                   // here we prepare the translation of the boat of the requested sink
      gp_Trsf translation;
-     gp_Vec vrt_displ(0.0,0.0,sink);
+     gp_Vec vrt_displ(translation_vect(0),translation_vect(1),translation_vect(2));
      translation.SetTranslation(vrt_displ);
 
-
-     gp_Trsf this_transf = translation;
+                                  // the rotation and translation are combined in a single transformation
+     gp_Trsf this_transf = translation*rotation;
 
                                   // the rotation and translation are combined in a single transformation
      gp_Trsf new_Transf = this_transf*prev_Transf;
@@ -867,6 +868,27 @@ gp_Trsf BoatModel::set_current_position(const double &sink)
                     ExcMessage("Transom edges don't possess a single intersection with horizontal plane: is transom not immersed any more?"));
          transomLeft = IntersectorLeft.Point(1);
          transomRight = IntersectorRight.Point(1);
+
+         double param;
+         double u;
+         double v;
+
+         IntersectorLeft.Parameters(1, u, v, param);
+         gp_Vec current_left_transom_tangent_vec = left_transom_bspline->DN(param,1);
+         gp_Dir current_left_transom_tangent_dir(current_left_transom_tangent_vec);
+         current_left_transom_tangent_dir.Transform(current_loc);
+         current_left_transom_tangent = Point<3>(current_left_transom_tangent_dir.X(),
+                                                 current_left_transom_tangent_dir.Y(),
+                                                 current_left_transom_tangent_dir.Z());
+
+         IntersectorRight.Parameters(1, u, v, param);
+         gp_Vec current_right_transom_tangent_vec = right_transom_bspline->DN(param,1);
+         gp_Dir current_right_transom_tangent_dir(current_right_transom_tangent_vec);
+         current_right_transom_tangent_dir.Transform(current_loc);
+         current_right_transom_tangent = Point<3>(current_right_transom_tangent_dir.X(),
+                                                  current_right_transom_tangent_dir.X(),
+                                                  current_right_transom_tangent_dir.Z());
+       
          transomLeft.Transform(current_loc);
          transomRight.Transform(current_loc);
          CurrentPointLeftTransom = Pnt(transomLeft);