*/
static void Rotate3DTriangle(const double* PP1,const double*PP2,const double*PP3, TranslationRotationMatrix& rotation_matrix)
{
- //initializes
- rotation_matrix.translate(PP1);
-
+ double P1w[3];
double P2w[3];
double P3w[3];
P2w[0]=PP2[0]; P2w[1]=PP2[1];P2w[2]=PP2[2];
// translating to set P1 at the origin
for (int i=0; i<3; i++)
{
- P2w[i]-=PP1[i];
- P3w[i]-=PP1[i];
+ P1w[i] = -PP1[i];
+ P2w[i] -= PP1[i];
+ P3w[i] -= PP1[i];
}
+ //initializes matrix
+ rotation_matrix.translate(P1w);
+
// rotating to set P2 on the Oxy plane
TranslationRotationMatrix A;
A.rotate_x(P2w);
*/
static void Rotate3DBipoint(const double* PP1,const double*PP2, TranslationRotationMatrix& rotation_matrix)
{
- //initializes
- rotation_matrix.translate(PP1);
-
+ double P1w[3];
double P2w[3];
P2w[0]=PP2[0]; P2w[1]=PP2[1];P2w[2]=PP2[2];
// translating to set P1 at the origin
for (int i=0; i<3; i++)
- P2w[i]-=PP1[i];
+ {
+ P1w[i] = -PP1[i];
+ P2w[i] -= PP1[i];
+ }
+
+ //initializes
+ rotation_matrix.translate(P1w);
// rotating to set P2 on the Oxy plane
TranslationRotationMatrix A;
}
}
+void MEDCouplingSkyLineArray::validSuperIndex(const std::string& func, int superIndex) const
+{
+ if(superIndex < 0 || superIndex >= _super_index->getNbOfElems())
+ {
+ std::ostringstream oss;
+ oss << "MEDCouplingSkyLineArray::" << func << ": invalid super index!";
+ throw INTERP_KERNEL::Exception(oss.str());
+ }
+}
+
+void MEDCouplingSkyLineArray::validIndex(const std::string& func, int idx) const
+{
+ if(idx < 0 || idx >= _index->getNbOfElems())
+ {
+ std::ostringstream oss;
+ oss << "MEDCouplingSkyLineArray::" << func << ": invalid index!";
+ throw INTERP_KERNEL::Exception(oss.str());
+ }
+}
+
+void MEDCouplingSkyLineArray::validSuperIndexAndIndex(const std::string& func, int superIndex, int index) const
+{
+ validSuperIndex(func, superIndex);
+ int idx = _super_index->begin()[superIndex] + index;
+ if(idx < 0 || idx >= _index->getNbOfElems())
+ {
+ std::ostringstream oss;
+ oss << "MEDCouplingSkyLineArray::" << func << ": invalid index!";
+ throw INTERP_KERNEL::Exception(oss.str());
+ }
+}
+
std::string MEDCouplingSkyLineArray::simpleRepr() const
{
std::ostringstream oss;
/**
* For a 2- or 3-level SkyLine array, return a copy of the absolute pack with given identifier.
*/
-void MEDCouplingSkyLineArray::getPackSafe(const int absolutePackId, std::vector<int> & pack) const
+void MEDCouplingSkyLineArray::getSimplePackSafe(const int absolutePackId, std::vector<int> & pack) const
{
if(absolutePackId < 0 || absolutePackId >= _index->getNbOfElems())
throw INTERP_KERNEL::Exception("MEDCouplingSkyLineArray::getPackSafe: invalid index!");
/**
* Same as getPackSafe, but directly returns a pointer to the internal data with the size of the pack.
*/
-const int * MEDCouplingSkyLineArray::getPackSafePtr(const int absolutePackId, int & packSize) const
+const int * MEDCouplingSkyLineArray::getSimplePackSafePtr(const int absolutePackId, int & packSize) const
{
if(absolutePackId < 0 || absolutePackId >= _index->getNbOfElems())
throw INTERP_KERNEL::Exception("MEDCouplingSkyLineArray::getPackSafe: invalid index!");
using namespace std;
checkSuperIndex("deletePack");
+ validSuperIndexAndIndex("deletePack", superIdx, idx);
int * vP = _values->getPointer();
int * siP(_super_index->getPointer()), *iP(_index->getPointer());
using namespace std;
checkSuperIndex("pushBackPack");
+ validSuperIndex("pushBackPack", superIdx);
int *siP(_super_index->getPointer()), *iP(_index->getPointer());
const int sz(distance(packBg, packEnd));
(siP[ii])++;
}
+/**
+ * Replace pack with absolute index 'idx' with the provided new pack. Function can be used either
+ * for 2-level SkyLine or 3-level SkyLine.
+ */
+void MEDCouplingSkyLineArray::replaceSimplePack(const int idx, const int * packBg, const int * packEnd)
+{
+ using namespace std;
+
+ validIndex("replaceSimplePack", idx);
+
+ int * iP(_index->getPointer());
+ int newSz = distance(packBg, packEnd);
+ const int start = iP[idx], end = iP[idx+1];
+
+ // _values
+ int initValSz = _values->getNbOfElems();
+ int deltaSz = newSz-(end-start); // can be negative
+ if (deltaSz)
+ {
+ if (deltaSz > 0)
+ _values->reAlloc(initValSz+deltaSz);
+ int *vP(_values->getPointer());
+ copy(vP+end, vP+initValSz, vP+end+deltaSz);
+ if (deltaSz < 0)
+ _values->reAlloc(initValSz+deltaSz);
+ }
+
+ // copy new pack
+ copy(packBg, packEnd, _values->getPointer()+start);
+
+ // _index
+ for(int ii = idx+1; ii < _index->getNbOfElems(); ii++)
+ iP[ii] += deltaSz;
+}
/**
* Replace pack with super index 'superIdx' and index 'idx' with the provided new pack.
+ * Function can be used only for 3-level SkyLine.
*/
void MEDCouplingSkyLineArray::replacePack(const int superIdx, const int idx, const int * packBg, const int * packEnd)
{
using namespace std;
checkSuperIndex("replacePack");
+ validSuperIndexAndIndex("replacePack", superIdx, idx);
- int * vP = _values->getPointer();
int * siP(_super_index->getPointer()), *iP(_index->getPointer());
int newSz = distance(packBg, packEnd);
const int start = iP[siP[superIdx]+idx], end = iP[siP[superIdx]+idx+1];
- // _values
- copy(vP+end, vP+_values->getNbOfElems(), vP+start);
- _values->reAlloc(_values->getNbOfElems() - (end-start));
- // _index
- int nt = _index->getNbOfElems();
- copy(iP+siP[superIdx]+idx+1, iP+nt, iP+siP[superIdx]+idx);
- _index->reAlloc(nt-1); iP = _index->getPointer(); // better not forget this ...
- for(int ii = siP[superIdx]+idx; ii < nt-1; ii++)
- iP[ii] -= (end-start);
+ // _values
+ int initValSz = _values->getNbOfElems();
+ int deltaSz = newSz-(end-start); // can be negative
+ if (deltaSz)
+ {
+ if (deltaSz > 0)
+ _values->reAlloc(initValSz+deltaSz);
+ int *vP(_values->getPointer());
+ copy(vP+end, vP+initValSz, vP+end+deltaSz);
+ if (deltaSz < 0)
+ _values->reAlloc(initValSz+deltaSz);
+ }
- // _super_index
- for(int ii = superIdx+1; ii < _super_index->getNbOfElems(); ii++)
- (siP[ii])--;
+ // copy new pack
+ copy(packBg, packEnd, _values->getPointer()+start);
+ // _index
+ for(int ii = siP[superIdx]+idx+1; ii < _index->getNbOfElems(); ii++)
+ iP[ii] += deltaSz;
}
std::string simpleRepr() const;
- void getPackSafe(const int absolutePackId, std::vector<int> & pack) const;
- const int * getPackSafePtr(const int absolutePackId, int & packSize) const;
+ void getSimplePackSafe(const int absolutePackId, std::vector<int> & pack) const;
+ const int * getSimplePackSafePtr(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 superIdx, const int idx);
void pushBackPack(const int superIdx, const int * packBg, const int * packEnd);
+ void replaceSimplePack(const int idx, const int * packBg, const int * packEnd);
void replacePack(const int superIdx, const int idx, const int * packBg, const int * packEnd);
void convertToPolyhedronConn( MCAuto<DataArrayInt>& c, MCAuto<DataArrayInt>& cI) const;
~MEDCouplingSkyLineArray();
void checkSuperIndex(const std::string& func) const;
+ void validSuperIndex(const std::string& func, int superIndex) const;
+ void validIndex(const std::string& func, int index) const;
+ void validSuperIndexAndIndex(const std::string& func, int superIndex, int index) const;
MCAuto<DataArrayInt> _super_index;
MCAuto<DataArrayInt> _index;
MCAuto<DataArrayInt> _values;
};
+
}
# endif
int split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI);
static bool Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords);
static void ComputeAllTypesInternal(std::set<INTERP_KERNEL::NormalizedCellType>& types, const DataArrayInt *nodalConnec, const DataArrayInt *nodalConnecIndex);
+ static bool OrderPointsAlongLine(const double * coo, int startNode, int endNode,
+ const int * c, const int * cI, const int *idsBg, const int *endBg,
+ std::vector<int> & pointIds);
+ static void ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode,
+ const std::vector<int>& insidePoints, std::vector<int>& modifiedFace);
public:
MEDCOUPLING_EXPORT static DataArrayInt *ComputeRangesFromTypeDistribution(const std::vector<int>& code);
MEDCOUPLING_EXPORT static const int N_MEDMEM_ORDER=24;
#include <cstring>
#include <limits>
#include <list>
+#include <assert.h>
using namespace MEDCoupling;
///@cond INTERNAL
/**
- *
+ * c, cI describe a wire mesh in 3D space.
+ *\return true if the segment is indeed split
*/
-bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, const double * startNode, const double *endNode,
+bool MEDCouplingUMesh::OrderPointsAlongLine(const double * coo, int startNode, int endNode,
const int * c, const int * cI, const int *idsBg, const int *endBg,
std::vector<int> & pointIds)
{
+ using namespace std;
+
+ const int SPACEDIM=3;
+ typedef pair<double, int> PairDI;
+ set< PairDI > x;
+ for (const int * it = idsBg; it != endBg; ++it)
+ {
+ assert(c[cI[*it]] == INTERP_KERNEL::NORM_SEG2);
+ int start = c[cI[*it]+1], end = c[cI[*it]+2];
+ x.insert(make_pair(coo[start*SPACEDIM], start)); // take only X coords
+ x.insert(make_pair(coo[end*SPACEDIM], end));
+ }
+ vector<PairDI> xx(x.begin(), x.end());
+ sort(xx.begin(),xx.end());
+ pointIds.reserve(xx.size());
+
+ // Keep what is inside [startNode, endNode]:
+ int go = 0;
+ for (vector<PairDI>::const_iterator it=xx.begin(); it != xx.end(); ++it)
+ {
+ const int idx = (*it).second;
+ if (!go)
+ {
+ if (idx == startNode)
+ go = 1;
+ if (idx == endNode)
+ go = 2;
+ continue;
+ }
+ if (idx == endNode || idx == startNode)
+ break;
+ pointIds.push_back(idx);
+ }
+ if (go == 2)
+ reverse(pointIds.begin(), pointIds.end());
+ return (pointIds.size() != 0);
+}
+
+void MEDCouplingUMesh::ReplaceEdgeInFace(const int * sIdxConn, const int * sIdxConnE, int startNode, int endNode,
+ const std::vector<int>& insidePoints, std::vector<int>& modifiedFace)
+{
+ using namespace std;
+ int dst = distance(sIdxConn, sIdxConnE);
+ modifiedFace.reserve(dst + insidePoints.size());
+ modifiedFace.resize(dst);
+ copy(sIdxConn, sIdxConnE, modifiedFace.data());
+
+ vector<int>::iterator shortEnd = modifiedFace.begin()+dst;
+ vector<int>::iterator startPos = find(modifiedFace.begin(), shortEnd , startNode);
+ if (startPos == shortEnd)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
+ vector<int>::iterator endPos = find(modifiedFace.begin(),shortEnd, endNode);
+ if (endPos == shortEnd)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ReplaceEdgeInFace: internal error, should never happen!");
+ int d = distance(startPos, endPos);
+ if (d > 0 || d == (dst-1))
+ modifiedFace.insert(++startPos, insidePoints.begin(), insidePoints.end());
+ else
+ modifiedFace.insert(++endPos, insidePoints.rbegin(), insidePoints.rend());
}
+
///@endcond
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for meshes with spaceDim=3 and meshDim=3!");
if(_types.size() != 1 || *(_types.begin()) != INTERP_KERNEL::NORM_POLYHED)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D : This method only works for polyhedrons! Call convertAllToPoly first.");
- MCAuto<DataArrayInt> desc(DataArrayInt::New()),descI(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
- MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc,descI,revDesc,revDescI));
- const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
+
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->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
const double * coo(_coords->begin());
+ MCAuto<DataArrayInt> ret(DataArrayInt::New());
+ std::cout << "BEGIN\n";
+ std::cout << connSla->simpleRepr();
- /*************************
- * STEP 1 -- faces
- *************************/
- // Build BBTree
- MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree());
- const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
- int nCell=getNumberOfCells(), nDescCell=mDesc->getNumberOfCells();
- BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
- // Surfaces - handle biggest first
- MCAuto<MEDCouplingFieldDouble> surfF = mDesc->getMeasureField(true);
- DataArrayDouble * surfs = surfF->getArray();
- // Normal field
- MCAuto<MEDCouplingFieldDouble> normalsF = mDesc->buildOrthogonalField();
- DataArrayDouble * normals = normalsF->getArray();
- const double * normalsP = normals->getConstPointer();
-
- // Sort faces by decreasing surface:
- vector<pair<double,int>> S;
- for(int i=0;i < surfs->getNumberOfTuples();i++){
- pair<double,int> p = make_pair(surfs->getIJ(i, 0), i);
- S.push_back(p);
- }
- sort(S.rbegin(),S.rend()); // reverse sort
- vector<bool> hit(nDescCell);
- fill(hit.begin(), hit.end(), false);
- vector<int> hitPoly; // the final result: which 3D cells have been modified.
+ {
+ /*************************
+ * STEP 1 -- faces
+ *************************/
+ MCAuto<DataArrayInt> descDNU(DataArrayInt::New()),descIDNU(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
+ MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(descDNU,descIDNU,revDesc,revDescI));
+ const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
+ const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
+ MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
+
+ // Build BBTree
+ MCAuto<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree());
+ const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
+ int nCell=getNumberOfCells(), nDescCell=mDesc->getNumberOfCells();
+ BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
+ // Surfaces - handle biggest first
+ MCAuto<MEDCouplingFieldDouble> surfF = mDesc->getMeasureField(true);
+ DataArrayDouble * surfs = surfF->getArray();
+ // Normal field
+ MCAuto<MEDCouplingFieldDouble> normalsF = mDesc->buildOrthogonalField();
+ DataArrayDouble * normals = normalsF->getArray();
+ const double * normalsP = normals->getConstPointer();
+
+ // Sort faces by decreasing surface:
+ vector<pair<double,int>> S;
+ for(int i=0;i < surfs->getNumberOfTuples();i++){
+ pair<double,int> p = make_pair(surfs->begin()[i], i);
+ S.push_back(p);
+ }
+ sort(S.rbegin(),S.rend()); // reverse sort
+ vector<bool> hit(nDescCell);
+ fill(hit.begin(), hit.end(), false);
+ vector<int> hitPoly; // the final result: which 3D cells have been modified.
- for( vector<pair<double,int>>::const_iterator it = S.begin(); it != S.end(); it++)
- {
- int sIdx = (*it).second;
- if (hit[sIdx]) continue;
+ for( vector<pair<double,int>>::const_iterator it = S.begin(); it != S.end(); it++)
+ {
+ int sIdx = (*it).second;
+ if (hit[sIdx]) continue;
+
+ vector<int> candidates, cands2;
+ myTree.getIntersectingElems(bbox+sIdx*2*SPACEDIM,candidates);
+ // 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 + sIdx*SPACEDIM, normalsP + *(it)*SPACEDIM, eps);
+ if (*it != sIdx && col)
+ cands2.push_back(*it);
+ }
+ if (!cands2.size())
+ continue;
+
+ // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later
+ INTERP_KERNEL::TranslationRotationMatrix rotation;
+ INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[sIdx]+1]),
+ coo+SPACEDIM*(cDesc[cIDesc[sIdx]+2]),
+ coo+SPACEDIM*(cDesc[cIDesc[sIdx]+3]), rotation);
+
+ 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);
+ for (int ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); 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);
+ MCAuto<DataArrayInt> idsGoodPlane = baryPartZ->findIdsInRange(ze_z-eps, ze_z+eps);
+ 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(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
+
+ if (!cc->getNumberOfTuples())
+ continue;
+ MCAuto<DataArrayInt> dsi = ccI->deltaShiftIndex();
- vector<int> candidates, cands2;
- myTree.getIntersectingElems(bbox+sIdx*2*SPACEDIM,candidates);
- // 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 + sIdx*SPACEDIM, normalsP + *(it)*SPACEDIM, eps);
- if (*it != sIdx && col)
- cands2.push_back(*it);
+ MCAuto<DataArrayInt> tmp = dsi->findIdsInRange(0, 2);
+ if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled.");
}
- if (!cands2.size())
- continue;
-
- // Now rotate, and match barycenters -- this is where we will bring Intersect2DMeshes later
- INTERP_KERNEL::TranslationRotationMatrix rotation;
- INTERP_KERNEL::TranslationRotationMatrix::Rotate3DTriangle(coo+SPACEDIM*(cDesc[cIDesc[sIdx]+1]),
- coo+SPACEDIM*(cDesc[cIDesc[sIdx]+2]),
- coo+SPACEDIM*(cDesc[cIDesc[sIdx]+3]), rotation);
-
- 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);
- for (int ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); 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);
- MCAuto<DataArrayInt> idsGoodPlane = baryPartZ->findIdsInRange(ze_z-eps, ze_z+eps);
- 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(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
-
- if (!cc->getNumberOfTuples())
- continue;
- MCAuto<DataArrayInt> dsi = ccI->deltaShiftIndex();
- {
- MCAuto<DataArrayInt> tmp = dsi->findIdsInRange(0, 2);
- if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
- throw INTERP_KERNEL::Exception("Non expected non-conformity. Only simple (=partition-like) non-conformities are handled.");
+ MCAuto<DataArrayInt> ids = dsi->findIdsEqual(1);
+ // Boundary face:
+ if (!ids->getNumberOfTuples())
+ continue;
+
+ double checkSurf=0.0;
+ for (const int * ii = ids->begin(); ii != ids->end(); ii++)
+ {
+ int faceIdx = cands2[idsGoodPlane->begin()[*ii]];
+ hit[faceIdx] = true;
+ checkSurf += surfs->begin()[faceIdx];
+ }
+ if (fabs(checkSurf - surfs->begin()[sIdx]) > eps)
+ {
+ ostringstream oss;
+ oss << "MEDCouplingUMesh::conformize3D: Non expected non-conformity. Only simple (=partition-like) non-conformities are handled. Face #" << sIdx << " violates this condition!";
+ throw INTERP_KERNEL::Exception(oss.str());
+ }
+
+ // For all polyhedrons using this face, replace connectivity:
+ vector<int> polyIndices, packsIds, facePack;
+ for (int ii=revDescIP[sIdx]; ii < revDescIP[sIdx+1]; ii++)
+ polyIndices.push_back(revDescP[ii]);
+ ret->pushBackValsSilent(polyIndices.data(),polyIndices.data()+polyIndices.size());
+
+ // 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)
+ {
+ if (packsIds[jj] == -1)
+ // The below should never happen - if a face is used several times, with a different layout of the nodes
+ // it means that is is already conform, so it is not hit by the algorithm. The algorithm only hits
+ // faces which are actually used only once, by a single cell. This is different for edges below.
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize3D: Could not find face in connectivity! Internal error.");
+ else
+ 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)];
+ int facePack2Sz;
+// vector<int> toto;
+// connSlaDesc->getSimplePackSafe(faceIdx, toto);
+ const int * facePack2 = connSlaDesc->getSimplePackSafePtr(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, facePack2+1, facePack2+facePack2Sz); // without the type
+ }
}
- // TODO check matching total surface!
-
- MCAuto<DataArrayInt> ids = dsi->findIdsEqual(1);
- // Boundary face:
- if (!ids->getNumberOfTuples())
- continue;
-
- for (const int * ii = ids->begin(); ii != ids->end(); ii++)
- hit[cands2[idsGoodPlane->getIJ(*ii,0)]] = true;
-
- // For all polyhedrons using this face, replace connectivity:
- vector<int> polyIndices, packsIds, facePack;
- for (int ii=*(revDescIP+sIdx); ii < *(revDescIP+sIdx+1); ii++)
- polyIndices.push_back(*(revDescP+ii));
- // 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, 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)];
- 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, 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
- /************************/
- // Now we have a face-conform mesh.
- mDesc = buildDescendingConnectivity(desc,descI,revDesc,revDescI);
- MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
- MCAuto<MEDCouplingUMesh> mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);
- const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
- bboxArr = mDesc2->getBoundingBoxForBBTree();
- const double *bbox2(bboxArr->begin());
- int nDesc2Cell=mDesc2->getNumberOfCells();
- BBTree<SPACEDIM,int> myTree(bbox2,0,0,nDesc2Cell,-eps);
-
- // Edges - handle longest first
- MCAuto<MEDCouplingFieldDouble> lenF = mDesc2->getMeasureField(true);
- DataArrayDouble * lens = lenF->getArray();
-
- // Sort edges by decreasing length:
- vector<pair<double,int>> S2;
- for(int i=0;i < lens->getNumberOfTuples();i++){
- pair<double,int> p = make_pair(lens->getIJ(i, 0), i);
- S2.push_back(p);
- }
- sort(S.rbegin(),S.rend()); // reverse sort
- hit.resize(nDesc2Cell);
- fill(hit.begin(), hit.end(), false);
+ std::cout << "BEFORE STEP2\n";
+ std::cout << connSla->simpleRepr();
- for( vector<pair<double,int>>::const_iterator it = S.begin(); it != S.end(); it++)
- {
- int eIdx = (*it).second;
- if (hit[eIdx]) continue;
-
- vector<int> candidates, cands2;
- myTree.getIntersectingElems(bbox+eIdx*2*SPACEDIM,candidates);
- // Keep only candidates colinear with current edge
- double vCurr[3];
- unsigned start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2];
- for (int i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar?
- vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3];
- for(vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
- {
- double vOther[3];
- unsigned start2 = cDesc2[cIDesc2[*it]+1], end2 = cDesc2[cIDesc2[*it]+2];
- for (int i3=0; i3 < 3; i3++)
- vCurr[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3];
- bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps);
- if (*it != eIdx && col)
- cands2.push_back(*it);
- }
- if (!cands2.size())
- continue;
-
- // Now rotate edges to bring them on Ox
- const double * startNode = coo+SPACEDIM*(cDesc[cIDesc[eIdx]+1]);
- const double * endNode = coo+SPACEDIM*(cDesc[cIDesc[eIdx]+2]);
- INTERP_KERNEL::TranslationRotationMatrix rotation;
- INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(startNode, endNode, rotation);
- MCAuto<MEDCouplingUMesh> mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false)); // false=zipCoords is called
- MCAuto<MEDCouplingUMesh> mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called
- MCAuto<DataArrayInt> nodeMap(mPartCand->zipCoordsTraducer());
- 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);
- for (int ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
- rotation.transform_vector(cooPartCand+SPACEDIM*ii);
- double ze_y = cooPartRef[1], ze_z = cooPartRef[2];
-
- // Eliminate all edges for which y or z is not null
- MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
- vector<int> compo; compo.push_back(1);
- MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
- compo[0] = 2;
- MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
- MCAuto<DataArrayInt> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
- MCAuto<DataArrayInt> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
- MCAuto<DataArrayInt> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
- if (!idsGoodLine->getNumberOfTuples())
- continue;
-
- // Now the ordering along the Ox axis:
- std::vector<int> insidePoints;
- bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), startNode, endNode,
- mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
- idsGoodLine->begin(), idsGoodLine->end(),
- /*out*/insidePoints);
- if (!isSplit) // current segment remains in one piece
- 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(baryPartXY->begin(), baryPartXY->getNumberOfTuples(), eps, cc, ccI);
-
- if (!cc->getNumberOfTuples())
- continue;
- MCAuto<DataArrayInt> dsi = ccI->deltaShiftIndex();
+ {
+ /************************
+ * STEP 2 -- edges
+ /************************/
+ // Now we have a face-conform mesh.
+
+ // Recompute descending
+ MCAuto<DataArrayInt> descDNU(DataArrayInt::New()),descIDNU(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescI(DataArrayInt::New());
+ // Rebuild desc connectivity with orientation this time!!
+ MCAuto<MEDCouplingUMesh> mDesc(buildDescendingConnectivity2(descDNU,descIDNU,revDesc,revDescI));
+ const int *revDescIP(revDescI->getConstPointer()), *revDescP(revDesc->getConstPointer());
+ const int *cDesc(mDesc->getNodalConnectivity()->begin()),*cIDesc(mDesc->getNodalConnectivityIndex()->begin());
+ MCAuto<MEDCouplingSkyLineArray> connSlaDesc(MEDCouplingSkyLineArray::New(mDesc->getNodalConnectivityIndex(), mDesc->getNodalConnectivity()));
+ std::cout << "writing VTK!\n";
+ mDesc->writeVTK("/tmp/toto_desc_confInter.vtu", true);
+ MCAuto<DataArrayInt> desc2(DataArrayInt::New()),descI2(DataArrayInt::New()),revDesc2(DataArrayInt::New()),revDescI2(DataArrayInt::New());
+ MCAuto<MEDCouplingUMesh> mDesc2 = mDesc->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);
+ std::cout << "writing VTK2!\n";
+ mDesc2->writeVTK("/tmp/toto_desc2_confInter.vtu", true);
+ const int *revDescIP2(revDescI2->getConstPointer()), *revDescP2(revDesc2->getConstPointer());
+ const int *cDesc2(mDesc2->getNodalConnectivity()->begin()),*cIDesc2(mDesc2->getNodalConnectivityIndex()->begin());
+ MCAuto<DataArrayDouble> bboxArr(mDesc2->getBoundingBoxForBBTree());
+ const double *bbox2(bboxArr->begin());
+ int nDesc2Cell=mDesc2->getNumberOfCells();
+ BBTree<SPACEDIM,int> myTree2(bbox2,0,0,nDesc2Cell,-eps);
+
+ // Edges - handle longest first
+ MCAuto<MEDCouplingFieldDouble> lenF = mDesc2->getMeasureField(true);
+ DataArrayDouble * lens = lenF->getArray();
+
+ // Sort edges by decreasing length:
+ vector<pair<double,int>> S;
+ for(int i=0;i < lens->getNumberOfTuples();i++){
+ pair<double,int> p = make_pair(lens->getIJ(i, 0), i);
+ S.push_back(p);
+ }
+ sort(S.rbegin(),S.rend()); // reverse sort
+ for( vector<pair<double,int>>::const_iterator it = S.begin(); it != S.end(); it++)
{
- MCAuto<DataArrayInt> tmp = dsi->findIdsInRange(0, 2);
- if (tmp->getNumberOfTuples() != dsi->getNumberOfTuples())
- throw INTERP_KERNEL::Exception("Non expected non-conformity. Only simple (=partition-like) non-conformities are handled.");
- }
- // TODO check matching total surface!
-
- MCAuto<DataArrayInt> ids = dsi->findIdsEqual(1);
- // Boundary face:
- if (!ids->getNumberOfTuples())
- continue;
-
- for (const int * ii = ids->begin(); ii != ids->end(); ii++)
- hit[cands2[idsGoodPlane->getIJ(*ii,0)]] = true;
-
- // For all polyhedrons using this face, replace connectivity:
- vector<int> polyIndices, packsIds, facePack;
- for (int ii=*(revDescIP+eIdx); ii < *(revDescIP+eIdx+1); ii++)
- polyIndices.push_back(*(revDescP+ii));
- // Current face connectivity
- const int * sIdxConn = cDesc + cIDesc[eIdx] + 1;
- const int * sIdxConnE = cDesc + cIDesc[eIdx+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, packsIds[jj]);
- // Insertion of new faces:
- for (const int * ii = ids->begin(); ii != ids->end(); ii++)
+ int eIdx = (*it).second;
+
+ vector<int> candidates, cands2;
+ myTree2.getIntersectingElems(bbox2+eIdx*2*SPACEDIM,candidates);
+ // Keep only candidates colinear with current edge
+ double vCurr[3];
+ unsigned start = cDesc2[cIDesc2[eIdx]+1], end = cDesc2[cIDesc2[eIdx]+2];
+ for (int i3=0; i3 < 3; i3++) // TODO: use fillSonCellNodalConnectivity2 or similar?
+ vCurr[i3] = coo[start*SPACEDIM+i3] - coo[end*SPACEDIM+i3];
+ for(vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
+ {
+ double vOther[3];
+ unsigned start2 = cDesc2[cIDesc2[*it]+1], end2 = cDesc2[cIDesc2[*it]+2];
+ for (int i3=0; i3 < 3; i3++)
+ vOther[i3] = coo[start2*SPACEDIM+i3] - coo[end2*SPACEDIM+i3];
+ bool col = INTERP_KERNEL::isColinear3D(vCurr, vOther, eps);
+ // Warning: different from faces: we need to keep eIdx in the final list of candidates because we need
+ // to have its nodes inside the sub mesh mPartCand below (needed in OrderPointsAlongLine())
+ if (col)
+ cands2.push_back(*it);
+ }
+ if (cands2.size() == 1 && cands2[0] == eIdx) // see warning above
+ continue;
+
+ // Now rotate edges to bring them on Ox
+ int startNode = cDesc2[cIDesc2[eIdx]+1];
+ int endNode = cDesc2[cIDesc2[eIdx]+2];
+ INTERP_KERNEL::TranslationRotationMatrix rotation;
+ INTERP_KERNEL::TranslationRotationMatrix::Rotate3DBipoint(coo+SPACEDIM*startNode, coo+SPACEDIM*endNode, rotation);
+ MCAuto<MEDCouplingUMesh> mPartRef(mDesc2->buildPartOfMySelfSlice(eIdx, eIdx+1,1,false)); // false=zipCoords is called
+ MCAuto<MEDCouplingUMesh> mPartCand(mDesc2->buildPartOfMySelf(&cands2[0], &cands2[0]+cands2.size(), true)); // true=zipCoords is called
+ MCAuto<DataArrayInt> nodeMap(mPartCand->zipCoordsTraducer());
+ int nbElemsNotM1;
{
- // Build pack from the face to insert:
- int faceIdx = cands2[idsGoodPlane->getIJ(*ii,0)];
- 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, facePack2+1, facePack2+facePack2Sz); // without the type
+ MCAuto<DataArrayInt> tmp(nodeMap->findIdsNotEqual(-1));
+ nbElemsNotM1 = tmp->getNbOfElems();
}
- }
-
+ MCAuto<DataArrayInt> nodeMapInv = nodeMap->invertArrayO2N2N2O(nbElemsNotM1);
+ 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);
+ for (int ii = 0; ii < mPartCand->_coords->getNumberOfTuples(); ii++)
+ rotation.transform_vector(cooPartCand+SPACEDIM*ii);
+ double ze_y = cooPartRef[1], ze_z = cooPartRef[2];
+
+ // Eliminate all edges for which y or z is not null
+ MCAuto<DataArrayDouble> baryPart = mPartCand->computeCellCenterOfMass();
+ vector<int> compo; compo.push_back(1);
+ MCAuto<DataArrayDouble> baryPartY = baryPart->keepSelectedComponents(compo);
+ compo[0] = 2;
+ MCAuto<DataArrayDouble> baryPartZ = baryPart->keepSelectedComponents(compo);
+ MCAuto<DataArrayInt> idsGoodLine1 = baryPartY->findIdsInRange(-eps, +eps);
+ MCAuto<DataArrayInt> idsGoodLine2 = baryPartZ->findIdsInRange(-eps, +eps);
+ MCAuto<DataArrayInt> idsGoodLine = idsGoodLine1->buildIntersection(idsGoodLine2);
+ if (!idsGoodLine->getNumberOfTuples())
+ continue;
+
+ // Now the ordering along the Ox axis:
+ std::vector<int> insidePoints;
+ bool isSplit = OrderPointsAlongLine(mPartCand->_coords->getConstPointer(), nodeMap->begin()[startNode], nodeMap->begin()[endNode],
+ mPartCand->getNodalConnectivity()->begin(), mPartCand->getNodalConnectivityIndex()->begin(),
+ idsGoodLine->begin(), idsGoodLine->end(),
+ /*out*/insidePoints);
+ std::cout << "edge " << eIdx << "\n";
+ for (int kk =0; kk < insidePoints.size(); kk++)
+ std::cout << " " << nodeMapInv->begin()[insidePoints[kk]] << "\n";
+
+ if (!isSplit) // current segment remains in one piece
+ continue;
+
+ // Get original node IDs in global coords array
+ for (std::vector<int>::iterator iit = insidePoints.begin(); iit!=insidePoints.end(); ++iit)
+ *iit = nodeMapInv->begin()[*iit];
+
+ vector<int> polyIndices, packsIds, facePack;
+ // For each face implying this edge
+// for (int ii=revDescIP2[eIdx]; ii < revDescIP2[eIdx+1]; ii++)
+// {
+// int faceIdx = abs(revDescP2[ii]) - 1;
+//
+// //connSlaDesc
+//
+// // Current face connectivity - find it in the original connectivity
+// const int * sIdxConn = cDesc + cIDesc[faceIdx] + 1;
+// const int * sIdxConnE = cDesc + cIDesc[faceIdx+1];
+//
+// vector<int> modifiedFace;
+// ReplaceEdgeInFace(sIdxConn, sIdxConnE, startNode, endNode, insidePoints, /*out*/modifiedFace);
+// std::cout << "modifiedFace WAS " << eIdx << "\n";
+// for (const int * kk =sIdxConn; kk < sIdxConnE; kk++)
+// std::cout << " " << *kk << "\n";
+// std::cout << "modifiedFace IS " << eIdx << "\n";
+// for (int kk =0; kk < modifiedFace.size(); kk++)
+// std::cout << " " << modifiedFace[kk] << "\n";
+//
+// int k3 = 0;
+// std::cout << "BEFORE\n";
+// std::cout << connSla->simpleRepr();
+// connSlaDesc->replacePack(faceIdx, modifiedFace);
+// std::cout << "AFTER\n";
+// std::cout << connSla->simpleRepr();
+// }
+ }
+ }
// Set back modified connectivity
- MCAuto<DataArrayInt> cAuto; cAuto.takeRef(_nodal_connec);
- MCAuto<DataArrayInt> cIAuto; cIAuto.takeRef(_nodal_connec_index);
connSla->convertToPolyhedronConn(cAuto, cIAuto);
-
-
- MCAuto<DataArrayInt> ret(DataArrayInt::New());
+ ret = ret->buildUniqueNotSorted();
return ret.retn();
}
pass
def testSwig2Hexa8HavingFacesWarped1(self):
- """ This test is bases on a "error" of interpolation detected. After investigation cell #3 of src is warped that leads to the fact that when trg is
+ """ This test is bases on a "error" of interpolation detected. After investigation cell #3 of src is warped that leads to the fact that when trg is
intersected with src the sum of intersection volume is greater than the volume of the trg cell.
A test that can be done is to split the cell #3 of src into tetrohedrons and by summing all the volumes it does not fit the volume computed of cell#3 unsplitted (expect for
GENERAL_24).
sla0.set3( superi.deepCopy(), index.deepCopy(), value.deepCopy() )
self.assertTrue( superi.isEqual( sla0.getSuperIndexArray() ))
- pack = sla0.getPackSafe(2)
+ pack = sla0.getSimplePackSafe(2)
self.assertEqual([9,10,1,12], pack)
ids = sla0.findPackIds([0,1], [9,10,1,12])
self.assertEqual([-1,1], ids)
self.assertEqual(siRef, si.getValues())
self.assertEqual(idxRef, idx.getValues())
self.assertEqual(cRef, val.getValues())
- pass
+ idxRef2 = [0,4,8,12,16,20,23,26,29,32,36,40,42,46,50,53,56, 60, 64, 68, 72, 76, 80 ]
+ cRef2 = [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, 300,300, 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.replacePack(1,2, [300,300])
+ si, idx, val = sla0.getSuperIndexArray(), sla0.getIndexArray(), sla0.getValuesArray()
+ self.assertEqual(siRef, si.getValues())
+ self.assertEqual(idxRef2, idx.getValues())
+ self.assertEqual(cRef2, val.getValues())
+
+ sla0.replacePack(1,2, [9,11,5,1])
+ si, idx, val = sla0.getSuperIndexArray(), sla0.getIndexArray(), sla0.getValuesArray()
+ self.assertEqual(siRef, si.getValues())
+ self.assertEqual(idxRef, idx.getValues())
+ self.assertEqual(cRef, val.getValues())
+
+ sla0.replaceSimplePack(11, [300,300]) # 11 is the abs index of pack (superIdx=1,idx=2)
+ si, idx, val = sla0.getSuperIndexArray(), sla0.getIndexArray(), sla0.getValuesArray()
+ self.assertEqual(siRef, si.getValues())
+ self.assertEqual(idxRef2, idx.getValues())
+ self.assertEqual(cRef2, val.getValues())
+
+ sla0.replaceSimplePack(11, [9,11,5,1]) # 11 is the abs index of pack (superIdx=1,idx=2)
+ 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
DataArrayInt* getIndexArray() const;
DataArrayInt* getValuesArray() const;
- void deletePack(const int i, const int j) throw(INTERP_KERNEL::Exception);
+ void deletePack(const int i, const int j) throw(INTERP_KERNEL::Exception);
%extend
{
return self->simpleRepr();
}
- PyObject *getPackSafe(int absolutePackId) const throw(INTERP_KERNEL::Exception)
+ PyObject *getSimplePackSafe(int absolutePackId) const throw(INTERP_KERNEL::Exception)
{
std::vector<int> ret;
- self->getPackSafe(absolutePackId,ret);
+ self->getSimplePackSafe(absolutePackId,ret);
return convertIntArrToPyList2(ret);
}
self->pushBackPack(i,vpack.data(), vpack.data()+vpack.size());
}
+ void replaceSimplePack(const int idx, PyObject *pack) throw(INTERP_KERNEL::Exception)
+ {
+ std::vector<int> vpack;
+ convertPyToNewIntArr3(pack,vpack);
+ self->replaceSimplePack(idx, vpack.data(), vpack.data()+vpack.size());
+ }
+
+ void replacePack(const int superIdx, const int idx, PyObject *pack) throw(INTERP_KERNEL::Exception)
+ {
+ std::vector<int> vpack;
+ convertPyToNewIntArr3(pack,vpack);
+ self->replacePack(superIdx, idx, vpack.data(), vpack.data()+vpack.size());
+ }
+
PyObject *convertToPolyhedronConn() const throw(INTERP_KERNEL::Exception)
{
MCAuto<DataArrayInt> d0=DataArrayInt::New();
def testSwig2Intersect2DMeshesQuadra1(self):
import cmath
def createDiagCircle(lX, lY, R, cells=[0,1]):
- """ A circle in a square box, cut along the diagonal.
+ """ A circle in a square box, cut along the diagonal.
"""
c = []
for i in range(8):
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()
+ mesh = MEDCouplingUMesh('merge', 3)
+ coo = DataArrayDouble([(0,0,0),(1,0,0),(0,1,0),(1,1,0),(0,0,1),(1,0,1),(0,1,1),(1,1,1),(1,0.5,0.5),
+ (2,0,0),(2,1,0),(2,0,1),(2,1,1),(1,0,2),(2,0,2),(1,1,2),(2,1,2)])
+ mesh.setCoords(coo)
+ c = DataArrayInt([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, 31, 9, 1, 3, 10, -1, 11, 12, 7, 5, -1,
+ 9, 11, 5, 1, -1, 1, 5, 7, 3, -1, 3, 7, 12, 10, -1, 10, 12, 11, 9, 31, 11, 5, 7, 12, -1,
+ 14, 16, 15, 13, -1, 11, 14, 13, 5, -1, 5, 13, 15, 7, -1, 7, 15, 16, 12, -1, 12, 16, 14, 11])
+ cI = DataArrayInt([0, 41, 71, 101])
+ mesh.setConnectivity(c, cI)
+ ret = mesh.conformize3D(1.0e-8)
+
+ mretDesc, _, _, _, _ = mesh.buildDescendingConnectivity()
+ mretDesc.writeVTK("/tmp/toto_conf_desc.vtu")
+ c, cI = mretDesc.getNodalConnectivity().getValues(), mretDesc.getNodalConnectivityIndex().getValues()
+ cRef = [5, 1, 0, 2, 3, 5, 5, 7, 6, 4, 5, 1, 5, 4, 0, 5, 0, 4, 6, 2, 5, 2, 6, 7, 3, 5, 3, 7, 8, 5,
+ 7, 5, 8, 5, 5, 1, 8, 5, 1, 3, 8, 5, 9, 1, 3, 10, 5, 11, 12, 7, 5, 5, 9, 11, 5, 1, 5, 3,
+ 7, 12, 10, 5, 10, 12, 11, 9, 5, 14, 16, 15, 13, 5, 11, 14, 13, 5, 5, 5, 13, 15, 7, 5, 7,
+ 15, 16, 12, 5, 12, 16, 14, 11]
+ cIRef = [0, 5, 10, 15, 20, 25, 29, 33, 37, 41, 46, 51, 56, 61, 66, 71, 76, 81, 86, 91]
+ self.assertEqual(19, mretDesc.getNumberOfCells())
+ self.assertEqual(cRef, c)
+ self.assertEqual(cIRef, cI)
+ self.assertEqual([1], ret.getValues())
+ pass
+
+ def testSwig2Conformize3D2(self):
+ """ More advanced test where edge merge is required. """
+ mesh = MEDCouplingUMesh('merge', 3)
+ coo = DataArrayDouble([(0,0,0),(1,0,0),(0,1,0),(1,1,0),(0,0,1),(1,0,1),(0,1,1),(1,1,1),(0,0,2),(1,0,2),(0,1,2),
+ (1,1,2),(1,0.6,0),(1,0.6,1),(1,0.3,1),(1,0.3,2),(2,0,0),(2,1,0),(2,0,1),(2,1,1),(2,0,2),(2,1,2)])
+ mesh.setCoords(coo)
+ c = DataArrayInt([31, 1, 0, 2, 3, 12, -1, 5, 13, 7, 6, 4, -1, 1, 5, 4, 0, -1, 0, 4, 6, 2, -1, 2, 6, 7, 3, -1, 3, 7, 13, 12, -1,
+ 5, 1, 12, 13, 31, 5, 4, 6, 7, 14, -1, 9, 15, 11, 10, 8, -1, 5, 9, 8, 4, -1, 4, 8, 10, 6, -1, 6, 10, 11, 7, -1,
+ 7, 11, 15, 14, -1, 9, 5, 14, 15, 31, 16, 1, 3, 17, -1, 18, 19, 7, 5, -1, 16, 18, 5, 1, -1, 1, 5, 7, 3, -1,
+ 3, 7, 19, 17, -1, 17, 19, 18, 16, 31, 18, 5, 7, 19, -1, 20, 21, 11, 9, -1, 18, 20, 9, 5, -1, 5, 9, 11, 7, -1,
+ 7, 11, 21, 19, -1, 19, 21, 20, 18])
+ cI = DataArrayInt([0, 37, 74, 104, 134])
+ mesh.setConnectivity(c, cI)
+
+ mretDesc, _, _, _, _ = mesh.buildDescendingConnectivity()
+ mretDesc2, _, _, _, _ = mretDesc.buildDescendingConnectivity()
mretDesc.writeVTK("/tmp/toto_desc.vtu")
+ mretDesc2.writeVTK("/tmp/toto_desc2.vtu")
- mret.conformize3D(1.0e-8)
+ ret = mesh.conformize3D(1.0e-8)
- mret.writeVTK("/tmp/toto_conf.vtu")
- mretDesc, _, _, _, _ = mret.buildDescendingConnectivity()
+ mretDesc, _, _, _, _ = mesh.buildDescendingConnectivity()
+ mretDesc2, _, _, _, _ = mretDesc.buildDescendingConnectivity()
mretDesc.writeVTK("/tmp/toto_conf_desc.vtu")
- return mret
+ mretDesc2.writeVTK("/tmp/toto_conf_desc2.vtu")
+ c, cI = mretDesc.getNodalConnectivity().getValues(), mretDesc.getNodalConnectivityIndex().getValues()
+ c2, cI2 = mretDesc.getNodalConnectivity().getValues(), mretDesc.getNodalConnectivityIndex().getValues()
+ cRef = []
+ cIRef = []
+ cRef2 = []
+ cIRef2 = []
+ print mretDesc2.getNumberOfCells()
+ return
+ self.assertEqual(22, mretDesc.getNumberOfCells())
+ self.assertEqual(40, mretDesc2.getNumberOfCells())
+ self.assertEqual(cRef, c)
+ self.assertEqual(cIRef, cI)
+ self.assertEqual(cRef2, c2)
+ self.assertEqual(cIRef2, cI2)
+ self.assertEqual([0,1,2,3], ret.getValues())
+ pass
+
if __name__ == '__main__':
unittest.main()
