-// Copyright (C) 2007-2013 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2014 CEA/DEN, EDF R&D
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
void MEDCouplingUMesh::checkCoherency() const
{
if(_mesh_dim<-1)
- throw INTERP_KERNEL::Exception("No mesh dimension specified !");
+ throw INTERP_KERNEL::Exception("No mesh dimension specified !");
if(_mesh_dim!=-1)
MEDCouplingPointSet::checkCoherency();
for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator iter=_types.begin();iter!=_types.end();iter++)
*
* \param [in] nbOfCells - estimation of the number of cell \a this mesh will contain.
*
+ * \if ENABLE_EXAMPLES
* \ref medcouplingcppexamplesUmeshStdBuild1 "Here is a C++ example".<br>
* \ref medcouplingpyexamplesUmeshStdBuild1 "Here is a Python example".
+ * \endif
*/
void MEDCouplingUMesh::allocateCells(int nbOfCells)
{
* \param [in] size - number of nodes constituting this cell.
* \param [in] nodalConnOfCell - the connectivity of the cell to add.
*
+ * \if ENABLE_EXAMPLES
* \ref medcouplingcppexamplesUmeshStdBuild1 "Here is a C++ example".<br>
* \ref medcouplingpyexamplesUmeshStdBuild1 "Here is a Python example".
+ * \endif
*/
void MEDCouplingUMesh::insertNextCell(INTERP_KERNEL::NormalizedCellType type, int size, const int *nodalConnOfCell)
{
* Compacts data arrays to release unused memory. This method is to be called after
* finishing cell insertion using \a this->insertNextCell().
*
+ * \if ENABLE_EXAMPLES
* \ref medcouplingcppexamplesUmeshStdBuild1 "Here is a C++ example".<br>
* \ref medcouplingpyexamplesUmeshStdBuild1 "Here is a Python example".
+ * \endif
*/
void MEDCouplingUMesh::finishInsertingCells()
{
*/
void MEDCouplingUMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const
{
- MEDCouplingPointSet::checkFastEquivalWith(other,prec);
- const MEDCouplingUMesh *otherC=dynamic_cast<const MEDCouplingUMesh *>(other);
+ MEDCouplingPointSet::checkFastEquivalWith(other,prec);
+ const MEDCouplingUMesh *otherC=dynamic_cast<const MEDCouplingUMesh *>(other);
if(!otherC)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::checkFastEquivalWith : Two meshes are not not unstructured !");
}
* \throw If the coordinates array is not set.
* \throw If the nodal connectivity of cells is not defined.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_getReverseNodalConnectivity "Here is a C++ example".<br>
* \ref py_mcumesh_getReverseNodalConnectivity "Here is a Python example".
+ * \endif
*/
void MEDCouplingUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const
{
* \throw If \a desc == NULL || \a descIndx == NULL || \a revDesc == NULL || \a
* revDescIndx == NULL.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_buildDescendingConnectivity "Here is a C++ example".<br>
* \ref py_mcumesh_buildDescendingConnectivity "Here is a Python example".
+ * \endif
* \sa buildDescendingConnectivity2()
*/
MEDCouplingUMesh *MEDCouplingUMesh::buildDescendingConnectivity(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx) const
* \throw If \a desc == NULL || \a descIndx == NULL || \a revDesc == NULL || \a
* revDescIndx == NULL.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_buildDescendingConnectivity2 "Here is a C++ example".<br>
* \ref py_mcumesh_buildDescendingConnectivity2 "Here is a Python example".
+ * \endif
* \sa buildDescendingConnectivity()
*/
MEDCouplingUMesh *MEDCouplingUMesh::buildDescendingConnectivity2(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx) const
* \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 MEDCouplingUMesh::buildDescendingConnectivity to compute the result.
* This method lists cell by cell in \b this which are its neighbors. To compute the result only connectivities are considered.
- * The a cell with id 'cellId' its neighbors are neighbors[neighborsIndx[cellId]:neighborsIndx[cellId+1]].
+ * The neighbor cells of cell having id 'cellId' are neighbors[neighborsIndx[cellId]:neighborsIndx[cellId+1]].
*
* \param [out] neighbors is an array storing all the neighbors of all cells in \b this. This array is newly allocated and should be dealt by the caller. \b neighborsIndx 2nd output
* parameter allows to select the right part in this array. The number of tuples is equal to the last values in \b neighborsIndx.
* excluding a set of meshdim-1 cells in input descending connectivity.
* Typically \b desc, \b descIndx, \b revDesc and \b revDescIndx input params are the result of MEDCouplingUMesh::buildDescendingConnectivity.
* This method lists cell by cell in \b this which are its neighbors. To compute the result only connectivities are considered.
- * The a cell with id 'cellId' its neighbors are neighbors[neighborsIndx[cellId]:neighborsIndx[cellId+1]].
+ * The neighbor cells of cell having id 'cellId' are neighbors[neighborsIndx[cellId]:neighborsIndx[cellId+1]].
*
* \param [in] desc descending connectivity array.
* \param [in] descIndx descending connectivity index array used to walk through \b desc.
* \param [out] neighborsIndx 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.
*/
void MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(const DataArrayInt *desc, const DataArrayInt *descIndx, const DataArrayInt *revDesc, const DataArrayInt *revDescIndx,
- DataArrayInt *&neighbors, DataArrayInt *&neighborsIndx) throw(INTERP_KERNEL::Exception)
+ DataArrayInt *&neighbors, DataArrayInt *&neighborsIndx)
{
if(!desc || !descIndx || !revDesc || !revDescIndx)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ComputeNeighborsOfCellsAdv some input array is empty !");
neighborsIndx=out1.retn();
}
+/*!
+ * \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 MEDCouplingUMesh::buildDescendingConnectivity to compute the result.
+ * This method lists node by node in \b this which are its neighbors. To compute the result only connectivities are considered.
+ * The neighbor nodes of node having id 'nodeId' are neighbors[neighborsIndx[cellId]:neighborsIndx[cellId+1]].
+ *
+ * \param [out] neighbors is an array storing all the neighbors of all nodes in \b this. This array is newly allocated and should be dealt by the caller. \b neighborsIndx 2nd output
+ * parameter allows to select the right part in this array. The number of tuples is equal to the last values in \b neighborsIndx.
+ * \param [out] neighborsIndx 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.
+ */
+void MEDCouplingUMesh::computeNeighborsOfNodes(DataArrayInt *&neighbors, DataArrayInt *&neighborsIdx) const
+{
+ checkFullyDefined();
+ int mdim(getMeshDimension()),nbNodes(getNumberOfNodes());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> desc(DataArrayInt::New()),descIndx(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescIndx(DataArrayInt::New());
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> mesh1D;
+ switch(mdim)
+ {
+ case 3:
+ {
+ mesh1D=explode3DMeshTo1D(desc,descIndx,revDesc,revDescIndx);
+ break;
+ }
+ case 2:
+ {
+ mesh1D=buildDescendingConnectivity(desc,descIndx,revDesc,revDescIndx);
+ break;
+ }
+ case 1:
+ {
+ mesh1D=const_cast<MEDCouplingUMesh *>(this);
+ mesh1D->incrRef();
+ break;
+ }
+ default:
+ {
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::computeNeighborsOfNodes : Mesh dimension supported are [3,2,1] !");
+ }
+ }
+ desc=DataArrayInt::New(); descIndx=DataArrayInt::New(); revDesc=0; revDescIndx=0;
+ mesh1D->getReverseNodalConnectivity(desc,descIndx);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret0(DataArrayInt::New());
+ ret0->alloc(desc->getNumberOfTuples(),1);
+ int *r0Pt(ret0->getPointer());
+ const int *c1DPtr(mesh1D->getNodalConnectivity()->begin()),*rn(desc->begin()),*rni(descIndx->begin());
+ for(int i=0;i<nbNodes;i++,rni++)
+ {
+ for(const int *oneDCellIt=rn+rni[0];oneDCellIt!=rn+rni[1];oneDCellIt++)
+ *r0Pt++=c1DPtr[3*(*oneDCellIt)+1]==i?c1DPtr[3*(*oneDCellIt)+2]:c1DPtr[3*(*oneDCellIt)+1];
+ }
+ neighbors=ret0.retn();
+ neighborsIdx=descIndx.retn();
+}
+
/// @cond INTERNAL
/*!
* \throw If the nodal connectivity of cells is node defined.
* \throw If dimension of \a this mesh is not either 2 or 3.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_convertToPolyTypes "Here is a C++ example".<br>
* \ref py_mcumesh_convertToPolyTypes "Here is a Python example".
+ * \endif
*/
void MEDCouplingUMesh::convertToPolyTypes(const int *cellIdsToConvertBg, const int *cellIdsToConvertEnd)
{
int dim=getMeshDimension();
if(dim<2 || dim>3)
throw INTERP_KERNEL::Exception("Invalid mesh dimension : must be 2 or 3 !");
- int nbOfCells=getNumberOfCells();
+ int nbOfCells(getNumberOfCells());
if(dim==2)
{
const int *connIndex=_nodal_connec_index->getConstPointer();
}
else
{
- int *connIndex=_nodal_connec_index->getPointer();
- int connIndexLgth=_nodal_connec_index->getNbOfElems();
- const int *connOld=_nodal_connec->getConstPointer();
- int connOldLgth=_nodal_connec->getNbOfElems();
- std::vector<int> connNew(connOld,connOld+connOldLgth);
+ int *connIndex(_nodal_connec_index->getPointer());
+ const int *connOld(_nodal_connec->getConstPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> connNew(DataArrayInt::New()),connNewI(DataArrayInt::New()); connNew->alloc(0,1); connNewI->alloc(1,1); connNewI->setIJ(0,0,0);
+ std::vector<bool> toBeDone(nbOfCells,false);
for(const int *iter=cellIdsToConvertBg;iter!=cellIdsToConvertEnd;iter++)
{
if(*iter>=0 && *iter<nbOfCells)
+ toBeDone[*iter]=true;
+ else
+ {
+ std::ostringstream oss; oss << "MEDCouplingUMesh::convertToPolyTypes : On rank #" << std::distance(cellIdsToConvertBg,iter) << " value is " << *iter << " which is not";
+ oss << " in range [0," << nbOfCells << ") !";
+ throw INTERP_KERNEL::Exception(oss.str().c_str());
+ }
+ }
+ for(int cellId=0;cellId<nbOfCells;cellId++)
+ {
+ int pos(connIndex[cellId]),posP1(connIndex[cellId+1]);
+ int lgthOld(posP1-pos-1);
+ if(toBeDone[cellId])
{
- int pos=connIndex[*iter];
- int posP1=connIndex[(*iter)+1];
- int lgthOld=posP1-pos-1;
- const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connNew[pos]);
- connNew[pos]=INTERP_KERNEL::NORM_POLYHED;
- unsigned nbOfFaces=cm.getNumberOfSons2(&connNew[pos+1],lgthOld);
- int *tmp=new int[nbOfFaces*lgthOld];
- int *work=tmp;
- for(int j=0;j<(int)nbOfFaces;j++)
+ const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connOld[pos]);
+ unsigned nbOfFaces(cm.getNumberOfSons2(connOld+pos+1,lgthOld));
+ int *tmp(new int[nbOfFaces*lgthOld+1]);
+ int *work=tmp; *work++=INTERP_KERNEL::NORM_POLYHED;
+ for(unsigned j=0;j<nbOfFaces;j++)
{
INTERP_KERNEL::NormalizedCellType type;
- unsigned offset=cm.fillSonCellNodalConnectivity2(j,&connNew[pos+1],lgthOld,work,type);
+ unsigned offset=cm.fillSonCellNodalConnectivity2(j,connOld+pos+1,lgthOld,work,type);
work+=offset;
*work++=-1;
}
- std::size_t newLgth=std::distance(tmp,work)-1;
- std::size_t delta=newLgth-lgthOld;
- std::transform(connIndex+(*iter)+1,connIndex+connIndexLgth,connIndex+(*iter)+1,std::bind2nd(std::plus<int>(),delta));
- connNew.insert(connNew.begin()+posP1,tmp+lgthOld,tmp+newLgth);
- std::copy(tmp,tmp+lgthOld,connNew.begin()+pos+1);
+ std::size_t newLgth(std::distance(tmp,work)-1);//-1 for last -1
+ connNew->pushBackValsSilent(tmp,tmp+newLgth);
+ connNewI->pushBackSilent(connNewI->back()+(int)newLgth);
delete [] tmp;
}
else
{
- std::ostringstream oss; oss << "MEDCouplingUMesh::convertToPolyTypes : On rank #" << std::distance(cellIdsToConvertBg,iter) << " value is " << *iter << " which is not";
- oss << " in range [0," << nbOfCells << ") !";
- throw INTERP_KERNEL::Exception(oss.str().c_str());
+ connNew->pushBackValsSilent(connOld+pos,connOld+posP1);
+ connNewI->pushBackSilent(connNewI->back()+posP1-pos);
}
}
- _nodal_connec->alloc((int)connNew.size(),1);
- int *newConnPtr=_nodal_connec->getPointer();
- std::copy(connNew.begin(),connNew.end(),newConnPtr);
+ setConnectivity(connNew,connNewI,false);//false because computeTypes called just behind.
}
computeTypes();
}
* \throw If \a this mesh contains polyhedrons with the valid connectivity.
* \throw If \a this mesh contains polyhedrons with odd number of nodes.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_arePolyhedronsNotCorrectlyOriented "Here is a C++ example".<br>
* \ref py_mcumesh_arePolyhedronsNotCorrectlyOriented "Here is a Python example".
+ * \endif
*/
void MEDCouplingUMesh::convertExtrudedPolyhedra()
{
if(cm.isDynamic())
{
switch(cm.getDimension())
- {
+ {
case 2:
{
INTERP_KERNEL::AutoPtr<int> tmp=new int[lgthOfCurCell-1];
newType=(lgthOfCurCell==3)?INTERP_KERNEL::NORM_SEG2:INTERP_KERNEL::NORM_POLYL;
break;
}
- }
+ }
ret=ret || (newType!=type);
conn[newPos]=newType;
newPos+=newLgth+1;
* \throw If the nodal connectivity of cells is not defined.
* \throw If the nodal connectivity includes an invalid id.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_getNodeIdsInUse "Here is a C++ example".<br>
* \ref py_mcumesh_getNodeIdsInUse "Here is a Python example".
+ * \endif
* \sa computeNodeIdsAlg()
*/
DataArrayInt *MEDCouplingUMesh::getNodeIdsInUse(int& nbrOfNodesInUse) const
* \throw If the nodal connectivity of cells is not defined.
* \throw If the nodal connectivity includes an invalid id.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_zipCoordsTraducer "Here is a C++ example".<br>
* \ref py_mcumesh_zipCoordsTraducer "Here is a Python example".
+ * \endif
*/
DataArrayInt *MEDCouplingUMesh::zipCoordsTraducer()
{
int MEDCouplingUMesh::AreCellsEqual(const int *conn, const int *connI, int cell1, int cell2, int compType)
{
switch(compType)
- {
+ {
case 0:
return AreCellsEqual0(conn,connI,cell1,cell2);
case 1:
return AreCellsEqual3(conn,connI,cell1,cell2);
case 7:
return AreCellsEqual7(conn,connI,cell1,cell2);
- }
+ }
throw INTERP_KERNEL::Exception("Unknown comparison asked ! Must be in 0,1,2,3 or 7.");
}
else
return 0;
}
-
+
return work!=tmp+sz1?1:0;
}
else
}
void MEDCouplingUMesh::FindCommonCellsAlg(int compType, int startCellId, const DataArrayInt *nodal, const DataArrayInt *nodalI, const DataArrayInt *revNodal, const DataArrayInt *revNodalI,
- DataArrayInt *& commonCellsArr, DataArrayInt *& commonCellsIArr) throw(INTERP_KERNEL::Exception)
+ DataArrayInt *& commonCellsArr, DataArrayInt *& commonCellsIArr)
{
MEDCouplingAutoRefCountObjectPtr<DataArrayInt> commonCells=DataArrayInt::New(),commonCellsI=DataArrayInt::New(); commonCells->alloc(0,1);
int nbOfCells=nodalI->getNumberOfTuples()-1;
* \return bool - \c true if all cells of \a other mesh are present in the \a this
* mesh.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_areCellsIncludedIn "Here is a C++ example".<br>
* \ref py_mcumesh_areCellsIncludedIn "Here is a Python example".
+ * \endif
* \sa checkDeepEquivalOnSameNodesWith()
* \sa checkGeoEquivalWith()
*/
* \throw If the nodal connectivity of cells is not defined.
* \throw If any cell id in the array \a begin is not valid.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_buildPartOfMySelf "Here is a C++ example".<br>
* \ref py_mcumesh_buildPartOfMySelf "Here is a Python example".
+ * \endif
*/
MEDCouplingPointSet *MEDCouplingUMesh::buildPartOfMySelf(const int *begin, const int *end, bool keepCoords) const
{
* \throw If the nodal connectivity of cells is not defined.
* \throw If any node id in \a begin is not valid.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_buildFacePartOfMySelfNode "Here is a C++ example".<br>
* \ref py_mcumesh_buildFacePartOfMySelfNode "Here is a Python example".
+ * \endif
*/
MEDCouplingPointSet *MEDCouplingUMesh::buildFacePartOfMySelfNode(const int *begin, const int *end, bool fullyIn) const
{
* \throw If the coordinates array is not set.
* \throw If the nodal connectivity of cells is not defined.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_buildBoundaryMesh "Here is a C++ example".<br>
* \ref py_mcumesh_buildBoundaryMesh "Here is a Python example".
+ * \endif
*/
MEDCouplingPointSet *MEDCouplingUMesh::buildBoundaryMesh(bool keepCoords) const
{
* \throw If the coordinates array is not set.
* \throw If the nodal connectivity of cells is node defined.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_findBoundaryNodes "Here is a C++ example".<br>
* \ref py_mcumesh_findBoundaryNodes "Here is a Python example".
+ * \endif
*/
DataArrayInt *MEDCouplingUMesh::findBoundaryNodes() const
{
* \warning This method modifies param \b otherDimM1OnSameCoords (for speed reasons).
*/
void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1OnSameCoords, DataArrayInt *& nodeIdsToDuplicate,
- DataArrayInt *& cellIdsNeededToBeRenum, DataArrayInt *& cellIdsNotModified) const throw(INTERP_KERNEL::Exception)
+ DataArrayInt *& cellIdsNeededToBeRenum, DataArrayInt *& cellIdsNotModified) const
{
checkFullyDefined();
otherDimM1OnSameCoords.checkFullyDefined();
* See \ref MEDCouplingArrayRenumbering for more info on renumbering modes.
* \throw If the nodal connectivity of cells is not defined.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_renumberNodesInConn "Here is a C++ example".<br>
* \ref py_mcumesh_renumberNodesInConn "Here is a Python example".
+ * \endif
*/
void MEDCouplingUMesh::renumberNodesInConn(const int *newNodeNumbersO2N)
{
* \throw If the coordinates array is not set.
* \throw If the nodal connectivity of cells is not defined.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_getCellsInBoundingBox "Here is a C++ example".<br>
* \ref py_mcumesh_getCellsInBoundingBox "Here is a Python example".
+ * \endif
*/
DataArrayInt *MEDCouplingUMesh::getCellsInBoundingBox(const double *bbox, double eps) const
{
*/
DataArrayInt *MEDCouplingUMesh::giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
{
-
+
MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
ret->alloc(0,1);
checkConnectivityFullyDefined();
* 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),
- _nodal_connec(0),_nodal_connec_index(0),
- _types(other._types)
+ _nodal_connec(0),_nodal_connec_index(0),
+ _types(other._types)
{
if(other._nodal_connec)
_nodal_connec=other._nodal_connec->performCpy(deepCopy);
* \return DataArrayDouble * - a new instance of DataArrayDouble. The caller is to
* delete this array using decrRef() as it is no more needed.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_getPartMeasureField "Here is a C++ example".<br>
* \ref py_mcumesh_getPartMeasureField "Here is a Python example".
+ * \endif
* \sa getMeasureField()
*/
DataArrayDouble *MEDCouplingUMesh::getPartMeasureField(bool isAbs, const int *begin, const int *end) const
* \throw If the mesh and space dimension is not as specified above.
* \sa buildOrthogonalField()
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_buildPartOrthogonalField "Here is a C++ example".<br>
* \ref py_mcumesh_buildPartOrthogonalField "Here is a Python example".
+ * \endif
*/
MEDCouplingFieldDouble *MEDCouplingUMesh::buildPartOrthogonalField(const int *begin, const int *end) const
{
*/
MEDCouplingFieldDouble *MEDCouplingUMesh::buildDirectionVectorField() const
{
- if(getMeshDimension()!=1)
+ if(getMeshDimension()!=1)
throw INTERP_KERNEL::Exception("Expected a umesh with meshDim == 1 for buildDirectionVectorField !");
- if(_types.size()!=1 || *(_types.begin())!=INTERP_KERNEL::NORM_SEG2)
- throw INTERP_KERNEL::Exception("Expected a umesh with only NORM_SEG2 type of elements for buildDirectionVectorField !");
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME);
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> array=DataArrayDouble::New();
- int nbOfCells=getNumberOfCells();
- int spaceDim=getSpaceDimension();
- array->alloc(nbOfCells,spaceDim);
- double *pt=array->getPointer();
- const double *coo=getCoords()->getConstPointer();
- std::vector<int> conn;
- conn.reserve(2);
- for(int i=0;i<nbOfCells;i++)
- {
- conn.resize(0);
- getNodeIdsOfCell(i,conn);
- pt=std::transform(coo+conn[1]*spaceDim,coo+(conn[1]+1)*spaceDim,coo+conn[0]*spaceDim,pt,std::minus<double>());
- }
- ret->setArray(array);
- ret->setMesh(this);
- ret->synchronizeTimeWithSupport();
- return ret.retn();
+ if(_types.size()!=1 || *(_types.begin())!=INTERP_KERNEL::NORM_SEG2)
+ throw INTERP_KERNEL::Exception("Expected a umesh with only NORM_SEG2 type of elements for buildDirectionVectorField !");
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> array=DataArrayDouble::New();
+ int nbOfCells=getNumberOfCells();
+ int spaceDim=getSpaceDimension();
+ array->alloc(nbOfCells,spaceDim);
+ double *pt=array->getPointer();
+ const double *coo=getCoords()->getConstPointer();
+ std::vector<int> conn;
+ conn.reserve(2);
+ for(int i=0;i<nbOfCells;i++)
+ {
+ conn.resize(0);
+ getNodeIdsOfCell(i,conn);
+ pt=std::transform(coo+conn[1]*spaceDim,coo+(conn[1]+1)*spaceDim,coo+conn[0]*spaceDim,pt,std::minus<double>());
+ }
+ ret->setArray(array);
+ ret->setMesh(this);
+ ret->synchronizeTimeWithSupport();
+ return ret.retn();
}
/*!
{
if(getMeshDimension()!=1)
throw INTERP_KERNEL::Exception("Expected a umesh with meshDim == 1 for project1D !");
- if(_types.size()!=1 || *(_types.begin())!=INTERP_KERNEL::NORM_SEG2)
- throw INTERP_KERNEL::Exception("Expected a umesh with only NORM_SEG2 type of elements for project1D !");
- if(getSpaceDimension()!=3)
- throw INTERP_KERNEL::Exception("Expected a umesh with spaceDim==3 for project1D !");
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> f=buildDirectionVectorField();
- const double *fPtr=f->getArray()->getConstPointer();
- double tmp[3];
- for(int i=0;i<getNumberOfCells();i++)
- {
- const double *tmp1=fPtr+3*i;
- tmp[0]=tmp1[1]*v[2]-tmp1[2]*v[1];
- tmp[1]=tmp1[2]*v[0]-tmp1[0]*v[2];
- tmp[2]=tmp1[0]*v[1]-tmp1[1]*v[0];
- double n1=INTERP_KERNEL::norm<3>(tmp);
- n1/=INTERP_KERNEL::norm<3>(tmp1);
- if(n1>eps)
- throw INTERP_KERNEL::Exception("UMesh::Projection 1D failed !");
- }
- const double *coo=getCoords()->getConstPointer();
- for(int i=0;i<getNumberOfNodes();i++)
- {
- std::transform(coo+i*3,coo+i*3+3,pt,tmp,std::minus<double>());
- std::transform(tmp,tmp+3,v,tmp,std::multiplies<double>());
- res[i]=std::accumulate(tmp,tmp+3,0.);
- }
+ if(_types.size()!=1 || *(_types.begin())!=INTERP_KERNEL::NORM_SEG2)
+ throw INTERP_KERNEL::Exception("Expected a umesh with only NORM_SEG2 type of elements for project1D !");
+ if(getSpaceDimension()!=3)
+ throw INTERP_KERNEL::Exception("Expected a umesh with spaceDim==3 for project1D !");
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> f=buildDirectionVectorField();
+ const double *fPtr=f->getArray()->getConstPointer();
+ double tmp[3];
+ for(int i=0;i<getNumberOfCells();i++)
+ {
+ const double *tmp1=fPtr+3*i;
+ tmp[0]=tmp1[1]*v[2]-tmp1[2]*v[1];
+ tmp[1]=tmp1[2]*v[0]-tmp1[0]*v[2];
+ tmp[2]=tmp1[0]*v[1]-tmp1[1]*v[0];
+ double n1=INTERP_KERNEL::norm<3>(tmp);
+ n1/=INTERP_KERNEL::norm<3>(tmp1);
+ if(n1>eps)
+ throw INTERP_KERNEL::Exception("UMesh::Projection 1D failed !");
+ }
+ const double *coo=getCoords()->getConstPointer();
+ for(int i=0;i<getNumberOfNodes();i++)
+ {
+ std::transform(coo+i*3,coo+i*3+3,pt,tmp,std::minus<double>());
+ std::transform(tmp,tmp+3,v,tmp,std::multiplies<double>());
+ res[i]=std::accumulate(tmp,tmp+3,0.);
+ }
}
/*!
* \a this is expected to be a mesh so that its space dimension is equal to its
* mesh dimension + 1. Furthermore only mesh dimension 1 and 2 are supported for the moment.
* Distance from \a ptBg to \a ptEnd is expected to be equal to the space dimension. \a this is also expected to be fully defined (connectivity and coordinates).
-
+ *
* WARNING, if there is some orphan nodes in \a this (nodes not fetched by any cells in \a this ( see MEDCouplingUMesh::zipCoords ) ) these nodes will ** not ** been taken
* into account in this method. Only cells and nodes lying on them are considered in the algorithm (even if one of these orphan nodes is closer than returned distance).
* A user that needs to consider orphan nodes should invoke DataArrayDouble::minimalDistanceTo method on the coordinates array of \a this.
MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bboxArr(getBoundingBoxForBBTree());
const double *bbox(bboxArr->begin());
switch(spaceDim)
- {
+ {
case 3:
{
BBTreeDst<3> myTree(bbox,0,0,nbCells);
}
default:
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoints : only spacedim 2 and 3 supported !");
- }
+ }
cellIds=ret1.retn();
return ret0.retn();
}
for(const int *zeCell=cellIdsBg;zeCell!=cellIdsEnd;zeCell++)
{
switch((INTERP_KERNEL::NormalizedCellType)nc[ncI[*zeCell]])
- {
+ {
case INTERP_KERNEL::NORM_TRI3:
{
double tmp=INTERP_KERNEL::DistanceFromPtToTriInSpaceDim3(pt,coords+3*nc[ncI[*zeCell]+1],coords+3*nc[ncI[*zeCell]+2],coords+3*nc[ncI[*zeCell]+3]);
}
default:
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoint3DSurfAlg : not managed cell type ! Supporting TRI3, QUAD4 and POLYGON !");
- }
+ }
}
}
ret0=std::numeric_limits<double>::max();
for(const int *zeCell=cellIdsBg;zeCell!=cellIdsEnd;zeCell++)
{
- switch((INTERP_KERNEL::NormalizedCellType)nc[ncI[*zeCell]])
- {
+ switch((INTERP_KERNEL::NormalizedCellType)nc[ncI[*zeCell]])
+ {
case INTERP_KERNEL::NORM_SEG2:
{
std::size_t uselessEntry=0;
}
default:
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoint2DCurveAlg : not managed cell type ! Supporting SEG2 !");
- }
+ }
}
}
* \throw If the coordinates array is not set.
* \throw If \a this->getMeshDimension() != \a this->getSpaceDimension().
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_getCellsContainingPoint "Here is a C++ example".<br>
* \ref py_mcumesh_getCellsContainingPoint "Here is a Python example".
+ * \endif
*/
void MEDCouplingUMesh::getCellsContainingPoint(const double *pos, double eps, std::vector<int>& elts) const
{
// end
};
+
+
+ /*!
+ * Warning the nodes in \a m should be decrRefed ! To avoid that Node * pointer be replaced by another instance.
+ */
+ INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge2(INTERP_KERNEL::NormalizedCellType typ, const int *bg, const double *coords2D, std::map<INTERP_KERNEL::Node *,int>& m)
+ {
+ INTERP_KERNEL::Edge *ret=0;
+ INTERP_KERNEL::Node *n0(new INTERP_KERNEL::Node(coords2D[2*bg[0]],coords2D[2*bg[0]+1])),*n1(new INTERP_KERNEL::Node(coords2D[2*bg[1]],coords2D[2*bg[1]+1]));
+ m[n0]=bg[0]; m[n1]=bg[1];
+ switch(typ)
+ {
+ case INTERP_KERNEL::NORM_SEG2:
+ {
+ ret=new INTERP_KERNEL::EdgeLin(n0,n1);
+ break;
+ }
+ case INTERP_KERNEL::NORM_SEG3:
+ {
+ INTERP_KERNEL::Node *n2(new INTERP_KERNEL::Node(coords2D[2*bg[2]],coords2D[2*bg[2]+1])); m[n2]=bg[2];
+ INTERP_KERNEL::EdgeLin *e1(new INTERP_KERNEL::EdgeLin(n0,n2)),*e2(new INTERP_KERNEL::EdgeLin(n2,n1));
+ INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
+ // is the SEG3 degenerated, and thus can be reduced to a SEG2?
+ bool colinearity(inters.areColinears());
+ delete e1; delete e2;
+ if(colinearity)
+ { ret=new INTERP_KERNEL::EdgeLin(n0,n1); }
+ else
+ { ret=new INTERP_KERNEL::EdgeArcCircle(n0,n2,n1); }
+ break;
+ }
+ default:
+ throw INTERP_KERNEL::Exception("MEDCouplingUMeshBuildQPFromEdge2 : Expecting a mesh with spaceDim==2 and meshDim==1 !");
+ }
+ return ret;
+ }
+
INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge(INTERP_KERNEL::NormalizedCellType typ, std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >& mapp2, const int *bg)
{
INTERP_KERNEL::Edge *ret=0;
switch(typ)
- {
+ {
case INTERP_KERNEL::NORM_SEG2:
{
ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[1]].first);
}
default:
throw INTERP_KERNEL::Exception("MEDCouplingUMeshBuildQPFromEdge : Expecting a mesh with spaceDim==2 and meshDim==1 !");
- }
+ }
return ret;
}
* 'mapp' returns a mapping between local numbering in submesh (represented by a Node*) and the global node numbering in 'mDesc'.
*/
INTERP_KERNEL::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMesh(const MEDCouplingUMesh *mDesc, const std::vector<int>& candidates,
- std::map<INTERP_KERNEL::Node *,int>& mapp)
- throw(INTERP_KERNEL::Exception)
+ std::map<INTERP_KERNEL::Node *,int>& mapp)
{
mapp.clear();
std::map<int, std::pair<INTERP_KERNEL::Node *,bool> > mapp2;//bool is for a flag specifying if node is boundary (true) or only a middle for SEG3.
* \throw If the coordinates array is not set.
* \throw If \a this->getMeshDimension() != \a this->getSpaceDimension().
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_getCellsContainingPoints "Here is a C++ example".<br>
* \ref py_mcumesh_getCellsContainingPoints "Here is a Python example".
+ * \endif
*/
void MEDCouplingUMesh::getCellsContainingPoints(const double *pos, int nbOfPoints, double eps,
MEDCouplingAutoRefCountObjectPtr<DataArrayInt>& elts, MEDCouplingAutoRefCountObjectPtr<DataArrayInt>& eltsIndex) const
}
/*else if(mDim==2)
{
-
+
}*/
else
throw INTERP_KERNEL::Exception("For spaceDim==3 only meshDim==3 implemented for getelementscontainingpoints !");
* its connectivity will remain unchanged. If the computation leads to a modification of nodal connectivity of a cell its geometric type will be modified to INTERP_KERNEL::NORM_POLYGON.
*
* \return a newly allocated array containing cellIds that have been modified if any. If no cells have been impacted by this method NULL is returned.
+ * \sa MEDCouplingUMesh::colinearize2D
*/
DataArrayInt *MEDCouplingUMesh::convexEnvelop2D()
{
int oldNbOfNodes=getNumberOfNodes();
MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords;
switch(policy)
- {
+ {
case 0:
{
newCoords=fillExtCoordsUsingTranslation(mesh1D,isQuad);
}
default:
throw INTERP_KERNEL::Exception("Not implemented extrusion policy : must be in (0) !");
- }
+ }
setCoords(newCoords);
MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret=buildExtrudedMeshFromThisLowLev(oldNbOfNodes,isQuad);
updateTime();
double cosangle=i+1<nbOfLevsInVec?(p0r[0]-tmp3[0])*(p1r[0]-tmp3[0])+(p0r[1]-tmp3[1])*(p1r[1]-tmp3[1]):(p2r[0]-tmp3[0])*(p1r[0]-tmp3[0])+(p2r[1]-tmp3[1])*(p1r[1]-tmp3[1]);
double angle=acos(cosangle/(radius*radius));
tmp->rotate(end,vecPlane,angle);
-
}
retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr);
}
MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coordsSafe;
int meshDim=getMeshDimension();
switch(conversionType)
- {
+ {
case 0:
switch(meshDim)
- {
+ {
case 1:
ret=convertLinearCellsToQuadratic1D0(conn,connI,coords,types);
connSafe=conn; connISafe=connI; coordsSafe=coords;
break;
default:
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::convertLinearCellsToQuadratic : conversion of type 0 mesh dimensions available are [1,2,3] !");
- }
+ }
break;
- case 1:
- {
- switch(meshDim)
- {
case 1:
- ret=convertLinearCellsToQuadratic1D0(conn,connI,coords,types);//it is not a bug. In 1D policy 0 and 1 are equals
- connSafe=conn; connISafe=connI; coordsSafe=coords;
- break;
- case 2:
- ret=convertLinearCellsToQuadratic2D1(conn,connI,coords,types);
- connSafe=conn; connISafe=connI; coordsSafe=coords;
- break;
- case 3:
- ret=convertLinearCellsToQuadratic3D1(conn,connI,coords,types);
- connSafe=conn; connISafe=connI; coordsSafe=coords;
- break;
+ {
+ switch(meshDim)
+ {
+ case 1:
+ ret=convertLinearCellsToQuadratic1D0(conn,connI,coords,types);//it is not a bug. In 1D policy 0 and 1 are equals
+ connSafe=conn; connISafe=connI; coordsSafe=coords;
+ break;
+ case 2:
+ ret=convertLinearCellsToQuadratic2D1(conn,connI,coords,types);
+ connSafe=conn; connISafe=connI; coordsSafe=coords;
+ break;
+ case 3:
+ ret=convertLinearCellsToQuadratic3D1(conn,connI,coords,types);
+ connSafe=conn; connISafe=connI; coordsSafe=coords;
+ break;
+ default:
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::convertLinearCellsToQuadratic : conversion of type 1 mesh dimensions available are [1,2,3] !");
+ }
+ break;
+ }
default:
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::convertLinearCellsToQuadratic : conversion of type 1 mesh dimensions available are [1,2,3] !");
- }
- break;
- }
- default:
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::convertLinearCellsToQuadratic : conversion type available are 0 (default, the simplest) and 1 (the most complex) !");
- }
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::convertLinearCellsToQuadratic : conversion type available are 0 (default, the simplest) and 1 (the most complex) !");
+ }
setConnectivity(connSafe,connISafe,false);
_types=types;
setCoords(coordsSafe);
*/
DataArrayInt *MEDCouplingUMesh::convertLinearCellsToQuadratic2D0(DataArrayInt *&conn, DataArrayInt *&connI, DataArrayDouble *& coords, std::set<INTERP_KERNEL::NormalizedCellType>& types) const
{
-
MEDCouplingAutoRefCountObjectPtr<DataArrayInt> desc(DataArrayInt::New()),descI(DataArrayInt::New()),tmp2(DataArrayInt::New()),tmp3(DataArrayInt::New());
MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m1D=buildDescendingConnectivity(desc,descI,tmp2,tmp3); tmp2=0; tmp3=0;
return convertLinearCellsToQuadratic2DAnd3D0(m1D,desc,descI,conn,connI,coords,types);
DataArrayInt *MEDCouplingUMesh::simplexize(int policy)
{
switch(policy)
- {
+ {
case 0:
return simplexizePol0();
case 1:
return simplexizePol1();
case (int) INTERP_KERNEL::PLANAR_FACE_5:
- return simplexizePlanarFace5();
+ return simplexizePlanarFace5();
case (int) INTERP_KERNEL::PLANAR_FACE_6:
- return simplexizePlanarFace6();
+ return simplexizePlanarFace6();
default:
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::simplexize : unrecognized policy ! Must be :\n - 0 or 1 (only available for meshdim=2) \n - PLANAR_FACE_5, PLANAR_FACE_6 (only for meshdim=3)");
- }
+ }
}
/*!
if((INTERP_KERNEL::NormalizedCellType)oldc[ci[0]]==INTERP_KERNEL::NORM_QUAD4)
{
const int tmp[8]={(int)INTERP_KERNEL::NORM_TRI3,oldc[ci[0]+1],oldc[ci[0]+2],oldc[ci[0]+3],
- (int)INTERP_KERNEL::NORM_TRI3,oldc[ci[0]+1],oldc[ci[0]+3],oldc[ci[0]+4]};
+ (int)INTERP_KERNEL::NORM_TRI3,oldc[ci[0]+1],oldc[ci[0]+3],oldc[ci[0]+4]};
pt=std::copy(tmp,tmp+8,pt);
ptI[1]=ptI[0]+4;
ptI[2]=ptI[0]+8;
if((INTERP_KERNEL::NormalizedCellType)oldc[ci[0]]==INTERP_KERNEL::NORM_QUAD4)
{
const int tmp[8]={(int)INTERP_KERNEL::NORM_TRI3,oldc[ci[0]+1],oldc[ci[0]+2],oldc[ci[0]+4],
- (int)INTERP_KERNEL::NORM_TRI3,oldc[ci[0]+2],oldc[ci[0]+3],oldc[ci[0]+4]};
+ (int)INTERP_KERNEL::NORM_TRI3,oldc[ci[0]+2],oldc[ci[0]+3],oldc[ci[0]+4]};
pt=std::copy(tmp,tmp+8,pt);
ptI[1]=ptI[0]+4;
ptI[2]=ptI[0]+8;
* \throw If \a this->getMeshDimension() != 2.
* \throw If \a this->getSpaceDimension() != 3.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_are2DCellsNotCorrectlyOriented "Here is a C++ example".<br>
* \ref py_mcumesh_are2DCellsNotCorrectlyOriented "Here is a Python example".
+ * \endif
*/
void MEDCouplingUMesh::are2DCellsNotCorrectlyOriented(const double *vec, bool polyOnly, std::vector<int>& cells) const
{
* \throw If \a this->getMeshDimension() != 2.
* \throw If \a this->getSpaceDimension() != 3.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_are2DCellsNotCorrectlyOriented "Here is a C++ example".<br>
* \ref py_mcumesh_are2DCellsNotCorrectlyOriented "Here is a Python example".
+ * \endif
*/
void MEDCouplingUMesh::orientCorrectly2DCells(const double *vec, bool polyOnly)
{
* \throw If the coordinates array is not set.
* \throw If the nodal connectivity of cells is not defined.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_arePolyhedronsNotCorrectlyOriented "Here is a C++ example".<br>
* \ref py_mcumesh_arePolyhedronsNotCorrectlyOriented "Here is a Python example".
+ * \endif
*/
void MEDCouplingUMesh::arePolyhedronsNotCorrectlyOriented(std::vector<int>& cells) const
{
* \throw If the nodal connectivity of cells is not defined.
* \throw If the reparation fails.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_arePolyhedronsNotCorrectlyOriented "Here is a C++ example".<br>
* \ref py_mcumesh_arePolyhedronsNotCorrectlyOriented "Here is a Python example".
+ * \endif
* \sa MEDCouplingUMesh::findAndCorrectBadOriented3DCells
*/
void MEDCouplingUMesh::orientCorrectlyPolyhedrons()
if(type==INTERP_KERNEL::NORM_POLYHED)
{
try
- {
+ {
if(!IsPolyhedronWellOriented(conn+connI[i]+1,conn+connI[i+1],coordsPtr))
TryToCorrectPolyhedronOrientation(conn+connI[i]+1,conn+connI[i+1],coordsPtr);
- }
+ }
catch(INTERP_KERNEL::Exception& e)
- {
+ {
std::ostringstream oss; oss << "Something wrong in polyhedron #" << i << " : " << e.what();
throw INTERP_KERNEL::Exception(oss.str().c_str());
- }
+ }
}
}
updateTime();
* \throw If the coordinates array is not set.
* \throw If the nodal connectivity of cells is not defined.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_findAndCorrectBadOriented3DExtrudedCells "Here is a C++ example".<br>
* \ref py_mcumesh_findAndCorrectBadOriented3DExtrudedCells "Here is a Python example".
+ * \endif
* \sa MEDCouplingUMesh::findAndCorrectBadOriented3DCells
*/
DataArrayInt *MEDCouplingUMesh::findAndCorrectBadOriented3DExtrudedCells()
{
INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)conn[connI[i]];
switch(type)
- {
+ {
case INTERP_KERNEL::NORM_TETRA4:
{
if(!IsTetra4WellOriented(conn+connI[i]+1,conn+connI[i+1],coordsPtr))
}
default:
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::orientCorrectly3DCells : Your mesh contains type of cell not supported yet ! send mail to anthony.geay@cea.fr to add it !");
- }
+ }
}
updateTime();
return ret.retn();
{
INTERP_KERNEL::NormalizedCellType t=(INTERP_KERNEL::NormalizedCellType)*conn;
switch(t)
- {
- case INTERP_KERNEL::NORM_TRI3:
- {
- FillInCompact3DMode(spaceDim,3,conn+1,coo,tmp);
- *pt=INTERP_KERNEL::triEdgeRatio(tmp);
- break;
- }
- case INTERP_KERNEL::NORM_QUAD4:
- {
- FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
- *pt=INTERP_KERNEL::quadEdgeRatio(tmp);
- break;
- }
- case INTERP_KERNEL::NORM_TETRA4:
- {
- FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
- *pt=INTERP_KERNEL::tetraEdgeRatio(tmp);
- break;
- }
+ {
+ case INTERP_KERNEL::NORM_TRI3:
+ {
+ FillInCompact3DMode(spaceDim,3,conn+1,coo,tmp);
+ *pt=INTERP_KERNEL::triEdgeRatio(tmp);
+ break;
+ }
+ case INTERP_KERNEL::NORM_QUAD4:
+ {
+ FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
+ *pt=INTERP_KERNEL::quadEdgeRatio(tmp);
+ break;
+ }
+ case INTERP_KERNEL::NORM_TETRA4:
+ {
+ FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
+ *pt=INTERP_KERNEL::tetraEdgeRatio(tmp);
+ break;
+ }
default:
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getEdgeRatioField : A cell with not manged type (NORM_TRI3, NORM_QUAD4 and NORM_TETRA4) has been detected !");
- }
+ }
conn+=connI[i+1]-connI[i];
}
ret->setName("EdgeRatio");
{
INTERP_KERNEL::NormalizedCellType t=(INTERP_KERNEL::NormalizedCellType)*conn;
switch(t)
- {
- case INTERP_KERNEL::NORM_TRI3:
- {
- FillInCompact3DMode(spaceDim,3,conn+1,coo,tmp);
- *pt=INTERP_KERNEL::triAspectRatio(tmp);
- break;
- }
- case INTERP_KERNEL::NORM_QUAD4:
- {
- FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
- *pt=INTERP_KERNEL::quadAspectRatio(tmp);
- break;
- }
- case INTERP_KERNEL::NORM_TETRA4:
- {
- FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
- *pt=INTERP_KERNEL::tetraAspectRatio(tmp);
- break;
- }
+ {
+ case INTERP_KERNEL::NORM_TRI3:
+ {
+ FillInCompact3DMode(spaceDim,3,conn+1,coo,tmp);
+ *pt=INTERP_KERNEL::triAspectRatio(tmp);
+ break;
+ }
+ case INTERP_KERNEL::NORM_QUAD4:
+ {
+ FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
+ *pt=INTERP_KERNEL::quadAspectRatio(tmp);
+ break;
+ }
+ case INTERP_KERNEL::NORM_TETRA4:
+ {
+ FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
+ *pt=INTERP_KERNEL::tetraAspectRatio(tmp);
+ break;
+ }
default:
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getAspectRatioField : A cell with not manged type (NORM_TRI3, NORM_QUAD4 and NORM_TETRA4) has been detected !");
- }
+ }
conn+=connI[i+1]-connI[i];
}
ret->setName("AspectRatio");
{
INTERP_KERNEL::NormalizedCellType t=(INTERP_KERNEL::NormalizedCellType)*conn;
switch(t)
- {
- case INTERP_KERNEL::NORM_QUAD4:
- {
- FillInCompact3DMode(3,4,conn+1,coo,tmp);
- *pt=INTERP_KERNEL::quadWarp(tmp);
- break;
- }
+ {
+ case INTERP_KERNEL::NORM_QUAD4:
+ {
+ FillInCompact3DMode(3,4,conn+1,coo,tmp);
+ *pt=INTERP_KERNEL::quadWarp(tmp);
+ break;
+ }
default:
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getWarpField : A cell with not manged type (NORM_QUAD4) has been detected !");
- }
+ }
conn+=connI[i+1]-connI[i];
}
ret->setName("Warp");
{
INTERP_KERNEL::NormalizedCellType t=(INTERP_KERNEL::NormalizedCellType)*conn;
switch(t)
- {
- case INTERP_KERNEL::NORM_QUAD4:
- {
- FillInCompact3DMode(3,4,conn+1,coo,tmp);
- *pt=INTERP_KERNEL::quadSkew(tmp);
- break;
- }
+ {
+ case INTERP_KERNEL::NORM_QUAD4:
+ {
+ FillInCompact3DMode(3,4,conn+1,coo,tmp);
+ *pt=INTERP_KERNEL::quadSkew(tmp);
+ break;
+ }
default:
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getSkewField : A cell with not manged type (NORM_QUAD4) has been detected !");
- }
+ }
conn+=connI[i+1]-connI[i];
}
ret->setName("Skew");
DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree(double arcDetEps) const
{
int mDim(getMeshDimension()),sDim(getSpaceDimension());
- if(mDim!=2 || sDim!=2)
+ if((mDim==3 && sDim==3) || (mDim==2 && sDim==3) || (mDim==1 && sDim==1) || ( mDim==1 && sDim==3)) // Compute refined boundary box for quadratic elements only in 2D.
return getBoundingBoxForBBTreeFast();
- else
+ if((mDim==2 && sDim==2) || (mDim==1 && sDim==2))
{
bool presenceOfQuadratic(false);
for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator it=_types.begin();it!=_types.end();it++)
if(cm.isQuadratic())
presenceOfQuadratic=true;
}
- if(presenceOfQuadratic)
+ if(!presenceOfQuadratic)
+ return getBoundingBoxForBBTreeFast();
+ if(mDim==2 && sDim==2)
return getBoundingBoxForBBTree2DQuadratic(arcDetEps);
else
- return getBoundingBoxForBBTreeFast();
+ return getBoundingBoxForBBTree1DQuadratic(arcDetEps);
}
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getBoundingBoxForBBTree : Managed dimensions are (mDim=1,sDim=1), (mDim=1,sDim=2), (mDim=1,sDim=3), (mDim=2,sDim=2), (mDim=2,sDim=3) and (mDim=3,sDim=3) !");
}
/*!
}
/*!
- * This method aggregate the bbox regarding foreach 2D cell in \a this the whole shape. So this method is particulary useful for 2D meshes having quadratic cells
- * because for this type of cells getBoundingBoxForBBTreeFast method may return invalid bounding boxes.
+ * This method aggregates the bbox of each 2D cell in \a this considering the whole shape. This method is particularly
+ * useful for 2D meshes having quadratic cells
+ * because for this type of cells getBoundingBoxForBBTreeFast method may return invalid bounding boxes (since it just considers
+ * the two extremities of the arc of circle).
*
* \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12)
* \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components.
* \throw If \a this is not fully defined.
* \throw If \a this is not a mesh with meshDimension equal to 2.
* \throw If \a this is not a mesh with spaceDimension equal to 2.
+ * \sa MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic
*/
DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arcDetEps) const
{
checkFullyDefined();
- int spaceDim(getSpaceDimension()),mDim(getMeshDimension()),nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
- if(mDim!=2 || spaceDim!=2)
+ 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!");
MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim);
double *bbox(ret->getPointer());
{
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=1e-12;
+ 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++)
pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
else
pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
- INTERP_KERNEL::Bounds b; pol->fillBounds(b); delete pol;
+ INTERP_KERNEL::Bounds b; b.prepareForAggregation(); pol->fillBounds(b); delete pol;
bbox[0]=b.getXMin(); bbox[1]=b.getXMax(); bbox[2]=b.getYMin(); bbox[3]=b.getYMax();
}
return ret.retn();
}
-/// @cond INTERNAL
-
-namespace ParaMEDMEMImpl
+/*!
+ * This method aggregates the bbox of each 1D cell in \a this considering the whole shape. This method is particularly
+ * useful for 2D meshes having quadratic cells
+ * because for this type of cells getBoundingBoxForBBTreeFast method may return invalid bounding boxes (since it just considers
+ * the two extremities of the arc of circle).
+ *
+ * \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12)
+ * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components.
+ * \throw If \a this is not fully defined.
+ * \throw If \a this is not a mesh with meshDimension equal to 1.
+ * \throw If \a this is not a mesh with spaceDimension equal to 2.
+ * \sa MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic
+ */
+DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic(double arcDetEps) const
{
- class ConnReader
- {
+ checkFullyDefined();
+ 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!");
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim);
+ double *bbox(ret->getPointer());
+ const double *coords(_coords->getConstPointer());
+ const int *conn(_nodal_connec->getConstPointer()),*connI(_nodal_connec_index->getConstPointer());
+ for(int i=0;i<nbOfCells;i++,bbox+=4,connI++)
+ {
+ 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++)
+ {
+ int nodeId(conn[*connI+1+j]);
+ nodes[j]=new INTERP_KERNEL::Node(coords[nodeId*2],coords[nodeId*2+1]);
+ }
+ if(!cm.isQuadratic())
+ edge=INTERP_KERNEL::QuadraticPolygon::BuildLinearEdge(nodes);
+ else
+ edge=INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(nodes);
+ const INTERP_KERNEL::Bounds& b(edge->getBounds());
+ bbox[0]=b.getXMin(); bbox[1]=b.getXMax(); bbox[2]=b.getYMin(); bbox[3]=b.getYMax(); edge->decrRef();
+ }
+ return ret.retn();
+}
+
+/// @cond INTERNAL
+
+namespace ParaMEDMEMImpl
+{
+ class ConnReader
+ {
public:
ConnReader(const int *c, int val):_conn(c),_val(val) { }
bool operator() (const int& pos) { return _conn[pos]!=_val; }
MEDCoupling1GTUMesh *MEDCouplingUMesh::convertIntoSingleGeoTypeMesh() const
{
checkConnectivityFullyDefined();
- if(_types.size()!=1)
+ if(_types.size()!=1)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::convertIntoSingleGeoTypeMesh : current mesh does not contain exactly one geometric type !");
INTERP_KERNEL::NormalizedCellType typ=*_types.begin();
MEDCouplingAutoRefCountObjectPtr<MEDCoupling1GTUMesh> ret=MEDCoupling1GTUMesh::New(getName(),typ);
DataArrayInt *MEDCouplingUMesh::convertNodalConnectivityToStaticGeoTypeMesh() const
{
checkConnectivityFullyDefined();
- if(_types.size()!=1)
+ if(_types.size()!=1)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::convertNodalConnectivityToStaticGeoTypeMesh : current mesh does not contain exactly one geometric type !");
INTERP_KERNEL::NormalizedCellType typ=*_types.begin();
const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
*/
MEDCouplingUMesh *MEDCouplingUMesh::AggregateSortedByTypeMeshesOnSameCoords(const std::vector<const MEDCouplingUMesh *>& ms,
DataArrayInt *&szOfCellGrpOfSameType,
- DataArrayInt *&idInMsOfCellGrpOfSameType) throw(INTERP_KERNEL::Exception)
+ DataArrayInt *&idInMsOfCellGrpOfSameType)
{
std::vector<const MEDCouplingUMesh *> ms2;
for(std::vector<const MEDCouplingUMesh *>::const_iterator it=ms.begin();it!=ms.end();it++)
* \throw If the coordinates array is not set.
* \throw If the nodal connectivity of cells is not defined.
*
+ * \if ENABLE_EXAMPLES
* \ref cpp_mcumesh_getPartBarycenterAndOwner "Here is a C++ example".<br>
* \ref py_mcumesh_getPartBarycenterAndOwner "Here is a Python example".
+ * \endif
*/
DataArrayDouble *MEDCouplingUMesh::getPartBarycenterAndOwner(const int *begin, const int *end) const
{
* \throw If the coordinates array is not set in none of the meshes.
* \throw If \a a[ *i* ]->getMeshDimension() < 0.
* \throw If the meshes in \a a are of different dimension (getMeshDimension()).
-*/
+ */
MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshes(std::vector<const MEDCouplingUMesh *>& a)
{
std::size_t sz=a.size();
ret.push_back(cm.getExtrudedType());
int deltaz=isQuad?2*nbOfNodesPerLev:nbOfNodesPerLev;
switch(flatType)
- {
+ {
case INTERP_KERNEL::NORM_POINT1:
{
ret.push_back(connBg[1]);
case INTERP_KERNEL::NORM_TRI6:
{
int conn[15]={connBg[1],connBg[2],connBg[3],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4],connBg[5],connBg[6],connBg[4]+deltaz,connBg[5]+deltaz,connBg[6]+deltaz,
- connBg[1]+nbOfNodesPerLev,connBg[2]+nbOfNodesPerLev,connBg[3]+nbOfNodesPerLev};
+ connBg[1]+nbOfNodesPerLev,connBg[2]+nbOfNodesPerLev,connBg[3]+nbOfNodesPerLev};
ret.insert(ret.end(),conn,conn+15);
break;
}
}
default:
throw INTERP_KERNEL::Exception("A flat type has been detected that has not its extruded representation !");
- }
+ }
}
/*!
}
}
-/*!
- * This method makes the assumption spacedimension == meshdimension == 2.
- * This method works only for linear cells.
- *
- * \return a newly allocated array containing the connectivity of a polygon type enum included (NORM_POLYGON in pos#0)
- */
-DataArrayInt *MEDCouplingUMesh::buildUnionOf2DMesh() const
+DataArrayInt *MEDCouplingUMesh::buildUnionOf2DMeshLinear(const MEDCouplingUMesh *skin, const DataArrayInt *n2o) const
{
- if(getMeshDimension()!=2 || getSpaceDimension()!=2)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMesh : meshdimension, spacedimension must be equal to 2 !");
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=computeSkin();
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> o2n=m->zipCoordsTraducer();
- int nbOfNodesExpected=m->getNumberOfNodes();
- if(m->getNumberOfCells()!=nbOfNodesExpected)
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMesh : the mesh 2D in input appears to be not in a single part or a quadratic 2D mesh !");
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> n2o=o2n->invertArrayO2N2N2O(m->getNumberOfNodes());
- const int *n2oPtr=n2o->getConstPointer();
+ int nbOfNodesExpected(skin->getNumberOfNodes());
+ const int *n2oPtr(n2o->getConstPointer());
MEDCouplingAutoRefCountObjectPtr<DataArrayInt> revNodal(DataArrayInt::New()),revNodalI(DataArrayInt::New());
- m->getReverseNodalConnectivity(revNodal,revNodalI);
- const int *revNodalPtr=revNodal->getConstPointer(),*revNodalIPtr=revNodalI->getConstPointer();
- const int *nodalPtr=m->getNodalConnectivity()->getConstPointer();
- const int *nodalIPtr=m->getNodalConnectivityIndex()->getConstPointer();
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New(); ret->alloc(nbOfNodesExpected+1,1);
- int *work=ret->getPointer(); *work++=INTERP_KERNEL::NORM_POLYGON;
+ skin->getReverseNodalConnectivity(revNodal,revNodalI);
+ const int *revNodalPtr(revNodal->getConstPointer()),*revNodalIPtr(revNodalI->getConstPointer());
+ const int *nodalPtr(skin->getNodalConnectivity()->getConstPointer());
+ const int *nodalIPtr(skin->getNodalConnectivityIndex()->getConstPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(nbOfNodesExpected+1,1);
+ int *work(ret->getPointer()); *work++=INTERP_KERNEL::NORM_POLYGON;
if(nbOfNodesExpected<1)
return ret.retn();
- int prevCell=0;
- int prevNode=nodalPtr[nodalIPtr[0]+1];
+ int prevCell(0),prevNode(nodalPtr[nodalIPtr[0]+1]);
*work++=n2oPtr[prevNode];
for(int i=1;i<nbOfNodesExpected;i++)
{
conn.erase(prevNode);
if(conn.size()==1)
{
- int curNode=*(conn.begin());
+ int curNode(*(conn.begin()));
*work++=n2oPtr[curNode];
std::set<int> shar(revNodalPtr+revNodalIPtr[curNode],revNodalPtr+revNodalIPtr[curNode+1]);
shar.erase(prevCell);
prevNode=curNode;
}
else
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMesh : presence of unexpected 2 !");
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected 2 !");
}
else
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMesh : presence of unexpected 1 !");
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected 1 !");
}
else
- throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMesh : presence of unexpected cell !");
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshLinear : presence of unexpected cell !");
+ }
+ return ret.retn();
+}
+
+DataArrayInt *MEDCouplingUMesh::buildUnionOf2DMeshQuadratic(const MEDCouplingUMesh *skin, const DataArrayInt *n2o) const
+{
+ int nbOfNodesExpected(skin->getNumberOfNodes());
+ int nbOfTurn(nbOfNodesExpected/2);
+ const int *n2oPtr(n2o->getConstPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> revNodal(DataArrayInt::New()),revNodalI(DataArrayInt::New());
+ skin->getReverseNodalConnectivity(revNodal,revNodalI);
+ const int *revNodalPtr(revNodal->getConstPointer()),*revNodalIPtr(revNodalI->getConstPointer());
+ const int *nodalPtr(skin->getNodalConnectivity()->getConstPointer());
+ const int *nodalIPtr(skin->getNodalConnectivityIndex()->getConstPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(nbOfNodesExpected+1,1);
+ int *work(ret->getPointer()); *work++=INTERP_KERNEL::NORM_QPOLYG;
+ if(nbOfNodesExpected<1)
+ return ret.retn();
+ int prevCell(0),prevNode(nodalPtr[nodalIPtr[0]+1]);
+ *work=n2oPtr[prevNode]; work[nbOfTurn]=n2oPtr[nodalPtr[nodalIPtr[0]+3]]; work++;
+ for(int i=1;i<nbOfTurn;i++)
+ {
+ if(nodalIPtr[prevCell+1]-nodalIPtr[prevCell]==4)
+ {
+ std::set<int> conn(nodalPtr+nodalIPtr[prevCell]+1,nodalPtr+nodalIPtr[prevCell]+3);
+ conn.erase(prevNode);
+ if(conn.size()==1)
+ {
+ int curNode(*(conn.begin()));
+ *work=n2oPtr[curNode];
+ std::set<int> shar(revNodalPtr+revNodalIPtr[curNode],revNodalPtr+revNodalIPtr[curNode+1]);
+ shar.erase(prevCell);
+ if(shar.size()==1)
+ {
+ int curCell(*(shar.begin()));
+ work[nbOfTurn]=n2oPtr[nodalPtr[nodalIPtr[curCell]+3]];
+ prevCell=curCell;
+ prevNode=curNode;
+ work++;
+ }
+ else
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected 2 !");
+ }
+ else
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected 1 !");
+ }
+ else
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMeshQuadratic : presence of unexpected cell !");
}
return ret.retn();
}
+/*!
+ * This method makes the assumption spacedimension == meshdimension == 2.
+ * This method works only for linear cells.
+ *
+ * \return a newly allocated array containing the connectivity of a polygon type enum included (NORM_POLYGON in pos#0)
+ */
+DataArrayInt *MEDCouplingUMesh::buildUnionOf2DMesh() const
+{
+ if(getMeshDimension()!=2 || getSpaceDimension()!=2)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMesh : meshdimension, spacedimension must be equal to 2 !");
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> skin(computeSkin());
+ int oldNbOfNodes(skin->getNumberOfNodes());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> o2n(skin->zipCoordsTraducer());
+ int nbOfNodesExpected(skin->getNumberOfNodes());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> n2o(o2n->invertArrayO2N2N2O(oldNbOfNodes));
+ int nbCells(skin->getNumberOfCells());
+ if(nbCells==nbOfNodesExpected)
+ return buildUnionOf2DMeshLinear(skin,n2o);
+ else if(2*nbCells==nbOfNodesExpected)
+ return buildUnionOf2DMeshQuadratic(skin,n2o);
+ else
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildUnionOf2DMesh : the mesh 2D in input appears to be not in a single part of a 2D mesh !");
+}
+
/*!
* This method makes the assumption spacedimension == meshdimension == 3.
* This method works only for linear cells.
w4=std::copy(c.begin(),c.end(),w4);
}
}
- types->transformWithIndArr(PARAMEDMEM2VTKTYPETRADUCER,PARAMEDMEM2VTKTYPETRADUCER+INTERP_KERNEL::NORM_MAXTYPE);
+ types->transformWithIndArr(PARAMEDMEM2VTKTYPETRADUCER,PARAMEDMEM2VTKTYPETRADUCER+INTERP_KERNEL::NORM_MAXTYPE+1);
types->writeVTK(ofs,8,"UInt8","types",byteData);
offsets->writeVTK(ofs,8,"Int32","offsets",byteData);
if(szFaceOffsets!=0)
return std::string("UnstructuredGrid");
}
+std::string MEDCouplingUMesh::getVTKFileExtension() const
+{
+ return std::string("vtu");
+}
+
/*!
* Partitions the first given 2D mesh using the second given 2D mesh as a tool, and
* returns a result mesh constituted by polygons.
MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2)
{
+ if(!m1 || !m2)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes : input meshes must be not NULL !");
m1->checkFullyDefined();
m2->checkFullyDefined();
if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
MEDCouplingUMesh *m1Desc=0,*m2Desc=0; // descending connec. meshes
DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0;
std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
- INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
- INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,
- m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
+ m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
addCoo, m2Desc,desc2,descIndx2,revDesc2,revDescIndx2);
revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef();
MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2);
/* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2);
// Step 4: Prepare final result:
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> addCooDa=DataArrayDouble::New();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> addCooDa(DataArrayDouble::New());
addCooDa->alloc((int)(addCoo.size())/2,2);
std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> addCoordsQuadraticDa=DataArrayDouble::New();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> addCoordsQuadraticDa(DataArrayDouble::New());
addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2);
std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
std::vector<const DataArrayDouble *> coordss(4);
coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
- MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo=DataArrayDouble::Aggregate(coordss);
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret=MEDCouplingUMesh::New("Intersect2D",2);
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn=DataArrayInt::New(); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> connI=DataArrayInt::New(); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c1=DataArrayInt::New(); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c2=DataArrayInt::New(); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo(DataArrayDouble::Aggregate(coordss));
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret(MEDCouplingUMesh::New("Intersect2D",2));
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New()); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> connI(DataArrayInt::New()); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c1(DataArrayInt::New()); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c2(DataArrayInt::New()); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
ret->setConnectivity(conn,connI,true);
ret->setCoords(coo);
cellNb1=c1.retn(); cellNb2=c2.retn();
return ret.retn();
}
+//tony to put in private of MEDCouplingUMesh
+MEDCouplingUMesh *BuildMesh1DCutFrom(const MEDCouplingUMesh *mesh1D, const std::vector< std::vector<int> >& intersectEdge2, const DataArrayDouble *coords1, const std::vector<double>& addCoo)
+{
+ int nCells(mesh1D->getNumberOfCells());
+ if(nCells!=(int)intersectEdge2.size())
+ throw INTERP_KERNEL::Exception("BuildMesh1DCutFrom : internal error # 1 !");
+ const DataArrayDouble *coo2(mesh1D->getCoords());
+ const int *c(mesh1D->getNodalConnectivity()->begin()),*ci(mesh1D->getNodalConnectivityIndex()->begin());
+ const double *coo2Ptr(coo2->begin());
+ int offset1(coords1->getNumberOfTuples());
+ int offset2(offset1+coo2->getNumberOfTuples());
+ int offset3(offset2+addCoo.size()/2);
+ std::vector<double> addCooQuad;
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cOut(DataArrayInt::New()),ciOut(DataArrayInt::New()); cOut->alloc(0,1); ciOut->alloc(1,1); ciOut->setIJ(0,0,0);
+ int tmp[4],cicnt(0);
+ for(int i=0;i<nCells;i++)
+ {
+ std::map<INTERP_KERNEL::Node *,int> m;
+ INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coo2Ptr,m));
+ const std::vector<int>& subEdges(intersectEdge2[i]);
+ int nbSubEdge(subEdges.size()/2);
+ for(int j=0;j<nbSubEdge;j++)
+ {
+ MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node> n1(MEDCouplingUMeshBuildQPNode(subEdges[2*j],coords1->begin(),offset1,coo2Ptr,offset2,addCoo));
+ MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node> n2(MEDCouplingUMeshBuildQPNode(subEdges[2*j+1],coords1->begin(),offset1,coo2Ptr,offset2,addCoo));
+ MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> e2(e->buildEdgeLyingOnMe(n1,n2));
+ INTERP_KERNEL::Edge *e2Ptr(e2);
+ if(dynamic_cast<INTERP_KERNEL::EdgeArcCircle *>(e2Ptr))
+ {
+ tmp[0]=INTERP_KERNEL::NORM_SEG3;
+ tmp[1]=subEdges[2*j]; tmp[2]=subEdges[2*j+1];
+ cicnt+=4;
+ cOut->insertAtTheEnd(tmp,tmp+4);
+ ciOut->pushBackSilent(cicnt);
+ }
+ else
+ {
+ tmp[0]=INTERP_KERNEL::NORM_SEG2;
+ tmp[1]=subEdges[2*j]; tmp[2]=subEdges[2*j+1]; tmp[3]=offset3+(int)addCooQuad.size()/2;
+ cicnt+=3;
+ cOut->insertAtTheEnd(tmp,tmp+3);
+ ciOut->pushBackSilent(cicnt);
+ }
+ }
+ //INTERP_KERNEL::Edge *e2(e->buildEdgeLyingOnMe());
+ for(std::map<INTERP_KERNEL::Node *,int>::const_iterator it2=m.begin();it2!=m.end();it2++)
+ (*it2).first->decrRef();
+ e->decrRef();
+ }
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret(MEDCouplingUMesh::New(mesh1D->getName(),1));
+ ret->setConnectivity(cOut,ciOut,true);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr3(DataArrayDouble::New());
+ arr3->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr4(DataArrayDouble::New()); arr4->useArray(&addCooQuad[0],false,C_DEALLOC,(int)addCooQuad.size()/2,2);
+ std::vector<const DataArrayDouble *> coordss(4);
+ coordss[0]=coords1; coordss[1]=mesh1D->getCoords(); coordss[2]=arr3; coordss[3]=arr4;
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr(DataArrayDouble::Aggregate(coordss));
+ ret->setCoords(arr);
+ return ret.retn();
+}
+
+/*!
+ * Partitions the first given 2D mesh using the second given 1D mesh as a tool.
+ * Thus the final result contains all nodes from m1 plus new nodes. However it doesn't necessarily contains
+ * all nodes from \a mesh1D.
+ * The meshes should be in 2D space. In addition, returns two arrays mapping cells of the resulting mesh to cells of the input.
+ *
+ * \param [in] mesh2D - the 2D mesh (spacedim=meshdim=2) to be intersected using \a mesh1D tool.
+ * \param [in] mesh1D - the 1D mesh (spacedim=2 meshdim=1) the is the tool that will be used to intersect \a mesh2D.
+ * \param [in] eps - precision used to perform intersections and localization operations.
+ * \param [out] splitMesh2D - the result of the split of \a mesh2D mesh.
+ * \param [out] splitMesh1D - the result of the split of \a mesh1D mesh.
+ * \param [out] cellIdInMesh2D - the array that gives for each cell id \a i in \a splitMesh2D the id in \a mesh2D it comes from.
+ * So this array has a number of tuples equal to the number of cells of \a splitMesh2D and a number of component equal to 1.
+ * \param [out] cellIdInMesh1D - the array that gives for each cell id \a i in \a splitMesh1D the 1 or 2 id(s) in \a splitMesh2D that \a i shares.
+ * So this array has a number of tuples equal to the number of cells of \a splitMesh1D and a number of components equal to 2.
+ */
+void MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *mesh2D, const MEDCouplingUMesh *mesh1D, double eps, MEDCouplingUMesh *&splitMesh2D, MEDCouplingUMesh *&splitMesh1D, DataArrayInt *&cellIdInMesh2D, DataArrayInt *&cellIdInMesh1D)
+{
+ if(!mesh2D || !mesh1D)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine : input meshes must be not NULL !");
+ mesh2D->checkFullyDefined();
+ mesh1D->checkFullyDefined();
+ if(mesh2D->getMeshDimension()!=2 || mesh2D->getSpaceDimension()!=2 || mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=2)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works with mesh2D with spacedim=meshdim=2 and mesh1D with meshdim=1 spaceDim=2 !");
+ // Step 1: compute all edge intersections (new nodes)
+ std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
+ std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
+ INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
+ INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
+ //
+ // Build desc connectivity
+ DataArrayInt *desc1(DataArrayInt::New()),*descIndx1(DataArrayInt::New()),*revDesc1(DataArrayInt::New()),*revDescIndx1(DataArrayInt::New());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m1Desc(mesh2D->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1));
+ Intersect1DMeshes(m1Desc,mesh1D,eps,intersectEdge1,colinear2,subDiv2,addCoo);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> addCooDa(DataArrayDouble::New());
+ addCooDa->useArray(&addCoo[0],false,C_DEALLOC,(int)addCoo.size()/2,2);
+ // Step 2: re-order newly created nodes according to the ordering found in m2
+ std::vector< std::vector<int> > intersectEdge2;
+ BuildIntersectEdges(m1Desc,mesh1D,addCoo,subDiv2,intersectEdge2);
+ subDiv2.clear();
+ //
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret1(BuildMesh1DCutFrom(mesh1D,intersectEdge2,mesh2D->getCoords(),addCoo));
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> baryRet1(ret1->getBarycenterAndOwner());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> elts,eltsIndex;
+ mesh2D->getCellsContainingPoints(baryRet1->begin(),baryRet1->getNumberOfTuples(),eps,elts,eltsIndex);
+ splitMesh1D=ret1.retn();
+}
/**
* Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the
std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI, std::vector<int>& cNb1, std::vector<int>& cNb2)
{
static const int SPACEDIM=2;
- const double *coo1=m1->getCoords()->getConstPointer();
- const int *conn1=m1->getNodalConnectivity()->getConstPointer();
- const int *connI1=m1->getNodalConnectivityIndex()->getConstPointer();
- int offset1=m1->getNumberOfNodes();
- const double *coo2=m2->getCoords()->getConstPointer();
- const int *conn2=m2->getNodalConnectivity()->getConstPointer();
- const int *connI2=m2->getNodalConnectivityIndex()->getConstPointer();
- int offset2=offset1+m2->getNumberOfNodes();
- int offset3=offset2+((int)addCoords.size())/2;
+ const double *coo1(m1->getCoords()->getConstPointer());
+ const int *conn1(m1->getNodalConnectivity()->getConstPointer()),*connI1(m1->getNodalConnectivityIndex()->getConstPointer());
+ int offset1(m1->getNumberOfNodes());
+ const double *coo2(m2->getCoords()->getConstPointer());
+ const int *conn2(m2->getNodalConnectivity()->getConstPointer()),*connI2(m2->getNodalConnectivityIndex()->getConstPointer());
+ int offset2(offset1+m2->getNumberOfNodes());
+ int offset3(offset2+((int)addCoords.size())/2);
MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree()),bbox2Arr(m2->getBoundingBoxForBBTree());
const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
// Here a BBTree on 2D-cells, not on segments:
BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
- int ncell1=m1->getNumberOfCells();
+ int ncell1(m1->getNumberOfCells());
crI.push_back(0);
for(int i=0;i<ncell1;i++)
{
MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
// pol1 is the full cell from mesh2, in QP format, with all the additional intersecting nodes.
pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
- desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
+ desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
//
std::set<INTERP_KERNEL::Edge *> edges1;// store all edges of pol1 that are NOT consumed by intersect cells. If any after iteration over candidates2 -> a part of pol1 should appear in result
std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
// pol2 is the new QP in the final merged result.
pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
- pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
+ pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
}
ii=0;
for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
{
- pol1.initLocationsWithOther(pol2s[ii]);
+ INTERP_KERNEL::ComposedEdge::InitLocationsWithOther(pol1,pol2s[ii]);
pol2s[ii].updateLocOfEdgeFromCrudeDataArray2(desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
//MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
pol1.buildPartitionsAbs(pol2s[ii],edges1,edgesBoundary2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
if(!edges1.empty())
{
try
- {
+ {
INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
- }
+ }
catch(INTERP_KERNEL::Exception& e)
- {
+ {
std::ostringstream oss; oss << "Error when computing residual of cell #" << i << " in source/m1 mesh ! Maybe the neighbours of this cell in mesh are not well connected !\n" << "The deep reason is the following : " << e.what();
throw INTERP_KERNEL::Exception(oss.str().c_str());
- }
+ }
}
for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
(*it).second->decrRef();
}
}
+void IKGeo2DInternalMapper2(INTERP_KERNEL::Node *n, const std::map<INTERP_KERNEL::Node *,int>& m, int forbVal0, int forbVal1, std::vector<int>& isect)
+{
+ std::map<INTERP_KERNEL::Node *,int>::const_iterator it(m.find(n));
+ if(it==m.end())
+ throw INTERP_KERNEL::Exception("Internal error in remapping !");
+ int v((*it).second);
+ if(v==forbVal0 || v==forbVal1)
+ return ;
+ if(std::find(isect.begin(),isect.end(),v)==isect.end())
+ isect.push_back(v);
+}
+
+bool IKGeo2DInternalMapper(const INTERP_KERNEL::ComposedEdge& c, const std::map<INTERP_KERNEL::Node *,int>& m, int forbVal0, int forbVal1, std::vector<int>& isect)
+{
+ int sz(c.size());
+ if(sz<=1)
+ return false;
+ bool presenceOfOn(false);
+ for(int i=0;i<sz;i++)
+ {
+ INTERP_KERNEL::ElementaryEdge *e(c[i]);
+ if(e->getLoc()!=INTERP_KERNEL::FULL_ON_1)
+ continue ;
+ IKGeo2DInternalMapper2(e->getStartNode(),m,forbVal0,forbVal1,isect);
+ IKGeo2DInternalMapper2(e->getEndNode(),m,forbVal0,forbVal1,isect);
+ }
+ return presenceOfOn;
+}
+
+/**
+ * This method split some of edges of 2D cells in \a this. The edges to be split are specified in \a subNodesInSeg and in \a subNodesInSegI using index storage mode.
+ * To do the work this method can optionnaly needs information about middle of subedges for quadratic cases if a minimal creation of new nodes is wanted.
+ * So this method try to reduce at most the number of new nodes. The only case that can lead this method to add nodes if a SEG3 is split without information of middle.
+ * \b WARNING : is returned value is different from 0 a call to MEDCouplingUMesh::mergeNodes is necessary to avoid to have a non conform mesh.
+ *
+ * \return int - the number of new nodes created (in most of cases 0).
+ *
+ * \throw If \a this is not coherent.
+ * \throw If \a this has not spaceDim equal to 2.
+ * \throw If \a this has not meshDim equal to 2.
+ * \throw If some subcells needed to be split are orphan.
+ * \sa MEDCouplingUMesh::conformize2D
+ */
+int MEDCouplingUMesh::split2DCells(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *midOpt, const DataArrayInt *midOptI)
+{
+ if(!desc || !descI || !subNodesInSeg || !subNodesInSegI)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCells : the 4 first arrays must be not null !");
+ desc->checkAllocated(); descI->checkAllocated(); subNodesInSeg->checkAllocated(); subNodesInSegI->checkAllocated();
+ if(getSpaceDimension()!=2 || getMeshDimension()!=2)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCells : This method only works for meshes with spaceDim=2 and meshDim=2 !");
+ if(midOpt==0 && midOptI==0)
+ {
+ split2DCellsLinear(desc,descI,subNodesInSeg,subNodesInSegI);
+ return 0;
+ }
+ else if(midOpt!=0 && midOptI!=0)
+ return split2DCellsQuadratic(desc,descI,subNodesInSeg,subNodesInSegI,midOpt,midOptI);
+ else
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCells : middle parameters must be set to null for all or not null for all.");
+}
+
/*!
- * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2).
- * It builds the descending connectivity of the two meshes, and then using a binary tree
- * it computes the edge intersections. This results in new points being created : they're stored in addCoo.
- * Documentation about parameters colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs().
+ * \b WARNING this method is \b potentially \b non \b const (if returned array is 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. So if a edge in \a this can be split into entire edges in \a this this method
+ * will suppress such edges to use sub edges in \a this. So this method does not add nodes in \a this if merged edges are both linear (INTERP_KERNEL::NORM_SEG2).
+ * In the other cases new nodes can be created. If any are created, they will be appended at the end of the coordinates object before the invokation of this method.
+ *
+ * Whatever the returned value, this method does not alter the order of cells in \a this neither the orientation of cells.
+ * The modified cells, if any, are systematically declared as NORM_POLYGON or NORM_QPOLYG depending on the initial quadraticness of geometric type.
+ *
+ * This method expects that \b this has a meshDim equal 2 and spaceDim equal to 2 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 2.
+ * \throw If \a this has not meshDim equal to 2.
+ * \sa MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::split2DCells
*/
-void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
- std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
- MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
- std::vector<double>& addCoo,
- MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2)
- throw(INTERP_KERNEL::Exception)
+DataArrayInt *MEDCouplingUMesh::conformize2D(double eps)
{
static const int SPACEDIM=2;
- // Build desc connectivity
- desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New();
- desc2=DataArrayInt::New();
- descIndx2=DataArrayInt::New();
- revDesc2=DataArrayInt::New();
- revDescIndx2=DataArrayInt::New();
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
- m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
- m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
- const int *c1=m1Desc->getNodalConnectivity()->getConstPointer();
- const int *ci1=m1Desc->getNodalConnectivityIndex()->getConstPointer();
+ checkCoherency();
+ if(getSpaceDimension()!=2 || getMeshDimension()!=2)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> desc1(DataArrayInt::New()),descIndx1(DataArrayInt::New()),revDesc1(DataArrayInt::New()),revDescIndx1(DataArrayInt::New());
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> mDesc(buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1));
+ const int *c(mDesc->getNodalConnectivity()->getConstPointer()),*ci(mDesc->getNodalConnectivityIndex()->getConstPointer()),*rd(revDesc1->getConstPointer()),*rdi(revDescIndx1->getConstPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bboxArr(mDesc->getBoundingBoxForBBTree());
+ const double *bbox(bboxArr->begin()),*coords(getCoords()->begin());
+ int nCell(getNumberOfCells()),nDescCell(mDesc->getNumberOfCells());
+ std::vector< std::vector<int> > intersectEdge(nDescCell),overlapEdge(nDescCell);
+ std::vector<double> addCoo;
+ BBTree<SPACEDIM,int> myTree(bbox,0,0,nDescCell,-eps);
+ INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
+ INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
+ for(int i=0;i<nDescCell;i++)
+ {
+ std::vector<int> candidates;
+ myTree.getIntersectingElems(bbox+i*2*SPACEDIM,candidates);
+ for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
+ if(*it>i)
+ {
+ std::map<INTERP_KERNEL::Node *,int> m;
+ INTERP_KERNEL::Edge *e1(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m)),
+ *e2(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[*it]],c+ci[*it]+1,coords,m));
+ INTERP_KERNEL::MergePoints merge;
+ INTERP_KERNEL::QuadraticPolygon c1,c2;
+ e1->intersectWith(e2,merge,c1,c2);
+ e1->decrRef(); e2->decrRef();
+ if(IKGeo2DInternalMapper(c1,m,c[ci[i]+1],c[ci[i]+2],intersectEdge[i]))
+ overlapEdge[i].push_back(*it);
+ if(IKGeo2DInternalMapper(c2,m,c[ci[*it]+1],c[ci[*it]+2],intersectEdge[*it]))
+ overlapEdge[*it].push_back(i);
+ for(std::map<INTERP_KERNEL::Node *,int>::const_iterator it2=m.begin();it2!=m.end();it2++)
+ (*it2).first->decrRef();
+ }
+ }
+ // splitting done. sort intersect point in intersectEdge.
+ std::vector< std::vector<int> > middle(nDescCell);
+ int nbOf2DCellsToBeSplit(0);
+ bool middleNeedsToBeUsed(false);
+ std::vector<bool> cells2DToTreat(nDescCell,false);
+ for(int i=0;i<nDescCell;i++)
+ {
+ std::vector<int>& isect(intersectEdge[i]);
+ int sz((int)isect.size());
+ if(sz>1)
+ {
+ std::map<INTERP_KERNEL::Node *,int> m;
+ INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2((INTERP_KERNEL::NormalizedCellType)c[ci[i]],c+ci[i]+1,coords,m));
+ e->sortSubNodesAbs(coords,isect);
+ e->decrRef();
+ for(std::map<INTERP_KERNEL::Node *,int>::const_iterator it2=m.begin();it2!=m.end();it2++)
+ (*it2).first->decrRef();
+ }
+ if(sz!=0)
+ {
+ int idx0(rdi[i]),idx1(rdi[i+1]);
+ if(idx1-idx0!=1)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::conformize2D : internal error #0 !");
+ if(!cells2DToTreat[rd[idx0]])
+ {
+ cells2DToTreat[rd[idx0]]=true;
+ nbOf2DCellsToBeSplit++;
+ }
+ // try to reuse at most eventual 'middle' of SEG3
+ std::vector<int>& mid(middle[i]);
+ mid.resize(sz+1,-1);
+ if((INTERP_KERNEL::NormalizedCellType)c[ci[i]]==INTERP_KERNEL::NORM_SEG3)
+ {
+ middleNeedsToBeUsed=true;
+ const std::vector<int>& candidates(overlapEdge[i]);
+ std::vector<int> trueCandidates;
+ for(std::vector<int>::const_iterator itc=candidates.begin();itc!=candidates.end();itc++)
+ if((INTERP_KERNEL::NormalizedCellType)c[ci[*itc]]==INTERP_KERNEL::NORM_SEG3)
+ trueCandidates.push_back(*itc);
+ int stNode(c[ci[i]+1]),endNode(isect[0]);
+ for(int j=0;j<sz+1;j++)
+ {
+ for(std::vector<int>::const_iterator itc=trueCandidates.begin();itc!=trueCandidates.end();itc++)
+ {
+ int tmpSt(c[ci[*itc]+1]),tmpEnd(c[ci[*itc]+2]);
+ if((tmpSt==stNode && tmpEnd==endNode) || (tmpSt==endNode && tmpEnd==stNode))
+ { mid[j]=*itc; break; }
+ }
+ stNode=endNode;
+ endNode=j<sz-1?isect[j+1]:c[ci[i]+2];
+ }
+ }
+ }
+ }
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New()),notRet(DataArrayInt::New()); ret->alloc(nbOf2DCellsToBeSplit,1);
+ if(nbOf2DCellsToBeSplit==0)
+ return ret.retn();
+ //
+ int *retPtr(ret->getPointer());
+ for(int i=0;i<nCell;i++)
+ if(cells2DToTreat[i])
+ *retPtr++=i;
+ //
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> mSafe,nSafe,oSafe,pSafe,qSafe,rSafe;
+ DataArrayInt *m(0),*n(0),*o(0),*p(0),*q(0),*r(0);
+ MEDCouplingUMesh::ExtractFromIndexedArrays(ret->begin(),ret->end(),desc1,descIndx1,m,n); mSafe=m; nSafe=n;
+ DataArrayInt::PutIntoToSkylineFrmt(intersectEdge,o,p); oSafe=o; pSafe=p;
+ if(middleNeedsToBeUsed)
+ { DataArrayInt::PutIntoToSkylineFrmt(middle,q,r); qSafe=q; rSafe=r; }
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> modif(static_cast<MEDCouplingUMesh *>(buildPartOfMySelf(ret->begin(),ret->end(),true)));
+ int nbOfNodesCreated(modif->split2DCells(mSafe,nSafe,oSafe,pSafe,qSafe,rSafe));
+ setCoords(modif->getCoords());//if nbOfNodesCreated==0 modif and this have the same coordinates pointer so this line has no effect. But for quadratic cases this line is important.
+ setPartOfMySelf(ret->begin(),ret->end(),*modif);
+ {
+ bool areNodesMerged; int newNbOfNodes;
+ if(nbOfNodesCreated!=0)
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp(mergeNodes(eps,areNodesMerged,newNbOfNodes));
+ }
+ return ret.retn();
+}
+
+/*!
+ * This non const method works on 2D mesh. This method scans every cell in \a this and look if each edge constituting this cell is not mergeable with neighbors edges of that cell.
+ * If yes, the cell is "repaired" to minimize at most its number of edges. So this method do not change the overall shape of cells in \a this (with eps precision).
+ * This method do not take care of shared edges between cells, so this method can lead to a non conform mesh (\a this). If a conform mesh is required you're expected
+ * to invoke MEDCouplingUMesh::mergeNodes and MEDCouplingUMesh::conformize2D right after this call.
+ * This method works on any 2D geometric types of cell (even static one). If a cell is touched its type becomes dynamic automaticaly. For 2D "repaired" quadratic cells
+ * new nodes for center of merged edges is are systematically created and appended at the end of the previously existing nodes.
+ *
+ * If the returned array is empty it means that nothing has changed in \a this (as if it were a const method). If the array is not empty the connectivity of \a this is modified
+ * using new instance, idem for coordinates.
+ *
+ * If \a this is constituted by only linear 2D cells, this method is close to the computation of the convex hull of each cells in \a this.
+ *
+ * \return DataArrayInt * - The list of cellIds in \a this that have at least one edge colinearized.
+ *
+ * \throw If \a this is not coherent.
+ * \throw If \a this has not spaceDim equal to 2.
+ * \throw If \a this has not meshDim equal to 2.
+ *
+ * \sa MEDCouplingUMesh::conformize2D, MEDCouplingUMesh::mergeNodes, MEDCouplingUMesh::convexEnvelop2D.
+ */
+DataArrayInt *MEDCouplingUMesh::colinearize2D(double eps)
+{
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
+ checkCoherency();
+ if(getSpaceDimension()!=2 || getMeshDimension()!=2)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::colinearize2D : This method only works for meshes with spaceDim=2 and meshDim=2 !");
+ INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
+ INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
+ int nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
+ const int *cptr(_nodal_connec->begin()),*ciptr(_nodal_connec_index->begin());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newc(DataArrayInt::New()),newci(DataArrayInt::New()); newci->alloc(nbOfCells+1,1); newc->alloc(0,1); newci->setIJ(0,0,0);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> appendedCoords(DataArrayDouble::New()); appendedCoords->alloc(0,1);//1 not 2 it is not a bug.
+ const double *coords(_coords->begin());
+ int *newciptr(newci->getPointer());
+ for(int i=0;i<nbOfCells;i++,newciptr++,ciptr++)
+ {
+ if(Colinearize2DCell(coords,cptr+ciptr[0],cptr+ciptr[1],nbOfNodes,newc,appendedCoords))
+ ret->pushBackSilent(i);
+ newciptr[1]=newc->getNumberOfTuples();
+ }
+ //
+ if(ret->empty())
+ return ret.retn();
+ if(!appendedCoords->empty())
+ {
+ appendedCoords->rearrange(2);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords(DataArrayDouble::Aggregate(getCoords(),appendedCoords));//treat info on components
+ //non const part
+ setCoords(newCoords);
+ }
+ //non const part
+ setConnectivity(newc,newci,true);
+ return ret.retn();
+}
+/*!
+ * \param [out] intersectEdge1 - for each cell in \a m1Desc returns the result of the split. The result is given using pair of int given resp start and stop.
+ * So for all edge \a i in \a m1Desc \a intersectEdge1[i] is of length 2*n where n is the number of sub edges.
+ * And for each j in [1,n) intersect[i][2*(j-1)+1]==intersect[i][2*j].
+ * \param [out] subDiv2 - for each cell in \a m2Desc returns nodes that split it using convention \a m1Desc first, then \a m2Desc, then addCoo
+ * \param [out] colinear2 - for each cell in \a m2Desc returns the edges in \a m1Desc that are colinear to it.
+ * \param [out] addCoo - nodes to be append at the end
+ */
+void MEDCouplingUMesh::Intersect1DMeshes(const MEDCouplingUMesh *m1Desc, const MEDCouplingUMesh *m2Desc, double eps,
+ std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2, std::vector<double>& addCoo)
+{
+ static const int SPACEDIM=2;
+ INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
+ INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
+ const int *c1(m1Desc->getNodalConnectivity()->getConstPointer()),*ci1(m1Desc->getNodalConnectivityIndex()->getConstPointer());
// Build BB tree of all edges in the tool mesh (second mesh)
MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree()),bbox2Arr(m2Desc->getBoundingBoxForBBTree());
const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
- int nDescCell1=m1Desc->getNumberOfCells();
- int nDescCell2=m2Desc->getNumberOfCells();
+ int nDescCell1(m1Desc->getNumberOfCells()),nDescCell2(m2Desc->getNumberOfCells());
intersectEdge1.resize(nDescCell1);
colinear2.resize(nDescCell2);
subDiv2.resize(nDescCell2);
BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps);
std::vector<int> candidates1(1);
- int offset1=m1->getNumberOfNodes();
- int offset2=offset1+m2->getNumberOfNodes();
+ int offset1(m1Desc->getNumberOfNodes());
+ int offset2(offset1+m2Desc->getNumberOfNodes());
for(int i=0;i<nDescCell1;i++) // for all edges in the first mesh
{
std::vector<int> candidates2; // edges of mesh2 candidate for intersection
else
intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i+1]);
}
+}
+
+/*!
+ * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2).
+ * It builds the descending connectivity of the two meshes, and then using a binary tree
+ * it computes the edge intersections. This results in new points being created : they're stored in addCoo.
+ * Documentation about parameters colinear2 and subDiv2 can be found in method QuadraticPolygon::splitAbs().
+ */
+void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
+ std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
+ MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
+ std::vector<double>& addCoo,
+ MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2)
+{
+ // Build desc connectivity
+ desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New();
+ desc2=DataArrayInt::New();
+ descIndx2=DataArrayInt::New();
+ revDesc2=DataArrayInt::New();
+ revDescIndx2=DataArrayInt::New();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
+ m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
+ m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
+ Intersect1DMeshes(m1Desc,m2Desc,eps,intersectEdge1,colinear2,subDiv2,addCoo);
m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef();
m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef();
}
* \param[out] intersectEdge the same content as subDiv, but correclty oriented.
*/
void MEDCouplingUMesh::BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
- const std::vector<double>& addCoo,
- const std::vector< std::vector<int> >& subDiv, std::vector< std::vector<int> >& intersectEdge)
+ const std::vector<double>& addCoo,
+ const std::vector< std::vector<int> >& subDiv, std::vector< std::vector<int> >& intersectEdge)
{
int offset1=m1->getNumberOfNodes();
int ncell=m2->getNumberOfCells();
void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector<int>& cut3DCurve, std::vector<int>& nodesOnPlane, const int *nodal3DSurf, const int *nodalIndx3DSurf,
const int *nodal3DCurve, const int *nodalIndx3DCurve,
const int *desc, const int *descIndx,
- std::vector< std::pair<int,int> >& cut3DSurf) throw(INTERP_KERNEL::Exception)
+ std::vector< std::pair<int,int> >& cut3DSurf)
{
std::set<int> nodesOnP(nodesOnPlane.begin(),nodesOnPlane.end());
int nbOf3DSurfCell=(int)cut3DSurf.size();
}
}
switch(res.size())
- {
+ {
case 2:
{
cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
else
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::AssemblyPointsFrom3DCurve : unexpected situation !");
}
- }
+ }
}
}
*/
void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair<int,int> >& cut3DSurf,
const int *desc, const int *descIndx,
- DataArrayInt *nodalRes, DataArrayInt *nodalResIndx, DataArrayInt *cellIds) const throw(INTERP_KERNEL::Exception)
+ DataArrayInt *nodalRes, DataArrayInt *nodalResIndx, DataArrayInt *cellIds) const
{
checkFullyDefined();
if(getMeshDimension()!=3 || getSpaceDimension()!=3)
* \sa MEDCouplingUMesh::ExtractFromIndexedArrays2
*/
void MEDCouplingUMesh::ExtractFromIndexedArrays(const int *idsOfSelectBg, const int *idsOfSelectEnd, const DataArrayInt *arrIn, const DataArrayInt *arrIndxIn,
- DataArrayInt* &arrOut, DataArrayInt* &arrIndexOut) throw(INTERP_KERNEL::Exception)
+ DataArrayInt* &arrOut, DataArrayInt* &arrIndexOut)
{
if(!arrIn || !arrIndxIn)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ExtractFromIndexedArrays : input pointer is NULL !");
* \sa MEDCouplingUMesh::ExtractFromIndexedArrays
*/
void MEDCouplingUMesh::ExtractFromIndexedArrays2(int idsOfSelectStart, int idsOfSelectStop, int idsOfSelectStep, const DataArrayInt *arrIn, const DataArrayInt *arrIndxIn,
- DataArrayInt* &arrOut, DataArrayInt* &arrIndexOut) throw(INTERP_KERNEL::Exception)
+ DataArrayInt* &arrOut, DataArrayInt* &arrIndexOut)
{
if(!arrIn || !arrIndxIn)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ExtractFromIndexedArrays2 : input pointer is NULL !");
*/
void MEDCouplingUMesh::SetPartOfIndexedArrays(const int *idsOfSelectBg, const int *idsOfSelectEnd, const DataArrayInt *arrIn, const DataArrayInt *arrIndxIn,
const DataArrayInt *srcArr, const DataArrayInt *srcArrIndex,
- DataArrayInt* &arrOut, DataArrayInt* &arrIndexOut) throw(INTERP_KERNEL::Exception)
+ DataArrayInt* &arrOut, DataArrayInt* &arrIndexOut)
{
if(arrIn==0 || arrIndxIn==0 || srcArr==0 || srcArrIndex==0)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::SetPartOfIndexedArrays : presence of null pointer in input parameter !");
* \sa MEDCouplingUMesh::SetPartOfIndexedArrays
*/
void MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx(const int *idsOfSelectBg, const int *idsOfSelectEnd, DataArrayInt *arrInOut, const DataArrayInt *arrIndxIn,
- const DataArrayInt *srcArr, const DataArrayInt *srcArrIndex) throw(INTERP_KERNEL::Exception)
+ const DataArrayInt *srcArr, const DataArrayInt *srcArrIndex)
{
if(arrInOut==0 || arrIndxIn==0 || srcArr==0 || srcArrIndex==0)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx : presence of null pointer in input parameter !");
*/
void MEDCouplingUMesh::SetPartOfIndexedArrays2(int start, int end, int step, const DataArrayInt *arrIn, const DataArrayInt *arrIndxIn,
const DataArrayInt *srcArr, const DataArrayInt *srcArrIndex,
- DataArrayInt* &arrOut, DataArrayInt* &arrIndexOut) throw(INTERP_KERNEL::Exception)
+ DataArrayInt* &arrOut, DataArrayInt* &arrIndexOut)
{
if(arrIn==0 || arrIndxIn==0 || srcArr==0 || srcArrIndex==0)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::SetPartOfIndexedArrays2 : presence of null pointer in input parameter !");
* \sa MEDCouplingUMesh::SetPartOfIndexedArrays2 MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx
*/
void MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx2(int start, int end, int step, DataArrayInt *arrInOut, const DataArrayInt *arrIndxIn,
- const DataArrayInt *srcArr, const DataArrayInt *srcArrIndex) throw(INTERP_KERNEL::Exception)
+ const DataArrayInt *srcArr, const DataArrayInt *srcArrIndex)
{
if(arrInOut==0 || arrIndxIn==0 || srcArr==0 || srcArrIndex==0)
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx2 : presence of null pointer in input parameter !");
MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> tmp=static_cast<MEDCouplingUMesh *>(buildPartOfMySelf((*it)->begin(),(*it)->end(),true));
MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cell;
switch(mdim)
- {
+ {
case 2:
cell=tmp->buildUnionOf2DMesh();
break;
break;
default:
throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSpreadZonesWithPoly : meshdimension supported are [2,3] ! Not implemented yet for others !");
- }
-
+ }
+
ret->insertNextCell((INTERP_KERNEL::NormalizedCellType)cell->getIJSafe(0,0),cell->getNumberOfTuples()-1,cell->getConstPointer()+1);
}
//
return ret0.retn();
}
+/*!
+ * It is the linear part of MEDCouplingUMesh::split2DCells. Here no additionnal nodes will be added in \b this. So coordinates pointer remain unchanged (is not even touch).
+ *
+ * \sa MEDCouplingUMesh::split2DCells
+ */
+void MEDCouplingUMesh::split2DCellsLinear(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI)
+{
+ checkConnectivityFullyDefined();
+ int ncells(getNumberOfCells()),lgthToReach(getMeshLength()+subNodesInSeg->getNumberOfTuples());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
+ const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
+ int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
+ int prevPosOfCi(ciPtr[0]);
+ for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
+ {
+ int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(0);
+ *cPtr++=(int)INTERP_KERNEL::NORM_POLYGON; *cPtr++=oldConn[prevPosOfCi+1];
+ for(int j=0;j<sz;j++)
+ {
+ int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]);
+ for(int k=0;k<sz2;k++)
+ *cPtr++=subPtr[offset2+k];
+ if(j!=sz-1)
+ *cPtr++=oldConn[prevPosOfCi+j+2];
+ deltaSz+=sz2;
+ }
+ prevPosOfCi=ciPtr[1];
+ ciPtr[1]=ciPtr[0]+1+sz+deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
+ }
+ if(c->end()!=cPtr)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsLinear : Some of edges to be split are orphan !");
+ _nodal_connec->decrRef();
+ _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_POLYGON);
+}
+
+int InternalAddPoint(const INTERP_KERNEL::Edge *e, int id, const double *coo, int startId, int endId, DataArrayDouble& addCoo, int& nodesCnter)
+{
+ if(id!=-1)
+ return id;
+ else
+ {
+ int ret(nodesCnter++);
+ double newPt[2];
+ e->getMiddleOfPoints(coo+2*startId,coo+2*endId,newPt);
+ addCoo.insertAtTheEnd(newPt,newPt+2);
+ return ret;
+ }
+}
+
+/// @cond INTERNAL
+
+void EnterTheResultOf2DCellFirst(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<int>& middles)
+{
+ int tmp[3];
+ int trueStart(start>=0?start:nbOfEdges+start);
+ tmp[0]=linOrArc?(int)INTERP_KERNEL::NORM_QPOLYG:(int)INTERP_KERNEL::NORM_POLYGON; tmp[1]=connBg[trueStart]; tmp[2]=connBg[stp];
+ newConnOfCell->insertAtTheEnd(tmp,tmp+3);
+ if(linOrArc)
+ {
+ if(stp-start>1)
+ {
+ int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
+ InternalAddPoint(e,-1,coords,tmp[1],tmp[2],*appendedCoords,tmp2);
+ middles.push_back(tmp3+offset);
+ }
+ else
+ middles.push_back(connBg[trueStart+nbOfEdges]);
+ }
+}
+
+void EnterTheResultOf2DCellMiddle(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<int>& middles)
+{
+ int tmpSrt(newConnOfCell->back()),tmpEnd(connBg[stp]);
+ newConnOfCell->pushBackSilent(tmpEnd);
+ if(linOrArc)
+ {
+ if(stp-start>1)
+ {
+ int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
+ InternalAddPoint(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2);
+ middles.push_back(tmp3+offset);
+ }
+ else
+ middles.push_back(connBg[start+nbOfEdges]);
+ }
+}
+
+void EnterTheResultOf2DCellEnd(const INTERP_KERNEL::Edge *e, int start, int stp, int nbOfEdges, bool linOrArc, const double *coords, const int *connBg, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords, std::vector<int>& middles)
+{
+ if(linOrArc)
+ {
+ if(stp-start>1)
+ {
+ int tmpSrt(connBg[start]),tmpEnd(connBg[stp]);
+ int tmp2(0),tmp3(appendedCoords->getNumberOfTuples()/2);
+ InternalAddPoint(e,-1,coords,tmpSrt,tmpEnd,*appendedCoords,tmp2);
+ middles.push_back(tmp3+offset);
+ }
+ else
+ middles.push_back(connBg[start+nbOfEdges]);
+ }
+}
+
+/// @cond INTERNAL
+
+/*!
+ * Returns true if a colinearization has been found in the given cell. If false is returned the content pushed in \a newConnOfCell is equal to [ \a connBg , \a connEnd ) .
+ * \a appendedCoords is a DataArrayDouble instance with number of components equal to one (even if the items are pushed by pair).
+ */
+bool MEDCouplingUMesh::Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords)
+{
+ std::size_t sz(std::distance(connBg,connEnd));
+ if(sz<3)//3 because 2+1(for the cell type) and 2 is the minimal number of edges of 2D cell.
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Colinearize2DCell : the input cell has invalid format !");
+ sz--;
+ INTERP_KERNEL::AutoPtr<int> tmpConn(new int[sz]);
+ const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connBg[0]));
+ unsigned nbs(cm.getNumberOfSons2(connBg+1,sz)),nbOfHit(0);
+ int posBaseElt(0),posEndElt(0),nbOfTurn(0);
+ INTERP_KERNEL::NormalizedCellType typeOfSon;
+ std::vector<int> middles;
+ bool ret(false);
+ for(;nbOfHit<nbs;nbOfTurn++)
+ {
+ cm.fillSonCellNodalConnectivity2(posBaseElt,connBg+1,sz,tmpConn,typeOfSon);
+ std::map<INTERP_KERNEL::Node *,int> m;
+ INTERP_KERNEL::Edge *e(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
+ posEndElt++;
+ nbOfHit++;
+ unsigned endI(nbs-nbOfHit);
+ for(unsigned i=0;i<endI;i++)
+ {
+ cm.fillSonCellNodalConnectivity2(posBaseElt+(int)i+1,connBg+1,sz,tmpConn,typeOfSon);
+ INTERP_KERNEL::Edge *eCand(MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m));
+ INTERP_KERNEL::EdgeIntersector *eint(INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand));
+ bool isColinear(eint->areColinears());
+ if(isColinear)
+ {
+ nbOfHit++;
+ posEndElt++;
+ ret=true;
+ }
+ delete eint;
+ eCand->decrRef();
+ if(!isColinear)
+ {
+ if(nbOfTurn==0)
+ {//look if the first edge of cell is not colinear with last edges in this case the start of nodal connectivity is shifted back
+ unsigned endII(nbs-nbOfHit-1);//warning nbOfHit can be modified, so put end condition in a variable.
+ for(unsigned ii=0;ii<endII;ii++)
+ {
+ cm.fillSonCellNodalConnectivity2(nbs-ii-1,connBg+1,sz,tmpConn,typeOfSon);
+ eCand=MEDCouplingUMeshBuildQPFromEdge2(typeOfSon,tmpConn,coords,m);
+ eint=INTERP_KERNEL::Edge::BuildIntersectorWith(e,eCand);
+ isColinear=eint->areColinears();
+ if(isColinear)
+ {
+ nbOfHit++;
+ posBaseElt--;
+ ret=true;
+ }
+ delete eint;
+ eCand->decrRef();
+ if(!isColinear)
+ break;
+ }
+ }
+ break;
+ }
+ }
+ //push [posBaseElt,posEndElt) in newConnOfCell using e
+ if(nbOfTurn==0)
+ EnterTheResultOf2DCellFirst(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
+ else if(nbOfHit!=nbs)
+ EnterTheResultOf2DCellMiddle(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
+ else
+ EnterTheResultOf2DCellEnd(e,posBaseElt,posEndElt,(int)nbs,cm.isQuadratic(),coords,connBg+1,offset,newConnOfCell,appendedCoords,middles);
+ posBaseElt=posEndElt;
+ for(std::map<INTERP_KERNEL::Node *,int>::const_iterator it=m.begin();it!=m.end();it++)
+ (*it).first->decrRef();
+ e->decrRef();
+ }
+ if(!middles.empty())
+ newConnOfCell->insertAtTheEnd(middles.begin(),middles.end());
+ return ret;
+}
+
+/*!
+ * It is the quadratic part of MEDCouplingUMesh::split2DCells. Here some additionnal nodes can be added at the end of coordinates array object.
+ *
+ * \return int - the number of new nodes created.
+ * \sa MEDCouplingUMesh::split2DCells
+ */
+int MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI)
+{
+ checkCoherency();
+ int ncells(getNumberOfCells()),lgthToReach(getMeshLength()+2*subNodesInSeg->getNumberOfTuples()),nodesCnt(getNumberOfNodes());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c(DataArrayInt::New()); c->alloc((std::size_t)lgthToReach);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> addCoo(DataArrayDouble::New()); addCoo->alloc(0,1);
+ const int *subPtr(subNodesInSeg->begin()),*subIPtr(subNodesInSegI->begin()),*descPtr(desc->begin()),*descIPtr(descI->begin()),*oldConn(getNodalConnectivity()->begin());
+ const int *midPtr(mid->begin()),*midIPtr(midI->begin());
+ const double *oldCoordsPtr(getCoords()->begin());
+ int *cPtr(c->getPointer()),*ciPtr(getNodalConnectivityIndex()->getPointer());
+ int prevPosOfCi(ciPtr[0]);
+ for(int i=0;i<ncells;i++,ciPtr++,descIPtr++)
+ {
+ int offset(descIPtr[0]),sz(descIPtr[1]-descIPtr[0]),deltaSz(sz);
+ for(int j=0;j<sz;j++)
+ { int sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]); deltaSz+=sz2; }
+ *cPtr++=(int)INTERP_KERNEL::NORM_QPOLYG; cPtr[0]=oldConn[prevPosOfCi+1];
+ for(int j=0;j<sz;j++)//loop over subedges of oldConn
+ {
+ int offset2(subIPtr[descPtr[offset+j]]),sz2(subIPtr[descPtr[offset+j]+1]-subIPtr[descPtr[offset+j]]),offset3(midIPtr[descPtr[offset+j]]);
+ if(sz2==0)
+ {
+ if(j<sz-1)
+ cPtr[1]=oldConn[prevPosOfCi+2+j];
+ cPtr[deltaSz]=oldConn[prevPosOfCi+1+j+sz]; cPtr++;
+ continue;
+ }
+ std::vector<INTERP_KERNEL::Node *> ns(3);
+ ns[0]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+j]+1]);
+ ns[1]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+(1+j)%sz]+1]);
+ ns[2]=new INTERP_KERNEL::Node(oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]],oldCoordsPtr[2*oldConn[prevPosOfCi+1+sz+j]+1]);
+ MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Edge> e(INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(ns));
+ for(int k=0;k<sz2;k++)//loop over subsplit of current subedge
+ {
+ cPtr[1]=subPtr[offset2+k];
+ cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+k],oldCoordsPtr,cPtr[0],cPtr[1],*addCoo,nodesCnt); cPtr++;
+ }
+ int tmpEnd(oldConn[prevPosOfCi+1+(j+1)%sz]);
+ if(j!=sz-1)
+ { cPtr[1]=tmpEnd; }
+ cPtr[deltaSz]=InternalAddPoint(e,midPtr[offset3+sz2],oldCoordsPtr,cPtr[0],tmpEnd,*addCoo,nodesCnt); cPtr++;
+ }
+ prevPosOfCi=ciPtr[1]; cPtr+=deltaSz;
+ ciPtr[1]=ciPtr[0]+1+2*deltaSz;//sz==old nb of nodes because (nb of subedges=nb of nodes for polygons)
+ }
+ if(c->end()!=cPtr)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split2DCellsQuadratic : Some of edges to be split are orphan !");
+ _nodal_connec->decrRef();
+ _nodal_connec=c.retn(); _types.clear(); _types.insert(INTERP_KERNEL::NORM_QPOLYG);
+ addCoo->rearrange(2);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo(DataArrayDouble::Aggregate(getCoords(),addCoo));//info are copied from getCoords() by using Aggregate
+ setCoords(coo);
+ return addCoo->getNumberOfTuples();
+}
+
MEDCouplingUMeshCellIterator::MEDCouplingUMeshCellIterator(MEDCouplingUMesh *mesh):_mesh(mesh),_cell(new MEDCouplingUMeshCell(mesh)),
- _own_cell(true),_cell_id(-1),_nb_cell(0)
+ _own_cell(true),_cell_id(-1),_nb_cell(0)
{
if(mesh)
{
}
MEDCouplingUMeshCellIterator::MEDCouplingUMeshCellIterator(MEDCouplingUMesh *mesh, MEDCouplingUMeshCell *itc, int bg, int end):_mesh(mesh),_cell(itc),
- _own_cell(false),_cell_id(bg-1),
- _nb_cell(end)
+ _own_cell(false),_cell_id(bg-1),
+ _nb_cell(end)
{
if(mesh)
mesh->incrRef();
}
MEDCouplingUMeshCellEntry::MEDCouplingUMeshCellEntry(MEDCouplingUMesh *mesh, INTERP_KERNEL::NormalizedCellType type, MEDCouplingUMeshCell *itc, int bg, int end):_mesh(mesh),_type(type),
- _itc(itc),
- _bg(bg),_end(end)
+ _itc(itc),
+ _bg(bg),_end(end)
{
if(_mesh)
_mesh->incrRef();
