From 70dec6e0f2ad3783253b1c94b4a7996e2c2a3ba3 Mon Sep 17 00:00:00 2001 From: imadhammani Date: Thu, 5 Sep 2024 18:54:28 +0200 Subject: [PATCH 1/2] XCore intersectMesh: added a function to extract a cell --- Cassiopee/XCore/XCore/PyTree.py | 12 +++ .../XCore/XCore/intersectMesh/extract.cpp | 94 +++++++++++++++++++ Cassiopee/XCore/XCore/xcore.cpp | 2 + Cassiopee/XCore/XCore/xcore.h | 2 + Cassiopee/XCore/srcs.py | 2 + 5 files changed, 112 insertions(+) create mode 100644 Cassiopee/XCore/XCore/intersectMesh/extract.cpp diff --git a/Cassiopee/XCore/XCore/PyTree.py b/Cassiopee/XCore/XCore/PyTree.py index 427db4836..9eea5c1a7 100644 --- a/Cassiopee/XCore/XCore/PyTree.py +++ b/Cassiopee/XCore/XCore/PyTree.py @@ -354,12 +354,14 @@ def prepareMeshesForIntersection(master, slave): tm = C.newPyTree(["M_adapted", zmo]) ts = C.newPyTree(["S_adapted", zso]) + ''' try: import Intersector.PyTree as XOR except: raise ImportError("XCore.PyTree: requires Intersector.PyTree module.") print("Closing meshes...", flush=True) ts = XOR.closeCells(ts) tm = XOR.closeCells(tm) + ''' return tm, ts @@ -387,3 +389,13 @@ def intersectMesh(master, slave): ts = XOR.closeCells(ts) return tm, ts + +def extractCell(a, cid): + z = I.getZones(a)[0] + + m = C.getFields(I.__GridCoordinates__, z, api=3)[0] + + zmo = I.createZoneNode("cell"+str(cid), xcore.extractCell(m, cid)) + + return zmo + diff --git a/Cassiopee/XCore/XCore/intersectMesh/extract.cpp b/Cassiopee/XCore/XCore/intersectMesh/extract.cpp new file mode 100644 index 000000000..ea70443d9 --- /dev/null +++ b/Cassiopee/XCore/XCore/intersectMesh/extract.cpp @@ -0,0 +1,94 @@ +#include "mesh.h" +#include "karray.h" + +PyObject *K_XCORE::extractCell(PyObject *self, PyObject *args) +{ + PyObject *MESH; + E_Int cid; + + if (!PYPARSETUPLE_(args, O_ I_, &MESH, &cid)) { + RAISE("Bad input."); + return NULL; + } + + Karray marray; + + E_Int ret; + + ret = Karray_parse_ngon(MESH, marray); + + if (ret != 0) return NULL; + + IMesh M(*marray.cn, marray.X, marray.Y, marray.Z, marray.npts); + + IMesh M_out; + + auto &C = M_out.C; + auto &F = M_out.F; + auto &X = M_out.X; + auto &Y = M_out.Y; + auto &Z = M_out.Z; + + E_Int NC = 1; + + E_Int NF = M.C[cid].size(); + + std::vector cell(NF); + + for (E_Int i = 0; i < NF; i++) cell[i] = i; + + C.push_back(cell); + + for (E_Int fid : M.C[cid]) { + F.push_back(M.F[fid]); + } + + std::map new_pids; + + E_Int NP = 0; + + for (auto &pn : F) { + for (E_Int &p : pn) { + auto it = new_pids.find(p); + + if (it == new_pids.end()) { + new_pids[p] = NP; + p = NP; + NP++; + } else { + p = it->second; + } + } + } + + X.resize(NP); + Y.resize(NP); + Z.resize(NP); + + for (const auto &pdat : new_pids) { + E_Int opid = pdat.first; + E_Int npid = pdat.second; + X[npid] = M.X[opid]; + Y[npid] = M.Y[opid]; + Z[npid] = M.Z[opid]; + } + + M_out.nc = NC; + M_out.nf = NF; + M_out.np = NP; + + /* + for (auto &cell : M_out.C) { + for (auto f : cell) printf("%d ", f); + puts(""); + } + + for (auto &face : M_out.F) { + for (auto f : cell) printf("%d ", f); + puts(""); + } + */ + + + return M_out.export_karray(); +} \ No newline at end of file diff --git a/Cassiopee/XCore/XCore/xcore.cpp b/Cassiopee/XCore/XCore/xcore.cpp index 38fa6fa3a..e041e9745 100644 --- a/Cassiopee/XCore/XCore/xcore.cpp +++ b/Cassiopee/XCore/XCore/xcore.cpp @@ -53,6 +53,8 @@ static PyMethodDef Pyxcore [] = {"IntersectMesh_TriangulateFaceSet", K_XCORE::IntersectMesh_TriangulateFaceSet, METH_VARARGS}, {"IntersectMesh_ExtractMesh", K_XCORE::IntersectMesh_ExtractMesh, METH_VARARGS}, + {"extractCell", K_XCORE::extractCell, METH_VARARGS}, + {NULL, NULL} }; diff --git a/Cassiopee/XCore/XCore/xcore.h b/Cassiopee/XCore/XCore/xcore.h index 48eb79811..1ebc97206 100644 --- a/Cassiopee/XCore/XCore/xcore.h +++ b/Cassiopee/XCore/XCore/xcore.h @@ -69,6 +69,8 @@ namespace K_XCORE PyObject *IntersectMesh_Init(PyObject *self, PyObject *args); PyObject *IntersectMesh_ExtractMesh(PyObject *self, PyObject *args); PyObject *IntersectMesh_TriangulateFaceSet(PyObject *self, PyObject *args); + + PyObject *extractCell(PyObject *self, PyObject *args); } #endif diff --git a/Cassiopee/XCore/srcs.py b/Cassiopee/XCore/srcs.py index b35bcf278..2f2f9c65f 100644 --- a/Cassiopee/XCore/srcs.py +++ b/Cassiopee/XCore/srcs.py @@ -20,6 +20,8 @@ 'XCore/common/mem.cpp', 'XCore/common/common.cpp', + 'XCore/intersectMesh/extract.cpp', + 'XCore/intersectMesh/IntersectMesh_Init.cpp', 'XCore/intersectMesh/IntersectMesh_TriangulateFaceSet.cpp', 'XCore/intersectMesh/IntersectMesh_ExtractMesh.cpp', From 12be5826b158185c7c792b4b2a417350115b2598 Mon Sep 17 00:00:00 2001 From: imadhammani Date: Thu, 5 Sep 2024 18:56:20 +0200 Subject: [PATCH 2/2] XCore intersectMesh: better closed edges extraction. OK with EPFL case. --- Cassiopee/XCore/XCore/intersectMesh/dcel.cpp | 25 ++++------- Cassiopee/XCore/XCore/intersectMesh/dcel.h | 3 +- Cassiopee/XCore/XCore/intersectMesh/hedge.h | 2 + .../XCore/intersectMesh/intersectMesh.cpp | 12 +++--- Cassiopee/XCore/XCore/intersectMesh/mesh.cpp | 42 +++++++++++++++++++ Cassiopee/XCore/XCore/intersectMesh/mesh.h | 3 ++ .../prepareMeshesForIntersection.cpp | 7 +--- Cassiopee/XCore/XCore/intersectMesh/smesh.cpp | 41 ++++++++++++------ Cassiopee/XCore/XCore/intersectMesh/smesh.h | 2 +- 9 files changed, 95 insertions(+), 42 deletions(-) diff --git a/Cassiopee/XCore/XCore/intersectMesh/dcel.cpp b/Cassiopee/XCore/XCore/intersectMesh/dcel.cpp index 045dc4706..d777f8980 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/dcel.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/dcel.cpp @@ -358,6 +358,7 @@ void Dcel::init_hedges_and_faces(Smesh &M, E_Int color) Vertex *P = xit->key; Hedge *h = new Hedge(P); + h->eid = i; list[p].push_back(h); @@ -367,6 +368,7 @@ void Dcel::init_hedges_and_faces(Smesh &M, E_Int color) Vertex *V = xit->key; Hedge *t = new Hedge(V); + t->eid = i; list[q].push_back(t); @@ -392,7 +394,7 @@ void Dcel::init_hedges_and_faces(Smesh &M, E_Int color) const E_Float *N = &pnormals[3*pid]; assert(Sign(K_MATH::norm(N, 3)-1) == 0); - sort_leaving_hedges(hedges, N); + sort_leaving_hedges(hedges, N, M); for (size_t i = 0; i < hedges.size(); i++) { Hedge *h = hedges[i]; @@ -419,23 +421,10 @@ void Dcel::init_hedges_and_faces(Smesh &M, E_Int color) Face *f = new Face; f->oid[color] = i; - //f->oid[color] = M.l2gf.at(i); - assert(M.E2F[first_edge][0] == (E_Int)i || M.E2F[first_edge][1] == E_Int(i)); Hedge *REP = (M.E2F[first_edge][0] == (E_Int)i) ? h : t; - /* - if (REP->left != NULL) { - Vertex *v = REP->orig; - point_write("point", v->x, v->y, v->z); - hedge_write("hedge", REP); - hedge_write("twin", REP->twin); - hedge_write("next", REP->next); - hedge_write("prev", REP->prev); - } - */ - assert(REP->left == NULL); f->rep = REP; @@ -1141,7 +1130,9 @@ void Dcel::find_intersections_3D(const Smesh &M, const Smesh &S) } } -void Dcel::sort_leaving_hedges(std::vector &leaving, const E_Float N[3]) const +void Dcel::sort_leaving_hedges(std::vector &leaving, + const E_Float N[3], + const Smesh &M) const { // Choose a vector that is not parallel to N @@ -1221,7 +1212,9 @@ void Dcel::sort_leaving_hedges(std::vector &leaving, const E_Float N[3] else { Hedge *h = leaving[i]; Hedge *w = leaving[j]; + assert(h->color != w->color); + Vertex *P = h->orig; Vertex *Q = h->twin->orig; if (cmp_vtx(P, Q) < 0) return true; @@ -1370,7 +1363,7 @@ void Dcel::resolve_hedges(const Smesh &M, const Smesh &S) E_Float NORM = K_MATH::norm(N, 3); assert(Sign(NORM -1) == 0); - sort_leaving_hedges(leaving, N); + sort_leaving_hedges(leaving, N, M); for (size_t i = 0; i < leaving.size(); i++) { Hedge *h = leaving[i]; diff --git a/Cassiopee/XCore/XCore/intersectMesh/dcel.h b/Cassiopee/XCore/XCore/intersectMesh/dcel.h index ab7b8391c..74a3d066a 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/dcel.h +++ b/Cassiopee/XCore/XCore/intersectMesh/dcel.h @@ -108,5 +108,6 @@ struct Dcel { void trace_hedge(Hedge *sh, const Smesh &M, const Smesh &S, E_Int hid); - void sort_leaving_hedges(std::vector &leaving, const E_Float N[3]) const; + void sort_leaving_hedges(std::vector &leaving, const E_Float N[3], + const Smesh &M) const; }; diff --git a/Cassiopee/XCore/XCore/intersectMesh/hedge.h b/Cassiopee/XCore/XCore/intersectMesh/hedge.h index e862e1f3a..760827f24 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/hedge.h +++ b/Cassiopee/XCore/XCore/intersectMesh/hedge.h @@ -36,6 +36,8 @@ struct Hedge { E_Int color; Cycle *cycle; + E_Int eid; + // Projection of orig E_Float proj_ox = -10000; E_Float proj_oy = -10000; diff --git a/Cassiopee/XCore/XCore/intersectMesh/intersectMesh.cpp b/Cassiopee/XCore/XCore/intersectMesh/intersectMesh.cpp index 35164dce3..433d27c13 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/intersectMesh.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/intersectMesh.cpp @@ -314,6 +314,11 @@ PyObject *K_XCORE::intersectMesh(PyObject *self, PyObject *args) IMesh S(*sarray.cn, sarray.X, sarray.Y, sarray.Z, sarray.npts); M.make_skin(); + S.make_skin(); + + M.orient_skin(OUT); + S.orient_skin(IN); + assert(M.patch.empty()); for (E_Int fid : M.skin) M.patch.insert(fid); @@ -334,15 +339,12 @@ PyObject *K_XCORE::intersectMesh(PyObject *self, PyObject *args) for (E_Int i = 0; i < spatch_size; i++) S.patch.insert(spatch[i]-1); - M.orient_skin(OUT); - S.orient_skin(IN); - // Extract surface meshes Smesh Mf(M); Smesh Sf(S); - //Mf.write_ngon("Mf"); - //Sf.write_ngon("Sf"); + Mf.write_ngon("Mf"); + Sf.write_ngon("Sf"); puts("Making point edges..."); Mf.make_point_edges(); diff --git a/Cassiopee/XCore/XCore/intersectMesh/mesh.cpp b/Cassiopee/XCore/XCore/intersectMesh/mesh.cpp index 90c09b8c7..62ccdfe01 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/mesh.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/mesh.cpp @@ -856,6 +856,37 @@ E_Int IMesh::face_contains_point(E_Int face, E_Float x, E_Float y, E_Float z) co return hit; } +void IMesh::extract_edge_points(E_Int a, E_Int b, std::list &points) +{ + E_Int ref = 0; + + points.clear(); + points.push_back(a); + points.push_back(b); + + do { + ref = 0; + + auto pos = points.begin(); + + assert(*std::prev(points.end()) == b); + + for (auto it = points.begin(); it != std::prev(points.end()); it++) { + E_Int a = *it; + E_Int b = *std::next(it); + + UEdge e(a, b); + + auto search = ecenter.find(e); + + if (search != ecenter.end()) { + points.insert(std::next(it), search->second); + ref = 1; + } + } + } while (ref); +} + IMesh IMesh::extract_conformized() { // Keep all the points @@ -880,6 +911,16 @@ IMesh IMesh::extract_conformized() E_Int p = pn[j]; E_Int q = pn[(j+1)%pn.size()]; + std::list epoints; + + extract_edge_points(p, q, epoints); + + epoints.pop_back(); + + for (auto it = epoints.begin(); it != epoints.end(); it++) + new_face.push_back(*it); + + /* UEdge e(p, q); auto it = ecenter.find(e); @@ -890,6 +931,7 @@ IMesh IMesh::extract_conformized() } else { new_face.push_back(p); } + */ } new_nf++; diff --git a/Cassiopee/XCore/XCore/intersectMesh/mesh.h b/Cassiopee/XCore/XCore/intersectMesh/mesh.h index bc1a38f9c..98ddc976e 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/mesh.h +++ b/Cassiopee/XCore/XCore/intersectMesh/mesh.h @@ -22,6 +22,7 @@ #include #include #include +#include #include "point.h" #include "xcore.h" @@ -112,6 +113,8 @@ struct IMesh { void make_edges(); + void extract_edge_points(E_Int a, E_Int b, std::list &pts); + inline bool face_is_quad(E_Int face) const { return F[face].size() == 4; } inline bool face_is_tri(E_Int face) const { return F[face].size() == 3; } diff --git a/Cassiopee/XCore/XCore/intersectMesh/prepareMeshesForIntersection.cpp b/Cassiopee/XCore/XCore/intersectMesh/prepareMeshesForIntersection.cpp index 1b4b4228d..042c5ac1a 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/prepareMeshesForIntersection.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/prepareMeshesForIntersection.cpp @@ -56,7 +56,7 @@ PyObject *K_XCORE::prepareMeshesForIntersection(PyObject *self, PyObject *args) puts("Preparing meshes for intersection..."); - // Init and orient master/slave meshes + // Init master/slave meshes IMesh M(*marray.cn, marray.X, marray.Y, marray.Z, marray.npts); IMesh S(*sarray.cn, sarray.X, sarray.Y, sarray.Z, sarray.npts); @@ -83,11 +83,6 @@ PyObject *K_XCORE::prepareMeshesForIntersection(PyObject *self, PyObject *args) //for (E_Int i = 0; i < S.nf; i++) { const auto &pn = S.F[fid]; size_t stride = pn.size(); - if (stride > 4) { - std::vector points; - for (auto p : pn) points.push_back(Point(S.X[p], S.Y[p], S.Z[p])); - point_write("face", points); - } assert(stride == 3 || stride == 4); E_Int keep = 1; diff --git a/Cassiopee/XCore/XCore/intersectMesh/smesh.cpp b/Cassiopee/XCore/XCore/intersectMesh/smesh.cpp index 38b52c567..85d7ea10d 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/smesh.cpp +++ b/Cassiopee/XCore/XCore/intersectMesh/smesh.cpp @@ -20,6 +20,7 @@ #include "primitives.h" #include "triangle.h" #include "mesh.h" +#include "io.h" #include #include @@ -51,6 +52,7 @@ bool Smesh::ccw_oriented(E_Int face) Smesh::Smesh(const IMesh &M) { F.resize(M.patch.size()); + //assert(M.patch.size() == M.skin.size()); nf = 0; np = 0; @@ -78,8 +80,11 @@ Smesh::Smesh(const IMesh &M) nf++; } + assert(np == g2lp.size()); + assert(np == l2gp.size()); + // Get the points - np = g2lp.size(); + //np = g2lp.size(); X.resize(np); Y.resize(np); Z.resize(np); @@ -115,41 +120,51 @@ void Smesh::make_edges() F2E.resize(F.size()); std::map edges; + ne = 0; + for (E_Int i = 0; i < nf; i++) { auto &face = F[i]; for (size_t j = 0; j < face.size(); j++) { E_Int p = face[j]; E_Int q = face[(j+1)%face.size()]; - o_edge e(p, q); - auto it = edges.find(e); + o_edge EDGE(p, q); + auto it = edges.find(EDGE); if (it == edges.end()) { - F2E[i].push_back(E.size()); - edges[e] = E.size(); - E.push_back(e); + F2E[i].push_back(ne); + edges[EDGE] = ne; + E.push_back(EDGE); + ne++; } else { F2E[i].push_back(it->second); } } } - ne = E.size(); - - //assert(F.size()+1 + X.size() == E.size() + 2); + assert(ne == E.size()); // Make edge_to_face E2F.resize(ne, {-1,-1}); + std::vector count(ne, 0); + for (E_Int i = 0; i < nf; i++) { - auto &face = F2E[i]; + const auto &face = F2E[i]; + const auto &pn = F[i]; + assert(face.size() == pn.size()); + for (size_t j = 0; j < face.size(); j++) { E_Int e = face[j]; count[e]++; - assert(count[e] <= 2); + if (E2F[e][0] == -1) E2F[e][0] = i; - else E2F[e][1] = i; + else { + assert(E2F[e][1] == -1); + E2F[e][1] = i; + } } } + // Check for (E_Int i = 0; i < ne; i++) { E_Int fi = E2F[i][0]; @@ -739,7 +754,7 @@ Smesh Smesh::extract_conformized() return ret; } -void Smesh::write_faces(const char *fname, const std::vector &faces) +void Smesh::write_faces(const char *fname, const std::vector &faces) const { std::map pmap; diff --git a/Cassiopee/XCore/XCore/intersectMesh/smesh.h b/Cassiopee/XCore/XCore/intersectMesh/smesh.h index b8547e426..8d4a7dbc2 100644 --- a/Cassiopee/XCore/XCore/intersectMesh/smesh.h +++ b/Cassiopee/XCore/XCore/intersectMesh/smesh.h @@ -120,7 +120,7 @@ struct Smesh { //std::vector locate(E_Float x, E_Float y, E_Float z) const; - void write_faces(const char *fname, const std::vector &faces); + void write_faces(const char *fname, const std::vector &faces) const; void write_ngon(const char *fname);