another step towards the moving hull model. now boat_model class has a reference (hydrostatic) TopLoc_Location and a current location. A method has been added to change the current location which accounts for the boat rigid motions. All the projectors and smoothers will now act on the current hull configuration.

andreamola · andreamola · commit afcc5d34c9fe · 2015-08-19T12:58:44.000Z
git-svn-id: svn://svn.sissa.it/openship@557 958231e9-8e66-0410-a4a0-f664109d2741
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -18,7 +18,7 @@ FIND_PACKAGE(OpenCASCADE MODULE REQUIRED)
 
 
 SET(TARGET "waveBem")
-set(CMAKE_BUILD_TYPE "Release")
+set(CMAKE_BUILD_TYPE "Debug")
 
 file(GLOB TARGET_SRC source/*cc)
 include_directories(
diff --git a/include/boat_model.h b/include/boat_model.h
@@ -85,8 +85,12 @@ class BoatModel{
     TopoDS_Shape left_undist_water_line;
 				     // this is the undisturbed left water surface
     TopoDS_Shape undisturbed_water_surface_face; 
-                                     // location of the translated curves (for roto-translated hull)
+                                     // location of the translated shapes (for hull roto-translated
+                                     // in HYDROSTATIC position)
     TopLoc_Location reference_loc;
+                                     // location of the translated curves (for hull roto-translated
+                                     //in CURRENT position)
+    TopLoc_Location current_loc;
 				     // surface normal projector on boat surface right side
     OpenCascade::NormalProjection<2> *boat_surface_right;
 				     // surface normal projector on boat surface left side
diff --git a/include/numerical_towing_tank.h b/include/numerical_towing_tank.h
@@ -255,7 +255,7 @@ std::vector<Handle(Geom_Curve)> curves;
 std::vector<bool> on_curve_option;
                                    // locations are needed when reference configuration
                                    // is roto-translated
-std::vector<TopLoc_Location> ref_locations;
+std::vector< TopLoc_Location *> smoothers_locations;
  				   // nodes flags on the scalar dof_handler
 std::vector<GeometryFlags> flags;
  				   // nodes flags on the vector dof_handler
diff --git a/source/boat_model.cc b/source/boat_model.cc
@@ -309,7 +309,7 @@ void BoatModel::start_iges_model(std::string igesFileName,
      sh = right_transf.Shape();
      BRepBuilderAPI_Transform left_transf(refl_sh, Tcomp);
      refl_sh = left_transf.Shape();
-   
+
 
 
      cout<<"The hull has been placed in the correct position"<<endl;
@@ -398,6 +398,7 @@ Standard_Boolean OK = ICW.Write ("transom.igs");
   TopoDS_Edge edge =  TopoDS::Edge(edgeExplorer.Current());
   edge.Location(keel_edge.Location());
   this->reference_loc = keel_edge.Location();
+  this->current_loc =  keel_edge.Location();
   BRepAdaptor_Curve gg_curve(edge);
   equiv_keel_bspline = BRep_Tool::Curve(edge,L,First,Last);
   gp_Trsf L_transformation = L.Transformation(); 
@@ -406,11 +407,6 @@ Standard_Boolean OK = ICW.Write ("transom.igs");
 
   xyPlane->Transform(L_inv.Transformation());
   
-  cout<<First<<" "<<Last<<endl;
-  cout<<"TTEESSTT1a: "<<Pnt(equiv_keel_bspline->Value(Last))<<endl;
-  cout<<"TTEESSTT1b: "<<Pnt(gg_curve.Value(gg_curve.LastParameter()))<<endl;
-  cout<<"TTEESSTT2a: "<<Pnt(equiv_keel_bspline->Value(First))<<endl;
-  cout<<"TTEESSTT2b: "<<Pnt(gg_curve.Value(gg_curve.FirstParameter()))<<endl;
   TopExp_Explorer edge3Explorer(left_transom_edge, TopAbs_EDGE);
   TopoDS_Edge edge3 =  TopoDS::Edge(edge3Explorer.Current());
   left_transom_bspline = BRep_Tool::Curve(edge3,L,First,Last);
@@ -811,6 +807,31 @@ cout<<"UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUU"<<endl;
 void BoatModel::set_current_position(const double &sink)
 {
 
+                                  //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);
+                                  //we first get the full transformation currently applied to the shape
+     TopLoc_Location prev_L = sh.Location(); 
+     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);
+     translation.SetTranslation(vrt_displ);
+                                  //the rotation and translation are combined in a single transformation
+     gp_Trsf new_Transf = translation*prev_Transf;
+     TopLoc_Location new_L(new_Transf);
+
+                                  //the transformation is applied to the two sides of the boat
+     sh.Location(new_L);
+     refl_sh.Location(new_L);
+     right_transom_edge.Location(new_L);
+     left_transom_edge.Location(new_L);
+     keel_edge.Location(new_L);
+     current_loc = new_L;
 
 
 }
diff --git a/source/free_surface.cc b/source/free_surface.cc
@@ -279,6 +279,8 @@ void FreeSurface<dim>::reinit() {
 template <int dim>
 void FreeSurface<dim>::initial_conditions(Vector<double> &dst) {
 
+  comp_dom.boat_model.set_current_position(0.03);
+
   // da trattare il caso di dumped solution (continuazione conto salvato)
   initial_time = 0.0;
   initial_wave_shape.set_time(initial_time);
@@ -295,8 +297,7 @@ void FreeSurface<dim>::initial_conditions(Vector<double> &dst) {
       Vinf(i) = instantWindValue(i);
   std::cout<<std::endl<<"Initial conditions: simulation time= "<<initial_time<<"   Vinf= ";
   instantWindValue.print(std::cout,4,false,true);
-  std::cout<<std::endl;
- cout<<"Check: "<<endl;
+
   std::vector<Point<dim> > support_points(comp_dom.dh.n_dofs());
   DoFTools::map_dofs_to_support_points<dim-1, dim>( *comp_dom.mapping, comp_dom.dh, support_points);
   std::vector<Point<dim> > vector_support_points(comp_dom.vector_dh.n_dofs());
@@ -322,7 +323,7 @@ void FreeSurface<dim>::initial_conditions(Vector<double> &dst) {
   for(unsigned int i=0; i<comp_dom.vector_dh.n_dofs(); ++i)
       {
       comp_dom.map_points(i) -= geom_res(i);
-      //std::cout<<i<<" "<<dst(i)<<std::endl;
+      //std::cout<<i<<" "<<geom_res(i)<<"  | "<<comp_dom.map_points(i)<<std::endl;
       } 
 
 
@@ -2673,7 +2674,7 @@ setup_jacobian_prec(t,y,yp,0.0);
          				   //this is the horizontal plane
             Handle(Geom_Plane) horPlane = new Geom_Plane(0.,0.,1.,-dP0(2));
             Handle(Geom_Curve) curve;
-            TopLoc_Location L = comp_dom.boat_model.reference_loc;
+            TopLoc_Location L = comp_dom.boat_model.current_loc;
             TopLoc_Location L_inv = L.Inverted();
             horPlane->Transform(L_inv.Transformation());
             if (comp_dom.boat_model.is_transom)
@@ -4192,14 +4193,18 @@ int FreeSurface<dim>::residual_and_jacobian(const double t,
         for (unsigned int k=3; k<7; ++k)
             { 
             unsigned int i = comp_dom.moving_point_ids[k];
+            //cout<<"Moving point id: "<<i<<endl;
             {
             Point <3> dP0 = comp_dom.support_points[i];
             Point <3> dP;
          				   //this is the horizontal plane
             Handle(Geom_Plane) horPlane = new Geom_Plane(0.,0.,1.,-dP0(2));
             Handle(Geom_Curve) curve;
-            TopLoc_Location L = comp_dom.boat_model.reference_loc;
+            TopLoc_Location L = comp_dom.boat_model.current_loc;
             TopLoc_Location L_inv = L.Inverted();
+            //gp_Pnt test_point(0.0,0.0,0.0);
+            //test_point.Transform(L_inv.Transformation());
+            //cout<<"####### "<<Pnt(test_point)<<endl;
             horPlane->Transform(L_inv.Transformation());
             if (comp_dom.boat_model.is_transom)
                {
diff --git a/source/numerical_towing_tank.cc b/source/numerical_towing_tank.cc
@@ -20,6 +20,7 @@
 #include "../include/boat_surface.h"
 #include <deal.II/numerics/matrix_tools.h>
 #include <deal.II/grid/grid_refinement.h>
+#include <BRepAdaptor_Curve.hxx>
 
 using namespace dealii;
 using namespace OpenCascade;
@@ -1651,15 +1652,15 @@ void NumericalTowingTank::initialize_smoother()
 
   curves.resize(7);
   on_curve_option.resize(7);
-  ref_locations.resize(7);
+  smoothers_locations.resize(7,NULL);
   
   // parameters for front keel smoothing
   base_points[0]   = boat_model.PointFrontBot;
   moving_points[0] = boat_model.PointFrontTop;
   curves[0] = boat_model.equiv_keel_bspline;
   boundary_ids[0] = 30;
   on_curve_option[0] = true;
-  ref_locations[0] = boat_model.reference_loc;
+  smoothers_locations[0] = &boat_model.current_loc;
 
   // parameters for rear keel/left transom smoothing
   if (boat_model.is_transom)
@@ -1669,7 +1670,7 @@ void NumericalTowingTank::initialize_smoother()
      curves[1] = boat_model.left_transom_bspline;
      boundary_ids[1] = 32;
      on_curve_option[1] = true;
-     ref_locations[1] = boat_model.reference_loc;
+     smoothers_locations[1] = &boat_model.current_loc;
      }
   else
      {
@@ -1678,7 +1679,7 @@ void NumericalTowingTank::initialize_smoother()
      curves[1] = boat_model.equiv_keel_bspline;
      boundary_ids[1] = 32;
      on_curve_option[1] = true;
-     ref_locations[1] = boat_model.reference_loc;
+     smoothers_locations[1] = &boat_model.current_loc;
      }
 
   // parameters for rear keel/right transom smoothing
@@ -1689,7 +1690,7 @@ void NumericalTowingTank::initialize_smoother()
      curves[2] = boat_model.right_transom_bspline;
      boundary_ids[2] = 37;
      on_curve_option[2] = true;
-     ref_locations[2] = boat_model.reference_loc;
+     smoothers_locations[2] = &boat_model.current_loc;
      }
   else
      {
@@ -1698,7 +1699,7 @@ void NumericalTowingTank::initialize_smoother()
      curves[2] = boat_model.equiv_keel_bspline;
      boundary_ids[2] = 32;
      on_curve_option[2] = true;
-     ref_locations[2] = boat_model.reference_loc;
+     smoothers_locations[2] = &boat_model.current_loc;
      }
 
   // parameters for front right water line smoothing
@@ -1764,7 +1765,7 @@ void NumericalTowingTank::initialize_smoother()
       double smoother_tolerance = boat_model.boatWetLength*1e-3;
       line_smoothers[i] = new LineSmoothing(smoothing_map_points,
 	  				    curves[i],
-                                            ref_locations[i],
+                                            smoothers_locations[i],
                                             vector_dh,
 					    boundary_dofs[i],
 					    base_point_ids[i],
@@ -1786,6 +1787,7 @@ void NumericalTowingTank::update_smoother()
   curves[0] = boat_model.equiv_keel_bspline;
   boundary_ids[0] = 30;
   on_curve_option[0] = true;
+  smoothers_locations[0] = &boat_model.current_loc;
 
   // parameters for rear keel/left transom smoothing
   if (boat_model.is_transom)
@@ -1795,6 +1797,7 @@ void NumericalTowingTank::update_smoother()
      curves[1] = boat_model.left_transom_bspline;
      boundary_ids[1] = 32;
      on_curve_option[1] = true;
+     smoothers_locations[1] = &boat_model.current_loc;
      }
   else
      {
@@ -1803,6 +1806,7 @@ void NumericalTowingTank::update_smoother()
      curves[1] = boat_model.equiv_keel_bspline;
      boundary_ids[1] = 32;
      on_curve_option[1] = true;
+     smoothers_locations[1] = &boat_model.current_loc;
      }
 
   // parameters for rear keel/right transom smoothing
@@ -1813,6 +1817,7 @@ void NumericalTowingTank::update_smoother()
      curves[2] = boat_model.right_transom_bspline;
      boundary_ids[2] = 37;
      on_curve_option[2] = true;
+     smoothers_locations[2] = &boat_model.current_loc;
      }
   else
      {
@@ -1821,6 +1826,7 @@ void NumericalTowingTank::update_smoother()
      curves[2] = boat_model.equiv_keel_bspline;
      boundary_ids[2] = 32;
      on_curve_option[2] = true;
+     smoothers_locations[2] = &boat_model.current_loc;
      }
 
   // parameters for front right water line smoothing
@@ -2121,48 +2127,60 @@ vector_constraints.distribute(smoothing_map_points);
 //this takes care of the bow and stern nodes 
 if (only_surf_smoothing == false)
 for (unsigned int k=3; k<7; ++k)
-    { 
-    unsigned int i = moving_point_ids[k];
-    //cout<<k<<" "<<i<<" "<<support_points[i]<<endl;
-         {
-         Point <3> dP0 = support_points[i];
-         Point <3> dP; 
-         				   //this is the horizontal plane
-         Handle(Geom_Plane) horPlane = new Geom_Plane(0.,0.,1.,-dP0(2));
-         Handle(Geom_Curve) curve;
-         if (boat_model.is_transom)
+            { 
+            unsigned int i = moving_point_ids[k];
             {
-            if (k==3 || k==4)
-               curve = boat_model.equiv_keel_bspline;
-            else if (k == 6)
-               curve = boat_model.left_transom_bspline;
+            Point <3> dP0 = support_points[i];
+            Point <3> dP;
+         				   //this is the horizontal plane
+            Handle(Geom_Plane) horPlane = new Geom_Plane(0.,0.,1.,-dP0(2));
+            Handle(Geom_Curve) curve;
+            TopLoc_Location L = boat_model.current_loc;
+            TopLoc_Location L_inv = L.Inverted();
+            
+            horPlane->Transform(L_inv.Transformation());
+            if (boat_model.is_transom)
+               {
+               if (k==3 || k==4)
+                  curve = boat_model.equiv_keel_bspline;
+               else if (k == 6)
+                  curve = boat_model.left_transom_bspline;
+               else
+                  curve = boat_model.right_transom_bspline;
+               }
             else
-               curve = boat_model.right_transom_bspline;
-            }
-         else
-            {
-            curve = boat_model.equiv_keel_bspline;
-            }
-         GeomAPI_IntCS Intersector(curve, horPlane);
-         int npoints = Intersector.NbPoints();
-         AssertThrow((npoints != 0), ExcMessage("Keel or transom curve is not intersecting with horizontal plane!"));
-         //cout<<"Number of intersections: "<<npoints<<endl;
-         double minDistance=1e7;
-         gp_Pnt P;
-         gp_Vec V1;
-         double t,u,v;
-         for (int j=0; j<npoints;++j)
-             {
-             Point<3> inters = Pnt(Intersector.Point(j+1));
-             Intersector.Parameters(j+1,u,v,t);
+               {
+               curve = boat_model.equiv_keel_bspline;
+               }
 
-             if (dP0.distance(inters) < minDistance)
+            TopoDS_Edge edge = BRepBuilderAPI_MakeEdge(curve);
+            edge.Location(L);
+            BRepAdaptor_Curve AC(edge);
+            gp_Pnt P;
+            gp_Vec V1;
+            GeomAPI_IntCS Intersector(curve, horPlane);
+            int npoints = Intersector.NbPoints();
+
+            AssertThrow((npoints != 0), ExcMessage("Keel or transom curve is not intersecting with horizontal plane!"));
+            double minDistance=1e7;
+            double t,u,v;
+            for (int j=0; j<npoints;++j)
                 {
-                minDistance = dP0.distance(inters);
-                dP = inters;
-                curve->D1(t,P,V1);
+                gp_Pnt int_point = Intersector.Point(j+1);
+                int_point.Transform(L.Transformation());
+                Point<3> inters = Pnt(int_point);
+                Intersector.Parameters(j+1,u,v,t);
+                if (dP0.distance(inters) < minDistance)
+                   {
+                   minDistance = dP0.distance(inters);
+                   dP = inters;
+                   AC.D1(t,P,V1);
+                   }
                 }
-             }
+            //cout<<"Check plane-curve intersection:"<<endl;
+            //cout<<"Origin: "<<dP0<<"   Proj: "<<dP<<"  dist: "<<minDistance<<endl;
+            //cout<<Pnt(P)<<endl;
+          /*
           // here temporarily for kcs hull tests
             if (minDistance > 0.5*boat_model.boatWetLength)
                {
@@ -2189,12 +2207,12 @@ for (unsigned int k=3; k<7; ++k)
                   }
                cout<<"New part two: "<<dP<<" | "<<V1.X()<<" "<<V1.Y()<<" "<<V1.Z()<<" | "<<dP0<<endl;
                }
-         
+         */
          std::set<unsigned int> duplicates = double_nodes_set[i];
           //duplicates.erase(i); 
           for (std::set<unsigned int>::iterator pos = duplicates.begin(); pos !=duplicates.end(); pos++)
               {
-	      smoothing_map_points(3*(*pos)) = dP(0)-ref_points[3*(*pos)](0);
+              smoothing_map_points(3*(*pos)) = dP(0)-ref_points[3*(*pos)](0);
               smoothing_map_points(3*(*pos)+1) = dP(1)-ref_points[3*(*pos)](1);
               smoothing_map_points(3*(*pos)+2) = dP(2)-ref_points[3*(*pos)](2);
               edges_tangents[3*(*pos)] = V1.X();
@@ -2225,7 +2243,6 @@ for (unsigned int i=0; i<vector_dh.n_dofs(); ++i)
     }
              
 //////////////////////////////////////////
-
                              // line smoothing on keel/transom is performed here and modifies smoothing_map_points
                               // it ONLY moves nodes on the LEFT side of the keel
 for(unsigned int i=0; i<3; ++i)
@@ -2249,7 +2266,6 @@ for(unsigned int i=0; i<3; ++i)
           
           }
    }
-
                               // line smoothing on water_lines is performed here and modifies smoothing_map_points
                               // it ONLY moves nodes on the water that are also near_boat
 if (only_surf_smoothing == false)
diff --git a/waveBem.prm b/waveBem.prm
