From 3112bcdfb22616504eea574d6caf538bb4ba8153 Mon Sep 17 00:00:00 2001 From: calderans Date: Fri, 5 May 2023 17:55:38 +0200 Subject: [PATCH 1/7] - adding new MorphMesh component --- CMakeLists.txt | 9 ++ morphMesh/CMakeLists.txt | 53 +++++++++ morphMesh/README.md | 3 + morphMesh/inc/gmds/morphMesh/MorphMesh.h | 60 ++++++++++ morphMesh/src/MorphMesh.cpp | 134 +++++++++++++++++++++++ morphMesh/src/main.cpp | 36 ++++++ 6 files changed, 295 insertions(+) create mode 100644 morphMesh/CMakeLists.txt create mode 100644 morphMesh/README.md create mode 100644 morphMesh/inc/gmds/morphMesh/MorphMesh.h create mode 100644 morphMesh/src/MorphMesh.cpp create mode 100644 morphMesh/src/main.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index e34e8a426..42d27e084 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -251,6 +251,15 @@ GMDS_ADD_COMPONENT( ON # must be covered ) +GMDS_ADD_COMPONENT( + MORPHMESH + morphMesh + GMDSmorphMesh + "morphing the mesh" + ON + OFF # must be covered +) + #============================================================================== #mettre les find package uniquement ou necessaire diff --git a/morphMesh/CMakeLists.txt b/morphMesh/CMakeLists.txt new file mode 100644 index 000000000..3c698f147 --- /dev/null +++ b/morphMesh/CMakeLists.txt @@ -0,0 +1,53 @@ +#============================================================================== +# LIBRARY DEFINTION (SOURCE FILES) +#============================================================================== +# Nommer tout en GMDS_MODULE_NAME, GMDS_SRC, ... dans les composants +set(GMDS_LIB ${LIB_GMDS_MORPHMESH}) +set(GMDS_LIB_PREFIX gmds/morphMesh) + +set(GMDS_INC + ${CMAKE_BINARY_DIR}/exports/${GMDS_LIB}_export.h + inc/gmds/morphMesh/MorphMesh.h + ) +set(GMDS_SRC + src/MorphMesh.cpp + src/main.cpp + ) +#============================================================================== +add_library(${GMDS_LIB} ${GMDS_INC} ${GMDS_SRC}) +#============================================================================== +include(GenerateExportHeader) +generate_export_header(${GMDS_LIB} + EXPORT_FILE_NAME ${CMAKE_BINARY_DIR}/exports/${GMDS_LIB}_export.h + EXPORT_MACRO_NAME ${GMDS_LIB}_API) +#============================================================================== +# TARGET DEFINITION +#============================================================================== +include(GNUInstallDirs) + +#LIBRARY TO INSTALL +target_link_libraries(${GMDS_LIB} PUBLIC + ${LIB_GMDS_IG} + ${LIB_GMDS_IG_ALGO} + ${LIB_GMDS_IO} ANN) + +target_compile_features(${GMDS_LIB} PUBLIC cxx_std_14) + +# INCLUDE TO INSTALL +target_include_directories(${GMDS_LIB} PUBLIC + $ + $ + ) +set_target_properties(${GMDS_LIB} PROPERTIES PUBLIC_HEADER "${GMDS_INC}") + +install(TARGETS ${GMDS_LIB} + EXPORT GMDS_SUITE + PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${GMDS_LIB_PREFIX} + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) +#============================================================================== +add_executable(morphing src/main.cpp) +target_link_libraries(morphing PRIVATE ${GMDS_LIB}) +target_compile_features(morphing PUBLIC cxx_std_14) +install(TARGETS morphing) \ No newline at end of file diff --git a/morphMesh/README.md b/morphMesh/README.md new file mode 100644 index 000000000..c494bd762 --- /dev/null +++ b/morphMesh/README.md @@ -0,0 +1,3 @@ +#morphMesh module + +This module propose algorithm to morph mesh using deformation function. \ No newline at end of file diff --git a/morphMesh/inc/gmds/morphMesh/MorphMesh.h b/morphMesh/inc/gmds/morphMesh/MorphMesh.h new file mode 100644 index 000000000..b42a7185d --- /dev/null +++ b/morphMesh/inc/gmds/morphMesh/MorphMesh.h @@ -0,0 +1,60 @@ +/*----------------------------------------------------------------------------*/ +#ifndef GMDS_MORPHMESH_H +#define GMDS_MORPHMESH_H + +/*----------------------------------------------------------------------------*/ +#include "LIB_GMDS_MORPHMESH_export.h" +#include "gmds/ig/Mesh.h" +/*----------------------------------------------------------------------------*/ +namespace gmds { +/*----------------------------------------------------------------------------*/ +namespace morphmesh { +/*----------------------------------------------------------------------------*/ +class LIB_GMDS_MORPHMESH_API MorphMesh +{ + + public: + /*------------------------------------------------------------------------*/ + /** @brief Constructor. + */ + MorphMesh(Mesh* AMesh, const std::vector &APoints, double ARadius); + + /*------------------------------------------------------------------------*/ + /** @brief Destructor. + */ + ~MorphMesh(); + + /*------------------------------------------------------------------------*/ + /** @brief Constructor. + */ + void execute(); + + private: + + /*Surement a retirer plus tard*/ + void markLockedCells(); + + math::Point computeLocalOrigin(const math::Point& AP); + + double findHomothetyRatio(const math::Point& AP, const Node& ANode); + + private: + + /* Locked nodes */ + int m_locked; + /* List of "homothetic points" that will define the mesh deformation */ + std::vector m_targets; + /* The radius of an area that will be modified by a homothetic point*/ + double m_radius; + /* Nodes on the surface of the mesh */ + std::vector m_surfNodes; + + /* The mesh to deform */ + Mesh* m_mesh; +}; +/*----------------------------------------------------------------------------*/ +}// namespace morphmesh +/*----------------------------------------------------------------------------*/ +}// namespace gmds +/*----------------------------------------------------------------------------*/ +#endif // GMDS_MORPHMESH_H diff --git a/morphMesh/src/MorphMesh.cpp b/morphMesh/src/MorphMesh.cpp new file mode 100644 index 000000000..e805ce495 --- /dev/null +++ b/morphMesh/src/MorphMesh.cpp @@ -0,0 +1,134 @@ +/*----------------------------------------------------------------------------*/ +#include "gmds/morphMesh/MorphMesh.h" +/*----------------------------------------------------------------------------*/ +namespace gmds { +/*----------------------------------------------------------------------------*/ +namespace morphmesh { +/*----------------------------------------------------------------------------*/ +MorphMesh::MorphMesh(Mesh *AMesh, const std::vector &APoints, double ARadius): m_mesh(AMesh) +{ + m_radius = ARadius; + m_targets = APoints; + m_locked = m_mesh->newMark(); +} +/*----------------------------------------------------------------------------*/ +MorphMesh::~MorphMesh() {} +/*----------------------------------------------------------------------------*/ +void +MorphMesh::execute() +{ + markLockedCells(); + + for(auto n : m_mesh->nodes()){ + Node node = m_mesh->get(n); + if(node.get().size() == 4) + m_surfNodes.push_back(n); + } + + for(const auto& p : m_targets){ + std::cout<<"Traitement point "<get(n).point(); + double dist = n_point.distance(p); + if(dist < min_dist){ + min_dist = dist; + n0 = n; + } + } + Node node0 = m_mesh->get(n0); + math::Point point0 = node0.point(); + // on calcul le ration d'homothétie pour ce point + double H = findHomothetyRatio(p,node0); + + + for(auto n : m_mesh->nodes()){ + if(!m_mesh->isMarked(n,m_locked)) { + Node node = m_mesh->get(n); + math::Point verifPoint = node.point(); + verifPoint.setY(point0.Y()); // On fait ça pour mettre les deux points sur le même plan + + double dist = verifPoint.distance(point0); + if (dist <= m_radius) { + std::cout << "Dans le radius" << std::endl; + // Si le point est couvert par la zone à modifier alors on continue le traitement + double w = dist / m_radius; + std::cout << "Poids = " << w << std::endl; + + math::Point n_origin = computeLocalOrigin(node.point()); + math::Vector3d height = math::Vector3d(node.point() - n_origin); + + std::cout << "before H = " << H << std::endl; + double localH = H + ((1 - H) * w); // On réduit le ratio d'homothétie par le poids + std::cout << "after H = " << localH << std::endl; + + math::Point result = n_origin + height * localH; + + std::cout << "before " << node << std::endl; + node.setPoint(result); + std::cout << "after " << node << std::endl; + } + } + } + } +} +/*----------------------------------------------------------------------------*/ +void +MorphMesh::markLockedCells() +{ + /** ================================== + * !! WIP !! + ================================= + + 1) this method should also lock nodes given by user + 2) maybe we should lock cells too, edge, faces or region + */ + for(auto n : m_mesh->nodes()){ + Node node = m_mesh->get(n); + if(node.Y() <= 0 ){ + m_mesh->mark(n,m_locked); + } + } +} +/*----------------------------------------------------------------------------*/ +math::Point +MorphMesh::computeLocalOrigin(const math::Point& AP) +{ + /* + * 1) on crée un segment entre de AP et (X,Y,0) + * 2) On intersecte le segment avec les faces lock + * - si on a un point alors on retourne ce point + * - sinon on retourne (X,Y,0) + */ + + return {AP.X(),0,AP.Z()}; +} +/*----------------------------------------------------------------------------*/ +double MorphMesh::findHomothetyRatio(const math::Point &AP, const Node& ANode) +{ + + math::Point p_node = ANode.point(); + // Le vecteur de déformation initial : La translation du noeud du maillage vers le point + math::Vector3d deform(AP - ANode.point()); + + + // On cherche l'origine du noeud initial + math::Point origin = computeLocalOrigin(ANode.point()); + // On fait un vecteur origin->ni pour calculer l'homothétie des deux vecteurs. + // Faire le vecteur dans ce sens là permet de facilement déduire le ration d'homothétie + math::Vector3d height(ANode.point() - origin); + + + double dist1 = origin.distance(ANode.point()); + double dist2 = origin.distance(AP); + + + return dist2/dist1; +} +/*----------------------------------------------------------------------------*/ +} // namespace MorphMesh +/*----------------------------------------------------------------------------*/ +} // namespace gmds +/*----------------------------------------------------------------------------*/ \ No newline at end of file diff --git a/morphMesh/src/main.cpp b/morphMesh/src/main.cpp new file mode 100644 index 000000000..5de49b0c3 --- /dev/null +++ b/morphMesh/src/main.cpp @@ -0,0 +1,36 @@ +/*------------------------------------------------------------------------*/ +#include "gmds/ig/MeshDoctor.h" +#include "gmds/io/IGMeshIOService.h" +#include "gmds/io/VTKReader.h" +#include "gmds/io/VTKWriter.h" +#include +/*------------------------------------------------------------------------*/ +using namespace gmds; +/*------------------------------------------------------------------------*/ +int main(int argc, char* argv[]){ + + std::string param_file(argv[1]); + + gmds::Mesh m_mesh(gmds::MeshModel(gmds::DIM3 | gmds::R | gmds::N | gmds::N2R | gmds::R2N)); + gmds::IGMeshIOService ioService(&m_mesh); + gmds::VTKReader vtkReader(&ioService); + vtkReader.setCellOptions(gmds::N|gmds::R); + vtkReader.setDataOptions(gmds::N|gmds::R); + vtkReader.read(param_file); + + gmds::MeshDoctor doc(&m_mesh); + doc.buildN2R(m_mesh.getModel()); + doc.updateUpwardConnectivity(); + + gmds::math::Point point(5,0.5,0); + std::vector points; + points.push_back(point); + + gmds::morphmesh::MorphMesh morph(&m_mesh,points,2.5); + morph.execute(); + + gmds::VTKWriter w(&ioService); + w.setCellOptions(gmds::N|gmds::R); + w.setDataOptions(gmds::N|gmds::R); + w.write("morphing.vtk"); +} \ No newline at end of file From 53c18f091b49629f3bead8fefae275983ca3dc3b Mon Sep 17 00:00:00 2001 From: calderans Date: Mon, 12 Jun 2023 15:33:37 +0200 Subject: [PATCH 2/7] - big update on morphMesh component --- morphMesh/inc/gmds/morphMesh/MorphMesh.h | 8 +- morphMesh/src/MorphMesh.cpp | 259 ++++++++++++++++++++--- morphMesh/src/main.cpp | 9 +- 3 files changed, 240 insertions(+), 36 deletions(-) diff --git a/morphMesh/inc/gmds/morphMesh/MorphMesh.h b/morphMesh/inc/gmds/morphMesh/MorphMesh.h index b42a7185d..09d76532f 100644 --- a/morphMesh/inc/gmds/morphMesh/MorphMesh.h +++ b/morphMesh/inc/gmds/morphMesh/MorphMesh.h @@ -34,7 +34,8 @@ class LIB_GMDS_MORPHMESH_API MorphMesh /*Surement a retirer plus tard*/ void markLockedCells(); - math::Point computeLocalOrigin(const math::Point& AP); + bool computeLocalOrigin(const math::Point& AP, math::Point& AResult); + bool computeLocalOriginZ(const math::Point& AP, math::Point& AResult); double findHomothetyRatio(const math::Point& AP, const Node& ANode); @@ -42,12 +43,17 @@ class LIB_GMDS_MORPHMESH_API MorphMesh /* Locked nodes */ int m_locked; + Variable* locked_faces; /* List of "homothetic points" that will define the mesh deformation */ std::vector m_targets; /* The radius of an area that will be modified by a homothetic point*/ double m_radius; /* Nodes on the surface of the mesh */ std::vector m_surfNodes; + std::vector m_locked_faces; + //std::map is_locked_node; + //std::map nodes_origins; + /* The mesh to deform */ Mesh* m_mesh; diff --git a/morphMesh/src/MorphMesh.cpp b/morphMesh/src/MorphMesh.cpp index e805ce495..02b1e5ae6 100644 --- a/morphMesh/src/MorphMesh.cpp +++ b/morphMesh/src/MorphMesh.cpp @@ -1,5 +1,9 @@ /*----------------------------------------------------------------------------*/ #include "gmds/morphMesh/MorphMesh.h" +#include "gmds/math/Plane.h" +#include "gmds/math/Triangle.h" +#include + /*----------------------------------------------------------------------------*/ namespace gmds { /*----------------------------------------------------------------------------*/ @@ -10,13 +14,15 @@ MorphMesh::MorphMesh(Mesh *AMesh, const std::vector &APoints, doubl m_radius = ARadius; m_targets = APoints; m_locked = m_mesh->newMark(); + + locked_faces = m_mesh->newVariable("locked_face"); } /*----------------------------------------------------------------------------*/ MorphMesh::~MorphMesh() {} /*----------------------------------------------------------------------------*/ -void -MorphMesh::execute() +void MorphMesh::execute() { + std::cout<<"Marking cells"<nodes()){ @@ -25,8 +31,14 @@ MorphMesh::execute() m_surfNodes.push_back(n); } + std::vector durations1; + std::vector durations2; + + std::cout<<"Start deformation"<get(n0); math::Point point0 = node0.point(); // on calcul le ration d'homothétie pour ce point - double H = findHomothetyRatio(p,node0); + //double H = findHomothetyRatio(p,node0); + double Y0 = point0.Y(); + double X0 = point0.X(); + double Z = (Y0 - p.Y())/2; + math::Point tmp; + computeLocalOrigin(point0, tmp); + double omega_height0 = Y0 - tmp.Y(); + math::Vector3d t; + t.setXYZ(0,p.Y()-Y0,0); + math::Vector3d tz; + tz.setXYZ(0,0,Z); - for(auto n : m_mesh->nodes()){ - if(!m_mesh->isMarked(n,m_locked)) { + for(auto n : m_mesh->nodes()) { + if (!m_mesh->isMarked(n, m_locked)) { Node node = m_mesh->get(n); - math::Point verifPoint = node.point(); - verifPoint.setY(point0.Y()); // On fait ça pour mettre les deux points sur le même plan - double dist = verifPoint.distance(point0); + double dist = node.point().distance(point0); if (dist <= m_radius) { - std::cout << "Dans le radius" << std::endl; + time_t start,end; + time(&start); // Si le point est couvert par la zone à modifier alors on continue le traitement - double w = dist / m_radius; - std::cout << "Poids = " << w << std::endl; + //double w = 1 - (dist / m_radius); + double sigmoid = 1/(1 + exp(-5*( ( (1 - (dist / m_radius) ) *2)-1) ) ); + + math::Point n_origin; - math::Point n_origin = computeLocalOrigin(node.point()); - math::Vector3d height = math::Vector3d(node.point() - n_origin); + bool origin_locked = computeLocalOrigin(node.point(), n_origin); - std::cout << "before H = " << H << std::endl; - double localH = H + ((1 - H) * w); // On réduit le ratio d'homothétie par le poids - std::cout << "after H = " << localH << std::endl; + double omega_height = Y0 - n_origin.Y(); + double n_height = node.point().Y() - n_origin.Y(); - math::Point result = n_origin + height * localH; + math::Vector3d final_t = (n_height / omega_height) * (sigmoid * t); + math::Point result = node.point() + final_t; - std::cout << "before " << node << std::endl; node.setPoint(result); - std::cout << "after " << node << std::endl; + time(&end); + durations1.emplace_back(double(end - start)); + } + + double distX = fabs(node.X() - X0); + if (distX <= m_radius) { + time_t start,end; + time(&start); + // double w = 0.5+2*cbrt((1 - (distX / m_radius)) - 0.5); + double sigmoid = 1 / (1 + exp(-5 * (((1 - (distX / m_radius)) * 2) - 1))); + + math::Point n_origin; + bool origin_locked = computeLocalOriginZ(node.point(), n_origin); + + double omega_height = fabs(Y0) - fabs(n_origin.Z()); + double n_height = fabs(node.point().Z()) - fabs(n_origin.Z()); + + int sens = node.Z() < 0 ? -1 : 1; + math::Vector3d final_t = (n_height / omega_height) * (sigmoid * (sens * tz)); + + math::Point result = node.point() + final_t; + node.setPoint(result); + time(&end); + durations1.emplace_back(double(end - start)); } } } + /*for(auto n : m_mesh->nodes()) { + if (!m_mesh->isMarked(n, m_locked)) { + Node node = m_mesh->get(n); + double distX = fabs(node.X() - X0); + if (distX <= m_radius) { + // double w = 0.5+2*cbrt((1 - (distX / m_radius)) - 0.5); + double sigmoid = 1 / (1 + exp(-5 * (((1 - (distX / m_radius)) * 2) - 1))); + + math::Point n_origin; + bool origin_locked = computeLocalOriginZ(node.point(), n_origin); + + double omega_height = fabs(Y0) - fabs(n_origin.Z()); + double n_height = fabs(node.point().Z()) - fabs(n_origin.Z()); + + int sens = node.Z() < 0 ? -1 : 1; + math::Vector3d final_t = (n_height / omega_height) * (sigmoid * (sens * tz)); + + math::Point result = node.point() + final_t; + node.setPoint(result); + } + } + }*/ + } + + time(&endglob); + + + double avg = durations1[0]; + for(int i = 1; iunmarkAll(m_locked); + m_mesh->freeMark(m_locked); } /*----------------------------------------------------------------------------*/ -void -MorphMesh::markLockedCells() +void MorphMesh::markLockedCells() { - /** ================================== - * !! WIP !! - ================================= + Variable* locked_regions = m_mesh->newVariable("m_locked"); + std::set regions; + std::set nodes; + nodes.insert(27006); + m_mesh->mark(27006,m_locked); + for (int i = 0; i < 15; ++i) { + for(auto n : nodes){ + Node node = m_mesh->get(n); + for(const auto& r : node.get()){ + regions.emplace(r); + } + } + nodes.clear(); + for(const auto& r : regions){ + locked_regions->set(r.id(),1); + for(auto n : r.getIDs()){ + if(!m_mesh->isMarked(n,m_locked)){ + nodes.emplace(n); + m_mesh->mark(n,m_locked); + } + } + } + regions.clear(); + } - 1) this method should also lock nodes given by user - 2) maybe we should lock cells too, edge, faces or region - */ for(auto n : m_mesh->nodes()){ Node node = m_mesh->get(n); if(node.Y() <= 0 ){ m_mesh->mark(n,m_locked); } } + + //Test origine locale décallée + /*for(auto r : m_mesh->regions()){ + //if((4440<=r && r<=4639) || (5140<=r && r<=5339) || (5840<=r && r<=6039) || (6540<=r && r<=6739)){ + if((37000<=r && r<=37999) || (39000<=r && r<=39999) || (41000<=r && r<=41999) || (43000<=r && r<=43999)){ + locked_regions->set(r,1); + Region region = m_mesh->get(r); + for(auto n : region.getIDs()){ + m_mesh->mark(n,m_locked); + } + } + }*/ + + + for(auto f : m_mesh->faces()){ + Face face = m_mesh->get(f); + + std::vector f_regions = face.getIDs(); + if(f_regions.size() == 2){ + //Une face entre un hexa lock et un non lock + if((locked_regions->value(f_regions[0]) == 1) && (locked_regions->value(f_regions[1]) == 0) || + (locked_regions->value(f_regions[0]) == 0) && (locked_regions->value(f_regions[1]) == 1)){ + locked_faces->set(f,1); + m_locked_faces.push_back(f); + } + } + } } /*----------------------------------------------------------------------------*/ -math::Point -MorphMesh::computeLocalOrigin(const math::Point& AP) +bool MorphMesh::computeLocalOrigin(const math::Point& AP, math::Point& AResult) { /* * 1) on crée un segment entre de AP et (X,Y,0) @@ -102,8 +226,76 @@ MorphMesh::computeLocalOrigin(const math::Point& AP) * - si on a un point alors on retourne ce point * - sinon on retourne (X,Y,0) */ + std::vector points; + points.emplace_back(AP.X(),0,AP.Z()); + points.emplace_back(AP.X()+0.001,0,AP.Z()+0.001); + points.emplace_back(AP.X()-0.001,0,AP.Z()-0.001); + + for (auto f : m_locked_faces) { + for(const auto& p : points) { + //if (locked_faces->value(f) == 1) { + math::Segment depth(AP, p); + Face face = m_mesh->get(f); + math::Plane plane(face.center(), face.normal()); + // on cherche si la profondeur a une intersection avec un plan d'une face locked + if (depth.intersect(plane, AResult)) { + // si oui on cherche si le point d'intersection est dans la face + // On fait ça en créant deux triangles qui couvrent entièrement la face + // et qui permettent d'utiliser la méthode isIn() + std::vector nodes = face.get(); + math::Triangle tri1(nodes[0].point(), nodes[1].point(), nodes[2].point()); + math::Triangle tri2(nodes[2].point(), nodes[3].point(), nodes[0].point()); + if (tri1.intersect(depth) || tri2.intersect(depth)) { + return true; + } + } + //} + } + } + AResult = {AP.X(),0,AP.Z()}; + + return false; +} +/*----------------------------------------------------------------------------*/ +bool MorphMesh::computeLocalOriginZ(const math::Point& AP, math::Point& AResult) +{ + /* + * 1) on crée un segment entre de AP et (X,Y,0) + * 2) On intersecte le segment avec les faces lock + * - si on a un point alors on retourne ce point + * - sinon on retourne (X,Y,0) + */ + std::vector points; + points.emplace_back(AP.X(),AP.Y(),0); + points.emplace_back(AP.X()+0.001,AP.Y()+0.001,0); + points.emplace_back(AP.X()-0.001,AP.Y()-0.001,0); + + for (auto f : m_locked_faces) { + for(const auto& p : points) { + //if (locked_faces->value(f) == 1) { + math::Segment depth(AP,math::Point(AP.X(),AP.Y(),0)); + Face face = m_mesh->get(f); + math::Plane plane(face.center(), face.normal()); + // on cherche si la profondeur a une intersection avec un plan d'une face locked + if (depth.intersect(plane, AResult)) { + // std::cout<<"test"< nodes = face.get(); + math::Triangle tri1(nodes[0].point(), nodes[1].point(), nodes[2].point()); + math::Triangle tri2(nodes[2].point(), nodes[3].point(), nodes[0].point()); + if (tri1.intersect(depth) || tri2.intersect(depth)) { + return true; + } + } + //} + } + } + + AResult = {AP.X(),AP.Y(),0}; - return {AP.X(),0,AP.Z()}; + return false; } /*----------------------------------------------------------------------------*/ double MorphMesh::findHomothetyRatio(const math::Point &AP, const Node& ANode) @@ -115,7 +307,8 @@ double MorphMesh::findHomothetyRatio(const math::Point &AP, const Node& ANode) // On cherche l'origine du noeud initial - math::Point origin = computeLocalOrigin(ANode.point()); + math::Point origin; + bool origin_locked = computeLocalOrigin(ANode.point(), origin); // On fait un vecteur origin->ni pour calculer l'homothétie des deux vecteurs. // Faire le vecteur dans ce sens là permet de facilement déduire le ration d'homothétie math::Vector3d height(ANode.point() - origin); @@ -124,6 +317,8 @@ double MorphMesh::findHomothetyRatio(const math::Point &AP, const Node& ANode) double dist1 = origin.distance(ANode.point()); double dist2 = origin.distance(AP); + std::cout<<"dist origin -> node "< point "< points; points.push_back(point); @@ -32,5 +33,7 @@ int main(int argc, char* argv[]){ gmds::VTKWriter w(&ioService); w.setCellOptions(gmds::N|gmds::R); w.setDataOptions(gmds::N|gmds::R); - w.write("morphing.vtk"); + w.write("morphingtest.vtk"); + + } \ No newline at end of file From 33330bc4c7c1d9bce998303d61f2c035aa206e0a Mon Sep 17 00:00:00 2001 From: calderans Date: Fri, 5 May 2023 17:55:38 +0200 Subject: [PATCH 3/7] - adding new MorphMesh component --- CMakeLists.txt | 10 ++ morphMesh/CMakeLists.txt | 53 +++++++++ morphMesh/README.md | 3 + morphMesh/inc/gmds/morphMesh/MorphMesh.h | 60 ++++++++++ morphMesh/src/MorphMesh.cpp | 134 +++++++++++++++++++++++ morphMesh/src/main.cpp | 36 ++++++ 6 files changed, 296 insertions(+) create mode 100644 morphMesh/CMakeLists.txt create mode 100644 morphMesh/README.md create mode 100644 morphMesh/inc/gmds/morphMesh/MorphMesh.h create mode 100644 morphMesh/src/MorphMesh.cpp create mode 100644 morphMesh/src/main.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 2c0a3bf48..fcfddddee 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -279,6 +279,16 @@ GMDS_ADD_COMPONENT( OFF # must be covered ) +GMDS_ADD_COMPONENT( + MORPHMESH + morphMesh + GMDSmorphMesh + "morphing the mesh" + OFF + OFF # must be covered +) + + #============================================================================== set (GMDS_INCLUDE_DIRS APPEND) diff --git a/morphMesh/CMakeLists.txt b/morphMesh/CMakeLists.txt new file mode 100644 index 000000000..3c698f147 --- /dev/null +++ b/morphMesh/CMakeLists.txt @@ -0,0 +1,53 @@ +#============================================================================== +# LIBRARY DEFINTION (SOURCE FILES) +#============================================================================== +# Nommer tout en GMDS_MODULE_NAME, GMDS_SRC, ... dans les composants +set(GMDS_LIB ${LIB_GMDS_MORPHMESH}) +set(GMDS_LIB_PREFIX gmds/morphMesh) + +set(GMDS_INC + ${CMAKE_BINARY_DIR}/exports/${GMDS_LIB}_export.h + inc/gmds/morphMesh/MorphMesh.h + ) +set(GMDS_SRC + src/MorphMesh.cpp + src/main.cpp + ) +#============================================================================== +add_library(${GMDS_LIB} ${GMDS_INC} ${GMDS_SRC}) +#============================================================================== +include(GenerateExportHeader) +generate_export_header(${GMDS_LIB} + EXPORT_FILE_NAME ${CMAKE_BINARY_DIR}/exports/${GMDS_LIB}_export.h + EXPORT_MACRO_NAME ${GMDS_LIB}_API) +#============================================================================== +# TARGET DEFINITION +#============================================================================== +include(GNUInstallDirs) + +#LIBRARY TO INSTALL +target_link_libraries(${GMDS_LIB} PUBLIC + ${LIB_GMDS_IG} + ${LIB_GMDS_IG_ALGO} + ${LIB_GMDS_IO} ANN) + +target_compile_features(${GMDS_LIB} PUBLIC cxx_std_14) + +# INCLUDE TO INSTALL +target_include_directories(${GMDS_LIB} PUBLIC + $ + $ + ) +set_target_properties(${GMDS_LIB} PROPERTIES PUBLIC_HEADER "${GMDS_INC}") + +install(TARGETS ${GMDS_LIB} + EXPORT GMDS_SUITE + PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${GMDS_LIB_PREFIX} + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) +#============================================================================== +add_executable(morphing src/main.cpp) +target_link_libraries(morphing PRIVATE ${GMDS_LIB}) +target_compile_features(morphing PUBLIC cxx_std_14) +install(TARGETS morphing) \ No newline at end of file diff --git a/morphMesh/README.md b/morphMesh/README.md new file mode 100644 index 000000000..c494bd762 --- /dev/null +++ b/morphMesh/README.md @@ -0,0 +1,3 @@ +#morphMesh module + +This module propose algorithm to morph mesh using deformation function. \ No newline at end of file diff --git a/morphMesh/inc/gmds/morphMesh/MorphMesh.h b/morphMesh/inc/gmds/morphMesh/MorphMesh.h new file mode 100644 index 000000000..b42a7185d --- /dev/null +++ b/morphMesh/inc/gmds/morphMesh/MorphMesh.h @@ -0,0 +1,60 @@ +/*----------------------------------------------------------------------------*/ +#ifndef GMDS_MORPHMESH_H +#define GMDS_MORPHMESH_H + +/*----------------------------------------------------------------------------*/ +#include "LIB_GMDS_MORPHMESH_export.h" +#include "gmds/ig/Mesh.h" +/*----------------------------------------------------------------------------*/ +namespace gmds { +/*----------------------------------------------------------------------------*/ +namespace morphmesh { +/*----------------------------------------------------------------------------*/ +class LIB_GMDS_MORPHMESH_API MorphMesh +{ + + public: + /*------------------------------------------------------------------------*/ + /** @brief Constructor. + */ + MorphMesh(Mesh* AMesh, const std::vector &APoints, double ARadius); + + /*------------------------------------------------------------------------*/ + /** @brief Destructor. + */ + ~MorphMesh(); + + /*------------------------------------------------------------------------*/ + /** @brief Constructor. + */ + void execute(); + + private: + + /*Surement a retirer plus tard*/ + void markLockedCells(); + + math::Point computeLocalOrigin(const math::Point& AP); + + double findHomothetyRatio(const math::Point& AP, const Node& ANode); + + private: + + /* Locked nodes */ + int m_locked; + /* List of "homothetic points" that will define the mesh deformation */ + std::vector m_targets; + /* The radius of an area that will be modified by a homothetic point*/ + double m_radius; + /* Nodes on the surface of the mesh */ + std::vector m_surfNodes; + + /* The mesh to deform */ + Mesh* m_mesh; +}; +/*----------------------------------------------------------------------------*/ +}// namespace morphmesh +/*----------------------------------------------------------------------------*/ +}// namespace gmds +/*----------------------------------------------------------------------------*/ +#endif // GMDS_MORPHMESH_H diff --git a/morphMesh/src/MorphMesh.cpp b/morphMesh/src/MorphMesh.cpp new file mode 100644 index 000000000..e805ce495 --- /dev/null +++ b/morphMesh/src/MorphMesh.cpp @@ -0,0 +1,134 @@ +/*----------------------------------------------------------------------------*/ +#include "gmds/morphMesh/MorphMesh.h" +/*----------------------------------------------------------------------------*/ +namespace gmds { +/*----------------------------------------------------------------------------*/ +namespace morphmesh { +/*----------------------------------------------------------------------------*/ +MorphMesh::MorphMesh(Mesh *AMesh, const std::vector &APoints, double ARadius): m_mesh(AMesh) +{ + m_radius = ARadius; + m_targets = APoints; + m_locked = m_mesh->newMark(); +} +/*----------------------------------------------------------------------------*/ +MorphMesh::~MorphMesh() {} +/*----------------------------------------------------------------------------*/ +void +MorphMesh::execute() +{ + markLockedCells(); + + for(auto n : m_mesh->nodes()){ + Node node = m_mesh->get(n); + if(node.get().size() == 4) + m_surfNodes.push_back(n); + } + + for(const auto& p : m_targets){ + std::cout<<"Traitement point "<get(n).point(); + double dist = n_point.distance(p); + if(dist < min_dist){ + min_dist = dist; + n0 = n; + } + } + Node node0 = m_mesh->get(n0); + math::Point point0 = node0.point(); + // on calcul le ration d'homothétie pour ce point + double H = findHomothetyRatio(p,node0); + + + for(auto n : m_mesh->nodes()){ + if(!m_mesh->isMarked(n,m_locked)) { + Node node = m_mesh->get(n); + math::Point verifPoint = node.point(); + verifPoint.setY(point0.Y()); // On fait ça pour mettre les deux points sur le même plan + + double dist = verifPoint.distance(point0); + if (dist <= m_radius) { + std::cout << "Dans le radius" << std::endl; + // Si le point est couvert par la zone à modifier alors on continue le traitement + double w = dist / m_radius; + std::cout << "Poids = " << w << std::endl; + + math::Point n_origin = computeLocalOrigin(node.point()); + math::Vector3d height = math::Vector3d(node.point() - n_origin); + + std::cout << "before H = " << H << std::endl; + double localH = H + ((1 - H) * w); // On réduit le ratio d'homothétie par le poids + std::cout << "after H = " << localH << std::endl; + + math::Point result = n_origin + height * localH; + + std::cout << "before " << node << std::endl; + node.setPoint(result); + std::cout << "after " << node << std::endl; + } + } + } + } +} +/*----------------------------------------------------------------------------*/ +void +MorphMesh::markLockedCells() +{ + /** ================================== + * !! WIP !! + ================================= + + 1) this method should also lock nodes given by user + 2) maybe we should lock cells too, edge, faces or region + */ + for(auto n : m_mesh->nodes()){ + Node node = m_mesh->get(n); + if(node.Y() <= 0 ){ + m_mesh->mark(n,m_locked); + } + } +} +/*----------------------------------------------------------------------------*/ +math::Point +MorphMesh::computeLocalOrigin(const math::Point& AP) +{ + /* + * 1) on crée un segment entre de AP et (X,Y,0) + * 2) On intersecte le segment avec les faces lock + * - si on a un point alors on retourne ce point + * - sinon on retourne (X,Y,0) + */ + + return {AP.X(),0,AP.Z()}; +} +/*----------------------------------------------------------------------------*/ +double MorphMesh::findHomothetyRatio(const math::Point &AP, const Node& ANode) +{ + + math::Point p_node = ANode.point(); + // Le vecteur de déformation initial : La translation du noeud du maillage vers le point + math::Vector3d deform(AP - ANode.point()); + + + // On cherche l'origine du noeud initial + math::Point origin = computeLocalOrigin(ANode.point()); + // On fait un vecteur origin->ni pour calculer l'homothétie des deux vecteurs. + // Faire le vecteur dans ce sens là permet de facilement déduire le ration d'homothétie + math::Vector3d height(ANode.point() - origin); + + + double dist1 = origin.distance(ANode.point()); + double dist2 = origin.distance(AP); + + + return dist2/dist1; +} +/*----------------------------------------------------------------------------*/ +} // namespace MorphMesh +/*----------------------------------------------------------------------------*/ +} // namespace gmds +/*----------------------------------------------------------------------------*/ \ No newline at end of file diff --git a/morphMesh/src/main.cpp b/morphMesh/src/main.cpp new file mode 100644 index 000000000..5de49b0c3 --- /dev/null +++ b/morphMesh/src/main.cpp @@ -0,0 +1,36 @@ +/*------------------------------------------------------------------------*/ +#include "gmds/ig/MeshDoctor.h" +#include "gmds/io/IGMeshIOService.h" +#include "gmds/io/VTKReader.h" +#include "gmds/io/VTKWriter.h" +#include +/*------------------------------------------------------------------------*/ +using namespace gmds; +/*------------------------------------------------------------------------*/ +int main(int argc, char* argv[]){ + + std::string param_file(argv[1]); + + gmds::Mesh m_mesh(gmds::MeshModel(gmds::DIM3 | gmds::R | gmds::N | gmds::N2R | gmds::R2N)); + gmds::IGMeshIOService ioService(&m_mesh); + gmds::VTKReader vtkReader(&ioService); + vtkReader.setCellOptions(gmds::N|gmds::R); + vtkReader.setDataOptions(gmds::N|gmds::R); + vtkReader.read(param_file); + + gmds::MeshDoctor doc(&m_mesh); + doc.buildN2R(m_mesh.getModel()); + doc.updateUpwardConnectivity(); + + gmds::math::Point point(5,0.5,0); + std::vector points; + points.push_back(point); + + gmds::morphmesh::MorphMesh morph(&m_mesh,points,2.5); + morph.execute(); + + gmds::VTKWriter w(&ioService); + w.setCellOptions(gmds::N|gmds::R); + w.setDataOptions(gmds::N|gmds::R); + w.write("morphing.vtk"); +} \ No newline at end of file From 0d4881765208a846c578162f5570af6b7b037d62 Mon Sep 17 00:00:00 2001 From: calderans Date: Mon, 12 Jun 2023 15:33:37 +0200 Subject: [PATCH 4/7] - big update on morphMesh component --- morphMesh/inc/gmds/morphMesh/MorphMesh.h | 8 +- morphMesh/src/MorphMesh.cpp | 259 ++++++++++++++++++++--- morphMesh/src/main.cpp | 9 +- 3 files changed, 240 insertions(+), 36 deletions(-) diff --git a/morphMesh/inc/gmds/morphMesh/MorphMesh.h b/morphMesh/inc/gmds/morphMesh/MorphMesh.h index b42a7185d..09d76532f 100644 --- a/morphMesh/inc/gmds/morphMesh/MorphMesh.h +++ b/morphMesh/inc/gmds/morphMesh/MorphMesh.h @@ -34,7 +34,8 @@ class LIB_GMDS_MORPHMESH_API MorphMesh /*Surement a retirer plus tard*/ void markLockedCells(); - math::Point computeLocalOrigin(const math::Point& AP); + bool computeLocalOrigin(const math::Point& AP, math::Point& AResult); + bool computeLocalOriginZ(const math::Point& AP, math::Point& AResult); double findHomothetyRatio(const math::Point& AP, const Node& ANode); @@ -42,12 +43,17 @@ class LIB_GMDS_MORPHMESH_API MorphMesh /* Locked nodes */ int m_locked; + Variable* locked_faces; /* List of "homothetic points" that will define the mesh deformation */ std::vector m_targets; /* The radius of an area that will be modified by a homothetic point*/ double m_radius; /* Nodes on the surface of the mesh */ std::vector m_surfNodes; + std::vector m_locked_faces; + //std::map is_locked_node; + //std::map nodes_origins; + /* The mesh to deform */ Mesh* m_mesh; diff --git a/morphMesh/src/MorphMesh.cpp b/morphMesh/src/MorphMesh.cpp index e805ce495..02b1e5ae6 100644 --- a/morphMesh/src/MorphMesh.cpp +++ b/morphMesh/src/MorphMesh.cpp @@ -1,5 +1,9 @@ /*----------------------------------------------------------------------------*/ #include "gmds/morphMesh/MorphMesh.h" +#include "gmds/math/Plane.h" +#include "gmds/math/Triangle.h" +#include + /*----------------------------------------------------------------------------*/ namespace gmds { /*----------------------------------------------------------------------------*/ @@ -10,13 +14,15 @@ MorphMesh::MorphMesh(Mesh *AMesh, const std::vector &APoints, doubl m_radius = ARadius; m_targets = APoints; m_locked = m_mesh->newMark(); + + locked_faces = m_mesh->newVariable("locked_face"); } /*----------------------------------------------------------------------------*/ MorphMesh::~MorphMesh() {} /*----------------------------------------------------------------------------*/ -void -MorphMesh::execute() +void MorphMesh::execute() { + std::cout<<"Marking cells"<nodes()){ @@ -25,8 +31,14 @@ MorphMesh::execute() m_surfNodes.push_back(n); } + std::vector durations1; + std::vector durations2; + + std::cout<<"Start deformation"<get(n0); math::Point point0 = node0.point(); // on calcul le ration d'homothétie pour ce point - double H = findHomothetyRatio(p,node0); + //double H = findHomothetyRatio(p,node0); + double Y0 = point0.Y(); + double X0 = point0.X(); + double Z = (Y0 - p.Y())/2; + math::Point tmp; + computeLocalOrigin(point0, tmp); + double omega_height0 = Y0 - tmp.Y(); + math::Vector3d t; + t.setXYZ(0,p.Y()-Y0,0); + math::Vector3d tz; + tz.setXYZ(0,0,Z); - for(auto n : m_mesh->nodes()){ - if(!m_mesh->isMarked(n,m_locked)) { + for(auto n : m_mesh->nodes()) { + if (!m_mesh->isMarked(n, m_locked)) { Node node = m_mesh->get(n); - math::Point verifPoint = node.point(); - verifPoint.setY(point0.Y()); // On fait ça pour mettre les deux points sur le même plan - double dist = verifPoint.distance(point0); + double dist = node.point().distance(point0); if (dist <= m_radius) { - std::cout << "Dans le radius" << std::endl; + time_t start,end; + time(&start); // Si le point est couvert par la zone à modifier alors on continue le traitement - double w = dist / m_radius; - std::cout << "Poids = " << w << std::endl; + //double w = 1 - (dist / m_radius); + double sigmoid = 1/(1 + exp(-5*( ( (1 - (dist / m_radius) ) *2)-1) ) ); + + math::Point n_origin; - math::Point n_origin = computeLocalOrigin(node.point()); - math::Vector3d height = math::Vector3d(node.point() - n_origin); + bool origin_locked = computeLocalOrigin(node.point(), n_origin); - std::cout << "before H = " << H << std::endl; - double localH = H + ((1 - H) * w); // On réduit le ratio d'homothétie par le poids - std::cout << "after H = " << localH << std::endl; + double omega_height = Y0 - n_origin.Y(); + double n_height = node.point().Y() - n_origin.Y(); - math::Point result = n_origin + height * localH; + math::Vector3d final_t = (n_height / omega_height) * (sigmoid * t); + math::Point result = node.point() + final_t; - std::cout << "before " << node << std::endl; node.setPoint(result); - std::cout << "after " << node << std::endl; + time(&end); + durations1.emplace_back(double(end - start)); + } + + double distX = fabs(node.X() - X0); + if (distX <= m_radius) { + time_t start,end; + time(&start); + // double w = 0.5+2*cbrt((1 - (distX / m_radius)) - 0.5); + double sigmoid = 1 / (1 + exp(-5 * (((1 - (distX / m_radius)) * 2) - 1))); + + math::Point n_origin; + bool origin_locked = computeLocalOriginZ(node.point(), n_origin); + + double omega_height = fabs(Y0) - fabs(n_origin.Z()); + double n_height = fabs(node.point().Z()) - fabs(n_origin.Z()); + + int sens = node.Z() < 0 ? -1 : 1; + math::Vector3d final_t = (n_height / omega_height) * (sigmoid * (sens * tz)); + + math::Point result = node.point() + final_t; + node.setPoint(result); + time(&end); + durations1.emplace_back(double(end - start)); } } } + /*for(auto n : m_mesh->nodes()) { + if (!m_mesh->isMarked(n, m_locked)) { + Node node = m_mesh->get(n); + double distX = fabs(node.X() - X0); + if (distX <= m_radius) { + // double w = 0.5+2*cbrt((1 - (distX / m_radius)) - 0.5); + double sigmoid = 1 / (1 + exp(-5 * (((1 - (distX / m_radius)) * 2) - 1))); + + math::Point n_origin; + bool origin_locked = computeLocalOriginZ(node.point(), n_origin); + + double omega_height = fabs(Y0) - fabs(n_origin.Z()); + double n_height = fabs(node.point().Z()) - fabs(n_origin.Z()); + + int sens = node.Z() < 0 ? -1 : 1; + math::Vector3d final_t = (n_height / omega_height) * (sigmoid * (sens * tz)); + + math::Point result = node.point() + final_t; + node.setPoint(result); + } + } + }*/ + } + + time(&endglob); + + + double avg = durations1[0]; + for(int i = 1; iunmarkAll(m_locked); + m_mesh->freeMark(m_locked); } /*----------------------------------------------------------------------------*/ -void -MorphMesh::markLockedCells() +void MorphMesh::markLockedCells() { - /** ================================== - * !! WIP !! - ================================= + Variable* locked_regions = m_mesh->newVariable("m_locked"); + std::set regions; + std::set nodes; + nodes.insert(27006); + m_mesh->mark(27006,m_locked); + for (int i = 0; i < 15; ++i) { + for(auto n : nodes){ + Node node = m_mesh->get(n); + for(const auto& r : node.get()){ + regions.emplace(r); + } + } + nodes.clear(); + for(const auto& r : regions){ + locked_regions->set(r.id(),1); + for(auto n : r.getIDs()){ + if(!m_mesh->isMarked(n,m_locked)){ + nodes.emplace(n); + m_mesh->mark(n,m_locked); + } + } + } + regions.clear(); + } - 1) this method should also lock nodes given by user - 2) maybe we should lock cells too, edge, faces or region - */ for(auto n : m_mesh->nodes()){ Node node = m_mesh->get(n); if(node.Y() <= 0 ){ m_mesh->mark(n,m_locked); } } + + //Test origine locale décallée + /*for(auto r : m_mesh->regions()){ + //if((4440<=r && r<=4639) || (5140<=r && r<=5339) || (5840<=r && r<=6039) || (6540<=r && r<=6739)){ + if((37000<=r && r<=37999) || (39000<=r && r<=39999) || (41000<=r && r<=41999) || (43000<=r && r<=43999)){ + locked_regions->set(r,1); + Region region = m_mesh->get(r); + for(auto n : region.getIDs()){ + m_mesh->mark(n,m_locked); + } + } + }*/ + + + for(auto f : m_mesh->faces()){ + Face face = m_mesh->get(f); + + std::vector f_regions = face.getIDs(); + if(f_regions.size() == 2){ + //Une face entre un hexa lock et un non lock + if((locked_regions->value(f_regions[0]) == 1) && (locked_regions->value(f_regions[1]) == 0) || + (locked_regions->value(f_regions[0]) == 0) && (locked_regions->value(f_regions[1]) == 1)){ + locked_faces->set(f,1); + m_locked_faces.push_back(f); + } + } + } } /*----------------------------------------------------------------------------*/ -math::Point -MorphMesh::computeLocalOrigin(const math::Point& AP) +bool MorphMesh::computeLocalOrigin(const math::Point& AP, math::Point& AResult) { /* * 1) on crée un segment entre de AP et (X,Y,0) @@ -102,8 +226,76 @@ MorphMesh::computeLocalOrigin(const math::Point& AP) * - si on a un point alors on retourne ce point * - sinon on retourne (X,Y,0) */ + std::vector points; + points.emplace_back(AP.X(),0,AP.Z()); + points.emplace_back(AP.X()+0.001,0,AP.Z()+0.001); + points.emplace_back(AP.X()-0.001,0,AP.Z()-0.001); + + for (auto f : m_locked_faces) { + for(const auto& p : points) { + //if (locked_faces->value(f) == 1) { + math::Segment depth(AP, p); + Face face = m_mesh->get(f); + math::Plane plane(face.center(), face.normal()); + // on cherche si la profondeur a une intersection avec un plan d'une face locked + if (depth.intersect(plane, AResult)) { + // si oui on cherche si le point d'intersection est dans la face + // On fait ça en créant deux triangles qui couvrent entièrement la face + // et qui permettent d'utiliser la méthode isIn() + std::vector nodes = face.get(); + math::Triangle tri1(nodes[0].point(), nodes[1].point(), nodes[2].point()); + math::Triangle tri2(nodes[2].point(), nodes[3].point(), nodes[0].point()); + if (tri1.intersect(depth) || tri2.intersect(depth)) { + return true; + } + } + //} + } + } + AResult = {AP.X(),0,AP.Z()}; + + return false; +} +/*----------------------------------------------------------------------------*/ +bool MorphMesh::computeLocalOriginZ(const math::Point& AP, math::Point& AResult) +{ + /* + * 1) on crée un segment entre de AP et (X,Y,0) + * 2) On intersecte le segment avec les faces lock + * - si on a un point alors on retourne ce point + * - sinon on retourne (X,Y,0) + */ + std::vector points; + points.emplace_back(AP.X(),AP.Y(),0); + points.emplace_back(AP.X()+0.001,AP.Y()+0.001,0); + points.emplace_back(AP.X()-0.001,AP.Y()-0.001,0); + + for (auto f : m_locked_faces) { + for(const auto& p : points) { + //if (locked_faces->value(f) == 1) { + math::Segment depth(AP,math::Point(AP.X(),AP.Y(),0)); + Face face = m_mesh->get(f); + math::Plane plane(face.center(), face.normal()); + // on cherche si la profondeur a une intersection avec un plan d'une face locked + if (depth.intersect(plane, AResult)) { + // std::cout<<"test"< nodes = face.get(); + math::Triangle tri1(nodes[0].point(), nodes[1].point(), nodes[2].point()); + math::Triangle tri2(nodes[2].point(), nodes[3].point(), nodes[0].point()); + if (tri1.intersect(depth) || tri2.intersect(depth)) { + return true; + } + } + //} + } + } + + AResult = {AP.X(),AP.Y(),0}; - return {AP.X(),0,AP.Z()}; + return false; } /*----------------------------------------------------------------------------*/ double MorphMesh::findHomothetyRatio(const math::Point &AP, const Node& ANode) @@ -115,7 +307,8 @@ double MorphMesh::findHomothetyRatio(const math::Point &AP, const Node& ANode) // On cherche l'origine du noeud initial - math::Point origin = computeLocalOrigin(ANode.point()); + math::Point origin; + bool origin_locked = computeLocalOrigin(ANode.point(), origin); // On fait un vecteur origin->ni pour calculer l'homothétie des deux vecteurs. // Faire le vecteur dans ce sens là permet de facilement déduire le ration d'homothétie math::Vector3d height(ANode.point() - origin); @@ -124,6 +317,8 @@ double MorphMesh::findHomothetyRatio(const math::Point &AP, const Node& ANode) double dist1 = origin.distance(ANode.point()); double dist2 = origin.distance(AP); + std::cout<<"dist origin -> node "< point "< points; points.push_back(point); @@ -32,5 +33,7 @@ int main(int argc, char* argv[]){ gmds::VTKWriter w(&ioService); w.setCellOptions(gmds::N|gmds::R); w.setDataOptions(gmds::N|gmds::R); - w.write("morphing.vtk"); + w.write("morphingtest.vtk"); + + } \ No newline at end of file From d0a507fdf699b70ebc1fc985f11f17a255230764 Mon Sep 17 00:00:00 2001 From: nicolas le goff Date: Tue, 29 Aug 2023 23:31:40 +0200 Subject: [PATCH 5/7] copied FastLocalize from the aero component into morphmesh --- morphMesh/CMakeLists.txt | 7 ++- morphMesh/inc/gmds/morphMesh/FastLocalize.h | 52 +++++++++++++++ morphMesh/src/FastLocalize.cpp | 70 +++++++++++++++++++++ morphMesh/tst/CMakeLists.txt | 13 ++++ morphMesh/tst/FastLocalizeTestSuite.h | 67 ++++++++++++++++++++ morphMesh/tst/main_test.cpp | 13 ++++ 6 files changed, 221 insertions(+), 1 deletion(-) create mode 100644 morphMesh/inc/gmds/morphMesh/FastLocalize.h create mode 100644 morphMesh/src/FastLocalize.cpp create mode 100644 morphMesh/tst/CMakeLists.txt create mode 100644 morphMesh/tst/FastLocalizeTestSuite.h create mode 100644 morphMesh/tst/main_test.cpp diff --git a/morphMesh/CMakeLists.txt b/morphMesh/CMakeLists.txt index 3c698f147..3ca77353e 100644 --- a/morphMesh/CMakeLists.txt +++ b/morphMesh/CMakeLists.txt @@ -7,11 +7,12 @@ set(GMDS_LIB_PREFIX gmds/morphMesh) set(GMDS_INC ${CMAKE_BINARY_DIR}/exports/${GMDS_LIB}_export.h + inc/gmds/morphMesh/FastLocalize.h inc/gmds/morphMesh/MorphMesh.h ) set(GMDS_SRC + src/FastLocalize.cpp src/MorphMesh.cpp - src/main.cpp ) #============================================================================== add_library(${GMDS_LIB} ${GMDS_INC} ${GMDS_SRC}) @@ -47,6 +48,10 @@ install(TARGETS ${GMDS_LIB} ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) #============================================================================== +if(WITH_TEST) + add_subdirectory(tst) +endif(WITH_TEST) +#============================================================================== add_executable(morphing src/main.cpp) target_link_libraries(morphing PRIVATE ${GMDS_LIB}) target_compile_features(morphing PUBLIC cxx_std_14) diff --git a/morphMesh/inc/gmds/morphMesh/FastLocalize.h b/morphMesh/inc/gmds/morphMesh/FastLocalize.h new file mode 100644 index 000000000..d44c808d5 --- /dev/null +++ b/morphMesh/inc/gmds/morphMesh/FastLocalize.h @@ -0,0 +1,52 @@ +/*----------------------------------------------------------------------------*/ +#ifndef GMDS_MORPHMESH_FASTLOCALIZE_H +#define GMDS_MORPHMESH_FASTLOCALIZE_H +/*----------------------------------------------------------------------------*/ +#include "LIB_GMDS_MORPHMESH_export.h" +# include +/*----------------------------------------------------------------------------*/ +#include +#include +/*----------------------------------------------------------------------------*/ +namespace gmds { +/*----------------------------------------------------------------------------*/ +class LIB_GMDS_MORPHMESH_API FastLocalize { + public: + + /*-------------------------------------------------------------------*/ + /** @brief Constructor. + * @param AMesh the mesh where we work on + */ + FastLocalize(Mesh *AMesh); + /*-------------------------------------------------------------------*/ + /** @brief Destructor. + */ + virtual ~FastLocalize(); + + /*-------------------------------------------------------------------*/ + /** @brief Given an input point \p APoint return the cell data that + * indicate what is the closest cell of \APoint (it can be a node or + * a face) + * @param APoint point to localize + * @return the cell data indicating where \p APoint is in this->m_mesh + */ + Cell::Data find(const math::Point& APoint); + /*-------------------------------------------------------------------*/ + + private: + void buildANNTree(); + private: + /** mesh we work on */ + Mesh *m_mesh; + std::map m_ann2gmds_id; + ANNkd_tree* m_kdTree; + ANNidxArray m_nnIdx; // near neighbor indices + ANNdistArray m_dists; // near neighbor distances + ANNpoint m_queryPt; // query point + ANNpointArray m_dataPts; // data points +}; +/*----------------------------------------------------------------------------*/ +} // namespace gmds +/*----------------------------------------------------------------------------*/ +#endif // GMDS_MORPHMESH_FASTLOCALIZE_H +/*----------------------------------------------------------------------------*/ diff --git a/morphMesh/src/FastLocalize.cpp b/morphMesh/src/FastLocalize.cpp new file mode 100644 index 000000000..12728ec8c --- /dev/null +++ b/morphMesh/src/FastLocalize.cpp @@ -0,0 +1,70 @@ +/*------------------------------------------------------------------------*/ +#include "gmds/morphMesh/FastLocalize.h" +//#include +/*------------------------------------------------------------------------*/ +using namespace gmds; +/*------------------------------------------------------------------------*/ +FastLocalize::FastLocalize(Mesh *AMesh): m_mesh(AMesh) { + int nb_pnts = m_mesh->getNbNodes(); + int k = 20; // max number of nearest neighbors + int dim = 3; // dimension + int maxPts = nb_pnts; // maximum number of data points + + int nPts; // actual number of data points + + m_queryPt = annAllocPt(3); + m_dataPts = annAllocPts(maxPts, dim); // allocate data points + m_nnIdx = new ANNidx[k]; // allocate near neigh indices + m_dists = new ANNdist[k]; // allocate near neighbor dists + + //======================================================== + // (1) Fill in the ANN structure for storing points + // + // Important: Points in APnts and dataPnts are stored with + // same index. + //======================================================== + nPts = 0; + for(auto n_id: m_mesh->nodes()){ + math::Point p = m_mesh->get(n_id).point(); + m_dataPts[nPts][0] = p.X(); + m_dataPts[nPts][1] = p.Y(); + m_dataPts[nPts][2] = p.Z(); + m_ann2gmds_id[nPts]=n_id; + nPts++; + }; + + //======================================================== + // (2) Build the search structure + //======================================================== + m_kdTree = new ANNkd_tree(m_dataPts, // the data points + nPts, // number of points + dim); // dimension of space +} +/*------------------------------------------------------------------------*/ +FastLocalize::~FastLocalize(){ + delete [] m_nnIdx; // clean things up + delete [] m_dists; + delete m_kdTree; + annDeallocPt(m_queryPt); + annDeallocPts(m_dataPts); + annClose(); // done with ANN + +} +/*------------------------------------------------------------------------*/ +Cell::Data +FastLocalize::find(const math::Point &APoint) +{ + int k = 5; // max number of nearest neighbors + m_queryPt[0]=APoint.X(); + m_queryPt[1]=APoint.Y(); + m_queryPt[2]=APoint.Z(); + + m_kdTree->annkSearch( // search + m_queryPt,// query point + k, + m_nnIdx, + m_dists, + 0.01); + return Cell::Data(0,m_ann2gmds_id[m_nnIdx[0]]); +} +/*------------------------------------------------------------------------*/ \ No newline at end of file diff --git a/morphMesh/tst/CMakeLists.txt b/morphMesh/tst/CMakeLists.txt new file mode 100644 index 000000000..ecab75d7e --- /dev/null +++ b/morphMesh/tst/CMakeLists.txt @@ -0,0 +1,13 @@ +add_executable(GMDS_MORPHMESH_TEST + FastLocalizeTestSuite.h + main_test.cpp + ) +#============================================================================== +target_link_libraries(GMDS_MORPHMESH_TEST PUBLIC + ${GMDS_LIB} + ${LIB_GMDS_IO} + GTest::gtest) +#============================================================================== +gtest_discover_tests(GMDS_MORPHMESH_TEST + WORKING_DIRECTORY ${CMAKE_BINARY_DIR}) +#============================================================================== diff --git a/morphMesh/tst/FastLocalizeTestSuite.h b/morphMesh/tst/FastLocalizeTestSuite.h new file mode 100644 index 000000000..7dbd3129a --- /dev/null +++ b/morphMesh/tst/FastLocalizeTestSuite.h @@ -0,0 +1,67 @@ +#include +#include +#include +#include +#include +#include +#include +#include +/*----------------------------------------------------------------------------*/ +/* TESTS UNITAIRES */ +/*----------------------------------------------------------------------------*/ + +TEST(MorphMesh_FastLocalizeTestSuite, test_FastLocalize_2D) +{ + // WE READ + gmds::Mesh m(gmds::MeshModel(gmds::DIM3 | gmds::F | gmds::N | gmds::E | gmds::N2E | gmds::N2F | gmds::F2N | gmds::E2N | gmds::F2E | gmds::E2F)); + + std::string dir(TEST_SAMPLES_DIR); + std::string vtk_file = dir + "/Aero/2D/C1_2D_0.1.vtk"; + + gmds::IGMeshIOService ioService(&m); + gmds::VTKReader vtkReader(&ioService); + vtkReader.setCellOptions(gmds::N | gmds::F); + vtkReader.read(vtk_file); + + gmds::MeshDoctor doc(&m); + doc.buildEdgesAndX2E(); + doc.updateUpwardConnectivity(); + + gmds::FastLocalize fl(&m); + gmds::Cell::Data data = fl.find(gmds::math::Point({-0.0338,0.74,0})); + ASSERT_EQ(data.dim, 0); + ASSERT_EQ(data.id, 734); + + data = fl.find(gmds::math::Point({-0.2,0.46,0})); + ASSERT_EQ(data.dim, 0); + ASSERT_EQ(data.id, 11); +} + + +TEST(MorphMesh_FastLocalizeTestSuite, test_FastLocalize_3D) +{ + // WE READ + gmds::Mesh m(gmds::MeshModel(gmds::DIM3 | gmds::F | gmds::N | gmds::E | gmds::N2E | gmds::N2F | gmds::F2N | gmds::E2N | gmds::F2E | gmds::E2F)); + + std::string dir(TEST_SAMPLES_DIR); + std::string vtk_file = dir + "/Aero/3D/C1_3D_0.1.vtk"; + + gmds::IGMeshIOService ioService(&m); + gmds::VTKReader vtkReader(&ioService); + vtkReader.setCellOptions(gmds::N | gmds::F); + vtkReader.read(vtk_file); + + gmds::MeshDoctor doc(&m); + doc.buildEdgesAndX2E(); + doc.updateUpwardConnectivity(); + + gmds::FastLocalize fl(&m); + + gmds::Cell::Data data = fl.find(gmds::math::Point({-1.0,0.0,0})); + ASSERT_EQ(data.dim, 0); + ASSERT_EQ(data.id, 14322); + + data = fl.find(gmds::math::Point({-1.33,0.678,-1.89})); + ASSERT_EQ(data.dim, 0); + ASSERT_EQ(data.id, 4048); +} \ No newline at end of file diff --git a/morphMesh/tst/main_test.cpp b/morphMesh/tst/main_test.cpp new file mode 100644 index 000000000..174d6c0cd --- /dev/null +++ b/morphMesh/tst/main_test.cpp @@ -0,0 +1,13 @@ +/*----------------------------------------------------------------------------*/ +#include +/*----------------------------------------------------------------------------*/ +// Files containing the different test suites to launch + +#include "FastLocalizeTestSuite.h" +/*----------------------------------------------------------------------------*/ +int main(int argc, char ** argv) { + ::testing::InitGoogleTest(&argc, argv); + return RUN_ALL_TESTS(); +} +/*----------------------------------------------------------------------------*/ + From 2bde8cbe91ac2e8f9cd1a8f17c5885f5248fcace Mon Sep 17 00:00:00 2001 From: nicolas le goff Date: Wed, 30 Aug 2023 08:32:26 +0200 Subject: [PATCH 6/7] check that fastlocalize works on a point cloud --- morphMesh/tst/FastLocalizeTestSuite.h | 36 +++++++++++++++++++++------ 1 file changed, 28 insertions(+), 8 deletions(-) diff --git a/morphMesh/tst/FastLocalizeTestSuite.h b/morphMesh/tst/FastLocalizeTestSuite.h index 7dbd3129a..8f294e9f9 100644 --- a/morphMesh/tst/FastLocalizeTestSuite.h +++ b/morphMesh/tst/FastLocalizeTestSuite.h @@ -13,7 +13,7 @@ TEST(MorphMesh_FastLocalizeTestSuite, test_FastLocalize_2D) { // WE READ - gmds::Mesh m(gmds::MeshModel(gmds::DIM3 | gmds::F | gmds::N | gmds::E | gmds::N2E | gmds::N2F | gmds::F2N | gmds::E2N | gmds::F2E | gmds::E2F)); + gmds::Mesh m(gmds::MeshModel(gmds::DIM3 | gmds::N | gmds::F | gmds::F2N)); std::string dir(TEST_SAMPLES_DIR); std::string vtk_file = dir + "/Aero/2D/C1_2D_0.1.vtk"; @@ -23,11 +23,13 @@ TEST(MorphMesh_FastLocalizeTestSuite, test_FastLocalize_2D) vtkReader.setCellOptions(gmds::N | gmds::F); vtkReader.read(vtk_file); - gmds::MeshDoctor doc(&m); - doc.buildEdgesAndX2E(); - doc.updateUpwardConnectivity(); + // fastlocalize works on a point cloud, no need for the faces + for(auto i: m.faces()) { + m.deleteFace(i); + } gmds::FastLocalize fl(&m); + gmds::Cell::Data data = fl.find(gmds::math::Point({-0.0338,0.74,0})); ASSERT_EQ(data.dim, 0); ASSERT_EQ(data.id, 734); @@ -41,7 +43,7 @@ TEST(MorphMesh_FastLocalizeTestSuite, test_FastLocalize_2D) TEST(MorphMesh_FastLocalizeTestSuite, test_FastLocalize_3D) { // WE READ - gmds::Mesh m(gmds::MeshModel(gmds::DIM3 | gmds::F | gmds::N | gmds::E | gmds::N2E | gmds::N2F | gmds::F2N | gmds::E2N | gmds::F2E | gmds::E2F)); + gmds::Mesh m(gmds::MeshModel(gmds::DIM3 | gmds::N | gmds::F | gmds::F2N)); std::string dir(TEST_SAMPLES_DIR); std::string vtk_file = dir + "/Aero/3D/C1_3D_0.1.vtk"; @@ -51,9 +53,10 @@ TEST(MorphMesh_FastLocalizeTestSuite, test_FastLocalize_3D) vtkReader.setCellOptions(gmds::N | gmds::F); vtkReader.read(vtk_file); - gmds::MeshDoctor doc(&m); - doc.buildEdgesAndX2E(); - doc.updateUpwardConnectivity(); + // fastlocalize works on a point cloud, no need for the faces + for(auto i: m.faces()) { + m.deleteFace(i); + } gmds::FastLocalize fl(&m); @@ -64,4 +67,21 @@ TEST(MorphMesh_FastLocalizeTestSuite, test_FastLocalize_3D) data = fl.find(gmds::math::Point({-1.33,0.678,-1.89})); ASSERT_EQ(data.dim, 0); ASSERT_EQ(data.id, 4048); +} + + +TEST(MorphMesh_FastLocalizeTestSuite, test_FastLocalize_nodes) +{ + // WE READ + gmds::Mesh m(gmds::MeshModel(gmds::DIM3 | gmds::N | gmds::F | gmds::F2N)); + + for(int i=0; i<50; i++) { + m.newNode(i,1,1); + } + + gmds::FastLocalize fl(&m); + + gmds::Cell::Data data = fl.find(gmds::math::Point({1.9,0.0,0})); + ASSERT_EQ(data.dim, 0); + ASSERT_EQ(data.id, 2); } \ No newline at end of file From 492890b580d221c89520b8826cb79dd47b591d88 Mon Sep 17 00:00:00 2001 From: calderans Date: Fri, 1 Sep 2023 11:31:20 +0200 Subject: [PATCH 7/7] - added elliptic morphing --- math/CMakeLists.txt | 6 +- morphMesh/CMakeLists.txt | 2 + morphMesh/inc/gmds/morphMesh/EllipticMorph.h | 55 ++++++ morphMesh/inc/gmds/morphMesh/FastLocalize.h | 6 +- morphMesh/src/EllipticMorph.cpp | 178 +++++++++++++++++++ morphMesh/src/FastLocalize.cpp | 16 +- morphMesh/src/main.cpp | 24 ++- 7 files changed, 269 insertions(+), 18 deletions(-) create mode 100644 morphMesh/inc/gmds/morphMesh/EllipticMorph.h create mode 100644 morphMesh/src/EllipticMorph.cpp diff --git a/math/CMakeLists.txt b/math/CMakeLists.txt index 7614a6751..de84c8500 100644 --- a/math/CMakeLists.txt +++ b/math/CMakeLists.txt @@ -37,6 +37,7 @@ set(GMDS_INC inc/gmds/math/Vector.h inc/gmds/math/VectorDyn.h inc/gmds/math/StreamComputation.h + inc/gmds/math/BezierSurface.h ) set(GMDS_SRC src/AxisAngleRotation.cpp @@ -67,10 +68,11 @@ set(GMDS_SRC src/VectorDyn.cpp src/Orientation.cpp src/Vector.cpp - src/TransfiniteInterpolation.cpp inc/gmds/math/BezierSurface.h src/BezierSurface.cpp) + src/TransfiniteInterpolation.cpp + src/BezierSurface.cpp) #============================================================================== find_package(Eigen3 3.3 REQUIRED) -#============================================================================== +#==============================================================================k add_library(${GMDS_LIB} ${GMDS_INC} ${GMDS_SRC}) #============================================================================== include(GenerateExportHeader) diff --git a/morphMesh/CMakeLists.txt b/morphMesh/CMakeLists.txt index 3ca77353e..cde2637d9 100644 --- a/morphMesh/CMakeLists.txt +++ b/morphMesh/CMakeLists.txt @@ -9,10 +9,12 @@ set(GMDS_INC ${CMAKE_BINARY_DIR}/exports/${GMDS_LIB}_export.h inc/gmds/morphMesh/FastLocalize.h inc/gmds/morphMesh/MorphMesh.h + inc/gmds/morphMesh/EllipticMorph.h ) set(GMDS_SRC src/FastLocalize.cpp src/MorphMesh.cpp + src/EllipticMorph.cpp ) #============================================================================== add_library(${GMDS_LIB} ${GMDS_INC} ${GMDS_SRC}) diff --git a/morphMesh/inc/gmds/morphMesh/EllipticMorph.h b/morphMesh/inc/gmds/morphMesh/EllipticMorph.h new file mode 100644 index 000000000..f56251d7c --- /dev/null +++ b/morphMesh/inc/gmds/morphMesh/EllipticMorph.h @@ -0,0 +1,55 @@ + +#ifndef GMDS_ELLIPTICMORPH_H +#define GMDS_ELLIPTICMORPH_H + +/*----------------------------------------------------------------------------*/ +#include "LIB_GMDS_MORPHMESH_export.h" +#include "gmds/ig/Mesh.h" +/*----------------------------------------------------------------------------*/ +namespace gmds { +/*----------------------------------------------------------------------------*/ +namespace morphmesh { +/*----------------------------------------------------------------------------*/ +class LIB_GMDS_MORPHMESH_API EllipticMorph +{ + + public: + /*------------------------------------------------------------------------*/ + /** @brief Constructor. + */ + EllipticMorph(Mesh* AMesh); + + /*------------------------------------------------------------------------*/ + /** @brief Destructor. + */ + ~EllipticMorph(); + + /*------------------------------------------------------------------------*/ + /** @brief Constructor. + */ + void execute(const std::vector>& AEllipses); + + private: + + /*Surement a retirer plus tard*/ + void markLockedCells(); + + + private: + + /* Locked nodes */ + int m_locked; + /* Nodes on the surface of the mesh */ + std::vector m_surfNodes; + std::vector m_locked_faces; + + + /* The mesh to deform */ + Mesh* m_mesh; +}; +/*----------------------------------------------------------------------------*/ +}// namespace morphmesh +/*----------------------------------------------------------------------------*/ +}// namespace gmds +/*----------------------------------------------------------------------------*/ +#endif // GMDS_ELLIPTICMORPH_H diff --git a/morphMesh/inc/gmds/morphMesh/FastLocalize.h b/morphMesh/inc/gmds/morphMesh/FastLocalize.h index d44c808d5..849218c42 100644 --- a/morphMesh/inc/gmds/morphMesh/FastLocalize.h +++ b/morphMesh/inc/gmds/morphMesh/FastLocalize.h @@ -17,7 +17,7 @@ class LIB_GMDS_MORPHMESH_API FastLocalize { /** @brief Constructor. * @param AMesh the mesh where we work on */ - FastLocalize(Mesh *AMesh); + FastLocalize(const std::vector &ANodes); /*-------------------------------------------------------------------*/ /** @brief Destructor. */ @@ -30,14 +30,14 @@ class LIB_GMDS_MORPHMESH_API FastLocalize { * @param APoint point to localize * @return the cell data indicating where \p APoint is in this->m_mesh */ - Cell::Data find(const math::Point& APoint); + TCellID find(const math::Point& APoint); /*-------------------------------------------------------------------*/ private: void buildANNTree(); private: /** mesh we work on */ - Mesh *m_mesh; + std::vector m_nodes; std::map m_ann2gmds_id; ANNkd_tree* m_kdTree; ANNidxArray m_nnIdx; // near neighbor indices diff --git a/morphMesh/src/EllipticMorph.cpp b/morphMesh/src/EllipticMorph.cpp new file mode 100644 index 000000000..f54472cd5 --- /dev/null +++ b/morphMesh/src/EllipticMorph.cpp @@ -0,0 +1,178 @@ +// +// Created by calderans on 30/08/23. +// + +#include "gmds/morphMesh/EllipticMorph.h" +#include "gmds/io/IGMeshIOService.h" +#include "gmds/morphMesh/FastLocalize.h" +/*----------------------------------------------------------------------------*/ + +using namespace gmds; +using namespace morphmesh; +/*----------------------------------------------------------------------------*/ +EllipticMorph::EllipticMorph(Mesh *AMesh): m_mesh(AMesh) { + + m_locked = m_mesh->newMark(); + markLockedCells(); + +} +/*----------------------------------------------------------------------------*/ +EllipticMorph::~EllipticMorph() {} +/*----------------------------------------------------------------------------*/ +void EllipticMorph::execute(const std::vector>& AEllipses) +{ + + std::vector nodes_int; + std::vector nodes_ext; + + + for(auto n : m_mesh->nodes()){ + Node node = m_mesh->get(n); + if(m_mesh->isMarked(n, m_locked)){ + nodes_int.push_back(node); + } + if((node.getIDs().size() == 4 || node.getIDs().size() == 2) && !m_mesh->isMarked(n,m_locked)){ + nodes_ext.push_back(node); + } + } + std::map vecs; + std::set modified_nodes; + + for(int i = 0; i < AEllipses.size()-1; i++) { + + double Xmin = AEllipses[i][0]; + double Xmax = AEllipses[i+1][0]; + + std::vector execution; + + for(auto n : m_mesh->nodes()){ + + Node node = m_mesh->get(n); + + if (node.X() >= Xmin && node.X() < Xmax && !m_mesh->isMarked(n, m_locked)) { + execution.push_back(node); + } + } + + double coefa_min = AEllipses[i][1]; + double coefb_min = AEllipses[i][2]; + double coefc_min = AEllipses[i][3]; + + double coefa_max = AEllipses[i+1][1]; + double coefb_max = AEllipses[i+1][2]; + double coefc_max = AEllipses[i+1][3]; + + FastLocalize fl_int(nodes_int); + FastLocalize fl_ext(nodes_ext); + + + for (auto n : execution) { + modified_nodes.insert(n.id()); + + math::Point p = n.point(); + + Node n_ext = m_mesh->get(fl_ext.find(p)); + math::Point p_int = m_mesh->get(fl_int.find(p)).point(); + math::Point p_ext = n_ext.point(); + math::Point p_axe(p.X(), 0, 0); + + double distXmin = p.X()-Xmin; + double distXmax = Xmax-Xmin; + double w = distXmin/distXmax; + + // math::Point p_origin = p.distance(p_int) < p.distance(p_axe) ? ; + + double coefy_min = p.Y() >= 0 ? coefa_min : coefc_min; + double coefy_max = p.Y() >= 0 ? coefa_max : coefc_max; + + double coefy_final = coefy_min * (1/(1 + exp(-5*( ( (1 - w ) *2)-1) ) )) + coefy_max * (1/(1 + exp(-5*( ( w *2)-1) ) )); + double coefz_final = coefb_min * (1/(1 + exp(-5*( ( (1 - w ) *2)-1) ) )) + coefb_max * (1/(1 + exp(-5*( ( w *2)-1) ) )); + + //double coefy_final = coefy_min * (1 - w) + coefy_max * (w); + //double coefz_final = coefb_min * (1 - w) + coefb_max * (w); + + + double theta = atan(n.Y() / n.Z()); + if (theta < 0) { + theta *= -1; + } + + double dist1 = p.distance(p_int); + double dist2 = p_ext.distance(p_int); + double omega = dist1 / dist2; + + // double z_new = omega*(coefa * cos(theta)); + // double y_new = omega*(coefb * sin(theta)); + double sinus = sin(theta); + double cosinus = cos(theta); + + if (n.Y() < 0) { + sinus *= -1; + } + if (n.Z() < 0) { + cosinus *= -1; + } + + math::Vector3d vec; + vec.setXYZ(0, (coefy_final * sinus) - p_ext.Y(), (coefz_final * cosinus) - p_ext.Z()); + + vec = omega * vec; + + vecs[n.id()] = vec; + } + } + + for(auto n : modified_nodes){ + + math::Vector3d vec(vecs[n]); + Node node = m_mesh->get(n); + double y = node.Y(); + double z = node.Z(); + + node.setY(y+vec.Y()); + node.setZ(z+vec.Z()); + } + +} +/*----------------------------------------------------------------------------*/ +void EllipticMorph::markLockedCells() +{ + Variable* locked_regions = m_mesh->newVariable("m_locked"); + + + std::set regions; + std::set nodes; + nodes.insert(27006); + m_mesh->mark(27006,m_locked); + for (int i = 0; i < 15; ++i) { + for(auto n : nodes){ + Node node = m_mesh->get(n); + + for(const auto& r : node.get()){ + regions.emplace(r); + } + } + nodes.clear(); + for(const auto& r : regions){ + locked_regions->set(r.id(),1); + for(auto n : r.getIDs()){ + if(!m_mesh->isMarked(n,m_locked)){ + nodes.emplace(n); + m_mesh->mark(n,m_locked); + + } + } + } + regions.clear(); + } + + + + for(auto n : m_mesh->nodes()){ + Node node = m_mesh->get(n); + math::Point p(node.point().X(),0,0); + if(node.point().distance(p) < 0.2){ + m_mesh->mark(n,m_locked); + } + } +} \ No newline at end of file diff --git a/morphMesh/src/FastLocalize.cpp b/morphMesh/src/FastLocalize.cpp index 12728ec8c..2a90de9b9 100644 --- a/morphMesh/src/FastLocalize.cpp +++ b/morphMesh/src/FastLocalize.cpp @@ -4,8 +4,8 @@ /*------------------------------------------------------------------------*/ using namespace gmds; /*------------------------------------------------------------------------*/ -FastLocalize::FastLocalize(Mesh *AMesh): m_mesh(AMesh) { - int nb_pnts = m_mesh->getNbNodes(); +FastLocalize::FastLocalize(const std::vector &ANodes): m_nodes(ANodes){ + int nb_pnts = m_nodes.size(); int k = 20; // max number of nearest neighbors int dim = 3; // dimension int maxPts = nb_pnts; // maximum number of data points @@ -24,14 +24,14 @@ FastLocalize::FastLocalize(Mesh *AMesh): m_mesh(AMesh) { // same index. //======================================================== nPts = 0; - for(auto n_id: m_mesh->nodes()){ - math::Point p = m_mesh->get(n_id).point(); + for(auto const &n: m_nodes){ + math::Point p = n.point(); m_dataPts[nPts][0] = p.X(); m_dataPts[nPts][1] = p.Y(); m_dataPts[nPts][2] = p.Z(); - m_ann2gmds_id[nPts]=n_id; + m_ann2gmds_id[nPts]=n.id(); nPts++; - }; + } //======================================================== // (2) Build the search structure @@ -51,7 +51,7 @@ FastLocalize::~FastLocalize(){ } /*------------------------------------------------------------------------*/ -Cell::Data +TCellID FastLocalize::find(const math::Point &APoint) { int k = 5; // max number of nearest neighbors @@ -65,6 +65,6 @@ FastLocalize::find(const math::Point &APoint) m_nnIdx, m_dists, 0.01); - return Cell::Data(0,m_ann2gmds_id[m_nnIdx[0]]); + return m_ann2gmds_id[m_nnIdx[0]]; } /*------------------------------------------------------------------------*/ \ No newline at end of file diff --git a/morphMesh/src/main.cpp b/morphMesh/src/main.cpp index f91395a84..845201539 100644 --- a/morphMesh/src/main.cpp +++ b/morphMesh/src/main.cpp @@ -4,6 +4,7 @@ #include "gmds/io/VTKReader.h" #include "gmds/io/VTKWriter.h" #include +#include /*------------------------------------------------------------------------*/ using namespace gmds; /*------------------------------------------------------------------------*/ @@ -23,12 +24,25 @@ int main(int argc, char* argv[]){ doc.buildFacesAndR2F(); doc.updateUpwardConnectivity(); - gmds::math::Point point(6,0.5,0); - std::vector points; - points.push_back(point); + //gmds::morphmesh::MorphMesh morph(&m_mesh,points,2.5); + //morph.execute(); - gmds::morphmesh::MorphMesh morph(&m_mesh,points,2.5); - morph.execute(); + morphmesh::EllipticMorph emorph(&m_mesh); + + std::vector> ellipses; + //ellipses.push_back({2,1,1,1}); + ellipses.push_back({2,1,1,1}); + ellipses.push_back({3,1.2,1.5,1}); + ellipses.push_back({4,1.2,1.5,1}); + ellipses.push_back({5,1.4,1.75,1.5}); + ellipses.push_back({7,1.4,1.75,1.5}); + ellipses.push_back({8,1,1,1}); + //ellipses.push_back({5.6,1.2,1.5,1}); + //ellipses.push_back({8.5,1.2,1.5,1}); + + emorph.execute(ellipses); + + std::cout<<"test"<