const int * siP(_super_index->begin()), * iP(_index->begin()), *vP(_values->begin());
int cnt = 0;
cI->alloc(_super_index->getNbOfElems(),1); // same number of super packs as number of cells
- std::vector<int> cVec, cIVec;
+ int * cIVecP(cI->getPointer());
MCAuto <DataArrayInt> dsi = _index->deltaShiftIndex();
int sz = dsi->accumulate(0) + dsi->getNbOfElems(); // think about it: one slot for the type, -1 at the end of each face of the cell
- cVec.reserve(sz);
- int * cVecP(cVec.data());
c->alloc(sz, 1);
+ int * cVecP(c->getPointer());
for (int i=0; i < _super_index->getNbOfElems()-1; i++)
{
- cIVec.push_back(cnt);
+ cIVecP[i]= cnt;
int endId = siP[i+1];
cVecP[cnt++] = INTERP_KERNEL::NORM_POLYHED;
for (int j=siP[i]; j < endId; j++)
cVecP[cnt++] = -1;
}
}
- cIVec.push_back(cnt);
- copy(cIVec.begin(), cIVec.end(), cI->getPointer());
- copy(cVec.begin(), cVec.begin()+cnt, c->getPointer());
+ cIVecP[_super_index->getNbOfElems()-1] = cnt;
}
std::size_t MEDCouplingSkyLineArray::getHeapMemorySizeWithoutChildren() const
oss << *i << " ";
oss << std::endl;
oss << " Values:" << std::endl;
- oss << " ";
+ oss << " ";
const int * v = _values->begin();
- int cnt = 0;
- for ( i = _index->begin(); v != _values->end(); ++v, ++cnt )
+ int cnt = 0, cntI = 0;
+ i = _index->begin();
+ for ( const int * si = super_index->begin()+1; v != _values->end(); ++v, ++cnt )
{
if ( cnt == *i )
{
+ if ( cntI == *si && cnt != 0)
+ {
+ oss << std::endl << " ";
+ ++si;
+ }
+
oss << "| ";
- ++i;
+ ++i; ++cntI;
}
oss << *v << " ";
}
std::copy(vP+iP[absolutePackId], vP+iP[absolutePackId+1],pack.begin());
}
+const int * MEDCouplingSkyLineArray::getPackSafePtr(const int absolutePackId, int & packSize) const
+{
+ if(absolutePackId < 0 || absolutePackId >= _index->getNbOfElems())
+ throw INTERP_KERNEL::Exception("MEDCouplingSkyLineArray::getPackSafe: invalid index!");
+ const int * iP(_index->begin()), *vP(_values->begin());
+ packSize = iP[absolutePackId+1]-iP[absolutePackId];
+ return vP+iP[absolutePackId];
+}
+
/**!
* For each given super-pack ID, provide the sub-index of the first matching pack. If no matching pack is found for the
* given super-pack -1 is returned.
* \param[in] superPackIndices the list of super-packs that should be inspected
- * \param[in] pack the pack that the function is looking for in each of the provided super-pack
+ * \param[in] packBg the pack that the function is looking for in each of the provided super-pack
+ * \param[in] packEnd the pack that the function is looking for in each of the provided super-pack
* \param[out] a vector of int, having the same size as superPackIndices and containing for each inspected super-pack
* the index of the first matching pack, or -1 if none found.
*/
-void MEDCouplingSkyLineArray::findPackIds(const std::vector<int> & superPackIndices, const std::vector<int> & pack,
- std::vector<int>& out) const
+void MEDCouplingSkyLineArray::findPackIds(const std::vector<int> & superPackIndices,
+ const int *packBg, const int *packEnd,
+ std::vector<int>& out) const
{
using namespace std;
checkSuperIndex("findPackIds");
+
+ int packSz = std::distance(packBg, packEnd);
+ if (!packSz)
+ throw INTERP_KERNEL::Exception("MEDCouplingSkyLineArray::findPackIds: void pack!");
+
out.resize(superPackIndices.size());
int i = 0;
const int * siP(_super_index->begin()), * iP(_index->begin()), *vP(_values->begin());
const int sPackIdx = *it;
// for each pack
for (int idx=siP[sPackIdx], j=0; idx < siP[sPackIdx+1]; idx++, j++)
- if (equal(&vP[iP[idx]], &vP[iP[idx+1]], pack.begin()))
- {
- out[i] = j;
- break;
- }
+ {
+ if (packSz == (iP[idx+1] - iP[idx]))
+ if (equal(&vP[iP[idx]], &vP[iP[idx+1]], packBg))
+ {
+ out[i] = j;
+ break;
+ }
+ }
}
}
int nt = _index->getNbOfElems();
copy(iP+siP[i]+j+1, iP+nt, iP+siP[i]+j);
_index->reAlloc(nt-1); iP = _index->getPointer(); // better not forget this ...
- for(int ii = siP[i]; ii < nt-1; ii++)
+ for(int ii = siP[i]+j; ii < nt-1; ii++)
iP[ii] -= (end-start);
// _super_index
/**!
* Insert a new pack in super-pack at index i. The pack is inserted at the end of the pack list of the chosen super-pack.
*/
-void MEDCouplingSkyLineArray::pushBackPack(const int i, const std::vector<int> & pack)
+void MEDCouplingSkyLineArray::pushBackPack(const int i, const int * packBg, const int * packEnd)
{
using namespace std;
checkSuperIndex("pushBackPack");
int *siP(_super_index->getPointer()), *iP(_index->getPointer());
- const int sz(pack.size());
+ const int sz(distance(packBg, packEnd));
// _values
_values->reAlloc(_values->getNbOfElems()+sz);
int *vP(_values->getPointer());
copy(vP+iP[siP[i+1]], vPE-sz, vP+iP[siP[i+1]]+sz);
// insert pack
- copy(pack.begin(), pack.end(), vP+iP[siP[i+1]]);
+ copy(packBg, packEnd, vP+iP[siP[i+1]]);
// _index
int nt = _index->getNbOfElems();
_index->reAlloc(nt+1); iP = _index->getPointer();
copy(iP+siP[i+1]+1, iP+nt, iP+siP[i+1]+2);
- iP[siP[i+1]+1] = iP[siP[i+1]] + pack.size();
+ iP[siP[i+1]+1] = iP[siP[i+1]] + sz;
for(int ii = siP[i+1]+2; ii < nt+1; ii++)
iP[ii] += sz;
std::string simpleRepr() const;
void getPackSafe(const int absolutePackId, std::vector<int> & pack) const;
- void findPackIds(const std::vector<int> & superPackIndices, const std::vector<int> & pack, std::vector<int>& out) const;
+ const int * getPackSafePtr(const int absolutePackId, int & packSize) const;
+ void findPackIds(const std::vector<int> & superPackIndices, const int *packBg, const int *packEnd,
+ std::vector<int>& out) const;
void deletePack(const int i, const int j);
- void pushBackPack(const int i, const std::vector<int> & pack);
+ void pushBackPack(const int i, const int * packBg, const int * packEnd);
void convertToPolyhedronConn( MCAuto<DataArrayInt>& c, MCAuto<DataArrayInt>& cI) const;
int *c(getNodalConnectivity()->getPointer()),*cI(getNodalConnectivityIndex()->getPointer());
MCAuto<MEDCouplingSkyLineArray> connSla(MEDCouplingSkyLineArray::BuildFromPolyhedronConn(getNodalConnectivity(), getNodalConnectivityIndex()));
const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
- MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivity(),
- mDesc->getNodalConnectivityIndex()));
+ MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
const double * cooDesc(mDesc->_coords->begin());
MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree());
const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
int nCell=getNumberOfCells(), nDescCell=mDesc->getNumberOfCells();
-
std::vector<double> addCoo;
BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
// Surfaces - handle biggest first
// Keep only candidates whose normal matches the normal of current face
for(vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
{
- bool col = INTERP_KERNEL::isColinear3D(normalsP + (*it)*sIdx, normalsP + *(it)*SPACEDIM, eps);
+ bool col = INTERP_KERNEL::isColinear3D(normalsP + sIdx*SPACEDIM, normalsP + *(it)*SPACEDIM, eps);
if (*it != sIdx && col)
cands2.push_back(*it);
}
cooDesc+SPACEDIM*(cDesc[cIDesc[sIdx]+2]),
cooDesc+SPACEDIM*(cDesc[cIDesc[sIdx]+3]), rotation);
- MCAuto<MEDCouplingUMesh> mPartRef(buildPartOfMySelf(&sIdx, &sIdx, false));
- MCAuto<MEDCouplingUMesh> mPartCand(buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false));
- double * cooPartRef(mPartCand->_coords->getPointer());
+ MCAuto<MEDCouplingUMesh> mPartRef(mDesc->buildPartOfMySelfSlice(sIdx, sIdx+1,1,false)); // false=zipCoords is called
+ MCAuto<MEDCouplingUMesh> mPartCand(mDesc->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), false)); // false=zipCoords is called
+ double * cooPartRef(mPartRef->_coords->getPointer());
double * cooPartCand(mPartCand->_coords->getPointer());
for (int ii = 0; ii < mPartRef->_coords->getNumberOfTuples(); ii++)
rotation.transform_vector(cooPartRef+SPACEDIM*ii);
rotation.transform_vector(cooPartCand+SPACEDIM*ii);
double ze_z = cooPartRef[2];
+ // Localize faces in 2D thanks to barycenters
MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
vector<int> compo; compo.push_back(2);
MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
if (!idsGoodPlane->getNumberOfTuples())
continue;
+ compo[0] = 0; compo.push_back(1);
+ MCAuto<DataArrayDouble> baryPartXY = baryPart->keepSelectedComponents(compo);
+ mPartRef->changeSpaceDimension(2,0.0);
MCAuto<DataArrayInt> cc(DataArrayInt::New()), ccI(DataArrayInt::New());
- mPartRef->getCellsContainingPoints(baryPart->begin(), baryPart->getNumberOfTuples(), eps, cc, ccI);
+ mPartRef->getCellsContainingPoints(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
- // Localize faces in 2D thanks to barycenters
if (!cc->getNumberOfTuples())
continue;
MCAuto<DataArrayInt> dsi = ccI->deltaShiftIndex();
for (const int * ii = ids->begin(); ii != ids->end(); ii++)
hit[cands2[idsGoodPlane->getIJ(*ii,0)]] = true;
- // Current face connectivity
- const int * sIdxConn = cDesc + cIDesc[sIdx] + 1;
-
// For all polyhedrons using this face, replace connectivity:
- vector<int> polyIndices, subPacksIds, facePack;
+ vector<int> polyIndices, packsIds, facePack;
for (int ii=*(revDescIP+sIdx); ii < *(revDescIP+sIdx+1); ii++)
polyIndices.push_back(*(revDescP+ii));
- connSla->findPackIds(polyIndices, facePack, subPacksIds);
+ // Current face connectivity
+ const int * sIdxConn = cDesc + cIDesc[sIdx] + 1;
+ const int * sIdxConnE = cDesc + cIDesc[sIdx+1];
+ connSla->findPackIds(polyIndices, sIdxConn, sIdxConnE, packsIds);
// Deletion of old faces
int jj=0;
for (vector<int>::const_iterator it=polyIndices.begin(); it!=polyIndices.end(); ++it, ++jj)
- connSla->deletePack(*it, subPacksIds[jj]);
+ connSla->deletePack(*it, packsIds[jj]);
// Insertion of new faces:
for (const int * ii = ids->begin(); ii != ids->end(); ii++)
{
// Build pack from the face to insert:
int faceIdx = cands2[idsGoodPlane->getIJ(*ii,0)];
- vector<int> facePack;
- connSlaDesc->getPackSafe(faceIdx, facePack);
+ int facePack2Sz;
+ vector<int> toto;
+ connSlaDesc->getPackSafe(faceIdx, toto);
+ const int * facePack2 = connSlaDesc->getPackSafePtr(faceIdx, facePack2Sz); // contains the type!
// Insert it in all hit polyhedrons:
for (vector<int>::const_iterator it=polyIndices.begin(); it!=polyIndices.end(); ++it)
- connSla->pushBackPack(*it, facePack);
+ connSla->pushBackPack(*it, facePack2+1, facePack2+facePack2Sz); // without the type
}
}
+ // Set back modified connectivity
+ MCAuto<DataArrayInt> cAuto; cAuto.takeRef(_nodal_connec);
+ MCAuto<DataArrayInt> cIAuto; cIAuto.takeRef(_nodal_connec_index);
+ connSla->convertToPolyhedronConn(cAuto, cIAuto);
+
// STEP 2 -- edges
// TODO
ids = sla0.findPackIds([0,1], [9,10,1,12])
self.assertEqual([-1,1], ids)
- sla0.deletePack(1, 0)
+ sla0.deletePack(1, 1)
si, idx, val = sla0.getSuperIndexArray(), sla0.getIndexArray(), sla0.getValuesArray()
- self.assertEqual([28,1,4,9,10,1,12], val.getValues())
- self.assertEqual([0,3,7], idx.getValues())
+ self.assertEqual([28,1,4,2,35,8], val.getValues())
+ self.assertEqual([0,3,6], idx.getValues())
self.assertEqual([0,1,2], si.getValues())
sla0.pushBackPack(0, [3,2,1,0])
si, idx, val = sla0.getSuperIndexArray(), sla0.getIndexArray(), sla0.getValuesArray()
self.assertEqual([0,2,3], si.getValues())
- self.assertEqual([0,3,7,11], idx.getValues())
- self.assertEqual([28,1,4,3,2,1,0, 9,10,1,12], val.getValues())
+ self.assertEqual([0,3,7,10], idx.getValues())
+ self.assertEqual([28,1,4,3,2,1,0, 2,35,8], val.getValues())
# Build connectivity from POLYHED connectivity
cI = [0,16,41]
self.assertEqual(cI1.getValues(), cI)
pass
+ def testMEDCouplingSkyLineArrayThreeLevels2(self):
+ si = [0, 9, 15, 21]
+ siRef = [0, 9, 16, 22]
+ idx = [0,4,8,12,16,20,23,26,29, 32,36,40,44,48,52, 56,60,64,68,72,76,80]
+ c = [1,0,2,3, 5,7,6,4, 1,5,4,0, 0,4,6,2, 2,6,7,3, 3,7,8, 7,5,8, 5,1,8, 1,3,8,
+ 9,1,3,10, 11,12,7,5, 9,11,5,1, 1,5,7,3, 3,7,12,10, 10,12,11,9,
+ 11,5,7,12, 14,16,15,13, 11,14,13,5, 5,13,15,7, 7,15,16,12, 12,16,14,11]
+ idxRef = [0,4,8,12,16,20,23,26,29,32,36,40,44,48,52,55,58, 62, 66, 70, 74, 78, 82 ]
+ cRef = [1,0,2,3, 5,7,6,4, 1,5,4,0, 0,4,6,2, 2,6,7,3, 3,7,8, 7,5,8, 5,1,8, 1,3,8,
+ 9,1,3,10, 11,12,7,5, 9,11,5,1, 3,7,12,10, 10,12,11,9, 3,7,8, 7,5,8,
+ 11,5,7,12, 14,16,15,13, 11,14,13,5, 5,13,15,7, 7,15,16,12, 12,16,14,11]
+ sla0 = MEDCouplingSkyLineArray()
+ sla0.set3( DataArrayInt(si), DataArrayInt(idx), DataArrayInt(c) )
+ ids = sla0.findPackIds([1], [1,5,7,3])
+ sla0.deletePack(1, ids[0])
+ sla0.pushBackPack(1, [3,7,8])
+ sla0.pushBackPack(1, [7,5,8])
+ si, idx, val = sla0.getSuperIndexArray(), sla0.getIndexArray(), sla0.getValuesArray()
+ self.assertEqual(siRef, si.getValues())
+ self.assertEqual(idxRef, idx.getValues())
+ self.assertEqual(cRef, val.getValues())
+ pass
+
+
def testMEDCouplingUMeshgenerateGraph(self):
# cartesian mesh 3x3
arr=DataArrayDouble(4) ; arr.iota()
%newobject MEDCoupling::MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeed;
%newobject MEDCoupling::MEDCouplingUMesh::buildNewNumberingFromCommNodesFrmt;
%newobject MEDCoupling::MEDCouplingUMesh::conformize2D;
+%newobject MEDCoupling::MEDCouplingUMesh::conformize3D;
%newobject MEDCoupling::MEDCouplingUMesh::colinearize2D;
%newobject MEDCoupling::MEDCouplingUMesh::rearrange2ConsecutiveCellTypes;
%newobject MEDCoupling::MEDCouplingUMesh::sortCellsInMEDFileFrmt;
PyObject *findPackIds(PyObject *superPackIndices, PyObject *pack) const throw(INTERP_KERNEL::Exception)
{
std::vector<int> vpack, vspIdx, out;
+
convertPyToNewIntArr3(superPackIndices,vspIdx);
convertPyToNewIntArr3(pack,vpack);
- self->findPackIds(vspIdx,vpack, out);
+
+ self->findPackIds(vspIdx, vpack.data(), vpack.data()+vpack.size(), out);
return convertIntArrToPyList2(out);
}
{
std::vector<int> vpack;
convertPyToNewIntArr3(pack,vpack);
- self->pushBackPack(i,vpack);
+ self->pushBackPack(i,vpack.data(), vpack.data()+vpack.size());
}
PyObject *convertToPolyhedronConn() const throw(INTERP_KERNEL::Exception)
MEDCouplingUMesh *buildSetInstanceFromThis(int spaceDim) const throw(INTERP_KERNEL::Exception);
//tools
DataArrayInt *conformize2D(double eps) throw(INTERP_KERNEL::Exception);
+ DataArrayInt *conformize3D(double eps) throw(INTERP_KERNEL::Exception);
DataArrayInt *colinearize2D(double eps) throw(INTERP_KERNEL::Exception);
void shiftNodeNumbersInConn(int delta) throw(INTERP_KERNEL::Exception);
std::vector<bool> getQuadraticStatus() const throw(INTERP_KERNEL::Exception);
from MEDCouplingDataForTest import MEDCouplingDataForTest
import rlcompleter,readline # this line has to be here, to ensure a usability of MEDCoupling/MEDLoader. B4 removing it please notify to anthony.geay@cea.fr
-class MEDCouplingBasicsTest5(unittest.TestCase):
+class MEDCouplingIntersectTest(unittest.TestCase):
def testSwig2NonRegressionBugIntersectMeshes1(self):
src=MEDCouplingUMesh("src",2)
src.setCoords(DataArrayDouble([-2.5,-3,-2.5,3,2.5,3],3,2))
self.assertTrue(m.getNodalConnectivityIndex().isEqual(DataArrayInt([0,13,20,27])))
pass
+ def testSwig2Conformize3D1(self):
+ """ Simple test where no edge merge is required, only face merging (first part of the algo) """
+ m = MEDCouplingCMesh("cube")
+ dac = DataArrayDouble([0.0, 1.0])
+ m.setCoordsAt(0, dac); m.setCoordsAt(1, dac); m.setCoordsAt(2, dac)
+ m = m.buildUnstructured()
+ m.convertToPolyTypes(range(m.getNumberOfCells()))
+ coo = m.getCoords().getValues()
+ coo.extend([1.0, 0.5, 0.5])
+ m.setCoords(DataArrayDouble(coo, len(coo)/3, 3))
+ #c, cI = m.getNodalConnectivity().getValues(), m.getNodalConnectivityIndex().getValues()
+ ## Initial:
+ #c = [31, 1, 0, 2, 3, -1, 5, 7, 6, 4, -1, 1, 5, 4, 0, -1, 0, 4, 6, 2, -1, 2, 6, 7, 3, -1, 3, 7, 5, 1]
+ c = [31, 1, 0, 2, 3, -1, 5, 7, 6, 4, -1, 1, 5, 4, 0, -1, 0, 4, 6, 2, -1, 2, 6, 7, 3,
+ -1, 3, 7, 8,
+ -1, 7, 5, 8,
+ -1, 5, 1, 8,
+ -1, 1, 3, 8]
+ cI = [0, len(c)]
+ m.setConnectivity(DataArrayInt(c), DataArrayInt(cI))
+ m2 = MEDCouplingCMesh("cubes")
+ m2.setCoordsAt(0, DataArrayDouble([1.0, 2.0]))
+ m2.setCoordsAt(1, dac)
+ m2.setCoordsAt(2, DataArrayDouble([0.0, 1.0, 2.0]))
+ m2 = m2.buildUnstructured(); m2.convertToPolyTypes(range(m2.getNumberOfCells()))
+ mret = MEDCouplingUMesh.MergeUMeshes([m, m2])
+ mret.mergeNodes(1.0e-8)
+ mret.writeVTK("/tmp/toto.vtu")
+ mretDesc, _, _, _, _ = mret.buildDescendingConnectivity()
+ mretDesc.writeVTK("/tmp/toto_desc.vtu")
+
+ mret.conformize3D(1.0e-8)
+
+ mret.writeVTK("/tmp/toto_conf.vtu")
+ mretDesc, _, _, _, _ = mret.buildDescendingConnectivity()
+ mretDesc.writeVTK("/tmp/toto_conf_desc.vtu")
+ return mret
if __name__ == '__main__':
unittest.main()
