Merge pull request #252 from imadhammani/main

XCore intersectMesh: update

Cassiopee/XCore/XCore/PyTree.py

@@ -482,7 +482,96 @@ def extractFacesFromPointTag(t, tag_name):
     arr = C.getFields(I.__GridCoordinates__, z, api=3)[0]
     return xcore.extractFacesFromPointTag(arr, tag[1])
 
-def icapsule_init(mp, sp):
+###############################################################################
+
+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 icapsuleSetMaster(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 icapsuleSetSlaves(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 icapsuleAdapt2(IC):
+    return xcore.icapsule_adapt2(IC)
+
+def icapsuleIntersect2(IC):
+    return xcore.icapsule_intersect2(IC)
+
+def icapsuleExtractMaster(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 icapsuleExtractSlaves(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 icapsuleInit(mp, sp):
     zm = I.getZones(mp)[0]
     marr = C.getFields(I.__GridCoordinates__, zm, api=3)[0]
 
@@ -500,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))
@@ -515,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]
 
@@ -545,23 +634,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):
+def icapsuleExtractSlave(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):

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,

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);

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 <http://www.gnu.org/licenses/>.
+*/
+#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<PyArrayObject *> IMesh::triangulate_skin(const std::vector<Py_BC> &bcs_in,
+    const std::unordered_map<E_Int, E_Int> &fid_to_bc)
+{
+    E_Int NF = nf;
+
+    owner.resize(nf + skin.size(), -1);
+    neigh.resize(nf + skin.size(), -1);
+
+    std::vector<PyArrayObject *> new_bcs(bcs_in.size());
+    std::vector<E_Int *> ptrs(bcs_in.size());
+    std::vector<E_Int> 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);
+}