class DiameterCalculator;
/*!
- * This class descibes all static elements (different from polygons and polyhedron) 3D, 2D and 1D.
+ * This class describes all static elements (different from polygons and polyhedron) 3D, 2D and 1D.
*/
class CellModel
{
void getRealSourceCoordinatesPermute(ConnType icellS, int offset, std::vector<double>& coordsS);
void getRealCoordinates(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS, std::vector<double>& coordsT, std::vector<double>& coordsS, int& orientation);
double getValueRegardingOption(double val) const;
- static void rotate3DTriangle( double* PP1, double*PP2, double*PP3,
+ static void Rotate3DTriangle( double* PP1, double*PP2, double*PP3,
TranslationRotationMatrix& rotation_matrix);
protected:
const ConnType *_connectT;
{
TranslationRotationMatrix rotation;
//rotate3DTriangle(Coords_A, &Coords_A[SPACEDIM*i_A1], &Coords_A[SPACEDIM*i_A2], rotation);
- rotate3DTriangle(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2, rotation);
+ Rotate3DTriangle(Coords_B, Coords_B+SPACEDIM*i_B1, Coords_B+SPACEDIM*i_B2, rotation);
for (int i=0; i<nb_NodesA; i++)
rotation.transform_vector(Coords_A+SPACEDIM*i);
for (int i=0; i<nb_NodesB; i++)
}
template<class MyMeshType, class MyMatrix>
- void PlanarIntersector<MyMeshType,MyMatrix>::rotate3DTriangle(double* PP1, double*PP2, double*PP3,
+ void PlanarIntersector<MyMeshType,MyMatrix>::Rotate3DTriangle(double* PP1, double*PP2, double*PP3,
TranslationRotationMatrix& rotation_matrix)
{
//initializes
rotate_vector(P);
}
-
+ static void Rotate3DTriangle(const double* PP1,const double*PP2,const double*PP3, TranslationRotationMatrix& rotation_matrix)
+ {
+ //initializes
+ rotation_matrix.translate(PP1);
+
+ double P2w[3];
+ double P3w[3];
+ P2w[0]=PP2[0]; P2w[1]=PP2[1];P2w[2]=PP2[2];
+ P3w[0]=PP3[0]; P3w[1]=PP3[1];P3w[2]=PP3[2];
+
+ // translating to set P1 at the origin
+ for (int i=0; i<3; i++)
+ {
+ P2w[i]-=PP1[i];
+ P3w[i]-=PP1[i];
+ }
+
+ // rotating to set P2 on the Oxy plane
+ TranslationRotationMatrix A;
+ A.rotate_x(P2w);
+ A.rotate_vector(P3w);
+ rotation_matrix.multiply(A);
+
+ //rotating to set P2 on the Ox axis
+ TranslationRotationMatrix B;
+ B.rotate_z(P2w);
+ B.rotate_vector(P3w);
+ rotation_matrix.multiply(B);
+
+ //rotating to set P3 on the Oxy plane
+ TranslationRotationMatrix C;
+ C.rotate_x(P3w);
+ rotation_matrix.multiply(C);
+ }
+
+
+
private:
static const double EPS;
static const unsigned ROT_SIZE=9;
// return std::fabs(x - y) < errTol;
}
+
+ /**
+ * Test whether two 3D vectors are colinear. The two vectors are expected to be of unit norm (not checked)
+ * Implemented by checking that the norm of the cross product is null.
+ */
+ inline bool isColinear3D(const double *v1, const double *v2, const double eps = DEFAULT_ABS_TOL)
+ {
+ double cros[3];
+ cross(v1, v2, cros);
+ return epsilonEqual(dot(cros, cros), 0.0, eps);
+ }
+
+
/**
* Compares doubles using a relative tolerance
* This is suitable mainly for comparing larger values to each other. Before performing the relative test,
using namespace MEDCoupling;
MEDCouplingSkyLineArray::MEDCouplingSkyLineArray():
- _index( DataArrayInt::New() ), _values( DataArrayInt::New() ), _sub_values( DataArrayInt::New() )
+ _index( DataArrayInt::New() ), _values( DataArrayInt::New() ), _super_index( DataArrayInt::New() )
{
}
MEDCouplingSkyLineArray* MEDCouplingSkyLineArray::New( const MEDCouplingSkyLineArray & other )
{
MEDCouplingSkyLineArray* ret = new MEDCouplingSkyLineArray();
+ ret->_super_index = other._super_index;
ret->_index = other._index;
ret->_values = other._values;
- ret->_sub_values = other._sub_values;
return ret;
}
+/**! Build a three level SkyLine array from the dynamic connectivity of a dynamic mesh (i.e. containing only
+ * polyhedrons or polygons).
+ * The input arrays are deep copied, contrary to the other ctors.
+ */
+MEDCouplingSkyLineArray * MEDCouplingSkyLineArray::BuildFromPolyhedronConn( const DataArrayInt* c, const DataArrayInt* cI )
+{
+ using namespace std;
+
+ MEDCouplingSkyLineArray* ret = new MEDCouplingSkyLineArray();
+
+ const int * cP(c->begin()), * cIP(cI->begin());
+ int prev = -1;
+ if (c->getNbOfElems() != *(cI->end()-1))
+ throw INTERP_KERNEL::Exception("MEDCouplingSkyLineArray::BuildFromDynamicConn: misformatted connectivity (wrong nb of tuples)!");
+ for (int i=0; i < cI->getNbOfElems(); i++)
+ {
+ int j = cIP[i];
+ if (cIP[i] < prev)
+ throw INTERP_KERNEL::Exception("MEDCouplingSkyLineArray::BuildFromDynamicConn: misformatted connectivity (indices not monotonic ascending)!");
+ prev = cIP[i];
+ if (i!=cI->getNbOfElems()-1)
+ if (cP[j] != INTERP_KERNEL::NORM_POLYHED)
+ throw INTERP_KERNEL::Exception("MEDCouplingSkyLineArray::BuildFromDynamicConn: connectivity containing other types than POLYHED!");
+ }
+
+ vector<int> superIdx, idx, vals;
+ int cnt = 0, cnt2 = 0;
+ superIdx.reserve(cI->getNbOfElems());
+ superIdx.push_back(0);
+ idx.push_back(0);
+ vals.resize(c->getNbOfElems()); // too much because of the type and the -1, but still better than push_back().
+ for (int i=0; i < cI->getNbOfElems()-1; i++)
+ {
+ int start = cIP[i]+1, end = cIP[i+1];
+ int * work = vals.data() + cnt;
+ const int * w = cP+start;
+ const int * w2 = find(w, cP+end, -1);
+ while (w2 != cP+end)
+ {
+ copy(w, w2, work);
+ int d = distance(w, w2);
+ cnt += d; work +=d;
+ idx.push_back(cnt); cnt2++;
+ w = w2+1; // skip the -1
+ w2 = find(w, cP+end, -1);
+ }
+ copy(w, cP+end, work);
+ cnt += distance(w, cP+end);
+ idx.push_back(cnt); cnt2++;
+ superIdx.push_back(cnt2);
+ }
+ ret->_super_index->alloc(superIdx.size(),1);
+ copy(superIdx.begin(), superIdx.end(), ret->_super_index->getPointer());
+ ret->_index->alloc(idx.size(),1);
+ copy(idx.begin(), idx.end(), ret->_index->getPointer());
+ ret->_values->alloc(cnt,1);
+ copy(vals.begin(), vals.begin()+cnt, ret->_values->getPointer());
+
+ return ret;
+}
+
+/**
+ * Convert a three-level SkyLineArray into a polyhedral connectivity.
+ * The super-packs are interpreted as cell description, and the packs represent the face connectivity.
+ */
+void MEDCouplingSkyLineArray::convertToPolyhedronConn( MCAuto<DataArrayInt>& c, MCAuto<DataArrayInt>& cI) const
+{
+ // TODO: in this case an iterator would be nice
+ using namespace std;
+
+ checkSuperIndex("convertToPolyhedronConn");
+
+ 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;
+ 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);
+
+ for (int i=0; i < _super_index->getNbOfElems()-1; i++)
+ {
+ cIVec.push_back(cnt);
+ int endId = siP[i+1];
+ cVecP[cnt++] = INTERP_KERNEL::NORM_POLYHED;
+ for (int j=siP[i]; j < endId; j++)
+ {
+ int startId2 = iP[j], endId2 = iP[j+1];
+ copy(vP+startId2, vP+endId2, cVecP+cnt);
+ cnt += endId2-startId2;
+ if(j != endId-1)
+ cVecP[cnt++] = -1;
+ }
+ }
+ cIVec.push_back(cnt);
+ copy(cIVec.begin(), cIVec.end(), cI->getPointer());
+ copy(cVec.begin(), cVec.begin()+cnt, c->getPointer());
+}
std::size_t MEDCouplingSkyLineArray::getHeapMemorySizeWithoutChildren() const
{
- return _index->getHeapMemorySizeWithoutChildren()+_values->getHeapMemorySizeWithoutChildren()+_sub_values->getHeapMemorySizeWithoutChildren();
+ return _index->getHeapMemorySizeWithoutChildren()+_values->getHeapMemorySizeWithoutChildren()+_super_index->getHeapMemorySizeWithoutChildren();
}
std::vector<const BigMemoryObject *> MEDCouplingSkyLineArray::getDirectChildrenWithNull() const
{
std::vector<const BigMemoryObject *> ret;
+ ret.push_back(_super_index);
ret.push_back(_index);
ret.push_back(_values);
- ret.push_back(_sub_values);
return ret;
}
else _values = DataArrayInt::New();
}
+void MEDCouplingSkyLineArray::set3( DataArrayInt* superIndex, DataArrayInt* index, DataArrayInt* value )
+{
+ _super_index=superIndex;
+ if ( (DataArrayInt*)_super_index ) _super_index->incrRef();
+ else _super_index = DataArrayInt::New();
+ set(index, value);
+}
+
+DataArrayInt* MEDCouplingSkyLineArray::getSuperIndexArray() const
+{
+ return ((MEDCouplingSkyLineArray*)this)->_super_index;
+}
+
+
DataArrayInt* MEDCouplingSkyLineArray::getIndexArray() const
{
return ((MEDCouplingSkyLineArray*)this)->_index;
return ((MEDCouplingSkyLineArray*)this)->_values;
}
+void MEDCouplingSkyLineArray::checkSuperIndex(const std::string& func) const
+{
+ if (!_super_index->getNbOfElems())
+ {
+ std::ostringstream oss;
+ oss << "MEDCouplingSkyLineArray::"<< func << ": not a three level SkyLineArray! Method is not available for two-level SkyLineArray.";
+ throw INTERP_KERNEL::Exception(oss.str());
+ }
+}
+
std::string MEDCouplingSkyLineArray::simpleRepr() const
{
std::ostringstream oss;
- oss << "MEDCouplingSkyLineArray" << std::endl;
- oss << " Nb of items: " << getNumberOf() << std::endl;
+ oss << "MEDCouplingSkyLineArray (" << this << ")" << std::endl;
+ MCAuto<DataArrayInt> super_index = _super_index->deepCopy();
+ if (_super_index->getNbOfElems())
+ oss << " Nb of super-packs: " << getSuperNumberOf() << std::endl;
+ else
+ {
+ super_index->alloc(2,1);
+ super_index->setIJSilent(0,0,0);
+ super_index->setIJSilent(1,0,_index->getNbOfElems()-1);
+ }
+ oss << " Nb of packs: " << getNumberOf() << std::endl;
oss << " Nb of values: " << getLength() << std::endl;
- oss << " Index:" << std::endl;
+
+ if (_super_index->getNbOfElems())
+ {
+ oss << " Super-indices:" << std::endl;
+ oss << " ";
+ const int * i = _super_index->begin();
+ for ( ; i != _super_index->end(); ++i )
+ oss << *i << " ";
+ oss << std::endl;
+ }
+
+ oss << " Indices:" << std::endl;
oss << " ";
const int * i = _index->begin();
for ( ; i != _index->end(); ++i )
oss << *i << " ";
oss << std::endl;
- oss << " Value:" << std::endl;
+ oss << " Values:" << std::endl;
oss << " ";
const int * v = _values->begin();
int cnt = 0;
return oss.str();
}
+
+void MEDCouplingSkyLineArray::getPackSafe(const int absolutePackId, std::vector<int> & pack) const
+{
+ if(absolutePackId < 0 || absolutePackId >= _index->getNbOfElems())
+ throw INTERP_KERNEL::Exception("MEDCouplingSkyLineArray::getPackSafe: invalid index!");
+ const int * iP(_index->begin()), *vP(_values->begin());
+ int sz = iP[absolutePackId+1]-iP[absolutePackId];
+ pack.resize(sz);
+ std::copy(vP+iP[absolutePackId], vP+iP[absolutePackId+1],pack.begin());
+}
+
+
+/**!
+ * 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[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
+{
+ using namespace std;
+
+ checkSuperIndex("findPackIds");
+ out.resize(superPackIndices.size());
+ int i = 0;
+ const int * siP(_super_index->begin()), * iP(_index->begin()), *vP(_values->begin());
+ for(vector<int>::const_iterator it=superPackIndices.begin(); it!=superPackIndices.end(); ++it, i++)
+ {
+ out[i] = -1;
+ 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;
+ }
+ }
+}
+
+/**!
+ * Delete pack number j in super-pack number i.
+ * \param[in] i is the super-pack number
+ * \param[in] j is the pack index inside the super-pack i.
+ */
+void MEDCouplingSkyLineArray::deletePack(const int i, const int j)
+{
+ using namespace std;
+
+ checkSuperIndex("deletePack");
+
+ int * vP = _values->getPointer();
+ int * siP(_super_index->getPointer()), *iP(_index->getPointer());
+ const int start = iP[siP[i]+j], end = iP[siP[i]+j+1];
+ // _values
+ copy(vP+end, vP+_values->getNbOfElems(), vP+start);
+ _values->reAlloc(_values->getNbOfElems() - (end-start));
+
+ // _index
+ 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++)
+ iP[ii] -= (end-start);
+
+ // _super_index
+ for(int ii = i+1; ii < _super_index->getNbOfElems(); ii++)
+ (siP[ii])--;
+}
+
+/**!
+ * 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)
+{
+ using namespace std;
+
+ checkSuperIndex("pushBackPack");
+
+ int *siP(_super_index->getPointer()), *iP(_index->getPointer());
+ const int sz(pack.size());
+
+ // _values
+ _values->reAlloc(_values->getNbOfElems()+sz);
+ int * vPE(_values->getPointer()+_values->getNbOfElems());
+ 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]]);
+
+ // _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();
+ for(int ii = siP[i+1]+2; ii < nt+1; ii++)
+ iP[ii] += sz;
+
+ // _super_index
+ for(int ii = i+1; ii < _super_index->getNbOfElems(); ii++)
+ (siP[ii])++;
+}
+
namespace MEDCoupling
{
/**!
- * Class allowing the easy manipulation of the indexed array format, where the first array is a set of offsets to
- * be used in the second array, to extract packs of values.
+ * Class allowing the easy manipulation of the indexed array format, where the first array (the "indices") is a set of offsets to
+ * be used in the second array (the "values"), to extract packs of values.
*
- * This class allows to pursuie this logic up to 3 levels, i.e. the first array points to packs in the second, which
- * itself points to identifiers in the third array.
+ * Example:
+ * index = [0,2,7,10]
+ * values = [1,24,2,33,6,7,10,11,9,28]
+ * which describes 3 packs of (integer) values : [1,24] and [2,33,6,7,10] and [11,9,28]
+ *
+ * Thus the index array is always monotic ascendant.
+ *
+ * This class allows to pursue this logic up to 3 levels, i.e. the first array (the "indices") points to packs in the
+ * second (the "values"), which itself points to identifiers in a third array (the "sub-values").
*
* This particularly useful for connectivity of pure polygonal/polyhedral meshes.
+ *
+ * Example:
+ * super-index = [0,1,3]
+ * index = [0,3,6,10]
+ * values = [28,1,4,2,35,8,9,10,1,12]
+ * which represent two 3 packs and two super-packs. The first super-pack is [[28,1,4]] and has only one pack [28,1,4].
+ * The second super-pack is [[2,35,8], [9,10,1,12]] and has two packs [2,35,8] and [9,10,1,12].
+ * Note that contrary to index, the integers in super-index are interpreted as being inclusive: the first super-pack
+ * goes from offset 0 to offset 1, inclusive. This is not the same for index, where the upper bound is exclusive.
*/
class MEDCOUPLING_EXPORT MEDCouplingSkyLineArray : public RefCountObject
{
static MEDCouplingSkyLineArray * New( DataArrayInt* index, DataArrayInt* value );
static MEDCouplingSkyLineArray * New( const MEDCouplingSkyLineArray & other );
+ static MEDCouplingSkyLineArray * BuildFromPolyhedronConn( const DataArrayInt* c, const DataArrayInt* cI );
+
std::size_t getHeapMemorySizeWithoutChildren() const;
std::vector<const BigMemoryObject *> getDirectChildrenWithNull() const;
void set( DataArrayInt* index, DataArrayInt* value );
+ void set3( DataArrayInt* superIndex, DataArrayInt* index, DataArrayInt* value );
+ int getSuperNumberOf() const { return _super_index->getNbOfElems()-1; }
int getNumberOf() const { return _index->getNbOfElems()-1; }
int getLength() const { return _values->getNbOfElems(); }
+
+ const int* getSuperIndex() const { return _super_index->begin(); }
const int* getIndex() const { return _index->begin(); }
const int* getValues() const { return _values->begin(); }
+ DataArrayInt* getSuperIndexArray() const;
DataArrayInt* getIndexArray() const;
DataArrayInt* getValuesArray() const;
std::string simpleRepr() const;
-// replaceWithPackFromOther()
+ 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;
+
+ void deletePack(const int i, const int j);
+ void pushBackPack(const int i, const std::vector<int> & pack);
+
+ void convertToPolyhedronConn( MCAuto<DataArrayInt>& c, MCAuto<DataArrayInt>& cI) const;
private:
MEDCouplingSkyLineArray();
~MEDCouplingSkyLineArray();
+ void checkSuperIndex(const std::string& func) const;
+
+ MCAuto<DataArrayInt> _super_index;
MCAuto<DataArrayInt> _index;
MCAuto<DataArrayInt> _values;
- MCAuto<DataArrayInt> _sub_values;
};
}
# endif
MEDCOUPLING_EXPORT DataArrayDouble *computePlaneEquationOf3DFaces() const;
MEDCOUPLING_EXPORT DataArrayInt *conformize2D(double eps);
MEDCOUPLING_EXPORT DataArrayInt *colinearize2D(double eps);
+ MEDCOUPLING_EXPORT DataArrayInt *conformize3D(double eps);
MEDCOUPLING_EXPORT int split2DCells(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *midOpt=0, const DataArrayInt *midOptI=0);
MEDCOUPLING_EXPORT static MEDCouplingUMesh *Build0DMeshFromCoords(DataArrayDouble *da);
MEDCOUPLING_EXPORT static MCAuto<MEDCouplingUMesh> Build1DMeshFromCoords(DataArrayDouble *da);
#include "MEDCouplingUMesh.hxx"
#include "MEDCoupling1GTUMesh.hxx"
+#include "MEDCouplingFieldDouble.hxx"
#include "CellModel.hxx"
#include "VolSurfUser.txx"
#include "InterpolationUtils.hxx"
#include "InterpKernelGeo2DEdgeLin.hxx"
#include "InterpKernelGeo2DEdgeArcCircle.hxx"
#include "InterpKernelGeo2DQuadraticPolygon.hxx"
+#include "TranslationRotationMatrix.hxx"
+#include "VectorUtils.hxx"
+#include "MEDCouplingSkyLineArray.hxx"
#include <sstream>
#include <fstream>
return ret.retn();
}
+/*!
+ * \b WARNING this method is \b potentially \b non \b const (if returned array is not empty).
+ * \b WARNING this method lead to have a non geometric type sorted mesh (for MED file users) !
+ * This method performs a conformization of \b this.
+ *
+ * Only polyhedron cells are supported. You can call convertAllToPoly()
+ *
+ * This method expects that \b this has a meshDim equal 3 and spaceDim equal to 3 too.
+ * This method expects that all nodes in \a this are not closer than \a eps.
+ * If it is not the case you can invoke MEDCouplingUMesh::mergeNodes before calling this method.
+ *
+ * \param [in] eps the relative error to detect merged edges.
+ * \return DataArrayInt * - The list of cellIds in \a this that have been subdivided. If empty, nothing changed in \a this (as if it were a const method). The array is a newly allocated array
+ * that the user is expected to deal with.
+ *
+ * \throw If \a this is not coherent.
+ * \throw If \a this has not spaceDim equal to 3.
+ * \throw If \a this has not meshDim equal to 3.
+ * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::convertAllToPoly()
+ */
+DataArrayInt *MEDCouplingUMesh::conformize3D(double eps)
+{
+ using namespace std;
+
+ static const int SPACEDIM=3;
+ checkConsistencyLight();
+ if(getSpaceDimension()!=3 || getMeshDimension()!=3)
+ 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->getNodalConnectivity(),
+ mDesc->getNodalConnectivityIndex()));
+ 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
+ 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
+ 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 + (*it)*sIdx, 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(cooDesc+SPACEDIM*(cDesc[cIDesc[sIdx]+1]),
+ 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());
+ 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];
+
+ 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;
+
+ MCAuto<DataArrayInt> cc(DataArrayInt::New()), ccI(DataArrayInt::New());
+ mPartRef->getCellsContainingPoints(baryPart->begin(), baryPart->getNumberOfTuples(), eps, cc, ccI);
+
+ // Localize faces in 2D thanks to barycenters
+ 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.");
+ }
+ // 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;
+
+ // Current face connectivity
+ const int * sIdxConn = cDesc + cIDesc[sIdx] + 1;
+
+ // For all polyhedrons using this face, replace connectivity:
+ vector<int> polyIndices, subPacksIds, facePack;
+ for (int ii=*(revDescIP+sIdx); ii < *(revDescIP+sIdx+1); ii++)
+ polyIndices.push_back(*(revDescP+ii));
+ connSla->findPackIds(polyIndices, facePack, subPacksIds);
+ // 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]);
+ // 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);
+ // Insert it in all hit polyhedrons:
+ for (vector<int>::const_iterator it=polyIndices.begin(); it!=polyIndices.end(); ++it)
+ connSla->pushBackPack(*it, facePack);
+ }
+ }
+
+ // STEP 2 -- edges
+ // TODO
+
+ MCAuto<DataArrayInt> ret(DataArrayInt::New());
+ return ret.retn();
+}
+
pass
+ def testMEDCouplingSkyLineArrayThreeLevels(self):
+ # [[28,1,4]] , [[2,35,8], [9,10,1,12]]
+ superi = DataArrayInt([ 0,1,3 ])
+ index = DataArrayInt ([ 0,3,6,10 ])
+ value = DataArrayInt ([ 28,1,4,2,35,8,9,10,1,12 ])
+
+ sla0 = MEDCouplingSkyLineArray()
+ self.assertEqual( -1, sla0.getSuperNumberOf() )
+ self.assertEqual( -1, sla0.getNumberOf() )
+ self.assertEqual( 0, sla0.getLength() )
+ sla0.set3( superi.deepCopy(), index.deepCopy(), value.deepCopy() )
+ self.assertTrue( superi.isEqual( sla0.getSuperIndexArray() ))
+
+ pack = sla0.getPackSafe(2)
+ self.assertEqual([9,10,1,12], pack)
+ ids = sla0.findPackIds([0,1], [9,10,1,12])
+ self.assertEqual([-1,1], ids)
+
+ sla0.deletePack(1, 0)
+ 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([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())
+
+ # Build connectivity from POLYHED connectivity
+ cI = [0,16,41]
+ c = [NORM_POLYHED, 1,2,3,-1, 2,3,4,-1, 3,4,5,-1, 4,5,6,
+ NORM_POLYHED, 7,8,9,10,-1, 9,10,11,12,-1, 3,4,5,6,-1, 5,6,7,8,-1, 9,10,11,12]
+ sla0 = MEDCouplingSkyLineArray.BuildFromPolyhedronConn(DataArrayInt(c), DataArrayInt(cI))
+ si, idx, val = sla0.getSuperIndexArray(), sla0.getIndexArray(), sla0.getValuesArray()
+ self.assertEqual([0,4,9], si.getValues())
+ self.assertEqual([0,3,6,9,12,16,20,24,28,32], idx.getValues())
+ self.assertEqual([1,2,3, 2,3,4, 3,4,5, 4,5,6,
+ 7,8,9,10, 9,10,11,12, 3,4,5,6, 5,6,7,8, 9,10,11,12], val.getValues())
+ c1, cI1 = sla0.convertToPolyhedronConn()
+ self.assertEqual(c1.getValues(), c)
+ self.assertEqual(cI1.getValues(), cI)
+ pass
+
def testMEDCouplingUMeshgenerateGraph(self):
# cartesian mesh 3x3
arr=DataArrayDouble(4) ; arr.iota()
%newobject MEDCoupling::DenseMatrix::__mul__;
%newobject MEDCoupling::MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell;
%newobject MEDCoupling::MEDCouplingGaussLocalization::buildRefCell;
+%newobject MEDCoupling::MEDCouplingSkyLineArray::BuildFromPolyhedronConn;
%feature("unref") MEDCouplingPointSet "$this->decrRef();"
%feature("unref") MEDCouplingMesh "$this->decrRef();"
class MEDCouplingSkyLineArray
{
public:
+ static MEDCouplingSkyLineArray *BuildFromPolyhedronConn( const DataArrayInt* c, const DataArrayInt* cI ) throw(INTERP_KERNEL::Exception);
+
void set( DataArrayInt* index, DataArrayInt* value );
+ void set3( DataArrayInt* superIndex, DataArrayInt* index, DataArrayInt* value );
+
+ int getSuperNumberOf() const;
int getNumberOf() const;
int getLength() const;
+
+ DataArrayInt* getSuperIndexArray() const;
DataArrayInt* getIndexArray() const;
DataArrayInt* getValuesArray() const;
+
+ void deletePack(const int i, const int j) throw(INTERP_KERNEL::Exception);
+
%extend
{
MEDCouplingSkyLineArray() throw(INTERP_KERNEL::Exception)
{
return self->simpleRepr();
}
-
+
+ PyObject *getPackSafe(int absolutePackId) const throw(INTERP_KERNEL::Exception)
+ {
+ std::vector<int> ret;
+ self->getPackSafe(absolutePackId,ret);
+ return convertIntArrToPyList2(ret);
+ }
+
+ 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);
+ return convertIntArrToPyList2(out);
+ }
+
+ void pushBackPack(const int i, PyObject *pack) throw(INTERP_KERNEL::Exception)
+ {
+ std::vector<int> vpack;
+ convertPyToNewIntArr3(pack,vpack);
+ self->pushBackPack(i,vpack);
+ }
+
+ PyObject *convertToPolyhedronConn() const throw(INTERP_KERNEL::Exception)
+ {
+ MCAuto<DataArrayInt> d0=DataArrayInt::New();
+ MCAuto<DataArrayInt> d1=DataArrayInt::New();
+ self->convertToPolyhedronConn(d0,d1);
+ PyObject *ret=PyTuple_New(2);
+ PyTuple_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(d0.retn()),SWIGTYPE_p_MEDCoupling__DataArrayInt, SWIG_POINTER_OWN | 0 ));
+ PyTuple_SetItem(ret,1,SWIG_NewPointerObj(SWIG_as_voidptr(d1.retn()),SWIGTYPE_p_MEDCoupling__DataArrayInt, SWIG_POINTER_OWN | 0 ));
+ return ret;
+ }
}
};
}
