From 8c6feadf41c96a1d82d87ed82882002d36ab1abf Mon Sep 17 00:00:00 2001 From: imadhammani Date: Thu, 12 Dec 2024 13:55:37 +0100 Subject: [PATCH 1/5] XCore intersectMesh: added pbrt ray-triangle intersection (WIP) --- Cassiopee/XCore/XCore/intersectMesh/AABB.h | 11 ++ Cassiopee/XCore/XCore/intersectMesh/point.h | 24 ++++ Cassiopee/XCore/XCore/intersectMesh/ray.cpp | 112 ++++++++++++++++++ Cassiopee/XCore/XCore/intersectMesh/ray.h | 13 +- .../removeIntersectingKPlanes.cpp | 45 ++++++- Cassiopee/XCore/XCore/intersectMesh/smesh.cpp | 2 + Cassiopee/XCore/XCore/intersectMesh/smesh.h | 19 ++- .../XCore/intersectMesh/smesh_extract.cpp | 2 +- .../XCore/XCore/intersectMesh/smesh_geom.cpp | 51 ++++++++ .../XCore/intersectMesh/smesh_locate.cpp | 23 ++-- 10 files changed, 272 insertions(+), 30 deletions(-) diff --git a/Cassiopee/XCore/XCore/intersectMesh/AABB.h b/Cassiopee/XCore/XCore/intersectMesh/AABB.h index 7e033faea..e02833137 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/AABB.h +++ b/Cassiopee/XCore/XCore/intersectMesh/AABB.h @@ -30,6 +30,17 @@ struct AABB { E_Float dx; E_Float dy; E_Float dz; + + bool is_point_inside(const E_Float p[3]) const + { + return (p[0] >= xmin && p[0] <= xmax && p[1] >= ymin && p[1] <= ymax && + p[2] >= zmin && p[2] <= zmax); + } + + void print() const + { + printf("[%f %f %f] - [%f %f %f]\n", xmin, ymin, zmin, xmax, ymax, zmax); + } }; const AABB AABB_HUGE = {EFLOATMIN, EFLOATMIN, EFLOATMIN, diff --git a/Cassiopee/XCore/XCore/intersectMesh/point.h b/Cassiopee/XCore/XCore/intersectMesh/point.h index a734d73c9..8f12239b2 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/point.h +++ b/Cassiopee/XCore/XCore/intersectMesh/point.h @@ -22,6 +22,30 @@ struct Point { E_Float x, y, z; + bool operator<(const Point &p) const + { + return (x < p.x) || (x == p.x && y < p.y) || + (x == p.x && y == p.y && z < p.z); + } + Point operator-(const Point &p) const + { + return {x-p.x, y-p.y, z-p.z}; + } + E_Float &operator[](E_Int idx) + { + return (idx == 0) ? x : ((idx == 1) ? y : z); + } + E_Float operator[](E_Int idx) const + { + return (idx == 0) ? x : ((idx == 1) ? y : z); + } +}; + +typedef Point Point3D; + +struct Point2D +{ + E_Float x, y; }; struct PointLoc { diff --git a/Cassiopee/XCore/XCore/intersectMesh/ray.cpp b/Cassiopee/XCore/XCore/intersectMesh/ray.cpp index b907e38fb..ff66ee598 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/ray.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/ray.cpp @@ -19,6 +19,118 @@ #include "ray.h" #include "primitives.h" #include "AABB.h" +#include "point.h" + +Ray::Ray(E_Float px, E_Float py, E_Float pz, E_Float dx, E_Float dy, E_Float dz) +{ + o[0] = px, o[1] = py, o[2] = pz; + d[0] = dx, d[1] = dy, d[2] = dz; + inv_dx = (d[0] != 0.0) ? 1.0/d[0] : EFLOATMAX; + inv_dy = (d[1] != 0.0) ? 1.0/d[1] : EFLOATMAX; + inv_dz = (d[2] != 0.0) ? 1.0/d[2] : EFLOATMAX; + + E_Float adx = fabs(d[0]); + E_Float ady = fabs(d[1]); + E_Float adz = fabs(d[2]); + + if (adx > ady && adx > adz) kz = 0; + else if (ady > adz) kz = 1; + else kz = 2; + + kx = kz+1; if (kx == 3) kx = 0; + ky = kx+1; if (ky == 3) ky = 0; + + if (d[kz] < 0.0) std::swap(kx, ky); + + Sx = d[kx] / d[kz]; + Sy = d[ky] / d[kz]; + Sz = 1.0 / d[kz]; +} + +bool Ray::intersect_triangle(const Point &tA, const Point &tB, + const Point &tC, E_Float &u, E_Float &v, E_Float &w, E_Float &t, + E_Float &x, E_Float &y, E_Float &z) const +{ + const Point O = (Point){o[0], o[1], o[2]}; + const Point A = tA - O; + const Point B = tB - O; + const Point C = tC - O; + + const E_Float Ax = A[kx] - Sx*A[kz]; + const E_Float Ay = A[ky] - Sy*A[kz]; + const E_Float Bx = B[kx] - Sx*B[kz]; + const E_Float By = B[ky] - Sy*B[kz]; + const E_Float Cx = C[kx] - Sx*C[kz]; + const E_Float Cy = C[ky] - Sy*C[kz]; + + u = Cx*By - Cy*Bx; + v = Ax*Cy - Ay*Cx; + w = Bx*Ay - By*Ax; + + //if (u == 0.0 || v == 0.0 || w == 0.0) { + // fprintf(stderr, "Warning: unsufficient floating-point precision. " + // "Results may be weird.\n"); + //} + + if ((u < 0.0 || v < 0.0 || w < 0.0) && (u > 0.0 || v > 0.0 || w > 0.0)) + return false; + + E_Float det = u + v + w; + + if (det == 0.0) return false; + + const E_Float Az = Sz*A[kz]; + const E_Float Bz = Sz*B[kz]; + const E_Float Cz = Sz*C[kz]; + t = u*Az + v*Bz + w*Cz; + + if (det < 0.0 && t >= 0.0) return false; + else if (det > 0.0 && t <= 0.0) return false; + + const E_Float inv_det = 1.0 / det; + u = u * inv_det; + v = v * inv_det; + w = w * inv_det; + t = t * inv_det; + + x = u * tA.x + v * tB.x + w * tC.x; + y = u * tA.y + v * tB.y + w * tC.y; + z = u * tA.z + v * tB.z + w * tC.z; + + return true; +} + +bool Ray::intersect_AABB(const AABB &box) const +{ + E_Float tmin = 0; + E_Float tmax = EFLOATMAX; + + E_Float boxMin[3] = { box.xmin, box.ymin, box.zmin }; + E_Float boxMax[3] = { box.xmax, box.ymax, box.zmax }; + + for (int i = 0; i < 3; i++) { + E_Float O = o[i]; + E_Float D = d[i]; + E_Float bmin = boxMin[i]; + E_Float bmax = boxMax[i]; + + if (d != 0) { + E_Float t1 = (bmin - O) / D; + E_Float t2 = (bmax - O) / D; + + if (t1 > t2) { E_Float temp = t1; t1 = t2; t2 = temp; } + + tmin = (t1 > tmin) ? t1 : tmin; + tmax = (t2 < tmax) ? t2 : tmax; + + if (tmin > tmax) return false; // No intersection + } else { + if (O < bmin || O > bmax) return false; // Parallel and outside slab + } + } + + return true; +} E_Int ray_point_orient(const E_Float o[3], const E_Float d[3], const E_Float fN[3], E_Float px, E_Float py, E_Float pz) diff --git a/Cassiopee/XCore/XCore/intersectMesh/ray.h b/Cassiopee/XCore/XCore/intersectMesh/ray.h index a80170a31..e5a9d5fdf 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/ray.h +++ b/Cassiopee/XCore/XCore/intersectMesh/ray.h @@ -28,15 +28,20 @@ struct Vec3 { }; struct Ray { - Point org; - Vec3 dir; + E_Float o[3]; + E_Float d[3]; E_Int kx, ky, kz; E_Float Sx, Sy, Sz; + E_Float inv_dx, inv_dy, inv_dz; Ray(Point O, Vec3 D); + Ray(E_Float px, E_Float py, E_Float pz, E_Float dx, E_Float dy, E_Float dz); - E_Int intersect_triangle(const Point &a, const Point &b, const Point &c, - TriangleIntersection &TI); + bool intersect_AABB(const AABB &box) const; + + bool intersect_triangle(const Point &a, const Point &b, const Point &c, + E_Float &u, E_Float &v, E_Float &w, E_Float &t, E_Float &x, + E_Float &y, E_Float &z) const; }; E_Int ray_point_orient(const E_Float o[3], const E_Float d[3], diff --git a/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp b/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp index 5e8c66231..c4aa454a5 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp @@ -55,6 +55,8 @@ void construct_faces(std::vector> &faces, NF++; } +#include + PyObject *handle_slave(const IMesh *M, const Smesh &Mf, Karray& sarray) { E_Int ni = sarray.ni; @@ -72,6 +74,11 @@ PyObject *handle_slave(const IMesh *M, const Smesh &Mf, Karray& sarray) // Indices of points to be projected std::vector proj_points; + //using std::chrono::high_resolution_clock; + //using std::chrono::duration_cast; + //using std::chrono::milliseconds; + //auto t1 = high_resolution_clock::now(); + for (E_Int j = 0; j < nj; j++) { for (E_Int i = 0; i < ni; i++) { @@ -86,7 +93,8 @@ PyObject *handle_slave(const IMesh *M, const Smesh &Mf, Karray& sarray) E_Float py = Ys[idx]; E_Float pz = Zs[idx]; - if (M->is_point_inside(px, py, pz)) { + if (Mf.is_point_inside(px, py, pz)) { + //if (M->is_point_inside(px, py, pz)) { kmax[base] = k-1; // Cache the point to be projected @@ -98,6 +106,9 @@ PyObject *handle_slave(const IMesh *M, const Smesh &Mf, Karray& sarray) } } + //assert(was_outside); + + if (!was_inside) { // i-j line completely outside of M // Projection points is the last point @@ -108,6 +119,10 @@ PyObject *handle_slave(const IMesh *M, const Smesh &Mf, Karray& sarray) } } + //auto t2 = high_resolution_clock::now(); + //auto ms_int = duration_cast(t2-t1); + //std::cout << ms_int.count() << "ms\n"; + assert(proj_points.size() == (size_t)nij); //point_write("proj_points.im", Xs, Ys, Zs, proj_points); @@ -544,6 +559,7 @@ PyObject *handle_slave2(IMesh *M, Karray& sarray, E_Int kmax) } #include "smesh.h" +#include PyObject *K_XCORE::removeIntersectingKPlanes(PyObject *self, PyObject *args) { @@ -562,13 +578,36 @@ PyObject *K_XCORE::removeIntersectingKPlanes(PyObject *self, PyObject *args) IMesh *M = (IMesh *)PyCapsule_GetPointer(MASTER, "IntersectMesh"); M->make_skin(); - M->make_bbox(); - M->hash_skin(); + //M->make_bbox(); + //M->hash_skin(); Smesh Mf = Smesh::Smesh_from_mesh_skin(*M, M->skin, false); printf("Mf: %d tris\n", Mf.nf); Mf.make_fcenters(); Mf.make_BVH(); + Mf.box.print(); + + /* + E_Float lo = -1.0; + E_Float up = 1.0; + std::uniform_real_distribution unif(lo, up); + std::default_random_engine re; + + E_Float px = 0, py = 0, pz = 0; + for (int i = 0; i < 1; i++) { + //E_Float dx = unif(re); + //E_Float dy = unif(re); + //E_Float dz = unif(re); + //E_Float n = sqrt(dx*dx + dy*dy + dz*dz); + //dx /= n; dy /= n; dz /= n; + E_Float dx = 1.0, dy = 0.0, dz = 0.0; + Ray ray(px, py, pz, dx, dy, dz); + bool inside = Mf.is_point_inside(ray); + assert(inside); + } + + return Py_None; + */ E_Int nslaves = PyList_Size(SLAVES); E_Int i, ret; diff --git a/Cassiopee/XCore/XCore/intersectMesh/smesh.cpp b/Cassiopee/XCore/XCore/intersectMesh/smesh.cpp index db3e74f1a..5c6690ec9 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/smesh.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/smesh.cpp @@ -194,6 +194,7 @@ Smesh::Smesh(const IMesh &M, const std::vector &fids, bool check_Euler_) make_edges(); make_point_faces(); make_point_edges(); + make_bbox(); } Smesh::Smesh(const Smesh &Mf, const std::set &fids, bool check_Euler_) @@ -253,6 +254,7 @@ Smesh::Smesh(const Smesh &Mf, const std::set &fids, bool check_Euler_) make_edges(); make_point_faces(); make_point_edges(); + make_bbox(); } o_edge::o_edge(E_Int P, E_Int Q) diff --git a/Cassiopee/XCore/XCore/intersectMesh/smesh.h b/Cassiopee/XCore/XCore/intersectMesh/smesh.h index 55c2832e3..3194fcd68 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/smesh.h +++ b/Cassiopee/XCore/XCore/intersectMesh/smesh.h @@ -27,22 +27,14 @@ #include "u_edge.h" #include "point.h" #include "BVH.h" +#include "AABB.h" #define LEFT 0 #define RIGHT 1 -struct Point3D -{ - E_Float x, y, z; -}; - -struct Point2D -{ - E_Float x, y; -}; - struct IMesh; -struct AABB; +struct Ray; +struct HitData; struct o_edge { E_Int p, q; @@ -111,6 +103,9 @@ struct Smesh { PointLoc &ploc, E_Float min_pdist) const; bool is_point_a_polygon_vertex(E_Float x, E_Float y, E_Float z, E_Int fid, PointLoc &ploc, E_Float min_pdist) const; + bool is_point_inside(E_Float px, E_Float py, E_Float pz) const; + bool is_point_inside(const Ray &ray) const; + void intersect_ray(const Ray &ray, E_Int node_idx, HitData &hit_data) const; void make_fcenters(); void make_fnormals(); @@ -130,7 +125,7 @@ struct Smesh { inline void compute_face_center(E_Int fid); // Hash - + AABB box; E_Int NX, NY, NZ, NXY, NXYZ; E_Float xmin, xmax, ymin, ymax, zmin, zmax; E_Float HX, HY, HZ; diff --git a/Cassiopee/XCore/XCore/intersectMesh/smesh_extract.cpp b/Cassiopee/XCore/XCore/intersectMesh/smesh_extract.cpp index f241e09c9..2c8783e7e 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/smesh_extract.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/smesh_extract.cpp @@ -169,7 +169,7 @@ std::set Smesh::extract_covering_faces(const Smesh &Sf, E_Float cur_pos[3] = {px, py, pz}; E_Int walk = 0; - E_Int max_walks = 20; + E_Int max_walks = 100; while (!found_tail && walk <= max_walks) { diff --git a/Cassiopee/XCore/XCore/intersectMesh/smesh_geom.cpp b/Cassiopee/XCore/XCore/intersectMesh/smesh_geom.cpp index 47b619122..588b2ff13 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/smesh_geom.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/smesh_geom.cpp @@ -19,6 +19,57 @@ #include "smesh.h" #include "primitives.h" #include "io.h" +#include "ray.h" + +struct HitData { + std::set tested; + std::set fids; + std::set hits; +}; + +void Smesh::intersect_ray(const Ray &ray, E_Int node_idx, HitData &hit_data) const +{ + const BVH_node &node = bvh_nodes[node_idx]; + if (!ray.intersect_AABB(node.box)) return; + if (!node.is_leaf()) { + intersect_ray(ray, node.left_node, hit_data); + intersect_ray(ray, node.left_node+1, hit_data); + return; + } + for (E_Int i = 0; i < node.tri_count; i++) { + E_Int tri = tri_idx[node.first_tri_idx+i]; + hit_data.tested.insert(tri); + const auto &pn = F[tri]; + //write_face("tri.im", tri); + E_Float u, v, w, t, x, y, z; + u = v = w = 2.0; + bool hit = ray.intersect_triangle({X[pn[0]], Y[pn[0]], Z[pn[0]]}, + {X[pn[1]], Y[pn[1]], Z[pn[1]]}, {X[pn[2]], Y[pn[2]], Z[pn[2]]}, + u, v, w, t, x, y, z); + + if (hit) { + hit_data.fids.insert(tri); + hit_data.hits.insert({x, y, z}); + } + } +} + +bool Smesh::is_point_inside(E_Float px, E_Float py, E_Float pz) const +{ + Ray ray(px, py, pz, 1, 0, 0); + if (!box.is_point_inside(ray.o)) return false; + HitData hit_data; + intersect_ray(ray, 0, hit_data); + return hit_data.hits.size() % 2 == 1; +} + +bool Smesh::is_point_inside(const Ray &ray) const +{ + if (!box.is_point_inside(ray.o)) return false; + HitData hit_data; + intersect_ray(ray, 0, hit_data); + return hit_data.hits.size() % 2 == 1; +} void Smesh::replace_by_projections(const std::vector &pids, const std::vector &plocs) diff --git a/Cassiopee/XCore/XCore/intersectMesh/smesh_locate.cpp b/Cassiopee/XCore/XCore/intersectMesh/smesh_locate.cpp index 7e3eb3eb4..05392b243 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/smesh_locate.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/smesh_locate.cpp @@ -417,16 +417,19 @@ void Smesh::make_bbox() E_Float dy = ymax - ymin; E_Float dz = zmax - zmin; - xmin = xmin - dx*0.01; - ymin = ymin - dy*0.01; - zmin = zmin - dz*0.01; - xmax = xmax + dx*0.01; - ymax = ymax + dy*0.01; - zmax = zmax + dz*0.01; - - HX = (dx != 0) ? (xmax - xmin) / NX : 1; - HY = (dy != 0) ? (ymax - ymin) / NY : 1; - HZ = (dz != 0) ? (zmax - zmin) / NZ : 1; + //xmin = xmin - dx*0.01; + //ymin = ymin - dy*0.01; + //zmin = zmin - dz*0.01; + //xmax = xmax + dx*0.01; + //ymax = ymax + dy*0.01; + //zmax = zmax + dz*0.01; + + //HX = (dx != 0) ? (xmax - xmin) / NX : 1; + //HY = (dy != 0) ? (ymax - ymin) / NY : 1; + //HZ = (dz != 0) ? (zmax - zmin) / NZ : 1; + + box.xmin = xmin, box.ymin = ymin, box.zmin = zmin; + box.xmax = xmax, box.ymax = ymax, box.zmax = zmax; } inline From d865f4c02de2b59fa0a5e02f5b7a613cae779f2e Mon Sep 17 00:00:00 2001 From: imadhammani Date: Tue, 17 Dec 2024 15:40:54 +0100 Subject: [PATCH 2/5] XCore intersectMesh: OK with relaxed Enoval --- Cassiopee/XCore/XCore/intersectMesh/mesh.cpp | 5 +- Cassiopee/XCore/XCore/intersectMesh/mesh.h | 3 +- Cassiopee/XCore/XCore/intersectMesh/o_edge.h | 40 +++++ .../removeIntersectingKPlanes.cpp | 4 +- .../XCore/XCore/intersectMesh/triangulate.cpp | 143 ++++++++++-------- 5 files changed, 126 insertions(+), 69 deletions(-) create mode 100644 Cassiopee/XCore/XCore/intersectMesh/o_edge.h diff --git a/Cassiopee/XCore/XCore/intersectMesh/mesh.cpp b/Cassiopee/XCore/XCore/intersectMesh/mesh.cpp index 73630bf8f..c63124557 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/mesh.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/mesh.cpp @@ -328,8 +328,6 @@ void IMesh::make_skin() owner.resize(nf, -1); neigh.resize(nf, -1); - std::vector count(nf, 0); - for (E_Int cid = 0; cid < nc; cid++) { const auto &pf = C[cid]; for (auto fid : pf) { @@ -339,8 +337,7 @@ void IMesh::make_skin() } for (E_Int fid = 0; fid < nf; fid++) { - if (neigh[fid] == -1) - skin.push_back(fid); + if (neigh[fid] == -1) skin.push_back(fid); } } diff --git a/Cassiopee/XCore/XCore/intersectMesh/mesh.h b/Cassiopee/XCore/XCore/intersectMesh/mesh.h index c0bd82082..f03bff7b2 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/mesh.h +++ b/Cassiopee/XCore/XCore/intersectMesh/mesh.h @@ -124,8 +124,7 @@ struct IMesh { void triangulate(const std::vector &fids); - std::vector triangulate_skin(const std::vector &bcs_in, - const std::unordered_map &fid_to_bc); + void triangulate_skin(std::unordered_map &fid_to_bc); void triangulate_skin(); diff --git a/Cassiopee/XCore/XCore/intersectMesh/o_edge.h b/Cassiopee/XCore/XCore/intersectMesh/o_edge.h new file mode 100644 index 000000000..facafdd85 --- /dev/null +++ b/Cassiopee/XCore/XCore/intersectMesh/o_edge.h @@ -0,0 +1,40 @@ +/* + Copyright 2013-2024 Onera. + + This file is part of Cassiopee. + + Cassiopee is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Cassiopee is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Cassiopee. If not, see . +*/ +#pragma once + +#include "xcore.h" + +struct o_edge { + E_Int p, q; + + o_edge(E_Int P, E_Int Q) + : p(P), q(Q) + {} +}; + +struct o_edge_cmp { + bool operator()(const o_edge &e, const o_edge &f) const + { + E_Int e_p = std::min(e.p, e.q); + E_Int e_q = std::max(e.p, e.q); + E_Int f_p = std::min(f.p, f.q); + E_Int f_q = std::max(f.p, f.q); + return (e_p < f_p) || (e_p == f_p && e_q < f_q); + } +}; diff --git a/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp b/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp index c4aa454a5..867a09776 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp @@ -95,7 +95,7 @@ PyObject *handle_slave(const IMesh *M, const Smesh &Mf, Karray& sarray) if (Mf.is_point_inside(px, py, pz)) { //if (M->is_point_inside(px, py, pz)) { - kmax[base] = k-1; + kmax[base] = k-2; // Cache the point to be projected E_Int proj_id = base + nij*kmax[base]; @@ -112,7 +112,7 @@ PyObject *handle_slave(const IMesh *M, const Smesh &Mf, Karray& sarray) if (!was_inside) { // i-j line completely outside of M // Projection points is the last point - kmax[base] = nk-1; + kmax[base] = nk-2; E_Int proj_id = base + nij*kmax[base]; proj_points.push_back(proj_id); } diff --git a/Cassiopee/XCore/XCore/intersectMesh/triangulate.cpp b/Cassiopee/XCore/XCore/intersectMesh/triangulate.cpp index e7d768418..2ac6133b3 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/triangulate.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/triangulate.cpp @@ -63,15 +63,28 @@ PyObject *K_XCORE::triangulate_skin(PyObject *self, PyObject *args) } } - std::vector new_bcs = M.triangulate_skin(bcs_in, fid_to_bc); + M.triangulate_skin(fid_to_bc); + + std::vector> new_bcs(nbcs); + for (const auto &it : fid_to_bc) { + new_bcs[it.second].push_back(it.first); + } PyObject *out = PyList_New(0); PyList_Append(out, M.export_karray()); PyObject *bcs_out = PyList_New(0); + for (const auto &new_bc : new_bcs) { - PyList_Append(bcs_out, (PyObject *)new_bc); - Py_DECREF(new_bc); + npy_intp dims[2]; + dims[0] = (npy_intp)new_bc.size(); + dims[1] = 1; + PyArrayObject *arr = (PyArrayObject *)PyArray_SimpleNew(1, dims, E_NPY_INT); + E_Int *ptr = (E_Int *)PyArray_DATA(arr); + memcpy(ptr, new_bc.data(), new_bc.size()*sizeof(E_Int)); + PyList_Append(bcs_out, (PyObject *)arr); + Py_DECREF(arr); } + PyList_Append(out, bcs_out); Py_DECREF(bcs_out); return out; @@ -81,98 +94,106 @@ PyObject *K_XCORE::triangulate_skin(PyObject *self, PyObject *args) } } -std::vector IMesh::triangulate_skin(const std::vector &bcs_in, - const std::unordered_map &fid_to_bc) +void IMesh::triangulate_skin(std::unordered_map &fid_to_bc) { - E_Int NF = nf; - - owner.resize(nf + skin.size(), -1); - neigh.resize(nf + skin.size(), -1); - - std::vector new_bcs(bcs_in.size()); - std::vector ptrs(bcs_in.size()); - std::vector counts(bcs_in.size(), 0); - - for (size_t i = 0; i < bcs_in.size(); i++) { - E_Int new_size = bcs_in[i].size * 2; - npy_intp dims[2]; - dims[0] = (npy_intp)new_size; - dims[1] = 1; - new_bcs[i] = (PyArrayObject *)PyArray_SimpleNew(1, dims, E_NPY_INT); - ptrs[i] = (E_Int *)PyArray_DATA(new_bcs[i]); + E_Int to_tri = 0; + + for (E_Int fid : skin) { + const auto &pn = F[fid]; + if (pn.size() == 3) continue; + assert(pn.size() == 4); + to_tri++; } - for (auto fid : skin) { - assert(neigh[fid] == -1); + F.resize(nf + 3*to_tri); + X.resize(np + to_tri, 0); + Y.resize(np + to_tri, 0); + Z.resize(np + to_tri, 0); + for (E_Int fid : skin) { const auto &pn = F[fid]; if (pn.size() == 3) continue; - - assert(pn.size() == 4); + for (E_Int pid : pn) { + X[np] += X[pid]; + Y[np] += Y[pid]; + Z[np] += Z[pid]; + } + X[np] *= 0.25; + Y[np] *= 0.25; + Z[np] *= 0.25; E_Int nodes[4] = {pn[0], pn[1], pn[2], pn[3]}; - F.push_back({nodes[2], nodes[3], nodes[0]}); - F[fid] = {nodes[0], nodes[1], nodes[2]}; - assert(F[fid].size() == 3); + F[fid] = {nodes[0], nodes[1], np}; + F[nf] = {nodes[1], nodes[2], np}; + F[nf+1] = {nodes[2], nodes[3], np}; + F[nf+2] = {nodes[3], nodes[0], np}; E_Int own = owner[fid]; auto &pf = C[own]; - pf.push_back(nf); - - owner[nf] = own; + size_t stride = pf.size(); + pf.resize(stride + 3); + for (int i = 0; i < 3; i++) pf[stride+i] = nf+i; - // Which bc (if any) does fid belong to? + // One-based auto it = fid_to_bc.find(fid+1); if (it != fid_to_bc.end()) { - E_Int bc_id = it->second; - ptrs[bc_id][counts[bc_id]++] = fid+1; - ptrs[bc_id][counts[bc_id]++] = nf+1; + for (int i = 0; i < 3; i++) fid_to_bc[nf+i+1] = it->second; } - nf++; + np += 1; + nf += 3; } - for (E_Int i = NF; i < nf; i++) - skin.push_back(i); - - for (size_t i = 0; i < counts.size(); i++) { - assert(counts[i] == 2*bcs_in[i].size); - } - - return new_bcs; + assert((size_t)nf == F.size()); + assert((size_t)np == X.size()); } void IMesh::triangulate_skin() { - E_Int NF = nf; + E_Int to_tri = 0; - owner.resize(nf + skin.size(), -1); - neigh.resize(nf + skin.size(), -1); + for (E_Int fid : skin) { + const auto &pn = F[fid]; + if (pn.size() == 3) continue; + assert(pn.size() == 4); + to_tri++; + } - for (auto fid : skin) { - assert(neigh[fid] == -1); + F.resize(nf + 3*to_tri); + X.resize(np + to_tri, 0); + Y.resize(np + to_tri, 0); + Z.resize(np + to_tri, 0); + for (E_Int fid : skin) { const auto &pn = F[fid]; if (pn.size() == 3) continue; - - assert(pn.size() == 4); + for (E_Int pid : pn) { + X[np] += X[pid]; + Y[np] += Y[pid]; + Z[np] += Z[pid]; + } + X[np] *= 0.25; + Y[np] *= 0.25; + Z[np] *= 0.25; E_Int nodes[4] = {pn[0], pn[1], pn[2], pn[3]}; - F.push_back({nodes[2], nodes[3], nodes[0]}); - F[fid] = {nodes[0], nodes[1], nodes[2]}; - assert(F[fid].size() == 3); - + F[fid] = {nodes[0], nodes[1], np}; + F[nf] = {nodes[1], nodes[2], np}; + F[nf+1] = {nodes[2], nodes[3], np}; + F[nf+2] = {nodes[3], nodes[0], np}; + E_Int own = owner[fid]; auto &pf = C[own]; - pf.push_back(nf); - - owner[nf] = own; + size_t stride = pf.size(); + pf.resize(stride + 3); + for (int i = 0; i < 3; i++) pf[stride+i] = nf+i; - nf++; + np += 1; + nf += 3; } - for (E_Int i = NF; i < nf; i++) - skin.push_back(i); + assert((size_t)nf == F.size()); + assert((size_t)np == X.size()); } \ No newline at end of file From 0d61bb6b4106e508f49eb7e3d8ce856ae8afb226 Mon Sep 17 00:00:00 2001 From: imadhammani Date: Tue, 7 Jan 2025 17:53:24 +0100 Subject: [PATCH 3/5] XCore intersectMesh: update --- Cassiopee/XCore/XCore/PyTree.py | 84 +++++++-- .../intersectMesh/IntersectMesh_Init.cpp | 9 + .../intersectMesh/backup_triangulate.cpp | 177 +++++++++++++++++ .../intersectMesh/check_manifoldness.cpp | 28 +++ .../XCore/XCore/intersectMesh/icapsule.h | 5 + .../XCore/intersectMesh/icapsule_adapt.cpp | 99 ++++++++++ .../XCore/intersectMesh/icapsule_extract.cpp | 112 +++++++++++ .../{icapsule.cpp => icapsule_init.cpp} | 171 ++++++++++++++--- .../intersectMesh/icapsule_intersect.cpp | 63 +++++++ .../XCore/XCore/intersectMesh/initSurf.cpp | 178 ++++++++++++++++++ Cassiopee/XCore/XCore/intersectMesh/mesh.h | 2 + .../XCore/XCore/intersectMesh/precise.cpp | 127 +++++++++++++ Cassiopee/XCore/XCore/intersectMesh/precise.h | 81 ++++++++ Cassiopee/XCore/XCore/intersectMesh/ray.cpp | 98 +--------- .../removeIntersectingKPlanes.cpp | 41 ++-- .../XCore/XCore/intersectMesh/smesh_geom.cpp | 7 +- Cassiopee/XCore/XCore/xcore.cpp | 6 + Cassiopee/XCore/XCore/xcore.h | 7 + Cassiopee/XCore/srcs.py | 8 +- 19 files changed, 1137 insertions(+), 166 deletions(-) create mode 100644 Cassiopee/XCore/XCore/intersectMesh/backup_triangulate.cpp create mode 100644 Cassiopee/XCore/XCore/intersectMesh/check_manifoldness.cpp create mode 100644 Cassiopee/XCore/XCore/intersectMesh/icapsule_adapt.cpp create mode 100644 Cassiopee/XCore/XCore/intersectMesh/icapsule_extract.cpp rename Cassiopee/XCore/XCore/intersectMesh/{icapsule.cpp => icapsule_init.cpp} (68%) create mode 100644 Cassiopee/XCore/XCore/intersectMesh/icapsule_intersect.cpp create mode 100644 Cassiopee/XCore/XCore/intersectMesh/initSurf.cpp create mode 100644 Cassiopee/XCore/XCore/intersectMesh/precise.cpp create mode 100644 Cassiopee/XCore/XCore/intersectMesh/precise.h diff --git a/Cassiopee/XCore/XCore/PyTree.py b/Cassiopee/XCore/XCore/PyTree.py index 389519fdf..22d4de6b4 100644 --- a/Cassiopee/XCore/XCore/PyTree.py +++ b/Cassiopee/XCore/XCore/PyTree.py @@ -482,6 +482,74 @@ def extractFacesFromPointTag(t, tag_name): arr = C.getFields(I.__GridCoordinates__, z, api=3)[0] return xcore.extractFacesFromPointTag(arr, tag[1]) +############################################################################### + +def icapsule_init2(): + return xcore.icapsule_init2() + +def icapsule_set_master(IC, m): + zones = I.getNodesFromType(m, 'Zone_t') + if len(zones) != 1: raise ValueError('Master should be one zone.') + zm = zones[0] + marr = C.getFields(I.__GridCoordinates__, zm, api=3)[0] + ctag = I.getNodeFromName(zm, 'keep') + if ctag is None: raise ValueError('Missing cell tags') + return xcore.icapsule_set_master(IC, marr, ctag[1]) + +def icapsule_set_slaves(IC, slaves): + sarrs = [] + ptags = [] + ctags = [] + + for slave in slaves: + + bases = I.getBases(slave) + + for base in bases: + zones = I.getZones(base) + for zone in zones: + sarr = C.getFields(I.__GridCoordinates__, zone, api=3)[0] + sarrs.append(sarr) + ptag = I.getNodeFromName(zone, 'tag') + if ptag is None: raise ValueError('Missing point tags.') + ptags.append(ptag[1]) + ctag = I.getNodeFromName(zone, 'keep') + if ctag is None: raise ValueError('Missing cell tags.') + ctags.append(ctag[1]) + + return xcore.icapsule_set_slaves(IC, sarrs, ptags, ctags) + +def icapsule_adapt2(IC): + return xcore.icapsule_adapt2(IC) + +def icapsule_intersect2(IC): + return xcore.icapsule_intersect2(IC) + +def icapsule_extract_master(IC): + marr, ctag = xcore.icapsule_extract_master(IC) + zm = I.createZoneNode("master", marr) + + cont = I.createUniqueChild(zm, I.__FlowSolutionCenters__, 'FlowSolution_t') + I._createUniqueChild(cont, 'GridLocation', 'GridLocation_t', value='CellCenter', ) + I.newDataArray("keep", value=ctag, parent=cont) + return zm + +def icapsule_extract_slaves(IC): + sarrs, ctags = xcore.icapsule_extract_slaves(IC) + assert(len(sarrs) == len(ctags)) + zones = [] + for i in range(len(sarrs)): + zs = I.createZoneNode("slave_"+str(i), sarrs[i]) + + cont = I.createUniqueChild(zs, I.__FlowSolutionCenters__, 'FlowSolution_t') + I._createUniqueChild(cont, 'GridLocation', 'GridLocation_t', value='CellCenter', ) + I.newDataArray("keep", value=ctags[i], parent=cont) + + zones.append(zs) + return zones + +############################################################################### + def icapsule_init(mp, sp): zm = I.getZones(mp)[0] marr = C.getFields(I.__GridCoordinates__, zm, api=3)[0] @@ -545,23 +613,9 @@ def icapsule_intersect(ma, sa): ts = C.newPyTree(['Base', slave_zones]) return tm, ts - -def icapsule_extract_master(IC): - marr = xcore.icapsule_extract_master(IC) - zm = I.createZoneNode("master", marr) - return zm - def icapsule_extract_slave(IC, index=0): sarr = xcore.icapsule_extract_slave(IC, index) - zs = I.createZoneNode("slave", sarr) - return zs - -def icapsule_extract_slaves(IC): - sarrs = xcore.icapsule_extract_slaves(IC) - zs = [] - for i in range(len(sarrs)): - z = I.createZoneNode("slave"+str(i), sarrs[i]) - zs.append(z) + zs = I.createZoneNode("slave_"+str(index), sarr) return zs def triangulateSkin(m): diff --git a/Cassiopee/XCore/XCore/intersectMesh/IntersectMesh_Init.cpp b/Cassiopee/XCore/XCore/intersectMesh/IntersectMesh_Init.cpp index 928228b8d..ec00bc039 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/IntersectMesh_Init.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/IntersectMesh_Init.cpp @@ -59,6 +59,15 @@ PyObject *K_XCORE::IntersectMesh_Init(PyObject *self, PyObject *args) for (E_Int i = 0; i < M->nc; i++) M->ctag[i] = (int)tags[i]; } + // Init surface mesh data + M->make_skin(); + M->Mf = Smesh::Smesh_from_mesh_skin(*M, M->skin, false); + printf("Mf: %d tris\n", M->Mf.nf); + M->Mf.make_fcenters(); + M->Mf.make_BVH(); + printf("Mesh bounds: "); + M->Mf.box.print(); + // Clean-up Karray_free_ngon(karray); diff --git a/Cassiopee/XCore/XCore/intersectMesh/backup_triangulate.cpp b/Cassiopee/XCore/XCore/intersectMesh/backup_triangulate.cpp new file mode 100644 index 000000000..97dac9ec1 --- /dev/null +++ b/Cassiopee/XCore/XCore/intersectMesh/backup_triangulate.cpp @@ -0,0 +1,177 @@ +/* + Copyright 2013-2024 Onera. + + This file is part of Cassiopee. + + Cassiopee is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Cassiopee is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Cassiopee. If not, see . +*/ +#include "mesh.h" +#include "common/Karray.h" + +PyObject *K_XCORE::triangulate_skin(PyObject *self, PyObject *args) +{ + PyObject *MESH; + if (!PYPARSETUPLE_(args, O_, &MESH)) { + RAISE("Bad input."); + return NULL; + } + + Karray array; + if (Karray_parse_ngon(MESH, array) != 0) { + RAISE("Mesh should be an NGon."); + return NULL; + } + + IMesh M(array); + M.make_skin(); + M.triangulate_skin(); + return M.export_karray(); +} + +/* +std::vector IMesh::triangulate_skin(const std::vector &bcs_in, + const std::unordered_map &fid_to_bc) +{ + E_Int NF = nf; + + owner.resize(nf + skin.size(), -1); + neigh.resize(nf + skin.size(), -1); + + std::vector new_bcs(bcs_in.size()); + std::vector ptrs(bcs_in.size()); + std::vector counts(bcs_in.size(), 0); + + for (size_t i = 0; i < bcs_in.size(); i++) { + E_Int new_size = bcs_in[i].size * 2; + npy_intp dims[2]; + dims[0] = (npy_intp)new_size; + dims[1] = 1; + new_bcs[i] = (PyArrayObject *)PyArray_SimpleNew(1, dims, E_NPY_INT); + ptrs[i] = (E_Int *)PyArray_DATA(new_bcs[i]); + } + + for (auto fid : skin) { + assert(neigh[fid] == -1); + + const auto &pn = F[fid]; + if (pn.size() == 3) continue; + + assert(pn.size() == 4); + + E_Int nodes[4] = {pn[0], pn[1], pn[2], pn[3]}; + + F.push_back({nodes[2], nodes[3], nodes[0]}); + F[fid] = {nodes[0], nodes[1], nodes[2]}; + assert(F[fid].size() == 3); + + E_Int own = owner[fid]; + auto &pf = C[own]; + pf.push_back(nf); + + owner[nf] = own; + + // Which bc (if any) does fid belong to? + auto it = fid_to_bc.find(fid+1); + if (it != fid_to_bc.end()) { + E_Int bc_id = it->second; + ptrs[bc_id][counts[bc_id]++] = fid+1; + ptrs[bc_id][counts[bc_id]++] = nf+1; + } + + nf++; + } + + for (E_Int i = NF; i < nf; i++) + skin.push_back(i); + + for (size_t i = 0; i < counts.size(); i++) { + assert(counts[i] == 2*bcs_in[i].size); + } + + return new_bcs; +} +*/ + +void IMesh::triangulate_skin() +{ + owner.clear(); + neigh.clear(); + owner.resize(nf, -1); + neigh.resize(nf, -1); + + for (size_t cid = 0; cid < nc; cid++) { + const auto &pf = C[cid]; + for (auto fid : pf) { + if (owner[fid] == -1) owner[fid] = cid; + else { + assert(neigh[fid] == -1); + neigh[fid] = cid; + } + } + } + + skin.clear(); + for (size_t fid = 0; fid < nf; fid++) { + if (neigh[fid] == -1) skin.push_back(fid); + } + + // One new point per triangulated face + E_Int NP = np + skin.size(); + + // Three new faces per triangulated face + E_Int NF = nf + 3*skin.size(); + + for (auto fid : skin) { + const auto &pn = F[fid]; + if (pn.size() == 3) continue; + + assert(pn.size() == 4); + + E_Int nodes[4] = {pn[0], pn[1], pn[2], pn[3]}; + + E_Float cx, cy, cz; + cx = cy = cz = 0; + for (int i = 0; i < 4; i++) { + cx += X[nodes[i]]; + cy += Y[nodes[i]]; + cz += Z[nodes[i]]; + } + cx *= 0.25; + cy *= 0.25; + cz *= 0.25; + + F[fid] = {nodes[0], nodes[1], np}; + F.push_back({nodes[1], nodes[2], np}); + F.push_back({nodes[2], nodes[3], np}); + F.push_back({nodes[3], nodes[0], np}); + + E_Int own = owner[fid]; + auto &pf = C[own]; + for (int i = 0; i < 3; i++) { + pf.push_back(nf+i); + owner.push_back(own); + neigh.push_back(-1); + } + + X.push_back(cx); + Y.push_back(cy); + Z.push_back(cz); + + nf += 3; + np += 1; + } + + assert(np == NP); + assert(nf == NF); +} diff --git a/Cassiopee/XCore/XCore/intersectMesh/check_manifoldness.cpp b/Cassiopee/XCore/XCore/intersectMesh/check_manifoldness.cpp new file mode 100644 index 000000000..ce026455e --- /dev/null +++ b/Cassiopee/XCore/XCore/intersectMesh/check_manifoldness.cpp @@ -0,0 +1,28 @@ +#include "mesh.h" +#include "common/Karray.h" + +PyObject *K_XCORE::check_manifoldness(PyObject *self, PyObject *args) +{ + PyObject *MESH; + if (!PYPARSETUPLE_(args, O_, &MESH)) { + RAISE("Bad input."); + return NULL; + } + + Karray array; + if (Karray_parse_ngon(MESH, array) != 0) { + RAISE("Input not NGon."); + return NULL; + } + + IMesh M(array); + + bool is_manifold = true; + + M.make_skin(); + Smesh Mf(M, M.skin, false); + + Karray_free_ngon(array); + + return PyBool_FromLong(is_manifold); +} \ No newline at end of file diff --git a/Cassiopee/XCore/XCore/intersectMesh/icapsule.h b/Cassiopee/XCore/XCore/intersectMesh/icapsule.h index e6c81a089..f2fc4c516 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/icapsule.h +++ b/Cassiopee/XCore/XCore/intersectMesh/icapsule.h @@ -24,6 +24,11 @@ struct ICapsule { IMesh M; std::vector Ss; + E_Float NEAR_EDGE_TOL = 1e-3; + E_Float NEAR_VERTEX_TOL = 1e-3; + + ICapsule(){} + ICapsule(const Karray &marray, const std::vector &sarrays, const std::vector &ptags); diff --git a/Cassiopee/XCore/XCore/intersectMesh/icapsule_adapt.cpp b/Cassiopee/XCore/XCore/intersectMesh/icapsule_adapt.cpp new file mode 100644 index 000000000..81d22fcea --- /dev/null +++ b/Cassiopee/XCore/XCore/intersectMesh/icapsule_adapt.cpp @@ -0,0 +1,99 @@ +#include "icapsule.h" + +PyObject *K_XCORE::icapsule_adapt2(PyObject *self, PyObject *args) +{ + PyObject *ICAPSULE; + + if (!PYPARSETUPLE_(args, O_, &ICAPSULE)) { + RAISE("Bad input"); + return NULL; + } + + if (!PyCapsule_IsValid(ICAPSULE, "ICAPSULE")) { + RAISE("Bad ICapsule hook."); + return NULL; + } + + ICapsule *icap = (ICapsule *)PyCapsule_GetPointer(ICAPSULE, "ICAPSULE"); + + auto &M = icap->M; + auto &Ss = icap->Ss; + + M.Mf = Smesh(M, M.skin, false); + auto &Mf = M.Mf; + + for (size_t i = 0; i < Ss.size(); i++) { + printf("Adapting for slave %zu\n", i); + + Mf.make_fcenters(); + Mf.make_fnormals(); + Mf.make_pnormals(); + Mf.make_point_faces(); + Mf.make_BVH(); + + auto &S = Ss[i]; + + Smesh Sf = Smesh::Smesh_from_point_tags(S, S.ptag.data(), true); + + Sf.make_fcenters(); + Sf.make_fnormals(); + Sf.make_pnormals(); + Sf.compute_min_distance_between_points(); + printf("Min dist: %f\n", Sf.min_pdist); + + //Sf.write_ngon("Sf_before_inter.im"); + + // Locate Sf points on Mf faces + auto plocs = Mf.locate2(Sf); + std::vector spids(Sf.np); + for (int i = 0; i < Sf.np; i++) spids[i] = i; + Sf.replace_by_projections(spids, plocs); + + // Extract the initial Mf faces that cover Sf + auto bfaces = Mf.extract_covering_faces(Sf, plocs); + + // Refinement loop + ICapsule::refine(Mf, bfaces, Sf); + + // Reconstruct S + Sf.reconstruct(S); + + // Reconstruct S skin + S.make_skin(); + + // Tag Sf faces + Sf.tag_faces(S); + } + + Mf.reconstruct(M); + + return Py_None; + + /* + PyObject *out = PyList_New(0); + PyList_Append(out, M.export_karray()); + + PyObject *slist = PyList_New(0); + for (const auto &S : Ss) { + PyList_Append(slist, S.export_karray()); + } + PyList_Append(out, slist); + Py_DECREF(slist); + + PyObject *tlist = PyList_New(0); + for (const auto &S : Ss) { + npy_intp dims[2]; + dims[0] = (npy_intp)S.ftag.size(); + dims[1] = 1; + PyArrayObject *arr = (PyArrayObject *)PyArray_SimpleNew(1, dims, E_NPY_INT); + E_Int *ptr = (E_Int *)PyArray_DATA(arr); + for (size_t i = 0; i < S.ftag.size(); i++) ptr[i] = S.ftag[i]+1; + PyList_Append(tlist, (PyObject *)arr); + Py_DECREF(arr); + } + PyList_Append(out, tlist); + Py_DECREF(tlist); + + return out; + */ +} \ No newline at end of file diff --git a/Cassiopee/XCore/XCore/intersectMesh/icapsule_extract.cpp b/Cassiopee/XCore/XCore/intersectMesh/icapsule_extract.cpp new file mode 100644 index 000000000..27e5be9fc --- /dev/null +++ b/Cassiopee/XCore/XCore/intersectMesh/icapsule_extract.cpp @@ -0,0 +1,112 @@ +#include "icapsule.h" + +PyObject *K_XCORE::icapsule_extract_master(PyObject *self, PyObject *args) +{ + PyObject *ICAPSULE; + if (!PYPARSETUPLE_(args, O_, &ICAPSULE)) { + RAISE("Bad input."); + return NULL; + } + + if (!PyCapsule_IsValid(ICAPSULE, "ICAPSULE")) { + RAISE("Bad ICapsule hook."); + return NULL; + } + + ICapsule *icap = (ICapsule *)PyCapsule_GetPointer(ICAPSULE, "ICAPSULE"); + const auto &M = icap->M; + + PyObject *out = PyList_New(0); + + PyObject *Mout = M.export_karray(); + PyList_Append(out, Mout); + Py_DECREF(Mout); + + // Extract cell tags + npy_intp dims[2]; + assert(M.ctag.size() == M.nc); + dims[0] = (npy_intp)M.nc; + dims[1] = 1; + PyArrayObject *arr = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE); + E_Float *ptr = (E_Float *)PyArray_DATA(arr); + for (E_Int i = 0; i < M.nc; i++) ptr[i] = M.ctag[i]; + PyList_Append(out, (PyObject *)arr); + + return out; +} + +PyObject *K_XCORE::icapsule_extract_slaves(PyObject *self, PyObject *args) +{ + PyObject *ICAPSULE; + + if (!PYPARSETUPLE_(args, O_, &ICAPSULE)) { + RAISE("Bad input."); + return NULL; + } + + if (!PyCapsule_IsValid(ICAPSULE, "ICAPSULE")) { + RAISE("Bad ICapsule hook."); + return NULL; + } + + ICapsule *icap = (ICapsule *)PyCapsule_GetPointer(ICAPSULE, "ICAPSULE"); + + PyObject *out = PyList_New(0); + + PyObject *arrays = PyList_New(0); + PyObject *ctags = PyList_New(0); + + for (size_t i = 0; i < icap->Ss.size(); i++) { + const auto &S = icap->Ss[i]; + PyObject *sarray = S.export_karray(); + PyList_Append(arrays, sarray); + Py_DECREF(sarray); + + npy_intp dims[2]; + assert(S.ctag.size() == S.nc); + dims[0] = (npy_intp)S.nc; + dims[1] = 1; + PyArrayObject *arr = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE); + E_Float *ptr = (E_Float *)PyArray_DATA(arr); + for (E_Int i = 0; i < S.nc; i++) ptr[i] = S.ctag[i]; + PyList_Append(ctags, (PyObject *)arr); + Py_DECREF(arr); + } + + PyList_Append(out, arrays); + Py_DECREF(arrays), + PyList_Append(out, ctags); + Py_DECREF(ctags); + + return out; +} + +PyObject *K_XCORE::icapsule_extract_slave(PyObject *self, PyObject *args) +{ + assert(0 && "Unimplemented"); + return Py_None; + /* + PyObject *ICAPSULE; + E_Int INDEX; + if (!PYPARSETUPLE_(args, O_ I_, &ICAPSULE, &INDEX)) { + RAISE("Bad input."); + return NULL; + } + + if (!PyCapsule_IsValid(ICAPSULE, "ICAPSULE")) { + RAISE("Bad ICapsule hook."); + return NULL; + } + + ICapsule *icap = (ICapsule *)PyCapsule_GetPointer(ICAPSULE, "ICAPSULE"); + + if (INDEX >= (E_Int)icap->Ss.size()) { + RAISE("Bad slave index."); + return NULL; + } + + auto Sout = icap->Ss[INDEX].export_karray(); + + return Sout; + */ +} \ No newline at end of file diff --git a/Cassiopee/XCore/XCore/intersectMesh/icapsule.cpp b/Cassiopee/XCore/XCore/intersectMesh/icapsule_init.cpp similarity index 68% rename from Cassiopee/XCore/XCore/intersectMesh/icapsule.cpp rename to Cassiopee/XCore/XCore/intersectMesh/icapsule_init.cpp index f531ca17e..96d79ecae 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/icapsule.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/icapsule_init.cpp @@ -6,10 +6,17 @@ #include "primitives.h" #include "dcel.h" -PyObject *K_XCORE::icapsule_extract_master(PyObject *self, PyObject *args) +PyObject *K_XCORE::icapsule_init2(PyObject *self, PyObject *args) { - PyObject *ICAPSULE; - if (!PYPARSETUPLE_(args, O_, &ICAPSULE)) { + ICapsule *ic = new ICapsule; + PyObject *hook = PyCapsule_New((void *)ic, "ICAPSULE", NULL); + return hook; +} + +PyObject *K_XCORE::icapsule_set_master(PyObject *self, PyObject *args) +{ + PyObject *ICAPSULE, *MASTER, *CTAG; + if (!PYPARSETUPLE_(args, OOO_, &ICAPSULE, &MASTER, &CTAG)) { RAISE("Bad input."); return NULL; } @@ -19,18 +26,37 @@ PyObject *K_XCORE::icapsule_extract_master(PyObject *self, PyObject *args) return NULL; } + Karray array; + if (Karray_parse_ngon(MASTER, array) != 0) { + RAISE("Failed to parse master NGon."); + return NULL; + } + + E_Float *ctag = NULL; + E_Int size; + E_Int ret = K_NUMPY::getFromNumpyArray(CTAG, ctag, size, true); + if (ret != 1 || size != array.ncells()) { + RAISE("Bad cell tag array."); + return NULL; + } + ICapsule *icap = (ICapsule *)PyCapsule_GetPointer(ICAPSULE, "ICAPSULE"); + icap->M = IMesh(array); - auto Mout = icap->M.export_karray(); + auto &M = icap->M; + M.set_tolerances(icap->NEAR_VERTEX_TOL, icap->NEAR_EDGE_TOL); + M.make_skin(); + M.orient_skin(OUT); + M.ctag.resize(M.nc); + memcpy(M.ctag.data(), ctag, M.nc * sizeof(E_Float)); - return Mout; + return Py_None; } -PyObject *K_XCORE::icapsule_extract_slave(PyObject *self, PyObject *args) +PyObject *K_XCORE::icapsule_set_slaves(PyObject *self, PyObject *args) { - PyObject *ICAPSULE; - E_Int INDEX; - if (!PYPARSETUPLE_(args, O_ I_, &ICAPSULE, &INDEX)) { + PyObject *ICAPSULE, *SLAVES, *PTAGS, *CTAGS; + if (!PYPARSETUPLE_(args, OOOO_, &ICAPSULE, &SLAVES, &PTAGS, &CTAGS)) { RAISE("Bad input."); return NULL; } @@ -42,41 +68,112 @@ PyObject *K_XCORE::icapsule_extract_slave(PyObject *self, PyObject *args) ICapsule *icap = (ICapsule *)PyCapsule_GetPointer(ICAPSULE, "ICAPSULE"); - if (INDEX >= (E_Int)icap->Ss.size()) { - RAISE("Bad slave index."); + if (!PyList_Check(SLAVES)) { + RAISE("Slaves should be a list of NGons."); return NULL; } - auto Sout = icap->Ss[INDEX].export_karray(); + if (!PyList_Check(PTAGS)) { + RAISE("Point tags should be a list of arrays."); + return NULL; + } - return Sout; -} + if (!PyList_Check(CTAGS)) { + RAISE("Cell tags should be a list of arrays."); + return NULL; + } -PyObject *K_XCORE::icapsule_extract_slaves(PyObject *self, PyObject *args) -{ - PyObject *ICAPSULE; + E_Int nslaves = PyList_Size(SLAVES); - if (!PYPARSETUPLE_(args, O_, &ICAPSULE)) { - RAISE("Bad input."); + // Meshes + + std::vector sarrays(nslaves); + E_Int ok = 0; + + for (E_Int i = 0; i < nslaves; i++) { + PyObject *SLAVE = PyList_GetItem(SLAVES, i); + E_Int ret = Karray_parse_ngon(SLAVE, sarrays[i]); + Py_DECREF(SLAVE); + if (ret != 0) { + RAISE("Slaves should be NGons."); + break; + } + ok++; + } + + if (ok != nslaves) { return NULL; } - if (!PyCapsule_IsValid(ICAPSULE, "ICAPSULE")) { - RAISE("Bad ICapsule hook."); + // Point tags + + std::vector ptags(nslaves, NULL); + + if (PyList_Size(PTAGS) != nslaves) { + RAISE("Ptags should be the same size as slaves."); return NULL; } - ICapsule *icap = (ICapsule *)PyCapsule_GetPointer(ICAPSULE, "ICAPSULE"); + for (E_Int i = 0; i < nslaves; i++) { + PyObject *PTAG = PyList_GetItem(PTAGS, i); + E_Int size = -1; + E_Int ret = K_NUMPY::getFromNumpyArray(PTAG, ptags[i], size, true); + //Py_DECREF(PTAG); + if (ret != 1 || size != sarrays[i].npoints()) { + RAISE("Ptags[i] should have size sarrays[i].npoints."); + return NULL; + } + } - PyObject *out = PyList_New(0); + // Cell tags - for (size_t i = 0; i < icap->Ss.size(); i++) { - PyObject *sarray = icap->Ss[i].export_karray(); - PyList_Append(out, sarray); - Py_DECREF(sarray); + std::vector ctags(nslaves, NULL); + + if (PyList_Size(CTAGS) != nslaves) { + RAISE("Cell tags list should be the same size as slaves."); + return NULL; } - return out; + for (E_Int i = 0; i < nslaves; i++) { + PyObject *CTAG = PyList_GetItem(CTAGS, i); + E_Int size = -1; + E_Int ret = K_NUMPY::getFromNumpyArray(CTAG, ctags[i], size, true); + //Py_DECREF(PTAG); + if (ret != 1 || size != sarrays[i].ncells()) { + RAISE("Ctags[i] should have size sarrays[i].npoints."); + return NULL; + } + } + + icap->Ss.reserve(sarrays.size()); + for (size_t i = 0; i < sarrays.size(); i++) { + icap->Ss.push_back(IMesh(sarrays[i])); + auto &S = icap->Ss[i]; + S.set_tolerances(icap->NEAR_VERTEX_TOL, icap->NEAR_EDGE_TOL); + S.make_skin(); + S.orient_skin(IN); + S.ptag.resize(S.np); + memcpy(S.ptag.data(), ptags[i], S.np*sizeof(E_Float)); + // Triangulate faces with tagged points + std::vector fids; + for (E_Int fid = 0; fid < S.nf; fid++) { + const auto &pn = S.F[fid]; + bool triangulate = true; + for (auto p : pn) { + if (S.ptag[p] != 1.0) { + triangulate = false; + break; + } + } + if (triangulate) fids.push_back(fid); + } + S.triangulate(fids); + // Cell tags + S.ctag.resize(S.nc); + memcpy(S.ctag.data(), ctags[i], S.nc*sizeof(E_Float)); + } + + return Py_None; } ICapsule::ICapsule(const Karray &marray, const std::vector &sarrays, @@ -354,13 +451,30 @@ PyObject *K_XCORE::icapsule_intersect(PyObject *self, PyObject *args) // Sf.write_ngon(fname); //} + auto plocs = Mf.locate2(Sf); + + /* + auto bfaces = Mf.extract_covering_faces(Sf, plocs); + Smesh Bf(Mf, bfaces, false); + Bf.write_ngon("Bf.im"); + Bf.make_fcenters(); + Bf.make_fnormals(); + Bf.make_pnormals(); + Bf.make_point_faces(); + Bf.make_BVH(); + plocs = Bf.locate2(Sf); + */ Dcel D = Dcel::intersect(Mf, Sf, plocs); + //Dcel D = Dcel::intersect(Bf, Sf, plocs); //E_Int nf_before_intersect = Sf.nf; D.reconstruct(Mf, Dcel::RED); + //D.reconstruct(Bf, Dcel::RED); + + //Mf.write_ngon("Mf_after_intersect.im"); D.reconstruct(Sf, Dcel::BLACK); @@ -385,6 +499,7 @@ PyObject *K_XCORE::icapsule_intersect(PyObject *self, PyObject *args) puts("Intersection OK\n"); fflush(stdout); } + //exit(0); Mf.reconstruct(M); diff --git a/Cassiopee/XCore/XCore/intersectMesh/icapsule_intersect.cpp b/Cassiopee/XCore/XCore/intersectMesh/icapsule_intersect.cpp new file mode 100644 index 000000000..0c7dc11a1 --- /dev/null +++ b/Cassiopee/XCore/XCore/intersectMesh/icapsule_intersect.cpp @@ -0,0 +1,63 @@ +#include "icapsule.h" +#include "dcel.h" + +PyObject *K_XCORE::icapsule_intersect2(PyObject *self, PyObject *args) +{ + PyObject *ICAPSULE; + + if (!PYPARSETUPLE_(args, O_, &ICAPSULE)) { + RAISE("Bad input"); + return NULL; + } + + if (!PyCapsule_IsValid(ICAPSULE, "ICAPSULE")) { + RAISE("Bad ICapsule hook."); + return NULL; + } + + ICapsule *icap = (ICapsule *)PyCapsule_GetPointer(ICAPSULE, "ICAPSULE"); + + auto &M = icap->M; + auto &Ss = icap->Ss; + + M.make_skin(); + M.orient_skin(OUT); + M.Mf = Smesh(M, M.skin, false); + auto &Mf = M.Mf; + + for (size_t i = 0; i < Ss.size(); i++) { + printf("Intersecting slave %zu ...\n", i); + + Mf.make_fcenters(); + Mf.make_fnormals(); + Mf.make_pnormals(); + Mf.make_point_faces(); + Mf.make_BVH(); + + auto &S = Ss[i]; + + Smesh Sf = Smesh::Smesh_from_tagged_faces(S, true); + Sf.make_fcenters(); + Sf.make_fnormals(); + Sf.make_pnormals(); + + auto plocs = Mf.locate2(Sf); + + Dcel D = Dcel::intersect(Mf, Sf, plocs); + + D.reconstruct(Mf, Dcel::RED); + D.reconstruct(Sf, Dcel::BLACK); + + Sf.reconstruct(S); + + // Tag Sf faces + Sf.tag_faces(S); + + printf("Done.\n"); + fflush(stdout); + } + + Mf.reconstruct(M); + + return Py_None; +} \ No newline at end of file diff --git a/Cassiopee/XCore/XCore/intersectMesh/initSurf.cpp b/Cassiopee/XCore/XCore/intersectMesh/initSurf.cpp new file mode 100644 index 000000000..014bc3a78 --- /dev/null +++ b/Cassiopee/XCore/XCore/intersectMesh/initSurf.cpp @@ -0,0 +1,178 @@ +#include "common/Karray.h" +#include "smesh.h" + +ierr Smesh_from_Karray(Smesh &Sf, const Karray &array) +{ + std::vector owner(array.nfaces(), -1); + std::vector neigh(array.nfaces(), -1); + E_Int nc = array.ncells(); + + for (E_Int cid = 0; cid < nc; cid++) { + E_Int nf = -1; + E_Int *pf = array.get_cell(cid, nf); + for (E_Int i = 0; i < nf; i++) { + E_Int fid = pf[i]-1; + if (owner[fid] == -1) { + owner[fid] = cid; + } else { + if (neigh[fid] != -1) { + IERROR("Face %d is shared between multiple cells.\n", fid); + return IERR_NON_MANIFOLD_FACE; + } + neigh[fid] = cid; + } + } + } + + std::vector skin; + E_Int nf = array.nfaces(); + for (E_Int fid = 0; fid < nf; fid++) { + if (neigh[fid] == -1) + skin.push_back(fid); + } + + Sf.np = 0; + std::map pmap; + Sf.F.resize(skin.size()); + Sf.nf = skin.size(); + + for (size_t i = 0; i < skin.size(); i++) { + E_Int gfid = skin[i]; + E_Int np = -1; + E_Int *pn = array.get_face(gfid, np); + auto &F = Sf.F[i]; + F.resize(np, -1); + for (E_Int j = 0; j < np; j++) { + E_Int gpid = pn[j]; + auto it = pmap.find(gpid); + if (it == pmap.end()) { + pmap[gpid] = Sf.np; + F[j] = Sf.np; + Sf.np++; + } else { + F[j] = it->second; + } + } + } + + Sf.P.resize(Sf.np); + Sf.X.resize(Sf.np); + Sf.Y.resize(Sf.np); + Sf.Z.resize(Sf.np); + for (const auto &pids : pmap) { + E_Int gpid = pids.first - 1; + E_Int lpid = pids.second; + auto &P = Sf.P[lpid]; + P.x() = array.x[gpid]; + P.y() = array.y[gpid]; + P.z() = array.z[gpid]; + Sf.X[lpid] = array.x[gpid]; + Sf.Y[lpid] = array.y[gpid]; + Sf.Z[lpid] = array.z[gpid]; + } + + Sf.Fc = Sf.F; + Sf.make_edge_data(); + + return IERR_OK; +} + +void Smesh::make_edge_data() +{ + F2E.resize(Fc.size()); + std::map edges; + ne = 0; + + for (E_Int fid = 0; fid < nf; fid++) { + auto &pn = Fc[fid]; + auto &pe = F2E[fid]; + pe.resize(pn.size()); + for (size_t i = 0; i < pn.size(); i++) { + E_Int p = pn[i]; + E_Int q = pn[(i+1)%pn.size()]; + o_edge e(p, q); + auto it = edges.find(e); + if (it == edges.end()) { + pe[i] = ne; + edges[e] = ne; + E.push_back(e); + ne++; + } else { + pe[i] = it->second; + } + } + } + + E2F.resize(ne); + + for (E_Int fid = 0; fid < nf; fid++) { + const auto &pe = F2E[fid]; + + for (E_Int eid : pe) { + if (E2F[eid][0] == -1) + E2F[eid][0] = fid; + else + E2F[eid][1] = fid; + } + } +} + +ierr Smesh::check_manifold() const +{ + std::vector count(ne, 0); + for (const auto &pe : F2E) { + for (auto eid : pe) { + count[eid]++; + if (count[eid] > 2) { + IERROR("Edge %d is shared between more than two faces.\n", eid); + return IERR_NON_MANIFOLD_EDGE; + } + } + } + return IERR_OK; +} + +PyObject *K_XCORE::init_surf(PyObject *self, PyObject *args) +{ + PyObject *MESH; + if (!PYPARSETUPLE_(args, O_, &MESH)) { + RAISE("Bad input."); + return NULL; + } + + Karray array; + if (Karray_parse_ngon(MESH, array) != 0) { + RAISE("Input mesh must be an NGon."); + return NULL; + } + + Smesh *Sf = new Smesh; + if (Smesh_from_Karray(*Sf, array) != IERR_OK) { + RAISE("Failed to create surface mesh from input NGon."); + Karray_free_ngon(array); + return NULL; + } + + PyObject *hook = PyCapsule_New((void *)Sf, NULL, NULL); + + return hook; +} + +PyObject *K_XCORE::check_surf_manifold(PyObject *self, PyObject *args) +{ + PyObject *SURF_HOOK; + if (!PYPARSETUPLE_(args, O_, &SURF_HOOK)) { + RAISE("Bad input."); + return NULL; + } + + Smesh *Sf = (Smesh *)PyCapsule_GetPointer(SURF_HOOK, NULL); + + bool is_manifold = Sf->check_manifold() != IERR_NON_MANIFOLD_EDGE; + + if (!is_manifold) { + IERROR("Surface is non-manifold.\n"); + } + + return PyBool_FromLong(is_manifold); +} \ No newline at end of file diff --git a/Cassiopee/XCore/XCore/intersectMesh/mesh.h b/Cassiopee/XCore/XCore/intersectMesh/mesh.h index f03bff7b2..6b93950bf 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/mesh.h +++ b/Cassiopee/XCore/XCore/intersectMesh/mesh.h @@ -99,6 +99,8 @@ struct IMesh { std::vector ptag; std::vector ctag; + Smesh Mf; + E_Float get_min_edge_length() const; void set_tolerances(E_Float near_vertex_tol, E_Float near_edge_tol) diff --git a/Cassiopee/XCore/XCore/intersectMesh/precise.cpp b/Cassiopee/XCore/XCore/intersectMesh/precise.cpp new file mode 100644 index 000000000..7eb930879 --- /dev/null +++ b/Cassiopee/XCore/XCore/intersectMesh/precise.cpp @@ -0,0 +1,127 @@ +/* + Copyright 2013-2024 Onera. + + This file is part of Cassiopee. + + Cassiopee is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Cassiopee is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Cassiopee. If not, see . +*/ + +#include "precise.h" + +static +std::vector merge_and_sort(const Expansion &e, const Expansion &f) +{ + std::vector sequence; + sequence.reserve(e.size() + f.size()); + for (size_t i = 0; i < e.size(); i++) sequence.push_back(e[i]); + for (size_t i = 0; i < f.size(); i++) sequence.push_back(f[i]); + std::sort(sequence.begin(), sequence.end(), + [&] (E_Float a, E_Float b) + { + return fabs(a) < fabs(b); + } + ); + return sequence; +} + +CompensatedSum fast_two_sum(E_Float a, E_Float b) +{ + assert(fabs(a) >= fabs(b)); + E_Float x = a + b; + E_Float bvirtual = x - a; + E_Float y = b - bvirtual; + return {x, y}; +} + +CompensatedSum two_sum(E_Float a, E_Float b) +{ + E_Float x = a + b; + E_Float bvirtual = x - a; + E_Float avirtual = x - bvirtual; + E_Float broundoff = b - bvirtual; + E_Float aroundoff = a - avirtual; + E_Float y = aroundoff + broundoff; + return {x, y}; +} + +CompensatedSum two_diff(E_Float a, E_Float b) +{ + E_Float x = a - b; + E_Float bvirtual = a - x; + E_Float avirtual = x + bvirtual; + E_Float broundoff = bvirtual - b; + E_Float aroundoff = a - avirtual; + E_Float y = aroundoff + broundoff; + return {x, y}; +} + +Expansion linear_sum(const Expansion &e, const Expansion &f) +{ + size_t m = e.size(), n = f.size(); + Expansion h(m + n); + auto g = merge_and_sort(e, f); + CompensatedSum Q = fast_two_sum(g[1], g[0]); + for (int i = 2; i < m + n; i++) { + CompensatedSum R = fast_two_sum(g[i], Q.y); + h[i-2] = R.y; + Q = two_sum(Q.x, R.x); + } + h[m+n-2] = Q.y; + h[m+n-1] = Q.x; + return h; +} + +// Note(Imad): for the 64bit double of C +constexpr E_Float splitter = (1 << 27) + 1; + +CompensatedSum split(E_Float a) +{ + E_Float c = splitter * a; + E_Float abig = c - a; + E_Float ahi = c - abig; + E_Float alo = a - ahi; + return {ahi, alo}; +} + +CompensatedSum two_product(E_Float a, E_Float b) +{ + E_Float x = a * b; + CompensatedSum A = split(a); + CompensatedSum B = split(b); + E_Float err1 = x - (A.x * B.x); + E_Float err2 = err1 - (A.y * B.x); + E_Float err3 = err2 - (A.x * B.y); + E_Float y = (A.y * B.y) - err3; + return {x, y}; +} + +Expansion difference_of_products(E_Float a, E_Float b, E_Float c, E_Float d) +{ + auto S1 = two_product(a, b); + auto S2 = -two_product(c, d); + Expansion E1(S1); + Expansion E2(S2); + Expansion E = linear_sum(E1, E2); + return E; +} + +int Expansion::sign() const +{ + for (E_Int i = m_terms.size()-1; i >= 0; i--) { + if (m_terms[i] == 0.0) continue; + if (m_terms[i] < 0.0) return -1; + return 1; + } + return 0; +} diff --git a/Cassiopee/XCore/XCore/intersectMesh/precise.h b/Cassiopee/XCore/XCore/intersectMesh/precise.h new file mode 100644 index 000000000..ed30df1a6 --- /dev/null +++ b/Cassiopee/XCore/XCore/intersectMesh/precise.h @@ -0,0 +1,81 @@ +/* + Copyright 2013-2024 Onera. + + This file is part of Cassiopee. + + Cassiopee is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Cassiopee is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Cassiopee. If not, see . +*/ + +#pragma once + +#include "xcore.h" +#include + +struct CompensatedSum { + E_Float x, y; + + CompensatedSum &operator-() + { + x = -x; + y = -y; + return *this; + } +}; + +struct Expansion { + std::vector m_terms; + + explicit Expansion(const CompensatedSum &s) + { + m_terms.resize(2); + m_terms[0] = s.y; + m_terms[1] = s.x; + } + + explicit Expansion(size_t m) { m_terms.resize(m); } + + size_t size() const { return m_terms.size(); } + + E_Float operator[](E_Int idx) const + { + assert(idx >= 0); + assert(idx < m_terms.size()); + return m_terms[idx]; + } + + E_Float &operator[](E_Int idx) + { + assert(idx >= 0); + assert(idx < m_terms.size()); + return m_terms[idx]; + } + + void print() const + { + for (auto term : m_terms) { + printf("%g + ", term); + } + printf("\n"); + } + + int sign() const; +}; + +CompensatedSum two_sum(E_Float a, E_Float b); +CompensatedSum fast_two_sum(E_Float a, E_Float b); +CompensatedSum two_diff(E_Float a, E_Float b); +CompensatedSum two_product(E_Float a, E_Float b); +Expansion linear_sum(const Expansion &e, const Expansion &f); +CompensatedSum split(E_Float a); +Expansion difference_of_products(E_Float a, E_Float b, E_Float c, E_Float d); \ No newline at end of file diff --git a/Cassiopee/XCore/XCore/intersectMesh/ray.cpp b/Cassiopee/XCore/XCore/intersectMesh/ray.cpp index ff66ee598..c681dea39 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/ray.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/ray.cpp @@ -118,7 +118,7 @@ bool Ray::intersect_AABB(const AABB &box) const E_Float t1 = (bmin - O) / D; E_Float t2 = (bmax - O) / D; - if (t1 > t2) { E_Float temp = t1; t1 = t2; t2 = temp; } + if (t1 > t2) { std::swap(t1, t2); } tmin = (t1 > tmin) ? t1 : tmin; tmax = (t2 < tmax) ? t2 : tmax; @@ -308,99 +308,3 @@ E_Int MollerTrumbore(E_Float px, E_Float py, E_Float pz, E_Float dx, E_Float dy, return 0; } - -/* -Ray::Ray(Point O, E_Float D[3]) -: org(O), dir(D) -{ - org.x = O.x; - org.y = O.y; - org.z = O.z; - - - E_Float adx = fabs(dir.x); - E_Float ady = fabs(dir.y); - E_Float adz = fabs(dir.z); - - if (adx > ady && adx > adz) kz = 0; - else if (ady > adz) kz = 1; - else kz = 2; - - kx = kz+1; if (kx == 3) kx = 0; - ky = kx+1; if (ky == 3) ky = 0; - - if (dir[kz] < 0.0) std::swap(kx, ky); - - Sx = dir[kx] / dir[kz]; - Sy = dir[ky] / dir[kz]; - Sz = 1.0 / dir[kz]; -} - -E_Int Ray::intersect_triangle(const Point &a, const Point &b, const Point &c, - TriangleIntersection &TI) -{ - // Calculate vertices relative to ray origin - const Point A = a - org; - const Point B = b - org; - const Point C = c - org; - - // Perform shear and scale of vertices - const E_Float Ax = A[kx] - Sx*A[kz]; - const E_Float Ay = A[ky] - Sy*A[kz]; - const E_Float Bx = B[kx] - Sx*B[kz]; - const E_Float By = B[ky] - Sy*B[kz]; - const E_Float Cx = C[kx] - Sx*C[kz]; - const E_Float Cy = C[ky] - Sy*C[kz]; - - // Calculate scale barycentric coordinates - E_Float U = Cx*By - Cy*Bx; - E_Float V = Ax*Cy - Ay*Cx; - E_Float W = Bx*Ay - By*Ax; - - // Fall back to test against edges using quad precision - if (U == 0.0 || V == 0.0 || W == 0.0) { - __float128 CxBy = (__float128)Cx * (__float128)By; - __float128 CyBx = (__float128)Cy * (__float128)Bx; - U = (E_Float)(CxBy - CyBx); - - __float128 AxCy = (__float128)Ax * (__float128)Cy; - __float128 AyCx = (__float128)Ay * (__float128)Cx; - V = (E_Float)(AxCy - AyCx); - - __float128 BxAy = (__float128)Bx * (__float128)Ay; - __float128 ByAx = (__float128)By * (__float128)Ax; - W = (E_Float)(BxAy - ByAx); - } - - // No backface culling - if ((U < 0.0 || V < 0.0 || W < 0.0) && - (U > 0.0 || V > 0.0 || W > 0.0)) - return 0; - - // Calculate determinant - E_Float det = U + V + W; - - if (det == 0.0) return 0; - - // Calculate scaled z-coordinates of vertices - // and use them to calculate the hit distance. - - const E_Float Az = Sz * A[kz]; - const E_Float Bz = Sz * B[kz]; - const E_Float Cz = Sz * C[kz]; - const E_Float T = U*Az + V*Bz + W*Cz; - - if (det < 0.0 && T >= 0) return 0; - else if (det > 0.0 && T <= 0) return 0; - - // Normalize U, V, W and T - const E_Float inv_det = 1.0 / det; - - TI.u = U * inv_det; - TI.v = V * inv_det; - TI.w = W * inv_det; - TI.t = T * inv_det; - - return 1; -} -*/ diff --git a/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp b/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp index 867a09776..ddb275c80 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp @@ -83,9 +83,9 @@ PyObject *handle_slave(const IMesh *M, const Smesh &Mf, Karray& sarray) for (E_Int i = 0; i < ni; i++) { E_Int base = i + ni*j; - bool was_inside = false; + //bool was_inside = true; - for (E_Int k = 0; k < nk; k++) { + for (E_Int k = nk-1; k >= 0; k--) { E_Int idx = base + nij*k; @@ -93,29 +93,28 @@ PyObject *handle_slave(const IMesh *M, const Smesh &Mf, Karray& sarray) E_Float py = Ys[idx]; E_Float pz = Zs[idx]; - if (Mf.is_point_inside(px, py, pz)) { + if (!Mf.is_point_inside(px, py, pz)) { //if (M->is_point_inside(px, py, pz)) { - kmax[base] = k-2; + kmax[base] = k-1; // Cache the point to be projected E_Int proj_id = base + nij*kmax[base]; proj_points.push_back(proj_id); - was_inside = true; + //was_outside = false; break; } } - - //assert(was_outside); - - if (!was_inside) { + /* + if (was_outside) { // i-j line completely outside of M // Projection points is the last point kmax[base] = nk-2; E_Int proj_id = base + nij*kmax[base]; proj_points.push_back(proj_id); } + */ } } @@ -560,6 +559,7 @@ PyObject *handle_slave2(IMesh *M, Karray& sarray, E_Int kmax) #include "smesh.h" #include +//#include "precise.h" PyObject *K_XCORE::removeIntersectingKPlanes(PyObject *self, PyObject *args) { @@ -575,17 +575,18 @@ PyObject *K_XCORE::removeIntersectingKPlanes(PyObject *self, PyObject *args) return NULL; } - IMesh *M = (IMesh *)PyCapsule_GetPointer(MASTER, "IntersectMesh"); - - M->make_skin(); - //M->make_bbox(); - //M->hash_skin(); + /* + E_Float a = 2.654351; + E_Float b = 8.321359; + E_Float c = 11.63157; + E_Float d = 29.68484; + auto E = difference_of_products(a, b, c, d); + E.print(); + printf("%d\n", E.sign()); + exit(0); + */ - Smesh Mf = Smesh::Smesh_from_mesh_skin(*M, M->skin, false); - printf("Mf: %d tris\n", Mf.nf); - Mf.make_fcenters(); - Mf.make_BVH(); - Mf.box.print(); + IMesh *M = (IMesh *)PyCapsule_GetPointer(MASTER, "IntersectMesh"); /* E_Float lo = -1.0; @@ -649,7 +650,7 @@ PyObject *K_XCORE::removeIntersectingKPlanes(PyObject *self, PyObject *args) for (E_Int i = 0; i < nslaves; i++) { printf("Projecting %d / %d\n", i+1, nslaves); //PyObject *st = handle_slave2(M, sarrays[i], kmax); - PyObject *st = handle_slave(M, Mf, sarrays[i]); + PyObject *st = handle_slave(M, M->Mf, sarrays[i]); PyList_Append(slaves_out, st); Py_DECREF(st); Karray_free_structured(sarrays[i]); diff --git a/Cassiopee/XCore/XCore/intersectMesh/smesh_geom.cpp b/Cassiopee/XCore/XCore/intersectMesh/smesh_geom.cpp index 588b2ff13..71f55cd33 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/smesh_geom.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/smesh_geom.cpp @@ -57,15 +57,12 @@ void Smesh::intersect_ray(const Ray &ray, E_Int node_idx, HitData &hit_data) con bool Smesh::is_point_inside(E_Float px, E_Float py, E_Float pz) const { Ray ray(px, py, pz, 1, 0, 0); - if (!box.is_point_inside(ray.o)) return false; - HitData hit_data; - intersect_ray(ray, 0, hit_data); - return hit_data.hits.size() % 2 == 1; + return is_point_inside(ray); } bool Smesh::is_point_inside(const Ray &ray) const { - if (!box.is_point_inside(ray.o)) return false; + //if (!box.is_point_inside(ray.o)) return false; HitData hit_data; intersect_ray(ray, 0, hit_data); return hit_data.hits.size() % 2 == 1; diff --git a/Cassiopee/XCore/XCore/xcore.cpp b/Cassiopee/XCore/XCore/xcore.cpp index 21f15091a..9810a1f3e 100644 --- a/Cassiopee/XCore/XCore/xcore.cpp +++ b/Cassiopee/XCore/XCore/xcore.cpp @@ -57,6 +57,12 @@ static PyMethodDef Pyxcore [] = {"IntersectMesh_Exit", K_XCORE::IntersectMesh_Exit, METH_VARARGS}, {"IntersectMesh_ExtractFaceSet", K_XCORE::IntersectMesh_ExtractFaceSet, METH_VARARGS}, + {"icapsule_init2", K_XCORE::icapsule_init2, METH_VARARGS}, + {"icapsule_set_master", K_XCORE::icapsule_set_master, METH_VARARGS}, + {"icapsule_set_slaves", K_XCORE::icapsule_set_slaves, METH_VARARGS}, + {"icapsule_adapt2", K_XCORE::icapsule_adapt2, METH_VARARGS}, + {"icapsule_intersect2", K_XCORE::icapsule_intersect2, METH_VARARGS}, + {"icapsule_init", K_XCORE::icapsule_init, METH_VARARGS}, {"icapsule_adapt", K_XCORE::icapsule_adapt, METH_VARARGS}, {"icapsule_intersect", K_XCORE::icapsule_intersect, METH_VARARGS}, diff --git a/Cassiopee/XCore/XCore/xcore.h b/Cassiopee/XCore/XCore/xcore.h index 12be22faa..4ce417213 100644 --- a/Cassiopee/XCore/XCore/xcore.h +++ b/Cassiopee/XCore/XCore/xcore.h @@ -69,6 +69,13 @@ namespace K_XCORE PyObject *IntersectMesh_Exit(PyObject *self, PyObject *args); PyObject *IntersectMesh_ExtractFaceSet(PyObject *self, PyObject *args); + PyObject *icapsule_init2(PyObject *self, PyObject *args); + PyObject *icapsule_set_master(PyObject *self, PyObject *args); + PyObject *icapsule_set_slaves(PyObject *self, PyObject *args); + PyObject *icapsule_adapt2(PyObject *self, PyObject *args); + PyObject *icapsule_intersect2(PyObject *self, PyObject *args); + + PyObject *icapsule_init(PyObject *self, PyObject *args); PyObject *icapsule_adapt(PyObject *self, PyObject *args); PyObject *icapsule_intersect(PyObject *self, PyObject *args); diff --git a/Cassiopee/XCore/srcs.py b/Cassiopee/XCore/srcs.py index 5374ecaf0..93617500f 100644 --- a/Cassiopee/XCore/srcs.py +++ b/Cassiopee/XCore/srcs.py @@ -23,11 +23,17 @@ 'XCore/intersectMesh/write.cpp', - 'XCore/intersectMesh/icapsule.cpp', + 'XCore/intersectMesh/icapsule_init.cpp', + 'XCore/intersectMesh/icapsule_adapt.cpp', + 'XCore/intersectMesh/icapsule_intersect.cpp', + 'XCore/intersectMesh/icapsule_extract.cpp', + 'XCore/intersectMesh/icapsule_refine.cpp', 'XCore/intersectMesh/triangulate.cpp', + 'XCore/intersectMesh/precise.cpp', + 'XCore/intersectMesh/mesh_io.cpp', 'XCore/intersectMesh/smesh.cpp', From 9781d4735f900f254e828fdeae52180a950ee523 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Tue, 7 Jan 2025 16:55:10 +0000 Subject: [PATCH 4/5] Apply autopep8 formatting --- Cassiopee/XCore/XCore/PyTree.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Cassiopee/XCore/XCore/PyTree.py b/Cassiopee/XCore/XCore/PyTree.py index 22d4de6b4..c957615cd 100644 --- a/Cassiopee/XCore/XCore/PyTree.py +++ b/Cassiopee/XCore/XCore/PyTree.py @@ -516,7 +516,7 @@ def icapsule_set_slaves(IC, slaves): ctag = I.getNodeFromName(zone, 'keep') if ctag is None: raise ValueError('Missing cell tags.') ctags.append(ctag[1]) - + return xcore.icapsule_set_slaves(IC, sarrs, ptags, ctags) def icapsule_adapt2(IC): From 5092e4c407b15b2931107685b2fc8c15c283456d Mon Sep 17 00:00:00 2001 From: imadhammani Date: Thu, 9 Jan 2025 17:50:36 +0100 Subject: [PATCH 5/5] XCore intersectMesh: split_connex + headers + camelCase --- Cassiopee/XCore/XCore/PyTree.py | 43 ++- .../intersectMesh/check_manifoldness.cpp | 18 ++ .../XCore/intersectMesh/icapsule_adapt.cpp | 18 ++ .../XCore/intersectMesh/icapsule_extract.cpp | 18 ++ .../XCore/intersectMesh/icapsule_init.cpp | 18 ++ .../intersectMesh/icapsule_intersect.cpp | 18 ++ .../XCore/XCore/intersectMesh/initSurf.cpp | 18 ++ .../removeIntersectingKPlanes.cpp | 6 +- Cassiopee/XCore/XCore/intersectMesh/split.cpp | 267 ++++++++++++++++++ Cassiopee/XCore/XCore/xcore.cpp | 2 + Cassiopee/XCore/XCore/xcore.h | 2 + Cassiopee/XCore/srcs.py | 2 + 12 files changed, 416 insertions(+), 14 deletions(-) create mode 100644 Cassiopee/XCore/XCore/intersectMesh/split.cpp diff --git a/Cassiopee/XCore/XCore/PyTree.py b/Cassiopee/XCore/XCore/PyTree.py index 07e2000ec..79c7ef7f2 100644 --- a/Cassiopee/XCore/XCore/PyTree.py +++ b/Cassiopee/XCore/XCore/PyTree.py @@ -484,10 +484,31 @@ def extractFacesFromPointTag(t, tag_name): ############################################################################### -def icapsule_init2(): +def splitConnex(m): + zones = I.getNodesFromType(m, 'Zone_t') + if len(zones) != 1: raise ValueError('Master should be one zone.') + zm = zones[0] + marr = C.getFields(I.__GridCoordinates__, zm, api=3)[0] + ptag = I.getNodeFromName(zm, 'tag') + if ptag is None: raise ValueError('Missing point tags') + ctag = I.getNodeFromName(zm, 'keep') + if ctag is None: raise ValueError('Missing cell tags') + new_arrs, new_ctags, new_ptags = xcore.split_connex(marr, ctag[1], ptag[1]) + zout = [] + for i in range(len(new_arrs)): + z = I.createZoneNode("zone"+str(i), new_arrs[i]) + cont = I.createUniqueChild(z, I.__FlowSolutionCenters__, 'FlowSolution_t') + I._createUniqueChild(cont, 'GridLocation', 'GridLocation_t', value='CellCenter', ) + I.newDataArray('keep', value=new_ctags[i], parent=cont) + cont = I.createUniqueChild(z, I.__FlowSolutionNodes__, 'FlowSolution_t') + I.newDataArray('tag', value=new_ptags[i], parent=cont) + zout.append(z) + return zout + +def icapsuleInit2(): return xcore.icapsule_init2() -def icapsule_set_master(IC, m): +def icapsuleSetMaster(IC, m): zones = I.getNodesFromType(m, 'Zone_t') if len(zones) != 1: raise ValueError('Master should be one zone.') zm = zones[0] @@ -496,7 +517,7 @@ def icapsule_set_master(IC, m): if ctag is None: raise ValueError('Missing cell tags') return xcore.icapsule_set_master(IC, marr, ctag[1]) -def icapsule_set_slaves(IC, slaves): +def icapsuleSetSlaves(IC, slaves): sarrs = [] ptags = [] ctags = [] @@ -519,13 +540,13 @@ def icapsule_set_slaves(IC, slaves): return xcore.icapsule_set_slaves(IC, sarrs, ptags, ctags) -def icapsule_adapt2(IC): +def icapsuleAdapt2(IC): return xcore.icapsule_adapt2(IC) -def icapsule_intersect2(IC): +def icapsuleIntersect2(IC): return xcore.icapsule_intersect2(IC) -def icapsule_extract_master(IC): +def icapsuleExtractMaster(IC): marr, ctag = xcore.icapsule_extract_master(IC) zm = I.createZoneNode("master", marr) @@ -534,7 +555,7 @@ def icapsule_extract_master(IC): I.newDataArray("keep", value=ctag, parent=cont) return zm -def icapsule_extract_slaves(IC): +def icapsuleExtractSlaves(IC): sarrs, ctags = xcore.icapsule_extract_slaves(IC) assert(len(sarrs) == len(ctags)) zones = [] @@ -550,7 +571,7 @@ def icapsule_extract_slaves(IC): ############################################################################### -def icapsule_init(mp, sp): +def icapsuleInit(mp, sp): zm = I.getZones(mp)[0] marr = C.getFields(I.__GridCoordinates__, zm, api=3)[0] @@ -568,7 +589,7 @@ def icapsule_init(mp, sp): return xcore.icapsule_init(marr, sarrs, tags) -def icapsule_adapt(IC): +def icapsuleAdapt(IC): marr, sarrs, stags = xcore.icapsule_adapt(IC) zm = I.createZoneNode("ma", marr) assert(len(sarrs) == len(stags)) @@ -583,7 +604,7 @@ def icapsule_adapt(IC): ts = C.newPyTree(['Base', slave_zones]) return tm, ts -def icapsule_intersect(ma, sa): +def icapsuleIntersect(ma, sa): zm = I.getZones(ma)[0] marr = C.getFields(I.__GridCoordinates__, zm, api=3)[0] @@ -613,7 +634,7 @@ def icapsule_intersect(ma, sa): ts = C.newPyTree(['Base', slave_zones]) return tm, ts -def icapsule_extract_slave(IC, index=0): +def icapsuleExtractSlave(IC, index=0): sarr = xcore.icapsule_extract_slave(IC, index) zs = I.createZoneNode("slave_"+str(index), sarr) return zs diff --git a/Cassiopee/XCore/XCore/intersectMesh/check_manifoldness.cpp b/Cassiopee/XCore/XCore/intersectMesh/check_manifoldness.cpp index ce026455e..b1f53f80a 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/check_manifoldness.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/check_manifoldness.cpp @@ -1,3 +1,21 @@ +/* + Copyright 2013-2024 Onera. + + This file is part of Cassiopee. + + Cassiopee is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Cassiopee is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Cassiopee. If not, see . +*/ #include "mesh.h" #include "common/Karray.h" diff --git a/Cassiopee/XCore/XCore/intersectMesh/icapsule_adapt.cpp b/Cassiopee/XCore/XCore/intersectMesh/icapsule_adapt.cpp index 81d22fcea..98090f55c 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/icapsule_adapt.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/icapsule_adapt.cpp @@ -1,3 +1,21 @@ +/* + Copyright 2013-2024 Onera. + + This file is part of Cassiopee. + + Cassiopee is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Cassiopee is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Cassiopee. If not, see . +*/ #include "icapsule.h" PyObject *K_XCORE::icapsule_adapt2(PyObject *self, PyObject *args) diff --git a/Cassiopee/XCore/XCore/intersectMesh/icapsule_extract.cpp b/Cassiopee/XCore/XCore/intersectMesh/icapsule_extract.cpp index 27e5be9fc..32e9b537b 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/icapsule_extract.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/icapsule_extract.cpp @@ -1,3 +1,21 @@ +/* + Copyright 2013-2024 Onera. + + This file is part of Cassiopee. + + Cassiopee is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Cassiopee is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Cassiopee. If not, see . +*/ #include "icapsule.h" PyObject *K_XCORE::icapsule_extract_master(PyObject *self, PyObject *args) diff --git a/Cassiopee/XCore/XCore/intersectMesh/icapsule_init.cpp b/Cassiopee/XCore/XCore/intersectMesh/icapsule_init.cpp index 96d79ecae..e154c44b9 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/icapsule_init.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/icapsule_init.cpp @@ -1,3 +1,21 @@ +/* + Copyright 2013-2024 Onera. + + This file is part of Cassiopee. + + Cassiopee is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Cassiopee is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Cassiopee. If not, see . +*/ #include #include "icapsule.h" #include "common/Karray.h" diff --git a/Cassiopee/XCore/XCore/intersectMesh/icapsule_intersect.cpp b/Cassiopee/XCore/XCore/intersectMesh/icapsule_intersect.cpp index 0c7dc11a1..e595710af 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/icapsule_intersect.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/icapsule_intersect.cpp @@ -1,3 +1,21 @@ +/* + Copyright 2013-2024 Onera. + + This file is part of Cassiopee. + + Cassiopee is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Cassiopee is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Cassiopee. If not, see . +*/ #include "icapsule.h" #include "dcel.h" diff --git a/Cassiopee/XCore/XCore/intersectMesh/initSurf.cpp b/Cassiopee/XCore/XCore/intersectMesh/initSurf.cpp index 014bc3a78..7922f71e2 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/initSurf.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/initSurf.cpp @@ -1,3 +1,21 @@ +/* + Copyright 2013-2024 Onera. + + This file is part of Cassiopee. + + Cassiopee is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Cassiopee is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Cassiopee. If not, see . +*/ #include "common/Karray.h" #include "smesh.h" diff --git a/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp b/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp index ddb275c80..4e8728aca 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/removeIntersectingKPlanes.cpp @@ -508,9 +508,9 @@ PyObject *handle_slave2(IMesh *M, Karray& sarray, E_Int kmax) char *varString, *eltType; K_ARRAY::getFromArray3(tpl, varString, f, ni, nj, nk, c, eltType); - E_Float* xt = f->begin(1); - E_Float* yt = f->begin(2); - E_Float* zt = f->begin(3); + E_Float *xt = f->begin(1); + E_Float *yt = f->begin(2); + E_Float *zt = f->begin(3); // Copy all the points up to kmax for (E_Int k = 0; k < kmax+1; k++) { diff --git a/Cassiopee/XCore/XCore/intersectMesh/split.cpp b/Cassiopee/XCore/XCore/intersectMesh/split.cpp new file mode 100644 index 000000000..26140a277 --- /dev/null +++ b/Cassiopee/XCore/XCore/intersectMesh/split.cpp @@ -0,0 +1,267 @@ +/* + Copyright 2013-2024 Onera. + + This file is part of Cassiopee. + + Cassiopee is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + Cassiopee is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with Cassiopee. If not, see . +*/ +#include "common/Karray.h" +#include +#include + +PyObject *K_XCORE::split_connex(PyObject *self, PyObject *args) +{ + PyObject *MESH, *CTAG, *PTAG; + if (!PYPARSETUPLE_(args, OOO_, &MESH, &CTAG, &PTAG)) { + RAISE("Bad input."); + return NULL; + } + + Karray array; + if (Karray_parse_ngon(MESH, array) != 0) { + RAISE("Failed to parse NGon."); + return NULL; + } + E_Int nc = array.ncells(); + E_Int np = array.npoints(); + + E_Float *ctag = NULL; + E_Int size = -1; + E_Int ret = K_NUMPY::getFromNumpyArray(CTAG, ctag, size, true); + if (ret != 1 || size != nc) { + RAISE("Bad ctag array."); + return NULL; + } + + E_Float *ptag = NULL; + ret = K_NUMPY::getFromNumpyArray(PTAG, ptag, size, true); + if (ret != 1 || size != np) { + RAISE("Bad ptag array."); + return NULL; + } + + // Cell-cell connectivity + E_Int nf = array.nfaces(); + std::vector owner(nf, -1), neigh(nf, -1); + for (E_Int cid = 0; cid < nc; cid++) { + E_Int stride = -1; + E_Int *pf = array.get_cell(cid, stride); + for (E_Int i = 0; i < stride; i++) { + E_Int fid = pf[i]-1; + if (owner[fid] == -1) owner[fid] = cid; + else neigh[fid] = cid; + } + } + + std::vector> clists; + + E_Int seed = 0; + std::vector visited(nc, 0); + + while (1) { + while (seed < nc && visited[seed]) seed++; + if (seed == nc) break; + + std::stack stk; + clists.push_back(std::vector()); + auto &cell_list = clists.back(); + assert(cell_list.empty()); + stk.push(seed); + + while (!stk.empty()) { + E_Int cid = stk.top(); + stk.pop(); + + if (visited[cid]) continue; + + visited[cid] = 1; + + cell_list.push_back(cid); + + E_Int stride = -1; + E_Int *pf = array.get_cell(cid, stride); + for (E_Int i = 0; i < stride; i++) { + E_Int fid = pf[i]-1; + E_Int nei = owner[fid] == cid ? neigh[fid] : owner[fid]; + if (nei == -1) continue; + if (!visited[nei]) stk.push(nei); + } + } + } + + assert(seed == nc); + printf("Connex parts: %lu\n", clists.size()); + { + E_Int size = 0; + for (const auto &clist : clists) size += clist.size(); + assert(size == nc); + } + + // Construct and export + PyObject *out = PyList_New(0); + + E_Float *X = array.x; + E_Float *Y = array.y; + E_Float *Z = array.z; + + PyObject *alist = PyList_New(0); + PyObject *ctlist = PyList_New(0); + PyObject *ptlist = PyList_New(0); + + for (const auto &clist : clists) { + // Build array from clist + E_Int sizeNGon = 0, sizeNFace = 0; + + std::map fmap; + std::map pmap; + nf = 0; + E_Int np = 0; + nc = clist.size(); + + for (E_Int cid : clist) { + E_Int stride = -1; + E_Int *pf = array.get_cell(cid, stride); + sizeNFace += stride; + + for (E_Int i = 0; i < stride; i++) { + E_Int fid = pf[i]; + + if (fmap.find(fid) == fmap.end()) { + fmap[fid] = nf; + E_Int NP = -1; + E_Int *pn = array.get_face(fid-1, NP); + sizeNGon += NP; + + for (E_Int j = 0; j < NP; j++) { + E_Int pid = pn[j]; + if (pmap.find(pid) == pmap.end()) { + pmap[pid] = np; + np++; + } + } + + nf++; + } + } + } + + const char *varString = "CoordinateX,CoordinateY,CoordinateZ"; + + PyObject *arr = K_ARRAY::buildArray3(3, varString, np, nc, nf, "NGON", + sizeNGon, sizeNFace, 3, false, 3); + + printf("np: %d\n", np); + printf("nf: %d\n", nf); + printf("nc: %d\n", nc); + + K_FLD::FldArrayF *f; + K_FLD::FldArrayI *cn; + K_ARRAY::getFromArray3(arr, f, cn); + + // Points + E_Float *px = f->begin(1); + E_Float *py = f->begin(2); + E_Float *pz = f->begin(3); + + for (const auto &pids : pmap) { + E_Int opid = pids.first-1; + E_Int npid = pids.second; + px[npid] = X[opid]; + py[npid] = Y[opid]; + pz[npid] = Z[opid]; + } + + // Faces + E_Int *indPG = cn->getIndPG(); + E_Int *ngon = cn->getNGon(); + indPG[0] = 0; + + for (const auto &fids : fmap) { + E_Int ofid = fids.first-1; + E_Int nfid = fids.second; + assert(nfid < nf); + E_Int NP = -1; + E_Int *pn = array.get_face(ofid, NP); + indPG[nfid+1] = NP; + } + for (E_Int fid = 0; fid < nf; fid++) indPG[fid+1] += indPG[fid]; + assert(indPG[nf] == sizeNGon); + + for (const auto &fids : fmap) { + E_Int ofid = fids.first-1; + E_Int nfid = fids.second; + E_Int NP = -1; + E_Int *pn = array.get_face(ofid, NP); + E_Int *ptr = &ngon[indPG[nfid]]; + for (E_Int i = 0; i < NP; i++) ptr[i] = pmap.at(pn[i])+1; + } + + // Cells + E_Int *indPH = cn->getIndPH(); + E_Int *nface = cn->getNFace(); + indPH[0] = 0; + nc = 0; + for (E_Int cid : clist) { + E_Int NF = -1; + E_Int *pf = array.get_cell(cid, NF); + indPH[nc+1] = NF; + nc++; + } + for (E_Int cid = 0; cid < nc; cid++) indPH[cid+1] += indPH[cid]; + assert(indPH[nc] == sizeNFace); + + nc = 0; + for (E_Int cid : clist) { + E_Int NF = -1; + E_Int *pf = array.get_cell(cid, NF); + E_Int *ptr = &nface[indPH[nc]]; + for (E_Int i = 0; i < NF; i++) ptr[i] = fmap.at(pf[i])+1 ; + nc++; + } + + npy_intp dims[2]; + dims[1] = 1; + + dims[0] = (npy_intp)nc; + PyArrayObject *ct = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE); + E_Float *ptr = (E_Float *)PyArray_DATA(ct); + for (size_t i = 0; i < nc; i++) ptr[i] = ctag[clist[i]]; + + dims[0] = (npy_intp)np; + PyArrayObject *pt = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE); + ptr = (E_Float *)PyArray_DATA(pt); + for (const auto &pids : pmap) { + E_Int npid = pids.second; + E_Int opid = pids.first-1; + ptr[npid] = ptag[opid]; + } + + + PyList_Append(alist, arr); + Py_DECREF(arr); + PyList_Append(ctlist, (PyObject *)ct); + Py_DECREF(ct); + PyList_Append(ptlist, (PyObject *)pt); + Py_DECREF(pt); + } + + PyList_Append(out, alist); + PyList_Append(out, ctlist); + PyList_Append(out, ptlist); + Py_DECREF(alist); + Py_DECREF(ctlist); + Py_DECREF(ptlist); + + return out; +} \ No newline at end of file diff --git a/Cassiopee/XCore/XCore/xcore.cpp b/Cassiopee/XCore/XCore/xcore.cpp index 9810a1f3e..e0aad1c24 100644 --- a/Cassiopee/XCore/XCore/xcore.cpp +++ b/Cassiopee/XCore/XCore/xcore.cpp @@ -63,6 +63,8 @@ static PyMethodDef Pyxcore [] = {"icapsule_adapt2", K_XCORE::icapsule_adapt2, METH_VARARGS}, {"icapsule_intersect2", K_XCORE::icapsule_intersect2, METH_VARARGS}, + {"split_connex", K_XCORE::split_connex, METH_VARARGS}, + {"icapsule_init", K_XCORE::icapsule_init, METH_VARARGS}, {"icapsule_adapt", K_XCORE::icapsule_adapt, METH_VARARGS}, {"icapsule_intersect", K_XCORE::icapsule_intersect, METH_VARARGS}, diff --git a/Cassiopee/XCore/XCore/xcore.h b/Cassiopee/XCore/XCore/xcore.h index 4ce417213..cb0be28d7 100644 --- a/Cassiopee/XCore/XCore/xcore.h +++ b/Cassiopee/XCore/XCore/xcore.h @@ -75,6 +75,8 @@ namespace K_XCORE PyObject *icapsule_adapt2(PyObject *self, PyObject *args); PyObject *icapsule_intersect2(PyObject *self, PyObject *args); + PyObject *split_connex(PyObject *self, PyObject *args); + PyObject *icapsule_init(PyObject *self, PyObject *args); PyObject *icapsule_adapt(PyObject *self, PyObject *args); diff --git a/Cassiopee/XCore/srcs.py b/Cassiopee/XCore/srcs.py index 0b6126e9d..10164898f 100644 --- a/Cassiopee/XCore/srcs.py +++ b/Cassiopee/XCore/srcs.py @@ -21,6 +21,8 @@ 'XCore/common/common.cpp', 'XCore/common/Karray.cpp', + 'XCore/intersectMesh/split.cpp', + 'XCore/intersectMesh/write.cpp', 'XCore/intersectMesh/icapsule_init.cpp',