Xiangyu/revise sycl diffusion

src/for_3D_build/bodies/mesh_helper_3d.cpp

@@ -403,111 +403,6 @@ void MeshFileHelpers::minimumDistance(StdVec<Real> &all_data_of_distance_between
     }
 }
 
-void MeshFileHelpers::vtuFileHeader(std::ofstream &out_file)
-{
-    out_file << "<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" byte_order=\"LittleEndian\" header_type=\"UInt64\">\n";
-    out_file << "<UnstructuredGrid>\n";
-    out_file << "<FieldData>\n";
-    out_file << "<DataArray type=\"Int32\" Name=\"ispatch\" NumberOfTuples=\"1\" format=\"ascii\" RangeMin=\"0\" RangeMax=\"0\">\n";
-    out_file << "0\n";
-    out_file << "</DataArray>\n";
-    out_file << "</FieldData>\n";
-}
-
-void MeshFileHelpers::vtuFileNodeCoordinates(std::ofstream &out_file, StdLargeVec<Vecd> &node_coordinates_,
-                                             StdLargeVec<StdVec<size_t>> &elements_nodes_connection_, SPHBody &bounds_, Real &rangemax)
-{
-    // Write point data
-    out_file << "<Piece NumberOfPoints=\"" << node_coordinates_.size() << "\" NumberOfCells=\"" << elements_nodes_connection_.size() << "\">\n";
-    out_file << "<PointData>\n";
-    out_file << "</PointData>\n";
-    out_file << "<CellData>\n";
-    out_file << "</CellData>\n";
-    out_file << "<Points>\n";
-    BoundingBox bounds = bounds_.getSPHSystemBounds();
-    Real first_max = bounds.first_.cwiseAbs().maxCoeff();
-    Real second_max = bounds.second_.cwiseAbs().maxCoeff();
-    rangemax = 1.03075 * (std::max(first_max, second_max));
-    out_file << "<DataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"3\" format=\"ascii\" RangeMin=\"0\" RangeMax=\"" << rangemax << "\">\n";
-
-    size_t total_nodes = node_coordinates_.size();
-    for (size_t node = 0; node != total_nodes; ++node)
-    {
-        out_file << node_coordinates_[node][0] << " " << node_coordinates_[node][1] << " " << node_coordinates_[node][2] << " \n";
-    }
-}
-
-void MeshFileHelpers::vtuFileInformationKey(std::ofstream &out_file, Real &rangemax)
-{
-    out_file << "<InformationKey name=\"L2_NORM_RANGE\" location=\"vtkDataArray\" length=\"2\">\n";
-    out_file << "<Value index=\"0\">\n";
-    out_file << "0\n";
-    out_file << "</Value>\n";
-    out_file << "<Value index=\"1\">\n";
-    out_file << rangemax << " \n";
-    out_file << "</Value>\n";
-    out_file << "</InformationKey>\n";
-    out_file << "<InformationKey name=\"L2_NORM_FINITE_RANGE\" location=\"vtkDataArray\" length=\"2\">\n";
-    out_file << "<Value index=\"0\">\n";
-    out_file << "0\n";
-    out_file << "</Value>\n";
-    out_file << "<Value index=\"1\">\n";
-    out_file << rangemax << " \n";
-    out_file << "</Value>\n";
-    out_file << "</InformationKey>\n";
-    out_file << "</DataArray>\n";
-    out_file << "</Points>\n";
-}
-
-void MeshFileHelpers::vtuFileCellConnectivity(std::ofstream &out_file, StdLargeVec<StdVec<size_t>> &elements_nodes_connection_, StdLargeVec<Vecd> &node_coordinates_)
-{
-    out_file << "<Cells>\n";
-    // Write Cell data
-    out_file << "<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\" RangeMin=\"0\" RangeMax=\"" << node_coordinates_.size() - 1 << "\">\n";
-
-    for (const auto &cell : elements_nodes_connection_)
-    {
-        for (const auto &vertex : cell)
-        {
-            out_file << vertex << " ";
-        }
-        out_file << "\n";
-    }
-
-    out_file << "</DataArray>\n";
-}
-
-void MeshFileHelpers::vtuFileOffsets(std::ofstream &out_file, StdLargeVec<StdVec<size_t>> &elements_nodes_connection_)
-{
-    out_file << "<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\" RangeMin=\"4\" RangeMax=\"" << 4 * elements_nodes_connection_.size() << "\">\n";
-
-    size_t offset = 0;
-    for (const auto &face : elements_nodes_connection_)
-    {
-        offset += face.size();
-        out_file << offset << " ";
-    }
-    out_file << "\n</DataArray>\n";
-}
-
-void MeshFileHelpers::vtuFileTypeOfCell(std::ofstream &out_file, StdLargeVec<StdVec<size_t>> &elements_nodes_connection_)
-{
-    size_t type = 10;
-    // Specifies type of cell 10 = tetrahedral
-    out_file << "<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\" RangeMin=\"10\" RangeMax=\"10\">\n";
-    for (const auto &types : elements_nodes_connection_)
-    {
-        for (size_t i = 0; i < types.size(); ++i)
-        {
-            out_file << type << " ";
-        }
-        out_file << "\n";
-    }
-    // Write face attribute data
-    out_file << "</DataArray>\n";
-    out_file << "</Cells>\n";
-}
-
 void MeshFileHelpers::numberOfNodesFluent(std::ifstream &mesh_file, size_t &number_of_points, std::string &text_line)
 {
 

src/for_3D_build/geometries/triangle_mesh_shape.cpp

@@ -3,22 +3,43 @@
 namespace SPH
 {
 //=================================================================================================//
-SimTK::ContactGeometry::TriangleMesh *TriangleMeshShape::generateTriangleMesh(const SimTK::PolygonalMesh &poly_mesh)
+TriangleMeshShape::TriangleMeshShape(const std::string &shape_name)
+    : Shape(shape_name), triangle_mesh_(nullptr) {}
+//=================================================================================================//
+void TriangleMeshShape::initializeFromPolygonalMesh(const SimTK::PolygonalMesh &poly_mesh)
 {
-    SimTK::ContactGeometry::TriangleMesh *triangle_mesh;
-    triangle_mesh = triangle_mesh_ptr_keeper_.createPtr<SimTK::ContactGeometry::TriangleMesh>(poly_mesh);
-    if (!SimTK::ContactGeometry::TriangleMesh::isInstance(*triangle_mesh))
+    triangle_mesh_ = triangle_mesh_ptr_keeper_.createPtr<TriangleMesh>(poly_mesh);
+    if (!TriangleMesh::isInstance(*triangle_mesh_))
     {
         std::cout << "\n Error the triangle mesh is not valid" << std::endl;
         std::cout << __FILE__ << ':' << __LINE__ << std::endl;
         throw;
     }
-    std::cout << "num of faces:" << triangle_mesh->getNumFaces() << std::endl;
+    std::cout << "TriangleMesh:" << name_ << std::endl;
+    std::cout << "num of vertices:" << triangle_mesh_->getNumVertices() << std::endl;
+    std::cout << "num of faces:" << triangle_mesh_->getNumFaces() << std::endl;
+
+    std::vector<std::array<Real, 3>> vertices;
+    vertices.reserve(triangle_mesh_->getNumVertices());
+    for (int i = 0; i < triangle_mesh_->getNumVertices(); i++)
+    {
+        const auto &p = triangle_mesh_->getVertexPosition(i);
+        vertices.push_back({Real(p[0]), Real(p[1]), Real(p[2])});
+    }
 
-    return triangle_mesh;
+    std::vector<std::array<int, 3>> faces;
+    faces.reserve(triangle_mesh_->getNumFaces());
+    for (int i = 0; i < triangle_mesh_->getNumFaces(); i++)
+    {
+        auto f1 = triangle_mesh_->getFaceVertex(i, 0);
+        auto f2 = triangle_mesh_->getFaceVertex(i, 1);
+        auto f3 = triangle_mesh_->getFaceVertex(i, 2);
+        faces.push_back({f1, f2, f3});
+    }
+    triangle_mesh_distance_.construct(vertices, faces);
 }
 //=================================================================================================//
-SimTK::ContactGeometry::TriangleMesh *TriangleMeshShape::getTriangleMesh()
+TriangleMesh *TriangleMeshShape::getTriangleMesh()
 {
     if (triangle_mesh_ == nullptr)
     {
@@ -31,51 +52,14 @@ SimTK::ContactGeometry::TriangleMesh *TriangleMeshShape::getTriangleMesh()
 //=================================================================================================//
 bool TriangleMeshShape::checkContain(const Vec3d &probe_point, bool BOUNDARY_INCLUDED)
 {
-    SimTKVec2 uv_coordinate;
-    bool inside = false; // note that direct prediction is not reliable sometime.
-    int face_id;
-
-    SimTKVec3 closest_pnt = triangle_mesh_->findNearestPoint(EigenToSimTK(probe_point), inside, face_id, uv_coordinate);
-    Vec3d from_face_to_pnt = probe_point - SimTKToEigen(closest_pnt);
-    Real distance_to_pnt = from_face_to_pnt.norm();
-    Vec3d direction_to_pnt = from_face_to_pnt / (distance_to_pnt + TinyReal);
-    Vec3d face_normal = SimTKToEigen(triangle_mesh_->getFaceNormal(face_id));
-    Real cosine_angle = face_normal.dot(direction_to_pnt);
-
-    int ite = 0;
-    while (fabs(cosine_angle) < Eps)
-    {
-        Vec3d jittered = probe_point; // jittering
-        for (int l = 0; l != probe_point.size(); ++l)
-            jittered[l] = probe_point[l] + rand_uniform(-0.5, 0.5) * (SqrtEps + distance_to_pnt * 0.1);
-        Vec3d from_face_to_jittered = jittered - SimTKToEigen(closest_pnt);
-        Vec3d direction_to_jittered = from_face_to_jittered / (from_face_to_jittered.norm() + TinyReal);
-        cosine_angle = face_normal.dot(direction_to_jittered);
-
-        ite++;
-        if (ite > 100)
-        {
-            std::cout << "\n Error: TriangleMeshShape::checkContain not able to achieve!  " << std::endl;
-            std::cout << __FILE__ << ':' << __LINE__ << std::endl;
-            exit(1);
-        }
-    }
+    Real distance = triangle_mesh_distance_.signed_distance(probe_point).distance;
 
-    return cosine_angle < 0.0 ? true : false;
+    return distance < 0.0 ? true : false;
 }
 //=================================================================================================//
 Vecd TriangleMeshShape::findClosestPoint(const Vecd &probe_point)
 {
-    bool inside = false;
-    int face_id;
-    SimTKVec2 norm;
-    SimTKVec3 closest_pnt = triangle_mesh_->findNearestPoint(SimTKVec3(probe_point[0], probe_point[1], probe_point[2]), inside, face_id, norm);
-    if (face_id < 0 && face_id > triangle_mesh_->getNumFaces())
-    {
-        std::cout << "\n Error the nearest point is not valid" << std::endl;
-        std::cout << __FILE__ << ':' << __LINE__ << std::endl;
-        throw;
-    }
+    auto closest_pnt = triangle_mesh_distance_.signed_distance(probe_point).nearest_point;
     return Vecd(closest_pnt[0], closest_pnt[1], closest_pnt[2]);
 }
 //=================================================================================================//
@@ -101,30 +85,35 @@ TriangleMeshShapeBrick::TriangleMeshShapeBrick(Vecd halfsize, int resolution, Ve
                                                const std::string &shape_name)
     : TriangleMeshShape(shape_name)
 {
-    SimTK::PolygonalMesh polymesh = SimTK::PolygonalMesh::createBrickMesh(SimTKVec3(halfsize[0], halfsize[1], halfsize[2]), resolution);
-    triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(SimTKVec3(translation[0], translation[1], translation[2])));
+    SimTK::PolygonalMesh poly_mesh = SimTK::PolygonalMesh::createBrickMesh(
+        SimTKVec3(halfsize[0], halfsize[1], halfsize[2]), resolution);
+    initializeFromPolygonalMesh(
+        poly_mesh.transformMesh(SimTKVec3(translation[0], translation[1], translation[2])));
 }
 //=================================================================================================//
-TriangleMeshShapeBrick::TriangleMeshShapeBrick(const TriangleMeshShapeBrick::ShapeParameters &shape_parameters,
-                                               const std::string &shape_name)
+TriangleMeshShapeBrick::TriangleMeshShapeBrick(
+    const TriangleMeshShapeBrick::ShapeParameters &shape_parameters, const std::string &shape_name)
     : TriangleMeshShapeBrick(shape_parameters.halfsize_, shape_parameters.resolution_,
                              shape_parameters.translation_, shape_name) {}
 //=================================================================================================//
 TriangleMeshShapeSphere::TriangleMeshShapeSphere(Real radius, int resolution, Vec3d translation,
                                                  const std::string &shape_name)
     : TriangleMeshShape(shape_name)
 {
-    SimTK::PolygonalMesh polymesh = SimTK::PolygonalMesh::createSphereMesh(radius, resolution);
-    triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(SimTKVec3(translation[0], translation[1], translation[2])));
+    SimTK::PolygonalMesh poly_mesh = SimTK::PolygonalMesh::createSphereMesh(radius, resolution);
+    initializeFromPolygonalMesh(
+        poly_mesh.transformMesh(SimTKVec3(translation[0], translation[1], translation[2])));
 }
 //=================================================================================================//
-TriangleMeshShapeCylinder::TriangleMeshShapeCylinder(SimTK::UnitVec3 axis, Real radius, Real halflength, int resolution,
-                                                     Vecd translation, const std::string &shape_name)
+TriangleMeshShapeCylinder::TriangleMeshShapeCylinder(
+    SimTK::UnitVec3 axis, Real radius, Real halflength, int resolution,
+    Vecd translation, const std::string &shape_name)
     : TriangleMeshShape(shape_name)
 {
-    SimTK::PolygonalMesh polymesh =
+    SimTK::PolygonalMesh poly_mesh =
         SimTK::PolygonalMesh::createCylinderMesh(axis, radius, halflength, resolution);
-    triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(SimTKVec3(translation[0], translation[1], translation[2])));
+    initializeFromPolygonalMesh(
+        poly_mesh.transformMesh(SimTKVec3(translation[0], translation[1], translation[2])));
 }
 //=================================================================================================//
 TriangleMeshShapeSTL::TriangleMeshShapeSTL(const std::string &filepathname, Vec3d translation,
@@ -137,43 +126,11 @@ TriangleMeshShapeSTL::TriangleMeshShapeSTL(const std::string &filepathname, Vec3
         std::cout << __FILE__ << ':' << __LINE__ << std::endl;
         throw;
     }
-    SimTK::PolygonalMesh polymesh;
-    polymesh.loadStlFile(filepathname);
-    polymesh.scaleMesh(scale_factor);
-    polymesh.transformMesh(SimTKVec3(translation[0], translation[1], translation[2]));
-    triangle_mesh_ = generateTriangleMesh(polymesh);
-
-    std::vector<std::array<Real, 3>> vertices;
-    vertices.reserve(polymesh.getNumVertices());
-    for (int i = 0; i < polymesh.getNumVertices(); i++)
-    {
-        const auto &p = polymesh.getVertexPosition(i);
-        vertices.push_back({Real(p[0]), Real(p[1]), Real(p[2])});
-    }
-
-    std::vector<std::array<int, 3>> faces;
-    faces.reserve(polymesh.getNumFaces());
-    for (int i = 0; i < polymesh.getNumFaces(); i++)
-    {
-        auto f1 = polymesh.getFaceVertex(i, 0);
-        auto f2 = polymesh.getFaceVertex(i, 1);
-        auto f3 = polymesh.getFaceVertex(i, 2);
-        faces.push_back({f1, f2, f3});
-    }
-    triangle_mesh_distance_.construct(vertices, faces);
-}
-//=================================================================================================//
-bool TriangleMeshShapeSTL::checkContain(const Vec3d &probe_point, bool BOUNDARY_INCLUDED)
-{
-    Real distance = triangle_mesh_distance_.signed_distance(probe_point).distance;
-
-    return distance < 0.0 ? true : false;
-}
-//=================================================================================================//
-Vecd TriangleMeshShapeSTL::findClosestPoint(const Vecd &probe_point)
-{
-    auto closest_pnt = triangle_mesh_distance_.signed_distance(probe_point).nearest_point;
-    return Vecd(closest_pnt[0], closest_pnt[1], closest_pnt[2]);
+    SimTK::PolygonalMesh poly_mesh;
+    poly_mesh.loadStlFile(filepathname);
+    poly_mesh.scaleMesh(scale_factor);
+    poly_mesh.transformMesh(SimTKVec3(translation[0], translation[1], translation[2]));
+    initializeFromPolygonalMesh(poly_mesh);
 }
 //=================================================================================================//
 } // namespace SPH