//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
-// Author : Anthony Geay (CEA/DEN)
+// Author : Anthony Geay (EDF R&D)
#include "MEDCouplingUMesh.hxx"
#include "MEDCouplingCMesh.hxx"
#include "InterpKernelGeo2DEdgeLin.hxx"
#include "InterpKernelGeo2DEdgeArcCircle.hxx"
#include "InterpKernelGeo2DQuadraticPolygon.hxx"
+#include "OrientationInverter.hxx"
#include "MEDCouplingUMesh_internal.hxx"
#include <sstream>
double MEDCouplingUMesh::EPS_FOR_POLYH_ORIENTATION=1.e-14;
/// @cond INTERNAL
-const INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::MEDMEM_ORDER[N_MEDMEM_ORDER] = { INTERP_KERNEL::NORM_POINT1, INTERP_KERNEL::NORM_SEG2, INTERP_KERNEL::NORM_SEG3, INTERP_KERNEL::NORM_SEG4, INTERP_KERNEL::NORM_POLYL, INTERP_KERNEL::NORM_TRI3, INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6, INTERP_KERNEL::NORM_TRI7, INTERP_KERNEL::NORM_QUAD8, INTERP_KERNEL::NORM_QUAD9, INTERP_KERNEL::NORM_POLYGON, INTERP_KERNEL::NORM_QPOLYG, INTERP_KERNEL::NORM_TETRA4, INTERP_KERNEL::NORM_PYRA5, INTERP_KERNEL::NORM_PENTA6, INTERP_KERNEL::NORM_HEXA8, INTERP_KERNEL::NORM_HEXGP12, INTERP_KERNEL::NORM_TETRA10, INTERP_KERNEL::NORM_PYRA13, INTERP_KERNEL::NORM_PENTA15, INTERP_KERNEL::NORM_HEXA20, INTERP_KERNEL::NORM_HEXA27, INTERP_KERNEL::NORM_POLYHED };
-const int MEDCouplingUMesh::MEDCOUPLING2VTKTYPETRADUCER[INTERP_KERNEL::NORM_MAXTYPE+1]={1,3,21,5,9,7,22,34,23,28,-1,-1,-1,-1,10,14,13,-1,12,-1,24,-1,16,27,-1,26,-1,29,-1,-1,25,42,36,4};
+const INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::MEDMEM_ORDER[N_MEDMEM_ORDER] = { INTERP_KERNEL::NORM_POINT1, INTERP_KERNEL::NORM_SEG2, INTERP_KERNEL::NORM_SEG3, INTERP_KERNEL::NORM_SEG4, INTERP_KERNEL::NORM_POLYL, INTERP_KERNEL::NORM_TRI3, INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6, INTERP_KERNEL::NORM_TRI7, INTERP_KERNEL::NORM_QUAD8, INTERP_KERNEL::NORM_QUAD9, INTERP_KERNEL::NORM_POLYGON, INTERP_KERNEL::NORM_QPOLYG, INTERP_KERNEL::NORM_TETRA4, INTERP_KERNEL::NORM_PYRA5, INTERP_KERNEL::NORM_PENTA6, INTERP_KERNEL::NORM_HEXA8, INTERP_KERNEL::NORM_HEXGP12, INTERP_KERNEL::NORM_TETRA10, INTERP_KERNEL::NORM_PYRA13, INTERP_KERNEL::NORM_PENTA15, INTERP_KERNEL::NORM_PENTA18, INTERP_KERNEL::NORM_HEXA20, INTERP_KERNEL::NORM_HEXA27, INTERP_KERNEL::NORM_POLYHED };
+const int MEDCouplingUMesh::MEDCOUPLING2VTKTYPETRADUCER[INTERP_KERNEL::NORM_MAXTYPE+1]={1,3,21,5,9,7,22,34,23,28,-1,-1,-1,-1,10,14,13,-1,12,-1,24,-1,16,27,-1,26,-1,29,32,-1,25,42,36,4};
/// @endcond
MEDCouplingUMesh *MEDCouplingUMesh::New()
void MEDCouplingUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const
{
checkFullyDefined();
- int nbOfNodes=getNumberOfNodes();
+ int nbOfNodes(getNumberOfNodes());
int *revNodalIndxPtr=(int *)malloc((nbOfNodes+1)*sizeof(int));
revNodalIndx->useArray(revNodalIndxPtr,true,C_DEALLOC,nbOfNodes+1,1);
std::fill(revNodalIndxPtr,revNodalIndxPtr+nbOfNodes+1,0);
- const int *conn=_nodal_connec->getConstPointer();
- const int *connIndex=_nodal_connec_index->getConstPointer();
- int nbOfCells=getNumberOfCells();
- int nbOfEltsInRevNodal=0;
+ const int *conn(_nodal_connec->begin()),*connIndex(_nodal_connec_index->begin());
+ int nbOfCells(getNumberOfCells()),nbOfEltsInRevNodal(0);
for(int eltId=0;eltId<nbOfCells;eltId++)
{
- const int *strtNdlConnOfCurCell=conn+connIndex[eltId]+1;
- const int *endNdlConnOfCurCell=conn+connIndex[eltId+1];
+ const int *strtNdlConnOfCurCell(conn+connIndex[eltId]+1),*endNdlConnOfCurCell(conn+connIndex[eltId+1]);
for(const int *iter=strtNdlConnOfCurCell;iter!=endNdlConnOfCurCell;iter++)
if(*iter>=0)//for polyhedrons
{
{
if(!desc || !descIndx || !revDesc || !revDescIndx)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ComputeNeighborsOfCellsAdv some input array is empty !");
- const int *descPtr=desc->getConstPointer();
- const int *descIPtr=descIndx->getConstPointer();
- const int *revDescPtr=revDesc->getConstPointer();
- const int *revDescIPtr=revDescIndx->getConstPointer();
+ const int *descPtr=desc->begin();
+ const int *descIPtr=descIndx->begin();
+ const int *revDescPtr=revDesc->begin();
+ const int *revDescIPtr=revDescIndx->begin();
//
int nbCells=descIndx->getNumberOfTuples()-1;
MCAuto<DataArrayInt> out0=DataArrayInt::New();
neighborsIndx=out1.retn();
}
+/*!
+ * Explodes \a this into edges whatever its dimension.
+ */
+MCAuto<MEDCouplingUMesh> MEDCouplingUMesh::explodeIntoEdges(MCAuto<DataArrayInt>& desc, MCAuto<DataArrayInt>& descIndex, MCAuto<DataArrayInt>& revDesc, MCAuto<DataArrayInt>& revDescIndx) const
+{
+ checkFullyDefined();
+ int mdim(getMeshDimension());
+ desc=DataArrayInt::New(); descIndex=DataArrayInt::New(); revDesc=DataArrayInt::New(); revDescIndx=DataArrayInt::New();
+ MCAuto<MEDCouplingUMesh> mesh1D;
+ switch(mdim)
+ {
+ case 3:
+ {
+ mesh1D=explode3DMeshTo1D(desc,descIndex,revDesc,revDescIndx);
+ break;
+ }
+ case 2:
+ {
+ mesh1D=buildDescendingConnectivity(desc,descIndex,revDesc,revDescIndx);
+ break;
+ }
+ default:
+ {
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::computeNeighborsOfNodes : Mesh dimension supported are [3,2] !");
+ }
+ }
+ return mesh1D;
+}
+
/*!
* \b WARNING this method do the assumption that connectivity lies on the coordinates set.
* For speed reasons no check of this will be done. This method calls
* The number of tuples is equal to the last values in \b neighborsIndx.
* \param [out] neighborsIdx is an array of size this->getNumberOfCells()+1 newly allocated and should
* be dealt by the caller. This arrays allow to use the first output parameter \b neighbors.
+ *
+ * \sa MEDCouplingUMesh::computeEnlargedNeighborsOfNodes
*/
void MEDCouplingUMesh::computeNeighborsOfNodes(DataArrayInt *&neighbors, DataArrayInt *&neighborsIdx) const
{
checkFullyDefined();
int mdim(getMeshDimension()),nbNodes(getNumberOfNodes());
MCAuto<DataArrayInt> desc(DataArrayInt::New()),descIndx(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescIndx(DataArrayInt::New());
- MCAuto<MEDCouplingUMesh> mesh1D;
+ MCConstAuto<MEDCouplingUMesh> mesh1D;
switch(mdim)
{
case 3:
}
case 1:
{
- mesh1D=const_cast<MEDCouplingUMesh *>(this);
- mesh1D->incrRef();
+ mesh1D.takeRef(this);
break;
}
default:
neighborsIdx=descIndx.retn();
}
+/*!
+ * Computes enlarged neighbors for each nodes in \a this. The behavior of this method is close to MEDCouplingUMesh::computeNeighborsOfNodes except that the neighborhood of each node is wider here.
+ * A node j is considered to be in the neighborhood of i if and only if there is a cell in \a this containing in its nodal connectivity both i and j.
+ * This method is useful to find ghost cells of a part of a mesh with a code based on fields on nodes.
+ *
+ * \sa MEDCouplingUMesh::computeNeighborsOfNodes
+ */
+void MEDCouplingUMesh::computeEnlargedNeighborsOfNodes(MCAuto<DataArrayInt> &neighbors, MCAuto<DataArrayInt>& neighborsIdx) const
+{
+ checkFullyDefined();
+ int nbOfNodes(getNumberOfNodes());
+ const int *conn(_nodal_connec->begin()),*connIndex(_nodal_connec_index->begin());
+ int nbOfCells(getNumberOfCells());
+ std::vector< std::set<int> > st0(nbOfNodes);
+ for(int eltId=0;eltId<nbOfCells;eltId++)
+ {
+ const int *strtNdlConnOfCurCell(conn+connIndex[eltId]+1),*endNdlConnOfCurCell(conn+connIndex[eltId+1]);
+ std::set<int> s(strtNdlConnOfCurCell,endNdlConnOfCurCell); s.erase(-1); //for polyhedrons
+ for(std::set<int>::const_iterator iter2=s.begin();iter2!=s.end();iter2++)
+ st0[*iter2].insert(s.begin(),s.end());
+ }
+ neighborsIdx=DataArrayInt::New(); neighborsIdx->alloc(nbOfNodes+1,1); neighborsIdx->setIJ(0,0,0);
+ {
+ int *neighIdx(neighborsIdx->getPointer());
+ for(std::vector< std::set<int> >::const_iterator it=st0.begin();it!=st0.end();it++,neighIdx++)
+ neighIdx[1]=neighIdx[0]+(*it).size()-1;
+ }
+ neighbors=DataArrayInt::New(); neighbors->alloc(neighborsIdx->back(),1);
+ {
+ const int *neighIdx(neighborsIdx->begin());
+ int *neigh(neighbors->getPointer()),nodeId(0);
+ for(std::vector< std::set<int> >::const_iterator it=st0.begin();it!=st0.end();it++,neighIdx++,nodeId++)
+ {
+ std::set<int> s(*it); s.erase(nodeId);
+ std::copy(s.begin(),s.end(),neigh+*neighIdx);
+ }
+ }
+}
+
/*!
* Converts specified cells to either polygons (if \a this is a 2D mesh) or
* polyhedrons (if \a this is a 3D mesh). The cells to convert are specified by an
int nbOfCells(getNumberOfCells());
if(dim==2)
{
- const int *connIndex=_nodal_connec_index->getConstPointer();
+ const int *connIndex=_nodal_connec_index->begin();
int *conn=_nodal_connec->getPointer();
for(const int *iter=cellIdsToConvertBg;iter!=cellIdsToConvertEnd;iter++)
{
oss << ", whereas other mesh dimension is set equal to " << otherOnSameCoordsThanThis.getMeshDimension() << " !";
throw INTERP_KERNEL::Exception(oss.str());
}
- int nbOfCellsToModify=(int)std::distance(cellIdsBg,cellIdsEnd);
+ std::size_t nbOfCellsToModify(std::distance(cellIdsBg,cellIdsEnd));
if(nbOfCellsToModify!=otherOnSameCoordsThanThis.getNumberOfCells())
{
std::ostringstream oss; oss << "MEDCouplingUMesh::setPartOfMySelf : cells ids length (" << nbOfCellsToModify << ") do not match the number of cells of other mesh (" << otherOnSameCoordsThanThis.getNumberOfCells() << ") !";
throw INTERP_KERNEL::Exception(oss.str());
}
- int nbOfCells=getNumberOfCells();
- bool easyAssign=true;
- const int *connI=_nodal_connec_index->getConstPointer();
- const int *connIOther=otherOnSameCoordsThanThis._nodal_connec_index->getConstPointer();
+ std::size_t nbOfCells(getNumberOfCells());
+ bool easyAssign(true);
+ const int *connI(_nodal_connec_index->begin());
+ const int *connIOther=otherOnSameCoordsThanThis._nodal_connec_index->begin();
for(const int *it=cellIdsBg;it!=cellIdsEnd && easyAssign;it++,connIOther++)
{
- if(*it>=0 && *it<nbOfCells)
+ if(*it>=0 && *it<(int)nbOfCells)
{
easyAssign=(connIOther[1]-connIOther[0])==(connI[*it+1]-connI[*it]);
}
throw INTERP_KERNEL::Exception(oss.str());
}
int nbOfCellsToModify=DataArray::GetNumberOfItemGivenBESRelative(start,end,step,"MEDCouplingUMesh::setPartOfMySelfSlice : ");
- if(nbOfCellsToModify!=otherOnSameCoordsThanThis.getNumberOfCells())
+ if(nbOfCellsToModify!=(int)otherOnSameCoordsThanThis.getNumberOfCells())
{
std::ostringstream oss; oss << "MEDCouplingUMesh::setPartOfMySelfSlice : cells ids length (" << nbOfCellsToModify << ") do not match the number of cells of other mesh (" << otherOnSameCoordsThanThis.getNumberOfCells() << ") !";
throw INTERP_KERNEL::Exception(oss.str());
DAInt neighIInit00(tmp11);
// Neighbor information of the mesh WITH the crack (some neighbors are removed):
DataArrayInt *idsTmp=0;
- bool b=m01->areCellsIncludedIn(&otherDimM1OnSameCoords,2,idsTmp);
+ m01->areCellsIncludedIn(&otherDimM1OnSameCoords,2,idsTmp);
DAInt ids(idsTmp);
// In the neighbor information remove the connection between high dimension cells and its low level constituents which are part
// of the frontier given in parameter (i.e. the cells of low dimension from the group delimiting the crack):
// Connex zone without the crack (to compute the next seed really)
int dnu;
DAInt connexCheck = MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeed(&seed, &seed+1, neighInit00,neighIInit00, -1, dnu);
- int cnt = 0;
+ std::size_t cnt(0);
for (int * ptr = connexCheck->getPointer(); cnt < connexCheck->getNumberOfTuples(); ptr++, cnt++)
hitCells->setIJ(*ptr,0,1);
// Connex zone WITH the crack (to identify cells lying on either part of the crack)
* \return INTERP_KERNEL::NormalizedCellType - enumeration item describing the cell type.
* \throw If \a cellId is invalid. Valid range is [0, \a this->getNumberOfCells() ).
*/
-INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::getTypeOfCell(int cellId) const
+INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::getTypeOfCell(std::size_t cellId) const
{
- const int *ptI=_nodal_connec_index->getConstPointer();
- const int *pt=_nodal_connec->getConstPointer();
- if(cellId>=0 && cellId<(int)_nodal_connec_index->getNbOfElems()-1)
+ const int *ptI(_nodal_connec_index->begin()),*pt(_nodal_connec->begin());
+ if(cellId<_nodal_connec_index->getNbOfElems()-1)
return (INTERP_KERNEL::NormalizedCellType) pt[ptI[cellId]];
else
{
/*!
* Returns nb of cells having the geometric type \a type. No throw if no cells in \a this has the geometric type \a type.
*/
-int MEDCouplingUMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
+std::size_t MEDCouplingUMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
{
- const int *ptI=_nodal_connec_index->getConstPointer();
- const int *pt=_nodal_connec->getConstPointer();
- int nbOfCells=getNumberOfCells();
- int ret=0;
- for(int i=0;i<nbOfCells;i++)
+ const int *ptI(_nodal_connec_index->begin()),*pt(_nodal_connec->begin());
+ std::size_t nbOfCells(getNumberOfCells()),ret(0);
+ for(std::size_t i=0;i<nbOfCells;i++)
if((INTERP_KERNEL::NormalizedCellType) pt[ptI[i]]==type)
ret++;
return ret;
* cleared before the appending.
* \throw If \a cellId is invalid. Valid range is [0, \a this->getNumberOfCells() ).
*/
-void MEDCouplingUMesh::getNodeIdsOfCell(int cellId, std::vector<int>& conn) const
+void MEDCouplingUMesh::getNodeIdsOfCell(std::size_t cellId, std::vector<int>& conn) const
{
- const int *ptI=_nodal_connec_index->getConstPointer();
- const int *pt=_nodal_connec->getConstPointer();
+ const int *ptI(_nodal_connec_index->begin()),*pt(_nodal_connec->begin());
for(const int *w=pt+ptI[cellId]+1;w!=pt+ptI[cellId+1];w++)
if(*w>=0)
conn.push_back(*w);
* Copy constructor. If 'deepCopy' is false \a this is a shallow copy of other.
* If 'deeCpy' is true all arrays (coordinates and connectivities) are deeply copied.
*/
-MEDCouplingUMesh::MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCopy):MEDCouplingPointSet(other,deepCopy),_mesh_dim(other._mesh_dim),
+MEDCouplingUMesh::MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCpy):MEDCouplingPointSet(other,deepCpy),_mesh_dim(other._mesh_dim),
_nodal_connec(0),_nodal_connec_index(0),
_types(other._types)
{
if(other._nodal_connec)
- _nodal_connec=other._nodal_connec->performCopyOrIncrRef(deepCopy);
+ _nodal_connec=other._nodal_connec->performCopyOrIncrRef(deepCpy);
if(other._nodal_connec_index)
- _nodal_connec_index=other._nodal_connec_index->performCopyOrIncrRef(deepCopy);
+ _nodal_connec_index=other._nodal_connec_index->performCopyOrIncrRef(deepCpy);
}
MEDCouplingUMesh::~MEDCouplingUMesh()
* \return int - the number of cells in \a this mesh.
* \throw If the nodal connectivity of cells is not defined.
*/
-int MEDCouplingUMesh::getNumberOfCells() const
+std::size_t MEDCouplingUMesh::getNumberOfCells() const
{
if(_nodal_connec_index)
return _nodal_connec_index->getNumberOfTuples()-1;
MCAuto<MEDCouplingFieldDouble> f=buildDirectionVectorField();
const double *fPtr=f->getArray()->getConstPointer();
double tmp[3];
- for(int i=0;i<getNumberOfCells();i++)
+ for(std::size_t i=0;i<getNumberOfCells();i++)
{
const double *tmp1=fPtr+3*i;
tmp[0]=tmp1[1]*v[2]-tmp1[2]*v[1];
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoints works only for spaceDim=meshDim+1 !");
if(meshDim!=2 && meshDim!=1)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoints : only mesh dimension 2 and 1 are implemented !");
- if(pts->getNumberOfComponents()!=spaceDim)
+ if((int)pts->getNumberOfComponents()!=spaceDim)
{
std::ostringstream oss; oss << "MEDCouplingUMesh::distanceToPoints : input pts DataArrayDouble has " << pts->getNumberOfComponents() << " components whereas it should be equal to " << spaceDim << " (mesh spaceDimension) !";
throw INTERP_KERNEL::Exception(oss.str());
updateTime();
}
+/*!
+ * This method invert orientation of all cells in \a this.
+ * After calling this method the absolute value of measure of cells in \a this are the same than before calling.
+ * This method only operates on the connectivity so coordinates are not touched at all.
+ */
+void MEDCouplingUMesh::invertOrientationOfAllCells()
+{
+ checkConnectivityFullyDefined();
+ std::set<INTERP_KERNEL::NormalizedCellType> gts(getAllGeoTypes());
+ int *conn(_nodal_connec->getPointer());
+ const int *conni(_nodal_connec_index->begin());
+ for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator gt=gts.begin();gt!=gts.end();gt++)
+ {
+ INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::OrientationInverter> oi(INTERP_KERNEL::OrientationInverter::BuildInstanceFrom(*gt));
+ MCAuto<DataArrayInt> cwt(giveCellsWithType(*gt));
+ for(const int *it=cwt->begin();it!=cwt->end();it++)
+ oi->operate(conn+conni[*it]+1,conn+conni[*it+1]);
+ }
+ updateTime();
+}
+
/*!
* Finds and fixes incorrectly oriented linear extruded volumes (INTERP_KERNEL::NORM_HEXA8,
* INTERP_KERNEL::NORM_PENTA6, INTERP_KERNEL::NORM_HEXGP12 etc) to respect the MED convention
DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arcDetEps) const
{
checkFullyDefined();
+ INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(arcDetEps);
+
int spaceDim(getSpaceDimension()),mDim(getMeshDimension()),nbOfCells(getNumberOfCells());
if(spaceDim!=2 || mDim!=2)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic : This method should be applied on mesh with mesh dimension equal to 2 and space dimension also equal to 2!");
{
const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[*connI]));
int sz(connI[1]-connI[0]-1);
- INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=arcDetEps;
std::vector<INTERP_KERNEL::Node *> nodes(sz);
INTERP_KERNEL::QuadraticPolygon *pol(0);
for(int j=0;j<sz;j++)
int spaceDim(getSpaceDimension()),mDim(getMeshDimension()),nbOfCells(getNumberOfCells());
if(spaceDim!=2 || mDim!=1)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic : This method should be applied on mesh with mesh dimension equal to 1 and space dimension also equal to 2!");
+ INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(arcDetEps);
MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim);
double *bbox(ret->getPointer());
const double *coords(_coords->begin());
{
const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[*connI]));
int sz(connI[1]-connI[0]-1);
- INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=arcDetEps;
std::vector<INTERP_KERNEL::Node *> nodes(sz);
INTERP_KERNEL::Edge *edge(0);
for(int j=0;j<sz;j++)
for(int i=0;i<nbOfCells;i++,nodalI++,retPtr+=4)
{
double matrix[16]={0,0,0,1,0,0,0,1,0,0,0,1,1,1,1,0},matrix2[16];
- if(nodalI[1]-nodalI[0]>=3)
+ if(nodalI[1]-nodalI[0]>=4)
{
+ double aa[3]={coor[nodal[nodalI[0]+1+1]*3+0]-coor[nodal[nodalI[0]+1+0]*3+0],
+ coor[nodal[nodalI[0]+1+1]*3+1]-coor[nodal[nodalI[0]+1+0]*3+1],
+ coor[nodal[nodalI[0]+1+1]*3+2]-coor[nodal[nodalI[0]+1+0]*3+2]}
+ ,bb[3]={coor[nodal[nodalI[0]+1+2]*3+0]-coor[nodal[nodalI[0]+1+0]*3+0],
+ coor[nodal[nodalI[0]+1+2]*3+1]-coor[nodal[nodalI[0]+1+0]*3+1],
+ coor[nodal[nodalI[0]+1+2]*3+2]-coor[nodal[nodalI[0]+1+0]*3+2]};
+ double cc[3]={aa[1]*bb[2]-aa[2]*bb[1],aa[2]*bb[0]-aa[0]*bb[2],aa[0]*bb[1]-aa[1]*bb[0]};
for(int j=0;j<3;j++)
{
int nodeId(nodal[nodalI[0]+1+j]);
throw INTERP_KERNEL::Exception(oss.str());
}
}
+ if(sqrt(cc[0]*cc[0]+cc[1]*cc[1]+cc[2]*cc[2])>1e-7)
+ {
+ INTERP_KERNEL::inverseMatrix(matrix,4,matrix2);
+ retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15];
+ }
+ else
+ {
+ if(nodalI[1]-nodalI[0]==4)
+ {
+ std::ostringstream oss; oss << "MEDCouplingUMesh::computePlaneEquationOf3DFaces : cell" << i << " : Presence of The 3 colinear points !";
+ throw INTERP_KERNEL::Exception(oss.str());
+ }
+ //
+ double dd[3]={0.,0.,0.};
+ for(int offset=nodalI[0]+1;offset<nodalI[1];offset++)
+ std::transform(coor+3*nodal[offset],coor+3*(nodal[offset]+1),dd,dd,std::plus<double>());
+ int nbOfNodesInCell(nodalI[1]-nodalI[0]-1);
+ std::transform(dd,dd+3,dd,std::bind2nd(std::multiplies<double>(),1./(double)nbOfNodesInCell));
+ std::copy(dd,dd+3,matrix+4*2);
+ INTERP_KERNEL::inverseMatrix(matrix,4,matrix2);
+ retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15];
+ }
}
else
{
std::ostringstream oss; oss << "MEDCouplingUMesh::computePlaneEquationOf3DFaces : invalid 2D cell #" << i << " ! Must be constitued by more than 3 nodes !";
throw INTERP_KERNEL::Exception(oss.str());
}
- INTERP_KERNEL::inverseMatrix(matrix,4,matrix2);
- retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15];
}
return ret.retn();
}
previousArrI=*arrIPtr;
*arrIPtr=(int)arrOut.size();
}
- if(arr->getNumberOfTuples()==(int)arrOut.size())
+ if(arr->getNumberOfTuples()==arrOut.size())
return false;
arr->alloc((int)arrOut.size(),1);
std::copy(arrOut.begin(),arrOut.end(),arr->getPointer());
