first 6dofs version that is advancing in time. tolerances in DAE solver have been invreased to 1e-2/1e-3

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

Author: andreamola

include/boat_model.h

@@ -43,7 +43,9 @@ class BoatModel{
 
     Point<3> compute_hydrostatic_moment(const double &sink);
 
-    gp_Trsf set_current_position(const Point<3> &translation_vect, const double &trim);
+    gp_Trsf set_current_position(const Point<3> &translation_vect,
+                                 const double &quaternion_scalar,
+                                 const Point<3> &quaternion_vect);
 
 //private: 
 				     // keel intersection with undisturbed free surface at bow

include/free_surface.h

@@ -461,6 +461,8 @@ public NewtonArgument
 
   TopLoc_Location restart_hull_location;
   Point<3> restart_hull_displacement;
+  Point<3> restart_hull_quat_vector;
+  double restart_hull_quat_scalar;
   Point<3> restart_transom_center_point;
   Point<3> restart_transom_left_point;
   Point<3> restart_transom_right_point;

include/numerical_towing_tank.h

@@ -328,6 +328,15 @@ double assigned_sink;
                                    // hydrostatic equilibrium position (which is assumed to be the angular
                                    // position assigned in the CAD file)
 double assigned_trim;
+
+                                   // hull moments of inertia
+double Ixx;
+double Ixy;
+double Ixz;
+double Iyy;
+double Iyz;
+double Izz;
+
                                    // here we have the maximum aspect ratio allowed for the cells
 double max_aspect_ratio;
                                    // the next is a factor determining the inclination of the mesh longitudinal

include/restart_nonlinear_problem_diff.h

@@ -93,7 +93,7 @@ public NewtonArgument
           }
 
         // we now add the rigid modes dofs
-  for (unsigned int d=0; d<6; ++d)
+  for (unsigned int d=0; d<13; ++d)
       {
       rigid_modes_indices.insert(d);
       }
@@ -163,7 +163,7 @@ public NewtonArgument
                                          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;
+      line_lengths[i] = 40;
   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();

source/boat_model.cc

@@ -84,6 +84,7 @@
 #include <GeomConvert.hxx>
 #include <Geom_BSplineSurface.hxx>
 #include <BRepAdaptor_Curve.hxx>
+#include <gp_Quaternion.hxx>
 
 #include <deal.II/grid/grid_reordering.h>
 #include <deal.II/grid/grid_tools.h>
@@ -805,15 +806,15 @@ cout<<"UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU"<<endl;
       return hydrostatic_force;
 }
 
-gp_Trsf BoatModel::set_current_position(const Point<3> &translation_vect, const double &trim)
+gp_Trsf BoatModel::set_current_position(const Point<3> &translation_vect,
+                                        const double &quaternion_scalar,
+                                        const Point<3> &quaternion_vect)
 {
 
-                                  //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);
+                                  //here we prepare the rotation of the boat of the requested trim angle
+     gp_Quaternion rot_quaternion(quaternion_vect(0), quaternion_vect(1), quaternion_vect(2), quaternion_scalar);    
      gp_Trsf rotation;
-     rotation.SetRotation(rot_axis,trim);
+     rotation.SetRotation(rot_quaternion);
                                   //we first get the full transformation currently applied to the shape
      TopLoc_Location prev_L = reference_loc; 
      gp_Trsf prev_Transf = prev_L.Transformation();
@@ -838,6 +839,14 @@ gp_Trsf BoatModel::set_current_position(const Point<3> &translation_vect, const
      keel_edge.Location(new_L);
      current_loc = new_L;
 
+     // SHOULD WE PRINT THE EULER ANGLES UP HERE? IS THE FACT WE HAVE A X AXIS
+     // DIRECTED BOW TO STERN GOING TO CHANGE SOMETHING?
+     double yaw_angle, pitch_angle, roll_angle;
+     rot_quaternion.GetEulerAngles(gp_YawPitchRoll, yaw_angle, pitch_angle, roll_angle);
+     cout<<"Current Yaw Angle: "<<yaw_angle<<endl;
+     cout<<"Current Pitch Angle: "<<pitch_angle<<endl;
+     cout<<"Current Roll Angle: "<<roll_angle<<endl;
+
 
      if (is_transom)
         {
@@ -904,7 +913,7 @@ gp_Trsf BoatModel::set_current_position(const Point<3> &translation_vect, const
 
 
 
-     return current_loc.Transformation();
+     return this_transf;//current_loc.Transformation();
 }