Back to index

salome-med  6.5.0
MEDCouplingUMesh.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "MEDCouplingUMesh.hxx"
00021 #include "MEDCouplingMemArray.txx"
00022 #include "MEDCouplingFieldDouble.hxx"
00023 #include "CellModel.hxx"
00024 #include "VolSurfUser.txx"
00025 #include "InterpolationUtils.hxx"
00026 #include "PointLocatorAlgos.txx"
00027 #include "BBTree.txx"
00028 #include "DirectedBoundingBox.hxx"
00029 #include "InterpKernelMeshQuality.hxx"
00030 #include "InterpKernelCellSimplify.hxx"
00031 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
00032 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
00033 #include "InterpKernelAutoPtr.hxx"
00034 #include "InterpKernelGeo2DNode.hxx"
00035 #include "InterpKernelGeo2DEdgeLin.hxx"
00036 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
00037 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
00038 
00039 #include <sstream>
00040 #include <fstream>
00041 #include <numeric>
00042 #include <cstring>
00043 #include <limits>
00044 #include <list>
00045 
00046 using namespace ParaMEDMEM;
00047 
00048 const char MEDCouplingUMesh::PART_OF_NAME[]="PartOf_";
00049 
00050 double MEDCouplingUMesh::EPS_FOR_POLYH_ORIENTATION=1.e-14;
00051 
00052 MEDCouplingUMesh *MEDCouplingUMesh::New()
00053 {
00054   return new MEDCouplingUMesh;
00055 }
00056 
00057 MEDCouplingUMesh *MEDCouplingUMesh::New(const char *meshName, int meshDim)
00058 {
00059   MEDCouplingUMesh *ret=new MEDCouplingUMesh;
00060   ret->setName(meshName);
00061   ret->setMeshDimension(meshDim);
00062   return ret;
00063 }
00064 
00065 MEDCouplingMesh *MEDCouplingUMesh::deepCpy() const
00066 {
00067   return clone(true);
00068 }
00069 
00070 MEDCouplingUMesh *MEDCouplingUMesh::clone(bool recDeepCpy) const
00071 {
00072   return new MEDCouplingUMesh(*this,recDeepCpy);
00073 }
00074 
00075 void MEDCouplingUMesh::updateTime() const
00076 {
00077   MEDCouplingPointSet::updateTime();
00078   if(_nodal_connec)
00079     {
00080       updateTimeWith(*_nodal_connec);
00081     }
00082   if(_nodal_connec_index)
00083     {
00084       updateTimeWith(*_nodal_connec_index);
00085     }
00086 }
00087 
00088 MEDCouplingUMesh::MEDCouplingUMesh():_iterator(-1),_mesh_dim(-2),
00089                                      _nodal_connec(0),_nodal_connec_index(0)
00090 {
00091 }
00092 
00099 void MEDCouplingUMesh::checkCoherency() const throw(INTERP_KERNEL::Exception)
00100 {
00101   if(_mesh_dim<-1)
00102     throw INTERP_KERNEL::Exception("No mesh dimension specified !");
00103   for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator iter=_types.begin();iter!=_types.end();iter++)
00104     {
00105       if((int)INTERP_KERNEL::CellModel::GetCellModel(*iter).getDimension()!=_mesh_dim)
00106         {
00107           std::ostringstream message;
00108           message << "Mesh invalid because dimension is " << _mesh_dim << " and there is presence of cell(s) with type " << (*iter);
00109           throw INTERP_KERNEL::Exception(message.str().c_str());
00110         }
00111     }
00112   if(_nodal_connec)
00113     {
00114       if(_nodal_connec->getNumberOfComponents()!=1)
00115         throw INTERP_KERNEL::Exception("Nodal connectivity array is expected to be with number of components set to one !");
00116       if(_nodal_connec->getInfoOnComponent(0)!="")
00117         throw INTERP_KERNEL::Exception("Nodal connectivity array is expected to have no info on its single component !");
00118     }
00119   if(_nodal_connec_index)
00120     {
00121       if(_nodal_connec_index->getNumberOfComponents()!=1)
00122         throw INTERP_KERNEL::Exception("Nodal connectivity index array is expected to be with number of components set to one !");
00123       if(_nodal_connec_index->getInfoOnComponent(0)!="")
00124         throw INTERP_KERNEL::Exception("Nodal connectivity index array is expected to have no info on its single component !");
00125     }
00126 }
00127 
00133 void MEDCouplingUMesh::checkCoherency1(double eps) const throw(INTERP_KERNEL::Exception)
00134 {
00135   checkCoherency();
00136   if(_mesh_dim==-1)
00137     return ;
00138   int meshDim=getMeshDimension();
00139   int nbOfNodes=getNumberOfNodes();
00140   int nbOfCells=getNumberOfCells();
00141   const int *ptr=_nodal_connec->getConstPointer();
00142   const int *ptrI=_nodal_connec_index->getConstPointer();
00143   for(int i=0;i<nbOfCells;i++)
00144     {
00145       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)ptr[ptrI[i]]);
00146       if((int)cm.getDimension()!=meshDim)
00147         {
00148           std::ostringstream oss;
00149           oss << "MEDCouplingUMesh::checkCoherency1 : cell << #" << i<< " with type Type " << cm.getRepr() << " in 'this' whereas meshdim == " << meshDim << " !";
00150           throw INTERP_KERNEL::Exception(oss.str().c_str());
00151         }
00152       int nbOfNodesInCell=ptrI[i+1]-ptrI[i]-1;
00153       if(!cm.isDynamic())
00154         if(nbOfNodesInCell!=(int)cm.getNumberOfNodes())
00155           {
00156             std::ostringstream oss;
00157             oss << "MEDCouplingUMesh::checkCoherency1 : cell #" << i << " with static Type '" << cm.getRepr() << "' has " <<  cm.getNumberOfNodes();
00158             oss << " nodes whereas in connectivity there is " << nbOfNodesInCell << " nodes ! Looks very bad !";
00159             throw INTERP_KERNEL::Exception(oss.str().c_str());
00160           }
00161       for(const int *w=ptr+ptrI[i]+1;w!=ptr+ptrI[i+1];w++)
00162         {
00163           int nodeId=*w;
00164           if(nodeId>=0)
00165             {
00166               if(nodeId>=nbOfNodes)
00167                 {
00168                   std::ostringstream oss; oss << "Cell #" << i << " is consituted of node #" << nodeId << " whereas there are only " << nbOfNodes << " nodes !";
00169                   throw INTERP_KERNEL::Exception(oss.str().c_str());
00170                 }
00171             }
00172           else if(nodeId<-1)
00173             {
00174               std::ostringstream oss; oss << "Cell #" << i << " is consituted of node #" << nodeId << " in connectivity ! sounds bad !";
00175               throw INTERP_KERNEL::Exception(oss.str().c_str());
00176             }
00177           else
00178             {
00179               if((INTERP_KERNEL::NormalizedCellType)(ptr[ptrI[i]])!=INTERP_KERNEL::NORM_POLYHED)
00180                 {
00181                   std::ostringstream oss; oss << "Cell #" << i << " is consituted of node #-1 in connectivity ! sounds bad !";
00182                   throw INTERP_KERNEL::Exception(oss.str().c_str());
00183                 }
00184             }
00185         }
00186     }
00187 }
00188 
00189 void MEDCouplingUMesh::checkCoherency2(double eps) const throw(INTERP_KERNEL::Exception)
00190 {
00191   checkCoherency1(eps);
00192 }
00193 
00194 void MEDCouplingUMesh::setMeshDimension(int meshDim)
00195 {
00196   if(meshDim<-1)
00197     throw INTERP_KERNEL::Exception("Invalid meshDim specified ! Must be greater or equal to -1 !");
00198   _mesh_dim=meshDim;
00199   declareAsNew();
00200 }
00201 
00202 void MEDCouplingUMesh::allocateCells(int nbOfCells)
00203 {
00204   if(_nodal_connec_index)
00205     {
00206       _nodal_connec_index->decrRef();
00207     }
00208   if(_nodal_connec)
00209     {
00210       _nodal_connec->decrRef();
00211     }
00212 
00213   _nodal_connec_index=DataArrayInt::New();
00214   _nodal_connec_index->alloc(nbOfCells+1,1);
00215   int *pt=_nodal_connec_index->getPointer();
00216   pt[0]=0;
00217   _nodal_connec=DataArrayInt::New();
00218   _nodal_connec->alloc(2*nbOfCells,1);
00219   _iterator=0;
00220   _types.clear();
00221   declareAsNew();
00222 }
00223 
00230 void MEDCouplingUMesh::insertNextCell(INTERP_KERNEL::NormalizedCellType type, int size, const int *nodalConnOfCell) throw(INTERP_KERNEL::Exception)
00231 {
00232   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
00233   if(_nodal_connec_index==0)
00234     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::insertNextCell : nodal connectivity not set ! invoke allocateCells before calling insertNextCell !");
00235   if((int)cm.getDimension()==_mesh_dim)
00236     {
00237       int nbOfElems=_nodal_connec_index->getNbOfElems()-1;
00238       if(_iterator>=nbOfElems)
00239         throw INTERP_KERNEL::Exception("MEDCouplingUMesh::insertNextCell : allocation of cells was wide enough ! Call insertNextCell with higher value or call finishInsertingCells !");
00240       int *pt=_nodal_connec_index->getPointer();
00241       int idx=pt[_iterator];
00242       
00243       _nodal_connec->writeOnPlace(idx,type,nodalConnOfCell,size);
00244       _types.insert(type);
00245       pt[++_iterator]=idx+size+1;
00246     }
00247   else
00248     {
00249       std::ostringstream oss; oss << "MEDCouplingUMesh::insertNextCell : cell type " << cm.getRepr() << " has a dimension " << cm.getDimension();
00250       oss << " whereas Mesh Dimension of current UMesh instance is set to " << _mesh_dim << " ! Please invoke \"setMeshDimension\" method before or invoke ";
00251       oss << "\"MEDCouplingUMesh::New\" static method with 2 parameters name and meshDimension !";
00252       throw INTERP_KERNEL::Exception(oss.str().c_str());
00253     }
00254 }
00255 
00259 void MEDCouplingUMesh::finishInsertingCells()
00260 {
00261   const int *pt=_nodal_connec_index->getConstPointer();
00262   int idx=pt[_iterator];
00263 
00264   _nodal_connec->reAlloc(idx);
00265   _nodal_connec_index->reAlloc(_iterator+1);
00266   _iterator=-1;
00267   _nodal_connec->declareAsNew();
00268   _nodal_connec_index->declareAsNew();
00269   updateTime();
00270 }
00271 
00276 MEDCouplingUMeshCellIterator *MEDCouplingUMesh::cellIterator()
00277 {
00278   return new MEDCouplingUMeshCellIterator(this);
00279 }
00280 
00287 MEDCouplingUMeshCellByTypeEntry *MEDCouplingUMesh::cellsByType() throw(INTERP_KERNEL::Exception)
00288 {
00289   if(!checkConsecutiveCellTypes())
00290     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::cellsByType : this mesh is not sorted by type !");
00291   return new MEDCouplingUMeshCellByTypeEntry(this);
00292 }
00293 
00294 std::set<INTERP_KERNEL::NormalizedCellType> MEDCouplingUMesh::getAllGeoTypes() const
00295 {
00296   return _types;
00297 }
00298 
00303 bool MEDCouplingUMesh::isEqual(const MEDCouplingMesh *other, double prec) const
00304 {
00305   const MEDCouplingUMesh *otherC=dynamic_cast<const MEDCouplingUMesh *>(other);
00306   if(!otherC)
00307     return false;
00308   if(!MEDCouplingPointSet::isEqual(other,prec))
00309     return false;
00310   if(_mesh_dim!=otherC->_mesh_dim)
00311     return false;
00312   if(_types!=otherC->_types)
00313     return false;
00314   if(_nodal_connec!=0 || otherC->_nodal_connec!=0)
00315     if(_nodal_connec==0 || otherC->_nodal_connec==0)
00316       return false;
00317   if(_nodal_connec!=otherC->_nodal_connec)
00318     if(!_nodal_connec->isEqual(*otherC->_nodal_connec))
00319       return false;
00320   if(_nodal_connec_index!=0 || otherC->_nodal_connec_index!=0)
00321     if(_nodal_connec_index==0 || otherC->_nodal_connec_index==0)
00322       return false;
00323   if(_nodal_connec_index!=otherC->_nodal_connec_index)
00324     if(!_nodal_connec_index->isEqual(*otherC->_nodal_connec_index))
00325       return false;
00326   return true;
00327 }
00328 
00329 bool MEDCouplingUMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
00330 {
00331   const MEDCouplingUMesh *otherC=dynamic_cast<const MEDCouplingUMesh *>(other);
00332   if(!otherC)
00333     return false;
00334   if(!MEDCouplingPointSet::isEqualWithoutConsideringStr(other,prec))
00335     return false;
00336   if(_mesh_dim!=otherC->_mesh_dim)
00337     return false;
00338   if(_types!=otherC->_types)
00339     return false;
00340   if(_nodal_connec!=0 || otherC->_nodal_connec!=0)
00341     if(_nodal_connec==0 || otherC->_nodal_connec==0)
00342       return false;
00343   if(_nodal_connec!=otherC->_nodal_connec)
00344     if(!_nodal_connec->isEqualWithoutConsideringStr(*otherC->_nodal_connec))
00345       return false;
00346   if(_nodal_connec_index!=0 || otherC->_nodal_connec_index!=0)
00347     if(_nodal_connec_index==0 || otherC->_nodal_connec_index==0)
00348       return false;
00349   if(_nodal_connec_index!=otherC->_nodal_connec_index)
00350     if(!_nodal_connec_index->isEqualWithoutConsideringStr(*otherC->_nodal_connec_index))
00351       return false;
00352   return true;
00353 }
00354 
00364 void MEDCouplingUMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
00365                                             DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception)
00366 {
00367   const MEDCouplingUMesh *otherC=dynamic_cast<const MEDCouplingUMesh *>(other);
00368   if(!otherC)
00369     throw INTERP_KERNEL::Exception("checkDeepEquivalWith : Two meshes are not not unstructured !");
00370   MEDCouplingMesh::checkFastEquivalWith(other,prec);
00371   if(_types!=otherC->_types)
00372     throw INTERP_KERNEL::Exception("checkDeepEquivalWith : Types are not equal !");
00373   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=MergeUMeshes(this,otherC);
00374   bool areNodesMerged;
00375   int newNbOfNodes;
00376   int oldNbOfNodes=getNumberOfNodes();
00377   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> da=m->buildPermArrayForMergeNode(prec,oldNbOfNodes,areNodesMerged,newNbOfNodes);
00378   //mergeNodes
00379   if(!areNodesMerged)
00380     throw INTERP_KERNEL::Exception("checkDeepEquivalWith : Nodes are incompatible ! ");
00381   const int *pt=std::find_if(da->getConstPointer()+oldNbOfNodes,da->getConstPointer()+da->getNbOfElems(),std::bind2nd(std::greater<int>(),oldNbOfNodes-1));
00382   if(pt!=da->getConstPointer()+da->getNbOfElems())
00383     throw INTERP_KERNEL::Exception("checkDeepEquivalWith : some nodes in other are not in this !");
00384   m->renumberNodes(da->getConstPointer(),newNbOfNodes);
00385   //
00386   nodeCor=da->substr(oldNbOfNodes);
00387   da=m->mergeNodes(prec,areNodesMerged,newNbOfNodes);
00388   if(nodeCor->isIdentity())
00389     {
00390       nodeCor->decrRef();
00391       nodeCor=0;
00392     }
00393   //
00394   da=m->zipConnectivityTraducer(cellCompPol);
00395   int maxId=*std::max_element(da->getConstPointer(),da->getConstPointer()+getNumberOfCells());
00396   pt=std::find_if(da->getConstPointer()+getNumberOfCells(),da->getConstPointer()+da->getNbOfElems(),std::bind2nd(std::greater<int>(),maxId));
00397   if(pt!=da->getConstPointer()+da->getNbOfElems())
00398     {
00399       nodeCor->decrRef(); nodeCor=0;
00400       throw INTERP_KERNEL::Exception("checkDeepEquivalWith : some cells in other are not in this !");
00401     }
00402   cellCor=DataArrayInt::New();
00403   cellCor->alloc(otherC->getNumberOfCells(),1);
00404   std::copy(da->getConstPointer()+getNumberOfCells(),da->getConstPointer()+da->getNbOfElems(),cellCor->getPointer());
00405   if(cellCor->isIdentity())
00406     {
00407       cellCor->decrRef();
00408       cellCor=0;
00409     }
00410 }
00411 
00422 void MEDCouplingUMesh::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
00423                                                        DataArrayInt *&cellCor) const throw(INTERP_KERNEL::Exception)
00424 {
00425   const MEDCouplingUMesh *otherC=dynamic_cast<const MEDCouplingUMesh *>(other);
00426   if(!otherC)
00427     throw INTERP_KERNEL::Exception("checkDeepEquivalOnSameNodesWith : Two meshes are not not unstructured !");
00428   MEDCouplingMesh::checkFastEquivalWith(other,prec);
00429   if(_types!=otherC->_types)
00430     throw INTERP_KERNEL::Exception("checkDeepEquivalOnSameNodesWith : Types are not equal !");
00431   if(_coords!=otherC->_coords)
00432     throw INTERP_KERNEL::Exception("checkDeepEquivalOnSameNodesWith : meshes do not share the same coordinates ! Use tryToShareSameCoordinates or call checkDeepEquivalWith !");
00433   std::vector<const MEDCouplingUMesh *> ms(2);
00434   ms[0]=this;
00435   ms[1]=otherC;
00436   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=MergeUMeshesOnSameCoords(ms);
00437   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> da=m->zipConnectivityTraducer(cellCompPol);
00438   int maxId=*std::max_element(da->getConstPointer(),da->getConstPointer()+getNumberOfCells());
00439   const int *pt=std::find_if(da->getConstPointer()+getNumberOfCells(),da->getConstPointer()+da->getNbOfElems(),std::bind2nd(std::greater<int>(),maxId));
00440   if(pt!=da->getConstPointer()+da->getNbOfElems())
00441     {
00442       throw INTERP_KERNEL::Exception("checkDeepEquivalOnSameNodesWith : some cells in other are not in this !");
00443     }
00444   cellCor=DataArrayInt::New();
00445   cellCor->alloc(otherC->getNumberOfCells(),1);
00446   std::copy(da->getConstPointer()+getNumberOfCells(),da->getConstPointer()+da->getNbOfElems(),cellCor->getPointer());
00447   if(cellCor->isIdentity())
00448     {
00449       cellCor->decrRef();
00450       cellCor=0;
00451     }
00452 }
00453 
00457 void MEDCouplingUMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const throw(INTERP_KERNEL::Exception)
00458 {
00459   const MEDCouplingUMesh *otherC=dynamic_cast<const MEDCouplingUMesh *>(other);
00460   if(!otherC)
00461     throw INTERP_KERNEL::Exception("checkFastEquivalWith : Two meshes are not not unstructured !");
00462   MEDCouplingPointSet::checkFastEquivalWith(other,prec);
00463   int nbOfCells=getNumberOfCells();
00464   if(nbOfCells<1)
00465     return ;
00466   bool status=true;
00467   status&=areCellsFrom2MeshEqual(otherC,0,prec);
00468   status&=areCellsFrom2MeshEqual(otherC,nbOfCells/2,prec);
00469   status&=areCellsFrom2MeshEqual(otherC,nbOfCells-1,prec);
00470   if(!status)
00471     throw INTERP_KERNEL::Exception("checkFastEquivalWith : Two meshes are not equal because on 3 test cells some difference have been detected !");
00472 }
00473 
00478 void MEDCouplingUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const throw(INTERP_KERNEL::Exception)
00479 {
00480   checkFullyDefined();
00481   int nbOfNodes=getNumberOfNodes();
00482   int *revNodalIndxPtr=new int[nbOfNodes+1];
00483   revNodalIndx->useArray(revNodalIndxPtr,true,CPP_DEALLOC,nbOfNodes+1,1);
00484   std::fill(revNodalIndxPtr,revNodalIndxPtr+nbOfNodes+1,0);
00485   const int *conn=_nodal_connec->getConstPointer();
00486   const int *connIndex=_nodal_connec_index->getConstPointer();
00487   int nbOfCells=getNumberOfCells();
00488   int nbOfEltsInRevNodal=0;
00489   for(int eltId=0;eltId<nbOfCells;eltId++)
00490     {
00491       const int *strtNdlConnOfCurCell=conn+connIndex[eltId]+1;
00492       const int *endNdlConnOfCurCell=conn+connIndex[eltId+1];
00493       for(const int *iter=strtNdlConnOfCurCell;iter!=endNdlConnOfCurCell;iter++)
00494         if(*iter>=0)//for polyhedrons
00495           {
00496             nbOfEltsInRevNodal++;
00497             revNodalIndxPtr[(*iter)+1]++;
00498           }
00499     }
00500   std::transform(revNodalIndxPtr+1,revNodalIndxPtr+nbOfNodes+1,revNodalIndxPtr,revNodalIndxPtr+1,std::plus<int>());
00501   int *revNodalPtr=new int[nbOfEltsInRevNodal];
00502   revNodal->useArray(revNodalPtr,true,CPP_DEALLOC,nbOfEltsInRevNodal,1);
00503   std::fill(revNodalPtr,revNodalPtr+nbOfEltsInRevNodal,-1);
00504   for(int eltId=0;eltId<nbOfCells;eltId++)
00505     {
00506       const int *strtNdlConnOfCurCell=conn+connIndex[eltId]+1;
00507       const int *endNdlConnOfCurCell=conn+connIndex[eltId+1];
00508       for(const int *iter=strtNdlConnOfCurCell;iter!=endNdlConnOfCurCell;iter++)
00509         if(*iter>=0)//for polyhedrons
00510           *std::find_if(revNodalPtr+revNodalIndxPtr[*iter],revNodalPtr+revNodalIndxPtr[*iter+1],std::bind2nd(std::equal_to<int>(),-1))=eltId;
00511     }
00512 }
00513 
00515 
00516 int MEDCouplingFastNbrer(int id, unsigned nb, const INTERP_KERNEL::CellModel& cm, bool compute, const int *conn1, const int *conn2)
00517 {
00518   return id;
00519 }
00520 
00521 int MEDCouplingOrientationSensitiveNbrer(int id, unsigned nb, const INTERP_KERNEL::CellModel& cm, bool compute, const int *conn1, const int *conn2)
00522 {
00523   if(!compute)
00524     return id+1;
00525   else
00526     {
00527       if(cm.getOrientationStatus(nb,conn1,conn2))
00528         return id+1;
00529       else
00530         return -(id+1);
00531     }
00532 }
00533 
00535 
00552 MEDCouplingUMesh *MEDCouplingUMesh::buildDescendingConnectivity(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx) const throw(INTERP_KERNEL::Exception)
00553 {
00554   return buildDescendingConnectivityGen(desc,descIndx,revDesc,revDescIndx,MEDCouplingFastNbrer);
00555 }
00556 
00565 MEDCouplingUMesh *MEDCouplingUMesh::buildDescendingConnectivity2(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx) const throw(INTERP_KERNEL::Exception)
00566 {
00567   return buildDescendingConnectivityGen(desc,descIndx,revDesc,revDescIndx,MEDCouplingOrientationSensitiveNbrer);
00568 }
00569 
00580 void MEDCouplingUMesh::computeNeighborsOfCells(DataArrayInt *&neighbors, DataArrayInt *&neighborsIdx) const throw(INTERP_KERNEL::Exception)
00581 {
00582   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> desc=DataArrayInt::New();
00583   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> descIndx=DataArrayInt::New();
00584   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> revDesc=DataArrayInt::New();
00585   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> revDescIndx=DataArrayInt::New();
00586   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> meshDM1=buildDescendingConnectivity(desc,descIndx,revDesc,revDescIndx);
00587   meshDM1=0;
00588   const int *descPtr=desc->getConstPointer();
00589   const int *descIPtr=descIndx->getConstPointer();
00590   const int *revDescPtr=revDesc->getConstPointer();
00591   const int *revDescIPtr=revDescIndx->getConstPointer();
00592   //
00593   int nbCells=getNumberOfCells();
00594   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> out0=DataArrayInt::New();
00595   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> out1=DataArrayInt::New(); out1->alloc(nbCells+1,1);
00596   int *out1Ptr=out1->getPointer();
00597   *out1Ptr++=0;
00598   std::vector<int> out0v;
00599   out0v.reserve(desc->getNumberOfTuples());
00600   for(int i=0;i<nbCells;i++,descIPtr++,out1Ptr++)
00601     {
00602       for(const int *w1=descPtr+descIPtr[0];w1!=descPtr+descIPtr[1];w1++)
00603         {
00604           std::set<int> s(revDescPtr+revDescIPtr[*w1],revDescPtr+revDescIPtr[(*w1)+1]);
00605           s.erase(i);
00606           out0v.insert(out0v.end(),s.begin(),s.end());
00607         }
00608       *out1Ptr=out0v.size();
00609     }
00610   out0->alloc((int)out0v.size(),1);
00611   std::copy(out0v.begin(),out0v.end(),out0->getPointer());
00612   neighbors=out0; out0->incrRef();
00613   neighborsIdx=out1; out1->incrRef();
00614 }
00615 
00617 
00622 MEDCouplingUMesh *MEDCouplingUMesh::buildDescendingConnectivityGen(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx, DimM1DescNbrer nbrer) const throw(INTERP_KERNEL::Exception)
00623 {
00624   checkFullyDefined();
00625   int nbOfCells=getNumberOfCells();
00626   int nbOfNodes=getNumberOfNodes();
00627   const int *conn=_nodal_connec->getConstPointer();
00628   const int *connIndex=_nodal_connec_index->getConstPointer();
00629   std::vector< std::vector<int> > descMeshConnB(nbOfCells);
00630   std::vector< std::vector<int> > revDescMeshConnB;
00631   std::vector< std::vector<int> > revNodalB(nbOfNodes);
00632   std::vector<int> meshDM1Conn;
00633   std::vector<int> meshDM1ConnIndex(1); meshDM1ConnIndex[0]=0;
00634   std::vector<int> meshDM1Type;
00635   for(int eltId=0;eltId<nbOfCells;eltId++)
00636     {
00637       int pos=connIndex[eltId];
00638       int posP1=connIndex[eltId+1];
00639       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[pos]);
00640       unsigned nbOfSons=cm.getNumberOfSons2(conn+pos+1,posP1-pos-1);
00641       int *tmp=new int[posP1-pos];
00642       for(unsigned i=0;i<nbOfSons;i++)
00643         {
00644           INTERP_KERNEL::NormalizedCellType cmsId;
00645           unsigned nbOfNodesSon=cm.fillSonCellNodalConnectivity2(i,conn+pos+1,posP1-pos-1,tmp,cmsId);
00646           const INTERP_KERNEL::CellModel& cms=INTERP_KERNEL::CellModel::GetCellModel(cmsId);
00647           std::set<int> shareableCells(revNodalB[tmp[0]].begin(),revNodalB[tmp[0]].end());
00648           for(unsigned j=1;j<nbOfNodesSon && !shareableCells.empty();j++)
00649             {
00650               std::set<int> tmp2(revNodalB[tmp[j]].begin(),revNodalB[tmp[j]].end());
00651               std::set<int> tmp3;
00652               std::set_intersection(tmp2.begin(),tmp2.end(),shareableCells.begin(),shareableCells.end(),inserter(tmp3,tmp3.begin()));
00653               shareableCells=tmp3;
00654             }
00655           std::list<int> shareableCellsL(shareableCells.begin(),shareableCells.end());
00656           std::set<int> ref(tmp,tmp+nbOfNodesSon);
00657           for(std::list<int>::iterator iter=shareableCellsL.begin();iter!=shareableCellsL.end();)
00658             {
00659               if(cms.isCompatibleWith((INTERP_KERNEL::NormalizedCellType)meshDM1Type[*iter]))
00660                 {
00661                   std::set<int> ref2(meshDM1Conn.begin()+meshDM1ConnIndex[*iter],meshDM1Conn.begin()+meshDM1ConnIndex[(*iter)+1]);
00662                   if(ref==ref2)
00663                     break;
00664                   else
00665                     iter=shareableCellsL.erase(iter);
00666                 }
00667               else
00668                 iter=shareableCellsL.erase(iter);
00669             }
00670           if(shareableCellsL.empty())
00671             {
00672               meshDM1Conn.insert(meshDM1Conn.end(),tmp,tmp+nbOfNodesSon);
00673               meshDM1ConnIndex.push_back(meshDM1ConnIndex.back()+nbOfNodesSon);
00674               int cellDM1Id=(int)meshDM1Type.size();
00675               meshDM1Type.push_back((int)cmsId);
00676               for(unsigned k=0;k<nbOfNodesSon;k++)
00677                 revNodalB[tmp[k]].push_back(cellDM1Id);
00678               revDescMeshConnB.resize(cellDM1Id+1);
00679               revDescMeshConnB.back().push_back(eltId);
00680               descMeshConnB[eltId].push_back(nbrer(cellDM1Id,0,cms,false,0,0));
00681             }
00682           else
00683             {
00684               int DM1cellId=shareableCellsL.front();
00685               revDescMeshConnB[DM1cellId].push_back(eltId);
00686               descMeshConnB[eltId].push_back(nbrer(DM1cellId,nbOfNodesSon,cms,true,tmp,&meshDM1Conn[meshDM1ConnIndex[DM1cellId]]));
00687             }
00688         }
00689       delete [] tmp;
00690     }
00691   revNodalB.clear();
00692   //
00693   std::string name="Mesh constituent of "; name+=getName();
00694   MEDCouplingUMesh *ret=MEDCouplingUMesh::New(name.c_str(),getMeshDimension()-1);
00695   ret->setCoords(getCoords());
00696   int nbOfCellsInConstituent=(int)meshDM1Type.size();
00697   ret->allocateCells(nbOfCellsInConstituent);
00698   revDescIndx->alloc(nbOfCellsInConstituent+1,1);
00699   int *tmp3=revDescIndx->getPointer(); tmp3[0]=0;
00700   for(int ii=0;ii<nbOfCellsInConstituent;ii++)
00701     {
00702       ret->insertNextCell((INTERP_KERNEL::NormalizedCellType)meshDM1Type[ii],meshDM1ConnIndex[ii+1]-meshDM1ConnIndex[ii],&meshDM1Conn[meshDM1ConnIndex[ii]]);
00703       tmp3[ii+1]=tmp3[ii]+((int)revDescMeshConnB[ii].size());
00704     }
00705   ret->finishInsertingCells();
00706   revDesc->alloc(tmp3[nbOfCellsInConstituent],1);
00707   tmp3=revDesc->getPointer();
00708   for(std::vector< std::vector<int> >::const_iterator iter2=revDescMeshConnB.begin();iter2!=revDescMeshConnB.end();iter2++)
00709     tmp3=std::copy((*iter2).begin(),(*iter2).end(),tmp3);
00710   meshDM1Type.clear(); meshDM1ConnIndex.clear(); meshDM1Conn.clear();
00711   descIndx->alloc(nbOfCells+1,1);
00712   tmp3=descIndx->getPointer(); tmp3[0]=0;
00713   for(int jj=0;jj<nbOfCells;jj++)
00714     tmp3[jj+1]=tmp3[jj]+((int)descMeshConnB[jj].size());
00715   desc->alloc(tmp3[nbOfCells],1);
00716   tmp3=desc->getPointer();
00717   for(std::vector< std::vector<int> >::const_iterator iter3=descMeshConnB.begin();iter3!=descMeshConnB.end();iter3++)
00718     tmp3=std::copy((*iter3).begin(),(*iter3).end(),tmp3);
00719   //
00720   return ret;
00721 }
00722 
00723 struct MEDCouplingAccVisit
00724 {
00725   MEDCouplingAccVisit():_new_nb_of_nodes(0) { }
00726   int operator()(int val) { if(val!=-1) return _new_nb_of_nodes++; else return -1; }
00727   int _new_nb_of_nodes;
00728 };
00729 
00731 
00732 
00744 void MEDCouplingUMesh::convertToPolyTypes(const int *cellIdsToConvertBg, const int *cellIdsToConvertEnd)
00745 {
00746   checkFullyDefined();
00747   int dim=getMeshDimension();
00748   if(dim<2 || dim>3)
00749     throw INTERP_KERNEL::Exception("Invalid mesh dimension : must be 2 or 3 !");
00750   int nbOfCells=getNumberOfCells();
00751   if(dim==2)
00752     {
00753       const int *connIndex=_nodal_connec_index->getConstPointer();
00754       int *conn=_nodal_connec->getPointer();
00755       for(const int *iter=cellIdsToConvertBg;iter!=cellIdsToConvertEnd;iter++)
00756         {
00757           if(*iter>=0 && *iter<nbOfCells)
00758             {
00759               const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[connIndex[*iter]]);
00760               if(!cm.isDynamic())
00761                 conn[connIndex[*iter]]=INTERP_KERNEL::NORM_POLYGON;
00762               else
00763                 conn[connIndex[*iter]]=INTERP_KERNEL::NORM_QPOLYG;
00764             }
00765           else
00766             {
00767               std::ostringstream oss; oss << "MEDCouplingUMesh::convertToPolyTypes : On rank #" << std::distance(cellIdsToConvertBg,iter) << " value is " << *iter << " which is not";
00768               oss << " in range [0," << nbOfCells << ") !";
00769               throw INTERP_KERNEL::Exception(oss.str().c_str());
00770             }
00771         }
00772     }
00773   else
00774     {
00775       int *connIndex=_nodal_connec_index->getPointer();
00776       int connIndexLgth=_nodal_connec_index->getNbOfElems();
00777       const int *connOld=_nodal_connec->getConstPointer();
00778       int connOldLgth=_nodal_connec->getNbOfElems();
00779       std::vector<int> connNew(connOld,connOld+connOldLgth);
00780       for(const int *iter=cellIdsToConvertBg;iter!=cellIdsToConvertEnd;iter++)
00781         {
00782           if(*iter>=0 && *iter<nbOfCells)
00783             {
00784               int pos=connIndex[*iter];
00785               int posP1=connIndex[(*iter)+1];
00786               int lgthOld=posP1-pos-1;
00787               const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connNew[pos]);
00788               connNew[pos]=INTERP_KERNEL::NORM_POLYHED;
00789               unsigned nbOfFaces=cm.getNumberOfSons2(&connNew[pos+1],lgthOld);
00790               int *tmp=new int[nbOfFaces*lgthOld];
00791               int *work=tmp;
00792               for(int j=0;j<(int)nbOfFaces;j++)
00793                 {
00794                   INTERP_KERNEL::NormalizedCellType type;
00795                   unsigned offset=cm.fillSonCellNodalConnectivity2(j,&connNew[pos+1],lgthOld,work,type);
00796                   work+=offset;
00797                   *work++=-1;
00798                 }
00799               std::size_t newLgth=std::distance(tmp,work)-1;
00800               std::size_t delta=newLgth-lgthOld;
00801               std::transform(connIndex+(*iter)+1,connIndex+connIndexLgth,connIndex+(*iter)+1,std::bind2nd(std::plus<int>(),delta));
00802               connNew.insert(connNew.begin()+posP1,tmp+lgthOld,tmp+newLgth);
00803               std::copy(tmp,tmp+lgthOld,connNew.begin()+pos+1);
00804               delete [] tmp;
00805             }
00806           else
00807             {
00808               std::ostringstream oss; oss << "MEDCouplingUMesh::convertToPolyTypes : On rank #" << std::distance(cellIdsToConvertBg,iter) << " value is " << *iter << " which is not";
00809               oss << " in range [0," << nbOfCells << ") !";
00810               throw INTERP_KERNEL::Exception(oss.str().c_str());
00811             }
00812         }
00813       _nodal_connec->alloc((int)connNew.size(),1);
00814       int *newConnPtr=_nodal_connec->getPointer();
00815       std::copy(connNew.begin(),connNew.end(),newConnPtr);
00816     }
00817   computeTypes();
00818 }
00819 
00825 void MEDCouplingUMesh::convertAllToPoly()
00826 {
00827   int nbOfCells=getNumberOfCells();
00828   std::vector<int> cellIds(nbOfCells);
00829   for(int i=0;i<nbOfCells;i++)
00830     cellIds[i]=i;
00831   convertToPolyTypes(&cellIds[0],&cellIds[0]+cellIds.size());
00832 }
00833 
00847 void MEDCouplingUMesh::convertExtrudedPolyhedra() throw(INTERP_KERNEL::Exception)
00848 {
00849   checkFullyDefined();
00850   if(getMeshDimension()!=3 || getSpaceDimension()!=3)
00851     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::convertExtrudedPolyhedra works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
00852   int nbOfCells=getNumberOfCells();
00853   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newCi=DataArrayInt::New();
00854   newCi->alloc(nbOfCells+1,1);
00855   int *newci=newCi->getPointer();
00856   const int *ci=_nodal_connec_index->getConstPointer();
00857   const int *c=_nodal_connec->getConstPointer();
00858   newci[0]=0;
00859   for(int i=0;i<nbOfCells;i++)
00860     {
00861       INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)c[ci[i]];
00862       if(type==INTERP_KERNEL::NORM_POLYHED)
00863         {
00864           if(std::count(c+ci[i]+1,c+ci[i+1],-1)!=0)
00865             {
00866               std::ostringstream oss; oss << "MEDCouplingUMesh::convertExtrudedPolyhedra : cell # " << i << " is a polhedron BUT it has NOT exactly 1 face !";
00867               throw INTERP_KERNEL::Exception(oss.str().c_str());
00868             }
00869           std::size_t n2=std::distance(c+ci[i]+1,c+ci[i+1]);
00870           if(n2%2!=0)
00871             {
00872               std::ostringstream oss; oss << "MEDCouplingUMesh::convertExtrudedPolyhedra : cell # " << i << " is a polhedron with 1 face but there is a mismatch of number of nodes in face should be even !";
00873               throw INTERP_KERNEL::Exception(oss.str().c_str());
00874             }
00875           int n1=(int)(n2/2);
00876           newci[i+1]=7*n1+2+newci[i];//6*n1 (nodal length) + n1+2 (number of faces) - 1 (number of '-1' separator is equal to number of faces -1) + 1 (for cell type)
00877         }
00878       else
00879         newci[i+1]=(ci[i+1]-ci[i])+newci[i];
00880     }
00881   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newC=DataArrayInt::New();
00882   newC->alloc(newci[nbOfCells],1);
00883   int *newc=newC->getPointer();
00884   for(int i=0;i<nbOfCells;i++)
00885     {
00886       INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)c[ci[i]];
00887       if(type==INTERP_KERNEL::NORM_POLYHED)
00888         {
00889           std::size_t n1=std::distance(c+ci[i]+1,c+ci[i+1])/2;
00890           newc=std::copy(c+ci[i],c+ci[i]+n1+1,newc);
00891           *newc++=-1;
00892           for(std::size_t j=0;j<n1;j++)
00893             {
00894               newc[j]=c[ci[i]+1+n1+(n1-j)%n1];
00895               newc[n1+5*j]=-1;
00896               newc[n1+5*j+1]=c[ci[i]+1+j];
00897               newc[n1+5*j+2]=c[ci[i]+1+(j+1)%n1];
00898               newc[n1+5*j+3]=c[ci[i]+1+(j+1)%n1+n1];
00899               newc[n1+5*j+4]=c[ci[i]+1+j+n1];
00900             }
00901           newc+=n1*6;
00902         }
00903       else
00904         newc=std::copy(c+ci[i],c+ci[i+1],newc);
00905     }
00906   _nodal_connec_index->decrRef(); _nodal_connec_index=newCi;
00907   _nodal_connec->decrRef(); _nodal_connec=newC;
00908   newC->incrRef(); newCi->incrRef();
00909 }
00910 
00916 void MEDCouplingUMesh::unPolyze()
00917 {
00918   checkFullyDefined();
00919   if(getMeshDimension()<=1)
00920     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::unPolyze works on umeshes with meshdim equals to 2 or 3 !");
00921   int nbOfCells=getNumberOfCells();
00922   if(nbOfCells<1)
00923     return ;
00924   int initMeshLgth=getMeshLength();
00925   int *conn=_nodal_connec->getPointer();
00926   int *index=_nodal_connec_index->getPointer();
00927   int posOfCurCell=0;
00928   int newPos=0;
00929   int lgthOfCurCell;
00930   for(int i=0;i<nbOfCells;i++)
00931     {
00932       lgthOfCurCell=index[i+1]-posOfCurCell;
00933       INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)conn[posOfCurCell];
00934       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
00935       INTERP_KERNEL::NormalizedCellType newType=INTERP_KERNEL::NORM_ERROR;
00936       int newLgth;
00937       if(cm.isDynamic())
00938         {
00939           if(cm.getDimension()==2)
00940             {
00941               INTERP_KERNEL::AutoPtr<int> tmp=new int[lgthOfCurCell-1];
00942               std::copy(conn+posOfCurCell+1,conn+posOfCurCell+lgthOfCurCell,(int *)tmp);
00943               newType=INTERP_KERNEL::CellSimplify::tryToUnPoly2D(cm.isQuadratic(),tmp,lgthOfCurCell-1,conn+newPos+1,newLgth);
00944             }
00945           if(cm.getDimension()==3)
00946             {
00947               int nbOfFaces,lgthOfPolyhConn;
00948               INTERP_KERNEL::AutoPtr<int> zipFullReprOfPolyh=INTERP_KERNEL::CellSimplify::getFullPolyh3DCell(type,conn+posOfCurCell+1,lgthOfCurCell-1,nbOfFaces,lgthOfPolyhConn);
00949               newType=INTERP_KERNEL::CellSimplify::tryToUnPoly3D(zipFullReprOfPolyh,nbOfFaces,lgthOfPolyhConn,conn+newPos+1,newLgth);
00950             }
00951           conn[newPos]=newType;
00952           newPos+=newLgth+1;
00953           posOfCurCell=index[i+1];
00954           index[i+1]=newPos;
00955         }
00956       else
00957         {
00958           std::copy(conn+posOfCurCell,conn+posOfCurCell+lgthOfCurCell,conn+newPos);
00959           newPos+=lgthOfCurCell;
00960           posOfCurCell+=lgthOfCurCell;
00961           index[i+1]=newPos;
00962         }
00963     }
00964   if(newPos!=initMeshLgth)
00965     _nodal_connec->reAlloc(newPos);
00966   computeTypes();
00967 }
00968 
00977 DataArrayInt *MEDCouplingUMesh::getNodeIdsInUse(int& nbrOfNodesInUse) const throw(INTERP_KERNEL::Exception)
00978 {
00979   nbrOfNodesInUse=-1;
00980   int nbOfNodes=getNumberOfNodes();
00981   DataArrayInt *ret=DataArrayInt::New();
00982   ret->alloc(nbOfNodes,1);
00983   int *traducer=ret->getPointer();
00984   std::fill(traducer,traducer+nbOfNodes,-1);
00985   int nbOfCells=getNumberOfCells();
00986   const int *connIndex=_nodal_connec_index->getConstPointer();
00987   int *conn=_nodal_connec->getPointer();
00988   for(int i=0;i<nbOfCells;i++)
00989     for(int j=connIndex[i]+1;j<connIndex[i+1];j++)
00990       if(conn[j]>=0)
00991         traducer[conn[j]]=1;
00992   nbrOfNodesInUse=(int)std::count(traducer,traducer+nbOfNodes,1);
00993   std::transform(traducer,traducer+nbOfNodes,traducer,MEDCouplingAccVisit());
00994   return ret;
00995 }
00996 
01003 DataArrayInt *MEDCouplingUMesh::zipCoordsTraducer() throw(INTERP_KERNEL::Exception)
01004 {
01005   int newNbOfNodes=-1;
01006   DataArrayInt *traducer=getNodeIdsInUse(newNbOfNodes);
01007   renumberNodes(traducer->getConstPointer(),newNbOfNodes);
01008   return traducer;
01009 }
01010 
01015 bool MEDCouplingUMesh::areCellsEqual(int cell1, int cell2, int compType) const
01016 {
01017   switch(compType)
01018     {
01019     case 0:
01020       return areCellsEqual0(cell1,cell2);
01021     case 1:
01022       return areCellsEqual1(cell1,cell2);
01023     case 2:
01024       return areCellsEqual2(cell1,cell2);
01025     }
01026   throw INTERP_KERNEL::Exception("Unknown comparison asked ! Must be in 0,1 or 2.");
01027 }
01028 
01032 bool MEDCouplingUMesh::areCellsEqual0(int cell1, int cell2) const
01033 {
01034   const int *conn=getNodalConnectivity()->getConstPointer();
01035   const int *connI=getNodalConnectivityIndex()->getConstPointer();
01036   if(connI[cell1+1]-connI[cell1]==connI[cell2+1]-connI[cell2])
01037     return std::equal(conn+connI[cell1]+1,conn+connI[cell1+1],conn+connI[cell2]+1);
01038   return false;
01039 }
01040 
01044 bool MEDCouplingUMesh::areCellsEqual1(int cell1, int cell2) const
01045 {
01046   const int *conn=getNodalConnectivity()->getConstPointer();
01047   const int *connI=getNodalConnectivityIndex()->getConstPointer();
01048   int sz=connI[cell1+1]-connI[cell1];
01049   if(sz==connI[cell2+1]-connI[cell2])
01050     {
01051       if(conn[connI[cell1]]==conn[connI[cell2]])
01052         {
01053           const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[connI[cell1]]);
01054           unsigned dim=cm.getDimension();
01055           if(dim!=3)
01056             {
01057               if(dim!=1)
01058                 {
01059                   int sz1=2*(sz-1);
01060                   int *tmp=new int[sz1];
01061                   int *work=std::copy(conn+connI[cell1]+1,conn+connI[cell1+1],tmp);
01062                   std::copy(conn+connI[cell1]+1,conn+connI[cell1+1],work);
01063                   work=std::search(tmp,tmp+sz1,conn+connI[cell2]+1,conn+connI[cell2+1]);
01064                   delete [] tmp;
01065                   return work!=tmp+sz1;
01066                 }
01067               else
01068                 return std::equal(conn+connI[cell1]+1,conn+connI[cell1+1],conn+connI[cell2]+1);//case of SEG2 and SEG3
01069             }
01070           else
01071             throw INTERP_KERNEL::Exception("MEDCouplingUMesh::areCellsEqual1 : not implemented yet for meshdim == 3 !");
01072         }
01073     }
01074   return false;
01075 }
01076 
01080 bool MEDCouplingUMesh::areCellsEqual2(int cell1, int cell2) const
01081 {
01082   const int *conn=getNodalConnectivity()->getConstPointer();
01083   const int *connI=getNodalConnectivityIndex()->getConstPointer();
01084   if(connI[cell1+1]-connI[cell1]==connI[cell2+1]-connI[cell2])
01085     {
01086       if(conn[connI[cell1]]==conn[connI[cell2]])
01087         {
01088           std::set<int> s1(conn+connI[cell1]+1,conn+connI[cell1+1]);
01089           std::set<int> s2(conn+connI[cell2]+1,conn+connI[cell2+1]);
01090           return s1==s2;
01091         }
01092     }
01093   return false;
01094 }
01095 
01100 bool MEDCouplingUMesh::areCellsFrom2MeshEqual(const MEDCouplingUMesh *other, int cellId, double prec) const
01101 {
01102   if(getTypeOfCell(cellId)!=other->getTypeOfCell(cellId))
01103     return false;
01104   std::vector<int> c1,c2;
01105   getNodeIdsOfCell(cellId,c1);
01106   other->getNodeIdsOfCell(cellId,c2);
01107   std::size_t sz=c1.size();
01108   if(sz!=c2.size())
01109     return false;
01110   for(std::size_t i=0;i<sz;i++)
01111     {
01112       std::vector<double> n1,n2;
01113       getCoordinatesOfNode(c1[0],n1);
01114       other->getCoordinatesOfNode(c2[0],n2);
01115       std::transform(n1.begin(),n1.end(),n2.begin(),n1.begin(),std::minus<double>());
01116       std::transform(n1.begin(),n1.end(),n1.begin(),std::ptr_fun<double,double>(fabs));
01117       if(*std::max_element(n1.begin(),n1.end())>prec)
01118         return false;
01119     }
01120   return true;
01121 }
01122 
01131 bool MEDCouplingUMesh::areCellsEqualInPool(const std::vector<int>& candidates, int compType, std::vector<int>& result) const
01132 {
01133   std::set<int> cand(candidates.begin(),candidates.end());
01134   cand.erase(-1);
01135   if(cand.size()<=1)
01136     return false;
01137   std::set<int>::const_iterator end=cand.end(); end--;
01138   bool ret=false;
01139   for(std::set<int>::const_iterator iter=cand.begin();iter!=end && !ret;iter++)
01140     {
01141       std::set<int>::const_iterator begin2=iter; begin2++;
01142       for(std::set<int>::const_iterator iter2=begin2;iter2!=cand.end();iter2++)
01143         {
01144           if(areCellsEqual(*iter,*iter2,compType))
01145             {
01146               if(!ret)
01147                 {
01148                   result.push_back(*iter);
01149                   ret=true;
01150                 }
01151               result.push_back(*iter2);
01152             }
01153         }
01154     }
01155   return ret;
01156 }
01157 
01166 template<int SPACEDIM>
01167 void MEDCouplingUMesh::findCommonCellsBase(int compType, std::vector<int>& res, std::vector<int>& resI) const
01168 {
01169   res.clear(); resI.clear();
01170   resI.push_back(0);
01171   std::vector<double> bbox;
01172   int nbOfCells=getNumberOfCells();
01173   getBoundingBoxForBBTree(bbox);
01174   double bb[2*SPACEDIM];
01175   double eps=getCaracteristicDimension();
01176   eps*=1.e-12;
01177   BBTree<SPACEDIM,int> myTree(&bbox[0],0,0,nbOfCells,-eps);
01178   const int *conn=getNodalConnectivity()->getConstPointer();
01179   const int *connI=getNodalConnectivityIndex()->getConstPointer();
01180   const double *coords=getCoords()->getConstPointer();
01181   std::vector<bool> isFetched(nbOfCells);
01182   for(int k=0;k<nbOfCells;k++)
01183     {
01184       if(!isFetched[k])
01185         {
01186           for(int j=0;j<SPACEDIM;j++)
01187             { bb[2*j]=std::numeric_limits<double>::max(); bb[2*j+1]=-std::numeric_limits<double>::max(); }
01188           for(const int *pt=conn+connI[k]+1;pt!=conn+connI[k+1];pt++)
01189             if(*pt>-1)
01190               {
01191                 for(int j=0;j<SPACEDIM;j++)
01192                   {
01193                     bb[2*j]=std::min(bb[2*j],coords[SPACEDIM*(*pt)+j]);
01194                     bb[2*j+1]=std::max(bb[2*j+1],coords[SPACEDIM*(*pt)+j]);
01195                   }
01196               }
01197           std::vector<int> candidates1;
01198           myTree.getIntersectingElems(bb,candidates1);
01199           std::vector<int> candidates;
01200           for(std::vector<int>::const_iterator iter=candidates1.begin();iter!=candidates1.end();iter++)
01201             if(!isFetched[*iter])
01202               candidates.push_back(*iter);
01203           if(areCellsEqualInPool(candidates,compType,res))
01204             {
01205               int pos=resI.back();
01206               resI.push_back((int)res.size());
01207               for(std::vector<int>::const_iterator it=res.begin()+pos;it!=res.end();it++)
01208                 isFetched[*it]=true;
01209             }
01210           isFetched[k]=true;
01211         }
01212     }
01213 }
01214 
01227 DataArrayInt *MEDCouplingUMesh::zipConnectivityTraducer(int compType) throw(INTERP_KERNEL::Exception)
01228 {
01229   int spaceDim=getSpaceDimension();
01230   int nbOfCells=getNumberOfCells();
01231   std::vector<int> commonCells;
01232   std::vector<int> commonCellsI;
01233   switch(spaceDim)
01234     {
01235     case 3:
01236       {
01237         findCommonCellsBase<3>(compType,commonCells,commonCellsI);
01238         break;
01239       }
01240     case 2:
01241       {
01242         findCommonCellsBase<2>(compType,commonCells,commonCellsI);
01243         break;
01244       }
01245     case 1:
01246       {
01247         findCommonCellsBase<1>(compType,commonCells,commonCellsI);
01248         break;
01249       }
01250     default:
01251       throw INTERP_KERNEL::Exception("Invalid spaceDimension : must be 1, 2 or 3.");
01252     }
01253   DataArrayInt *ret=DataArrayInt::New();
01254   ret->alloc(nbOfCells,1);
01255   int *retPtr=ret->getPointer();
01256   std::fill(retPtr,retPtr+nbOfCells,0);
01257   const std::size_t nbOfTupleSmCells=commonCellsI.size()-1;
01258   int id=-1;
01259   std::vector<int> cellsToKeep;
01260   for(std::size_t i=0;i<nbOfTupleSmCells;i++)
01261     {
01262       for(std::vector<int>::const_iterator it=commonCells.begin()+commonCellsI[i];it!=commonCells.begin()+commonCellsI[i+1];it++)
01263         retPtr[*it]=id;
01264       id--;
01265     }
01266   id=0;
01267   std::map<int,int> m;
01268   for(int i=0;i<nbOfCells;i++)
01269     {
01270       int val=retPtr[i];
01271       if(val==0)
01272         {
01273           retPtr[i]=id++;
01274           cellsToKeep.push_back(i);
01275         }
01276       else
01277         {
01278           std::map<int,int>::const_iterator iter=m.find(val);
01279           if(iter==m.end())
01280             {
01281               m[val]=id;
01282               retPtr[i]=id++;
01283               cellsToKeep.push_back(i);
01284             }
01285           else
01286             retPtr[i]=(*iter).second;
01287         }
01288     }
01289   MEDCouplingUMesh *self=(MEDCouplingUMesh *)buildPartOfMySelf(&cellsToKeep[0],&cellsToKeep[0]+cellsToKeep.size(),true);
01290   setConnectivity(self->getNodalConnectivity(),self->getNodalConnectivityIndex(),true);
01291   self->decrRef();
01292   return ret;
01293 }
01294 
01306 bool MEDCouplingUMesh::areCellsIncludedIn(const MEDCouplingUMesh *other, int compType, DataArrayInt *& arr) const throw(INTERP_KERNEL::Exception)
01307 {
01308   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> mesh=MergeUMeshesOnSameCoords(this,other);
01309   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> o2n=mesh->zipConnectivityTraducer(compType);
01310   int nbOfCells=getNumberOfCells();
01311   arr=o2n->substr(nbOfCells);
01312   arr->setName(other->getName());
01313   int tmp;
01314   if(other->getNumberOfCells()==0)
01315     return true;
01316   return arr->getMaxValue(tmp)<nbOfCells;
01317 }
01318 
01323 DataArrayInt *MEDCouplingUMesh::mergeNodes(double precision, bool& areNodesMerged, int& newNbOfNodes)
01324 {
01325   DataArrayInt *ret=buildPermArrayForMergeNode(precision,-1,areNodesMerged,newNbOfNodes);
01326   if(areNodesMerged)
01327     renumberNodes(ret->getConstPointer(),newNbOfNodes);
01328   return ret;
01329 }
01330 
01334 DataArrayInt *MEDCouplingUMesh::mergeNodes2(double precision, bool& areNodesMerged, int& newNbOfNodes)
01335 {
01336   DataArrayInt *ret=buildPermArrayForMergeNode(precision,-1,areNodesMerged,newNbOfNodes);
01337   if(areNodesMerged)
01338     renumberNodes2(ret->getConstPointer(),newNbOfNodes);
01339   return ret;
01340 }
01341 
01348 void MEDCouplingUMesh::tryToShareSameCoordsPermute(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception)
01349 {
01350   const DataArrayDouble *coords=other.getCoords();
01351   if(!coords)
01352     throw INTERP_KERNEL::Exception("tryToShareSameCoordsPermute : No coords specified in other !");
01353   if(!_coords)
01354     throw INTERP_KERNEL::Exception("tryToShareSameCoordsPermute : No coords specified in this whereas there is any in other !");
01355   int otherNbOfNodes=other.getNumberOfNodes();
01356   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=MergeNodesArray(&other,this);
01357   _coords->incrRef();
01358   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> oldCoords=_coords;
01359   setCoords(newCoords);
01360   bool areNodesMerged;
01361   int newNbOfNodes;
01362   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> da=buildPermArrayForMergeNode(epsilon,otherNbOfNodes,areNodesMerged,newNbOfNodes);
01363   if(!areNodesMerged)
01364     {
01365       setCoords(oldCoords);
01366       throw INTERP_KERNEL::Exception("tryToShareSameCoordsPermute fails : no nodes are mergeable with specified given epsilon !");
01367     }
01368   int maxId=*std::max_element(da->getConstPointer(),da->getConstPointer()+otherNbOfNodes);
01369   const int *pt=std::find_if(da->getConstPointer()+otherNbOfNodes,da->getConstPointer()+da->getNbOfElems(),std::bind2nd(std::greater<int>(),maxId));
01370   if(pt!=da->getConstPointer()+da->getNbOfElems())
01371     {
01372       setCoords(oldCoords);
01373       throw INTERP_KERNEL::Exception("tryToShareSameCoordsPermute fails : some nodes in this are not in other !");
01374     }
01375   setCoords(oldCoords);
01376   renumberNodesInConn(da->getConstPointer()+otherNbOfNodes);
01377   setCoords(coords);
01378 }
01379 
01386 MEDCouplingPointSet *MEDCouplingUMesh::buildPartOfMySelf(const int *begin, const int *end, bool keepCoords) const
01387 {
01388   if(getMeshDimension()!=-1)
01389     {
01390       MEDCouplingUMesh *ret=buildPartOfMySelfKeepCoords(begin,end);
01391       if(!keepCoords)
01392         ret->zipCoords();
01393       return ret;
01394     }
01395   else
01396     {
01397       if(end-begin!=1)
01398         throw INTERP_KERNEL::Exception("-1D mesh has only one cell !");
01399       if(begin[0]!=0)
01400         throw INTERP_KERNEL::Exception("-1D mesh has only one cell : 0 !");
01401       incrRef();
01402       return const_cast<MEDCouplingUMesh *>(this);
01403     }
01404 }
01405 
01406 DataArrayInt *MEDCouplingUMesh::getCellIdsFullyIncludedInNodeIds(const int *partBg, const int *partEnd) const
01407 {
01408   std::vector<int> cellIdsKept;
01409   fillCellIdsToKeepFromNodeIds(partBg,partEnd,true,cellIdsKept);
01410   DataArrayInt *ret=DataArrayInt::New();
01411   ret->alloc((int)cellIdsKept.size(),1);
01412   std::copy(cellIdsKept.begin(),cellIdsKept.end(),ret->getPointer());
01413   return ret;
01414 }
01415 
01427 void MEDCouplingUMesh::fillCellIdsToKeepFromNodeIds(const int *begin, const int *end, bool fullyIn, std::vector<int>& cellIdsKept) const
01428 {
01429   std::set<int> fastFinder(begin,end);
01430   int nbOfCells=getNumberOfCells();
01431   const int *conn=getNodalConnectivity()->getConstPointer();
01432   const int *connIndex=getNodalConnectivityIndex()->getConstPointer();
01433   for(int i=0;i<nbOfCells;i++)
01434     {
01435       std::set<int> connOfCell(conn+connIndex[i]+1,conn+connIndex[i+1]);
01436       connOfCell.erase(-1);//polyhedron separator
01437       int refLgth=(int)connOfCell.size();
01438       std::set<int> locMerge;
01439       std::insert_iterator< std::set<int> > it(locMerge,locMerge.begin());
01440       std::set_intersection(connOfCell.begin(),connOfCell.end(),fastFinder.begin(),fastFinder.end(),it);
01441       if(((int)locMerge.size()==refLgth && fullyIn) || (locMerge.size()!=0 && !fullyIn))
01442         cellIdsKept.push_back(i);
01443     }
01444 }
01445 
01449 DataArrayInt *MEDCouplingUMesh::getCellIdsLyingOnNodes(const int *begin, const int *end, bool fullyIn) const
01450 {
01451   std::vector<int> cellIdsKept;
01452   fillCellIdsToKeepFromNodeIds(begin,end,fullyIn,cellIdsKept);
01453   DataArrayInt *ret=DataArrayInt::New();
01454   ret->alloc((int)cellIdsKept.size(),1);
01455   std::copy(cellIdsKept.begin(),cellIdsKept.end(),ret->getPointer());
01456   ret->setName(getName());
01457   return ret;
01458 }
01459 
01466 MEDCouplingPointSet *MEDCouplingUMesh::buildPartOfMySelfNode(const int *begin, const int *end, bool fullyIn) const
01467 {
01468   std::vector<int> cellIdsKept;
01469   fillCellIdsToKeepFromNodeIds(begin,end,fullyIn,cellIdsKept);
01470   return buildPartOfMySelf(&cellIdsKept[0],&cellIdsKept[0]+cellIdsKept.size(),true);
01471 }
01472 
01479 MEDCouplingPointSet *MEDCouplingUMesh::buildFacePartOfMySelfNode(const int *begin, const int *end, bool fullyIn) const
01480 {
01481   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> desc,descIndx,revDesc,revDescIndx;
01482   desc=DataArrayInt::New(); descIndx=DataArrayInt::New(); revDesc=DataArrayInt::New(); revDescIndx=DataArrayInt::New();
01483   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> subMesh=buildDescendingConnectivity(desc,descIndx,revDesc,revDescIndx);
01484   desc=0; descIndx=0; revDesc=0; revDescIndx=0;
01485   return subMesh->buildPartOfMySelfNode(begin,end,fullyIn);
01486 }
01487 
01494 MEDCouplingPointSet *MEDCouplingUMesh::buildBoundaryMesh(bool keepCoords) const
01495 {
01496   DataArrayInt *desc=DataArrayInt::New();
01497   DataArrayInt *descIndx=DataArrayInt::New();
01498   DataArrayInt *revDesc=DataArrayInt::New();
01499   DataArrayInt *revDescIndx=DataArrayInt::New();
01500   //
01501   MEDCouplingUMesh *meshDM1=buildDescendingConnectivity(desc,descIndx,revDesc,revDescIndx);
01502   revDesc->decrRef();
01503   desc->decrRef();
01504   descIndx->decrRef();
01505   int nbOfCells=meshDM1->getNumberOfCells();
01506   const int *revDescIndxC=revDescIndx->getConstPointer();
01507   std::vector<int> boundaryCells;
01508   for(int i=0;i<nbOfCells;i++)
01509     if(revDescIndxC[i+1]-revDescIndxC[i]==1)
01510       boundaryCells.push_back(i);
01511   revDescIndx->decrRef();
01512   MEDCouplingPointSet *ret=meshDM1->buildPartOfMySelf(&boundaryCells[0],&boundaryCells[0]+boundaryCells.size(),keepCoords);
01513   meshDM1->decrRef();
01514   return ret;
01515 }
01516 
01522 DataArrayInt *MEDCouplingUMesh::findCellsIdsOnBoundary() const throw(INTERP_KERNEL::Exception)
01523 {
01524   checkFullyDefined();
01525   DataArrayInt *desc=DataArrayInt::New();
01526   DataArrayInt *descIndx=DataArrayInt::New();
01527   DataArrayInt *revDesc=DataArrayInt::New();
01528   DataArrayInt *revDescIndx=DataArrayInt::New();
01529   //
01530   MEDCouplingUMesh *meshDM1=buildDescendingConnectivity(desc,descIndx,revDesc,revDescIndx);
01531   meshDM1->decrRef();
01532   desc->decrRef();
01533   descIndx->decrRef();
01534   //
01535   DataArrayInt *tmp=revDescIndx->deltaShiftIndex();
01536   DataArrayInt *faceIds=tmp->getIdsEqual(1);
01537   tmp->decrRef();
01538   int nbOfFaces=faceIds->getNumberOfTuples();
01539   const int *faces=faceIds->getConstPointer();
01540   std::set<int> ret;
01541   for(const int *w=faces;w!=faces+nbOfFaces;w++)
01542     ret.insert(revDesc->getIJ(revDescIndx->getIJ(*w,0),0));
01543   faceIds->decrRef();
01544   //
01545   revDescIndx->decrRef();
01546   revDesc->decrRef();
01547   //
01548   DataArrayInt *ret2=DataArrayInt::New();
01549   ret2->alloc((int)ret.size(),1);
01550   std::copy(ret.begin(),ret.end(),ret2->getPointer());
01551   ret2->setName("BoundaryCells");
01552   return ret2;
01553 }
01554 
01558 void MEDCouplingUMesh::findBoundaryNodes(std::vector<int>& nodes) const
01559 {
01560   DataArrayInt *desc=DataArrayInt::New();
01561   DataArrayInt *descIndx=DataArrayInt::New();
01562   DataArrayInt *revDesc=DataArrayInt::New();
01563   DataArrayInt *revDescIndx=DataArrayInt::New();
01564   //
01565   MEDCouplingUMesh *meshDM1=buildDescendingConnectivity(desc,descIndx,revDesc,revDescIndx);
01566   revDesc->decrRef();
01567   desc->decrRef();
01568   descIndx->decrRef();
01569   std::set<int> ret;
01570   int nbOfCells=meshDM1->getNumberOfCells();
01571   const int *revDescIndxC=revDescIndx->getConstPointer();
01572   std::vector<int> boundaryCells;
01573   for(int i=0;i<nbOfCells;i++)
01574     if(revDescIndxC[i+1]-revDescIndxC[i]==1)
01575       boundaryCells.push_back(i);
01576   revDescIndx->decrRef();
01577   const int *conn=meshDM1->getNodalConnectivity()->getConstPointer();
01578   const int *connIndx=meshDM1->getNodalConnectivityIndex()->getConstPointer();
01579   for(std::vector<int>::const_iterator iter=boundaryCells.begin();iter!=boundaryCells.end();iter++)
01580     for(int k=connIndx[*iter]+1;k<connIndx[(*iter)+1];k++)
01581       ret.insert(conn[k]);
01582   nodes.resize(ret.size());
01583   std::copy(ret.begin(),ret.end(),nodes.begin());
01584   //
01585   meshDM1->decrRef();
01586 }
01587 
01588 MEDCouplingUMesh *MEDCouplingUMesh::buildUnstructured() const throw(INTERP_KERNEL::Exception)
01589 {
01590   incrRef();
01591   return const_cast<MEDCouplingUMesh *>(this);
01592 }
01593 
01594 /*
01595  * This method renumber 'this' using 'newNodeNumbers' array of size this->getNumberOfNodes.
01596  * newNbOfNodes specifies the *std::max_element(newNodeNumbers,newNodeNumbers+this->getNumberOfNodes())
01597  * This value is asked because often known by the caller of this method.
01598  * This method, contrary to MEDCouplingMesh::renumberCells does NOT conserve the number of nodes before and after.
01599  *
01600  * @param newNodeNumbers array specifying the new numbering in old2New convention.
01601  * @param newNbOfNodes the new number of nodes.
01602  */
01603 void MEDCouplingUMesh::renumberNodes(const int *newNodeNumbers, int newNbOfNodes)
01604 {
01605   MEDCouplingPointSet::renumberNodes(newNodeNumbers,newNbOfNodes);
01606   renumberNodesInConn(newNodeNumbers);
01607 }
01608 
01609 /*
01610  * This method renumber 'this' using 'newNodeNumbers' array of size this->getNumberOfNodes.
01611  * newNbOfNodes specifies the *std::max_element(newNodeNumbers,newNodeNumbers+this->getNumberOfNodes())
01612  * This value is asked because often known by the caller of this method.
01613  * This method, contrary to MEDCouplingMesh::renumberCells does NOT conserve the number of nodes before and after.
01614  * The difference with ParaMEDMEM::MEDCouplingUMesh::renumberNodes method is in the fact that the barycenter of merged nodes is computed here.
01615  *
01616  * @param newNodeNumbers array specifying the new numbering.
01617  * @param newNbOfNodes the new number of nodes.
01618  */
01619 void MEDCouplingUMesh::renumberNodes2(const int *newNodeNumbers, int newNbOfNodes)
01620 {
01621   MEDCouplingPointSet::renumberNodes2(newNodeNumbers,newNbOfNodes);
01622   renumberNodesInConn(newNodeNumbers);
01623 }
01624 
01631 void MEDCouplingUMesh::renumberNodesInConn(const int *newNodeNumbersO2N)
01632 {
01633   checkConnectivityFullyDefined();
01634   int *conn=getNodalConnectivity()->getPointer();
01635   const int *connIndex=getNodalConnectivityIndex()->getConstPointer();
01636   int nbOfCells=getNumberOfCells();
01637   for(int i=0;i<nbOfCells;i++)
01638     for(int iconn=connIndex[i]+1;iconn!=connIndex[i+1];iconn++)
01639       {
01640         int& node=conn[iconn];
01641         if(node>=0)//avoid polyhedron separator
01642           {
01643             node=newNodeNumbersO2N[node];
01644           }
01645       }
01646   _nodal_connec->declareAsNew();
01647   updateTime();
01648 }
01649 
01657 void MEDCouplingUMesh::shiftNodeNumbersInConn(int delta) throw(INTERP_KERNEL::Exception)
01658 {
01659   checkConnectivityFullyDefined();
01660   int *conn=getNodalConnectivity()->getPointer();
01661   const int *connIndex=getNodalConnectivityIndex()->getConstPointer();
01662   int nbOfCells=getNumberOfCells();
01663   for(int i=0;i<nbOfCells;i++)
01664     for(int iconn=connIndex[i]+1;iconn!=connIndex[i+1];iconn++)
01665       {
01666         int& node=conn[iconn];
01667         if(node>=0)//avoid polyhedron separator
01668           {
01669             node+=delta;
01670           }
01671       }
01672   _nodal_connec->declareAsNew();
01673   updateTime();
01674 }
01675 
01690 void MEDCouplingUMesh::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
01691 {
01692   checkConnectivityFullyDefined();
01693   int nbCells=getNumberOfCells();
01694   const int *array=old2NewBg;
01695   if(check)
01696     array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+nbCells);
01697   //
01698   const int *conn=_nodal_connec->getConstPointer();
01699   const int *connI=_nodal_connec_index->getConstPointer();
01700   DataArrayInt *newConn=DataArrayInt::New();
01701   newConn->alloc(_nodal_connec->getNumberOfTuples(),_nodal_connec->getNumberOfComponents());
01702   newConn->copyStringInfoFrom(*_nodal_connec);
01703   DataArrayInt *newConnI=DataArrayInt::New();
01704   newConnI->alloc(_nodal_connec_index->getNumberOfTuples(),_nodal_connec_index->getNumberOfComponents());
01705   newConnI->copyStringInfoFrom(*_nodal_connec_index);
01706   //
01707   int *newC=newConn->getPointer();
01708   int *newCI=newConnI->getPointer();
01709   int loc=0;
01710   newCI[0]=loc;
01711   for(int i=0;i<nbCells;i++)
01712     {
01713       std::size_t pos=std::distance(array,std::find(array,array+nbCells,i));
01714       int nbOfElts=connI[pos+1]-connI[pos];
01715       newC=std::copy(conn+connI[pos],conn+connI[pos+1],newC);
01716       loc+=nbOfElts;
01717       newCI[i+1]=loc;
01718     }
01719   //
01720   setConnectivity(newConn,newConnI);
01721   //
01722   newConn->decrRef();
01723   newConnI->decrRef();
01724   if(check)
01725     delete [] const_cast<int *>(array);
01726 }
01727 
01733 void MEDCouplingUMesh::getCellsInBoundingBox(const double *bbox, double eps, std::vector<int>& elems) const
01734 {
01735   if(getMeshDimension()==-1)
01736     {
01737       elems.push_back(0);
01738       return;
01739     }
01740   int dim=getSpaceDimension();
01741   double* elem_bb=new double[2*dim];
01742   const int* conn      = getNodalConnectivity()->getConstPointer();
01743   const int* conn_index= getNodalConnectivityIndex()->getConstPointer();
01744   const double* coords = getCoords()->getConstPointer();
01745   int nbOfCells=getNumberOfCells();
01746   for ( int ielem=0; ielem<nbOfCells;ielem++ )
01747     {
01748       for (int i=0; i<dim; i++)
01749         {
01750           elem_bb[i*2]=std::numeric_limits<double>::max();
01751           elem_bb[i*2+1]=-std::numeric_limits<double>::max();
01752         }
01753 
01754       for (int inode=conn_index[ielem]+1; inode<conn_index[ielem+1]; inode++)//+1 due to offset of cell type.
01755         {
01756           int node= conn[inode];
01757           if(node>=0)//avoid polyhedron separator
01758             {
01759               for (int idim=0; idim<dim; idim++)
01760                 {
01761                   if ( coords[node*dim+idim] < elem_bb[idim*2] )
01762                     {
01763                       elem_bb[idim*2] = coords[node*dim+idim] ;
01764                     }
01765                   if ( coords[node*dim+idim] > elem_bb[idim*2+1] )
01766                     {
01767                       elem_bb[idim*2+1] = coords[node*dim+idim] ;
01768                     }
01769                 }
01770             }
01771         }
01772       if (intersectsBoundingBox(elem_bb, bbox, dim, eps))
01773         {
01774           elems.push_back(ielem);
01775         }
01776     }
01777   delete [] elem_bb;
01778 }
01779 
01785 void MEDCouplingUMesh::getCellsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bbox, double eps, std::vector<int>& elems)
01786 {
01787   if(getMeshDimension()==-1)
01788     {
01789       elems.push_back(0);
01790       return;
01791     }
01792   int dim=getSpaceDimension();
01793   double* elem_bb=new double[2*dim];
01794   const int* conn      = getNodalConnectivity()->getConstPointer();
01795   const int* conn_index= getNodalConnectivityIndex()->getConstPointer();
01796   const double* coords = getCoords()->getConstPointer();
01797   int nbOfCells=getNumberOfCells();
01798   for ( int ielem=0; ielem<nbOfCells;ielem++ )
01799     {
01800       for (int i=0; i<dim; i++)
01801         {
01802           elem_bb[i*2]=std::numeric_limits<double>::max();
01803           elem_bb[i*2+1]=-std::numeric_limits<double>::max();
01804         }
01805 
01806       for (int inode=conn_index[ielem]+1; inode<conn_index[ielem+1]; inode++)//+1 due to offset of cell type.
01807         {
01808           int node= conn[inode];
01809           if(node>=0)//avoid polyhedron separator
01810             {
01811               for (int idim=0; idim<dim; idim++)
01812                 {
01813                   if ( coords[node*dim+idim] < elem_bb[idim*2] )
01814                     {
01815                       elem_bb[idim*2] = coords[node*dim+idim] ;
01816                     }
01817                   if ( coords[node*dim+idim] > elem_bb[idim*2+1] )
01818                     {
01819                       elem_bb[idim*2+1] = coords[node*dim+idim] ;
01820                     }
01821                 }
01822             }
01823         }
01824       if (intersectsBoundingBox(bbox, elem_bb, dim, eps))
01825         {
01826           elems.push_back(ielem);
01827         }
01828     }
01829   delete [] elem_bb;
01830 }
01831 
01835 INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::getTypeOfCell(int cellId) const
01836 {
01837   const int *ptI=_nodal_connec_index->getConstPointer();
01838   const int *pt=_nodal_connec->getConstPointer();
01839   return (INTERP_KERNEL::NormalizedCellType) pt[ptI[cellId]];
01840 }
01841 
01845 int MEDCouplingUMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
01846 {
01847   const int *ptI=_nodal_connec_index->getConstPointer();
01848   const int *pt=_nodal_connec->getConstPointer();
01849   int nbOfCells=getNumberOfCells();
01850   int ret=0;
01851   for(int i=0;i<nbOfCells;i++)
01852     if((INTERP_KERNEL::NormalizedCellType) pt[ptI[i]]==type)
01853       ret++;
01854   return ret;
01855 }
01856 
01862 void MEDCouplingUMesh::getNodeIdsOfCell(int cellId, std::vector<int>& conn) const
01863 {
01864   const int *ptI=_nodal_connec_index->getConstPointer();
01865   const int *pt=_nodal_connec->getConstPointer();
01866   for(const int *w=pt+ptI[cellId]+1;w!=pt+ptI[cellId+1];w++)
01867     if(*w>=0)
01868       conn.push_back(*w);
01869 }
01870 
01871 std::string MEDCouplingUMesh::simpleRepr() const
01872 {
01873   static const char msg0[]="No coordinates specified !";
01874   std::ostringstream ret;
01875   ret << "Unstructured mesh with name : \"" << getName() << "\"\n";
01876   ret << "Description of mesh : \"" << getDescription() << "\"\n";
01877   int tmpp1,tmpp2;
01878   double tt=getTime(tmpp1,tmpp2);
01879   ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
01880   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
01881   ret << "Mesh dimension : " << _mesh_dim << "\nSpace dimension : ";
01882   if(_coords!=0)
01883     {
01884       const int spaceDim=getSpaceDimension();
01885       ret << spaceDim << "\nInfo attached on space dimension : ";
01886       for(int i=0;i<spaceDim;i++)
01887         ret << "\"" << _coords->getInfoOnComponent(i) << "\" ";
01888       ret << "\n";
01889     }
01890   else
01891     ret << msg0 << "\n";
01892   ret << "Number of nodes : ";
01893   if(_coords!=0)
01894     ret << getNumberOfNodes() << "\n";
01895   else
01896     ret << msg0 << "\n";
01897   ret << "Number of cells : ";
01898   if(_nodal_connec!=0 && _nodal_connec_index!=0)
01899     ret << getNumberOfCells() << "\n";
01900   else
01901     ret << "No connectivity specified !" << "\n";
01902   ret << "Cell types present : ";
01903   for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator iter=_types.begin();iter!=_types.end();iter++)
01904     {
01905       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(*iter);
01906       ret << cm.getRepr() << " ";
01907     }
01908   ret << "\n";
01909   return ret.str();
01910 }
01911 
01912 std::string MEDCouplingUMesh::advancedRepr() const
01913 {
01914   std::ostringstream ret;
01915   ret << simpleRepr();
01916   ret << "\nCoordinates array : \n___________________\n\n";
01917   if(_coords)
01918     _coords->reprWithoutNameStream(ret);
01919   else
01920     ret << "No array set !\n";
01921   ret << "\n\nConnectivity arrays : \n_____________________\n\n";
01922   reprConnectivityOfThisLL(ret);
01923   return ret.str();
01924 }
01925 
01926 std::string MEDCouplingUMesh::reprConnectivityOfThis() const
01927 {
01928   std::ostringstream ret;
01929   reprConnectivityOfThisLL(ret);
01930   return ret.str();
01931 }
01932 
01943 MEDCouplingUMesh *MEDCouplingUMesh::buildSetInstanceFromThis(int spaceDim) const throw(INTERP_KERNEL::Exception)
01944 {
01945   int mdim=getMeshDimension();
01946   if(mdim<0)
01947     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSetInstanceFromThis : invalid mesh dimension ! Should be >= 0 !");
01948   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret=MEDCouplingUMesh::New(getName(),mdim);
01949   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp1,tmp2;
01950   bool needToCpyCT=true;
01951   if(!_nodal_connec)
01952     {
01953       tmp1=DataArrayInt::New(); tmp1->alloc(0,1);
01954       needToCpyCT=false;
01955     }
01956   else
01957     {
01958       tmp1=_nodal_connec;
01959       tmp1->incrRef();
01960     }
01961   if(!_nodal_connec_index)
01962     {
01963       tmp2=DataArrayInt::New(); tmp2->alloc(1,1); tmp2->setIJ(0,0,0);
01964       needToCpyCT=false;
01965     }
01966   else
01967     {
01968       tmp2=_nodal_connec_index;
01969       tmp2->incrRef();
01970     }
01971   ret->setConnectivity(tmp1,tmp2,false);
01972   if(needToCpyCT)
01973     ret->_types=_types;
01974   if(!_coords)
01975     {
01976       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=DataArrayDouble::New(); coords->alloc(0,spaceDim);
01977       ret->setCoords(coords);
01978     }
01979   else
01980     ret->setCoords(_coords);
01981   ret->incrRef();
01982   return ret;
01983 }
01984 
01985 void MEDCouplingUMesh::reprConnectivityOfThisLL(std::ostringstream& stream) const
01986 {
01987   if(_nodal_connec!=0 && _nodal_connec_index!=0)
01988     {
01989       int nbOfCells=getNumberOfCells();
01990       const int *c=_nodal_connec->getConstPointer();
01991       const int *ci=_nodal_connec_index->getConstPointer();
01992       for(int i=0;i<nbOfCells;i++)
01993         {
01994           const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)c[ci[i]]);
01995           stream << "Cell #" << i << " " << cm.getRepr() << " : ";
01996           std::copy(c+ci[i]+1,c+ci[i+1],std::ostream_iterator<int>(stream," "));
01997           stream << "\n";
01998         }
01999     }
02000   else
02001     stream << "Connectivity not defined !\n";
02002 }
02003 
02004 int MEDCouplingUMesh::getNumberOfNodesInCell(int cellId) const
02005 {
02006   const int *ptI=_nodal_connec_index->getConstPointer();
02007   const int *pt=_nodal_connec->getConstPointer();
02008   if(pt[ptI[cellId]]!=INTERP_KERNEL::NORM_POLYHED)
02009     return ptI[cellId+1]-ptI[cellId]-1;
02010   else
02011     return (int)std::count_if(pt+ptI[cellId]+1,pt+ptI[cellId+1],std::bind2nd(std::not_equal_to<int>(),-1));
02012 }
02013 
02018 std::set<INTERP_KERNEL::NormalizedCellType> MEDCouplingUMesh::getTypesOfPart(const int *begin, const int *end) const throw(INTERP_KERNEL::Exception)
02019 {
02020   checkFullyDefined();
02021   std::set<INTERP_KERNEL::NormalizedCellType> ret;
02022   const int *conn=_nodal_connec->getConstPointer();
02023   const int *connIndex=_nodal_connec_index->getConstPointer();
02024   for(const int *w=begin;w!=end;w++)
02025     ret.insert((INTERP_KERNEL::NormalizedCellType)conn[connIndex[*w]]);
02026   return ret;
02027 }
02028 
02033 void MEDCouplingUMesh::setConnectivity(DataArrayInt *conn, DataArrayInt *connIndex, bool isComputingTypes)
02034 {
02035   DataArrayInt::SetArrayIn(conn,_nodal_connec);
02036   DataArrayInt::SetArrayIn(connIndex,_nodal_connec_index);
02037   if(isComputingTypes)
02038     computeTypes();
02039   declareAsNew();
02040 }
02041 
02046 MEDCouplingUMesh::MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCopy):MEDCouplingPointSet(other,deepCopy),_iterator(-1),_mesh_dim(other._mesh_dim),
02047                                                                                  _nodal_connec(0),_nodal_connec_index(0),
02048                                                                                 _types(other._types)
02049 {
02050   if(other._nodal_connec)
02051     _nodal_connec=other._nodal_connec->performCpy(deepCopy);
02052   if(other._nodal_connec_index)
02053     _nodal_connec_index=other._nodal_connec_index->performCpy(deepCopy);
02054 }
02055 
02056 MEDCouplingUMesh::~MEDCouplingUMesh()
02057 {
02058   if(_nodal_connec)
02059     _nodal_connec->decrRef();
02060   if(_nodal_connec_index)
02061     _nodal_connec_index->decrRef();
02062 }
02063 
02067 void MEDCouplingUMesh::computeTypes()
02068 {
02069   if(_nodal_connec && _nodal_connec_index)
02070     {
02071       _types.clear();
02072       const int *conn=_nodal_connec->getConstPointer();
02073       const int *connIndex=_nodal_connec_index->getConstPointer();
02074       int nbOfElem=_nodal_connec_index->getNbOfElems()-1;
02075       for(const int *pt=connIndex;pt!=connIndex+nbOfElem;pt++)
02076         _types.insert((INTERP_KERNEL::NormalizedCellType)conn[*pt]);
02077     }
02078 }
02079 
02083 void MEDCouplingUMesh::checkFullyDefined() const throw(INTERP_KERNEL::Exception)
02084 {
02085   if(!_nodal_connec_index || !_nodal_connec || !_coords)
02086     throw INTERP_KERNEL::Exception("Reverse nodal connectivity computation requires full connectivity and coordinates set in unstructured mesh.");
02087 }
02088 
02092 void MEDCouplingUMesh::checkConnectivityFullyDefined() const throw(INTERP_KERNEL::Exception)
02093 {
02094   if(!_nodal_connec_index || !_nodal_connec)
02095     throw INTERP_KERNEL::Exception("Reverse nodal connectivity computation requires full connectivity set in unstructured mesh.");
02096 }
02097 
02098 int MEDCouplingUMesh::getNumberOfCells() const
02099 { 
02100   if(_nodal_connec_index)
02101     if(_iterator==-1)
02102       return _nodal_connec_index->getNumberOfTuples()-1;
02103     else
02104       return _iterator;
02105   else
02106     if(_mesh_dim==-1)
02107       return 1;
02108     else
02109       throw INTERP_KERNEL::Exception("Unable to get number of cells because no connectivity specified !");
02110 }
02111 
02112 int MEDCouplingUMesh::getMeshDimension() const
02113 {
02114   if(_mesh_dim<-1)
02115     throw INTERP_KERNEL::Exception("No mesh dimension specified !");
02116   return _mesh_dim;
02117 }
02118 
02122 int MEDCouplingUMesh::getMeshLength() const
02123 {
02124   return _nodal_connec->getNbOfElems();
02125 }
02126 
02130 void MEDCouplingUMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
02131 {
02132   MEDCouplingPointSet::getTinySerializationInformation(tinyInfoD,tinyInfo,littleStrings);
02133   tinyInfo.push_back(getMeshDimension());
02134   tinyInfo.push_back(getNumberOfCells());
02135   if(_nodal_connec)
02136     tinyInfo.push_back(getMeshLength());
02137   else
02138     tinyInfo.push_back(-1);
02139 }
02140 
02144 bool MEDCouplingUMesh::isEmptyMesh(const std::vector<int>& tinyInfo) const
02145 {
02146   return tinyInfo[6]<=0;
02147 }
02148 
02153 void MEDCouplingUMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
02154 {
02155   MEDCouplingPointSet::resizeForUnserialization(tinyInfo,a1,a2,littleStrings);
02156   if(tinyInfo[5]!=-1)
02157     a1->alloc(tinyInfo[7]+tinyInfo[6]+1,1);
02158 }
02159 
02163 void MEDCouplingUMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
02164 {
02165   MEDCouplingPointSet::serialize(a1,a2);
02166   if(getMeshDimension()>-1)
02167     {
02168       a1=DataArrayInt::New();
02169       a1->alloc(getMeshLength()+getNumberOfCells()+1,1);
02170       int *ptA1=a1->getPointer();
02171       const int *conn=getNodalConnectivity()->getConstPointer();
02172       const int *index=getNodalConnectivityIndex()->getConstPointer();
02173       ptA1=std::copy(index,index+getNumberOfCells()+1,ptA1);
02174       std::copy(conn,conn+getMeshLength(),ptA1);
02175     }
02176   else
02177     a1=0;
02178 }
02179 
02184 void MEDCouplingUMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector<std::string>& littleStrings)
02185 {
02186   MEDCouplingPointSet::unserialization(tinyInfoD,tinyInfo,a1,a2,littleStrings);
02187   setMeshDimension(tinyInfo[5]);
02188   if(tinyInfo[7]!=-1)
02189     {
02190       // Connectivity
02191       const int *recvBuffer=a1->getConstPointer();
02192       DataArrayInt* myConnecIndex=DataArrayInt::New();
02193       myConnecIndex->alloc(tinyInfo[6]+1,1);
02194       std::copy(recvBuffer,recvBuffer+tinyInfo[6]+1,myConnecIndex->getPointer());
02195       DataArrayInt* myConnec=DataArrayInt::New();
02196       myConnec->alloc(tinyInfo[7],1);
02197       std::copy(recvBuffer+tinyInfo[6]+1,recvBuffer+tinyInfo[6]+1+tinyInfo[7],myConnec->getPointer());
02198       setConnectivity(myConnec, myConnecIndex) ;
02199       myConnec->decrRef();
02200       myConnecIndex->decrRef();
02201     }
02202 }
02203 
02209 MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoords(const int *begin, const int *end) const
02210 {
02211   checkFullyDefined();
02212   int ncell=getNumberOfCells();
02213   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret=MEDCouplingUMesh::New();
02214   ret->_mesh_dim=_mesh_dim;
02215   ret->setCoords(_coords);
02216   std::size_t nbOfElemsRet=std::distance(begin,end);
02217   int *connIndexRet=new int[nbOfElemsRet+1];
02218   connIndexRet[0]=0;
02219   const int *conn=_nodal_connec->getConstPointer();
02220   const int *connIndex=_nodal_connec_index->getConstPointer();
02221   int newNbring=0;
02222   for(const int *work=begin;work!=end;work++,newNbring++)
02223     {
02224       if(*work>=0 && *work<ncell)
02225         connIndexRet[newNbring+1]=connIndexRet[newNbring]+connIndex[*work+1]-connIndex[*work];
02226       else
02227         {
02228           delete [] connIndexRet;
02229           std::ostringstream oss; oss << "MEDCouplingUMesh::buildPartOfMySelfKeepCoords : On pos #" << std::distance(begin,work) << " input cell id =" << *work << " should be in [0," << ncell << ") !";
02230           throw INTERP_KERNEL::Exception(oss.str().c_str());
02231         }
02232     }
02233   int *connRet=new int[connIndexRet[nbOfElemsRet]];
02234   int *connRetWork=connRet;
02235   std::set<INTERP_KERNEL::NormalizedCellType> types;
02236   for(const int *work=begin;work!=end;work++)
02237     {
02238       types.insert((INTERP_KERNEL::NormalizedCellType)conn[connIndex[*work]]);
02239       connRetWork=std::copy(conn+connIndex[*work],conn+connIndex[*work+1],connRetWork);
02240     }
02241   DataArrayInt *connRetArr=DataArrayInt::New();
02242   connRetArr->useArray(connRet,true,CPP_DEALLOC,connIndexRet[nbOfElemsRet],1);
02243   DataArrayInt *connIndexRetArr=DataArrayInt::New();
02244   connIndexRetArr->useArray(connIndexRet,true,CPP_DEALLOC,(int)nbOfElemsRet+1,1);
02245   ret->setConnectivity(connRetArr,connIndexRetArr,false);
02246   ret->_types=types;
02247   connRetArr->decrRef();
02248   connIndexRetArr->decrRef();
02249   ret->copyTinyInfoFrom(this);
02250   std::string name(getName());
02251   std::size_t sz=strlen(PART_OF_NAME);
02252   if(name.length()>=sz)
02253     name=name.substr(0,sz);
02254   if(name!=PART_OF_NAME)
02255     {
02256       std::ostringstream stream; stream << PART_OF_NAME << getName();
02257       ret->setName(stream.str().c_str());
02258     }
02259   else
02260     ret->setName(getName());
02261   ret->incrRef();
02262   return ret;
02263 }
02264 
02274 MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureField(bool isAbs) const
02275 {
02276   std::string name="MeasureOfMesh_";
02277   name+=getName();
02278   int nbelem=getNumberOfCells();
02279   MEDCouplingFieldDouble *field=MEDCouplingFieldDouble::New(ON_CELLS);
02280   field->setName(name.c_str());
02281   DataArrayDouble* array=DataArrayDouble::New();
02282   array->alloc(nbelem,1);
02283   double *area_vol=array->getPointer();
02284   field->setArray(array) ;
02285   array->decrRef();
02286   field->setMesh(const_cast<MEDCouplingUMesh *>(this));
02287   if(getMeshDimension()!=-1)
02288     {
02289       int ipt;
02290       INTERP_KERNEL::NormalizedCellType type;
02291       int dim_space=getSpaceDimension();
02292       const double *coords=getCoords()->getConstPointer();
02293       const int *connec=getNodalConnectivity()->getConstPointer();
02294       const int *connec_index=getNodalConnectivityIndex()->getConstPointer();
02295       for(int iel=0;iel<nbelem;iel++)
02296         {
02297           ipt=connec_index[iel];
02298           type=(INTERP_KERNEL::NormalizedCellType)connec[ipt];
02299           area_vol[iel]=INTERP_KERNEL::computeVolSurfOfCell2<int,INTERP_KERNEL::ALL_C_MODE>(type,connec+ipt+1,connec_index[iel+1]-ipt-1,coords,dim_space);
02300         }
02301       if(isAbs)
02302         std::transform(area_vol,area_vol+nbelem,area_vol,std::ptr_fun<double,double>(fabs));
02303     }
02304   else
02305     {
02306       area_vol[0]=std::numeric_limits<double>::max();
02307     }
02308   return field;
02309 }
02310 
02315 DataArrayDouble *MEDCouplingUMesh::getPartMeasureField(bool isAbs, const int *begin, const int *end) const
02316 {
02317   std::string name="PartMeasureOfMesh_";
02318   name+=getName();
02319   int nbelem=(int)std::distance(begin,end);
02320   DataArrayDouble* array=DataArrayDouble::New();
02321   array->setName(name.c_str());
02322   array->alloc(nbelem,1);
02323   double *area_vol=array->getPointer();
02324   if(getMeshDimension()!=-1)
02325     {
02326       int ipt;
02327       INTERP_KERNEL::NormalizedCellType type;
02328       int dim_space=getSpaceDimension();
02329       const double *coords=getCoords()->getConstPointer();
02330       const int *connec=getNodalConnectivity()->getConstPointer();
02331       const int *connec_index=getNodalConnectivityIndex()->getConstPointer();
02332       for(const int *iel=begin;iel!=end;iel++)
02333         {
02334           ipt=connec_index[*iel];
02335           type=(INTERP_KERNEL::NormalizedCellType)connec[ipt];
02336           *area_vol++=INTERP_KERNEL::computeVolSurfOfCell2<int,INTERP_KERNEL::ALL_C_MODE>(type,connec+ipt+1,connec_index[*iel+1]-ipt-1,coords,dim_space);
02337         }
02338       if(isAbs)
02339         std::transform(array->getPointer(),area_vol,array->getPointer(),std::ptr_fun<double,double>(fabs));
02340     }
02341   else
02342     {
02343       area_vol[0]=std::numeric_limits<double>::max();
02344     }
02345   return array;
02346 }
02347 
02352 MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureFieldOnNode(bool isAbs) const
02353 {
02354   MEDCouplingFieldDouble *tmp=getMeasureField(isAbs);
02355   std::string name="MeasureOnNodeOfMesh_";
02356   name+=getName();
02357   int nbNodes=getNumberOfNodes();
02358   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_NODES);
02359   double cst=1./((double)getMeshDimension()+1.);
02360   DataArrayDouble* array=DataArrayDouble::New();
02361   array->alloc(nbNodes,1);
02362   double *valsToFill=array->getPointer();
02363   std::fill(valsToFill,valsToFill+nbNodes,0.);
02364   const double *values=tmp->getArray()->getConstPointer();
02365   DataArrayInt *da=DataArrayInt::New();
02366   DataArrayInt *daInd=DataArrayInt::New();
02367   getReverseNodalConnectivity(da,daInd);
02368   const int *daPtr=da->getConstPointer();
02369   const int *daIPtr=daInd->getConstPointer();
02370   for(int i=0;i<nbNodes;i++)
02371     for(const int *cell=daPtr+daIPtr[i];cell!=daPtr+daIPtr[i+1];cell++)
02372       valsToFill[i]+=cst*values[*cell];
02373   ret->setMesh(this);
02374   da->decrRef();
02375   daInd->decrRef();
02376   ret->setArray(array);
02377   array->decrRef();
02378   tmp->decrRef();
02379   return ret;
02380 }
02381 
02386 MEDCouplingFieldDouble *MEDCouplingUMesh::buildOrthogonalField() const
02387 {
02388   if((getMeshDimension()!=2) && (getMeshDimension()!=1 || getSpaceDimension()!=2))
02389     throw INTERP_KERNEL::Exception("Expected a umesh with ( meshDim == 2 spaceDim == 2 or 3 ) or ( meshDim == 1 spaceDim == 2 ) !");
02390   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
02391   DataArrayDouble *array=DataArrayDouble::New();
02392   int nbOfCells=getNumberOfCells();
02393   int nbComp=getMeshDimension()+1;
02394   array->alloc(nbOfCells,nbComp);
02395   double *vals=array->getPointer();
02396   const int *connI=_nodal_connec_index->getConstPointer();
02397   const int *conn=_nodal_connec->getConstPointer();
02398   const double *coords=_coords->getConstPointer();
02399   if(getMeshDimension()==2)
02400     {
02401       if(getSpaceDimension()==3)
02402         {
02403           DataArrayDouble *loc=getBarycenterAndOwner();
02404           const double *locPtr=loc->getConstPointer();
02405           for(int i=0;i<nbOfCells;i++,vals+=3)
02406             {
02407               int offset=connI[i];
02408               INTERP_KERNEL::crossprod<3>(locPtr+3*i,coords+3*conn[offset+1],coords+3*conn[offset+2],vals);
02409               double n=INTERP_KERNEL::norm<3>(vals);
02410               std::transform(vals,vals+3,vals,std::bind2nd(std::multiplies<double>(),1./n));
02411             }
02412           loc->decrRef();
02413         }
02414       else
02415         {
02416           for(int i=0;i<nbOfCells;i++)
02417             { vals[3*i]=0.; vals[3*i+1]=0.; vals[3*i+2]=1.; }
02418         }
02419     }
02420   else//meshdimension==1
02421     {
02422       double tmp[2];
02423       for(int i=0;i<nbOfCells;i++)
02424         {
02425           int offset=connI[i];
02426           std::transform(coords+2*conn[offset+2],coords+2*conn[offset+2]+2,coords+2*conn[offset+1],tmp,std::minus<double>());
02427           double n=INTERP_KERNEL::norm<2>(tmp);
02428           std::transform(tmp,tmp+2,tmp,std::bind2nd(std::multiplies<double>(),1./n));
02429           *vals++=-tmp[1];
02430           *vals++=tmp[0];
02431         }
02432     }
02433   ret->setArray(array);
02434   array->decrRef();
02435   ret->setMesh(this);
02436   return ret;
02437 }
02438 
02443 MEDCouplingFieldDouble *MEDCouplingUMesh::buildPartOrthogonalField(const int *begin, const int *end) const
02444 {
02445   if((getMeshDimension()!=2) && (getMeshDimension()!=1 || getSpaceDimension()!=2))
02446     throw INTERP_KERNEL::Exception("Expected a umesh with ( meshDim == 2 spaceDim == 2 or 3 ) or ( meshDim == 1 spaceDim == 2 ) !");
02447   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
02448   DataArrayDouble *array=DataArrayDouble::New();
02449   std::size_t nbelems=std::distance(begin,end);
02450   int nbComp=getMeshDimension()+1;
02451   array->alloc((int)nbelems,nbComp);
02452   double *vals=array->getPointer();
02453   const int *connI=_nodal_connec_index->getConstPointer();
02454   const int *conn=_nodal_connec->getConstPointer();
02455   const double *coords=_coords->getConstPointer();
02456   if(getMeshDimension()==2)
02457     {
02458       if(getSpaceDimension()==3)
02459         {
02460           DataArrayDouble *loc=getPartBarycenterAndOwner(begin,end);
02461           const double *locPtr=loc->getConstPointer();
02462           for(const int *i=begin;i!=end;i++,vals+=3,locPtr+=3)
02463             {
02464               int offset=connI[*i];
02465               INTERP_KERNEL::crossprod<3>(locPtr,coords+3*conn[offset+1],coords+3*conn[offset+2],vals);
02466               double n=INTERP_KERNEL::norm<3>(vals);
02467               std::transform(vals,vals+3,vals,std::bind2nd(std::multiplies<double>(),1./n));
02468             }
02469           loc->decrRef();
02470         }
02471       else
02472         {
02473           for(std::size_t i=0;i<nbelems;i++)
02474             { vals[3*i]=0.; vals[3*i+1]=0.; vals[3*i+2]=1.; }
02475         }
02476     }
02477   else//meshdimension==1
02478     {
02479       double tmp[2];
02480       for(const int *i=begin;i!=end;i++)
02481         {
02482           int offset=connI[*i];
02483           std::transform(coords+2*conn[offset+2],coords+2*conn[offset+2]+2,coords+2*conn[offset+1],tmp,std::minus<double>());
02484           double n=INTERP_KERNEL::norm<2>(tmp);
02485           std::transform(tmp,tmp+2,tmp,std::bind2nd(std::multiplies<double>(),1./n));
02486           *vals++=-tmp[1];
02487           *vals++=tmp[0];
02488         }
02489     }
02490   ret->setArray(array);
02491   array->decrRef();
02492   ret->setMesh(this);
02493   return ret;
02494 }
02495 
02500 MEDCouplingFieldDouble *MEDCouplingUMesh::buildDirectionVectorField() const
02501 {
02502    if(getMeshDimension()!=1)
02503     throw INTERP_KERNEL::Exception("Expected a umesh with meshDim == 1 for buildDirectionVectorField !");
02504    if(_types.size()!=1 || *(_types.begin())!=INTERP_KERNEL::NORM_SEG2)
02505      throw INTERP_KERNEL::Exception("Expected a umesh with only NORM_SEG2 type of elements for buildDirectionVectorField !");
02506    MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
02507    DataArrayDouble *array=DataArrayDouble::New();
02508    int nbOfCells=getNumberOfCells();
02509    int spaceDim=getSpaceDimension();
02510    array->alloc(nbOfCells,spaceDim);
02511    double *pt=array->getPointer();
02512    const double *coo=getCoords()->getConstPointer();
02513    std::vector<int> conn;
02514    conn.reserve(2);
02515    for(int i=0;i<nbOfCells;i++)
02516      {
02517        conn.resize(0);
02518        getNodeIdsOfCell(i,conn);
02519        pt=std::transform(coo+conn[1]*spaceDim,coo+(conn[1]+1)*spaceDim,coo+conn[0]*spaceDim,pt,std::minus<double>());
02520      }
02521    ret->setArray(array);
02522    array->decrRef();
02523    ret->setMesh(this);
02524    return ret;   
02525 }
02526 
02541 MEDCouplingUMesh *MEDCouplingUMesh::buildSlice3D(const double *origin, const double *vec, double eps, DataArrayInt *&cellIds) const throw(INTERP_KERNEL::Exception)
02542 {
02543   checkFullyDefined();
02544   if(getMeshDimension()!=3 || getSpaceDimension()!=3)
02545     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
02546   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> candidates=getCellIdsCrossingPlane(origin,vec,eps);
02547   if(candidates->empty())
02548     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D : No 3D cells in this intercepts the specified plane considering bounding boxes !");
02549   std::vector<int> nodes;
02550   std::vector<int> cellIds2D,cellIds1D;
02551   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> subMesh=static_cast<MEDCouplingUMesh*>(buildPartOfMySelf(candidates->begin(),candidates->end(),false));
02552   subMesh->findNodesOnPlane(origin,vec,eps,nodes);
02553   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> desc1=DataArrayInt::New(),desc2=DataArrayInt::New();
02554   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> descIndx1=DataArrayInt::New(),descIndx2=DataArrayInt::New();
02555   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> revDesc1=DataArrayInt::New(),revDesc2=DataArrayInt::New();
02556   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> revDescIndx1=DataArrayInt::New(),revDescIndx2=DataArrayInt::New();
02557   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> mDesc2=subMesh->buildDescendingConnectivity(desc2,descIndx2,revDesc2,revDescIndx2);//meshDim==2 spaceDim==3
02558   revDesc2=0; revDescIndx2=0;
02559   mDesc2->fillCellIdsToKeepFromNodeIds(&nodes[0],&nodes[0]+nodes.size(),true,cellIds2D);
02560   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> mDesc1=mDesc2->buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1);//meshDim==1 spaceDim==3
02561   revDesc1=0; revDescIndx1=0;
02562   mDesc1->fillCellIdsToKeepFromNodeIds(&nodes[0],&nodes[0]+nodes.size(),true,cellIds1D);
02563   //
02564   std::vector<int> cut3DCurve(mDesc1->getNumberOfCells(),-2);
02565   for(std::vector<int>::const_iterator it=cellIds1D.begin();it!=cellIds1D.end();it++)
02566     cut3DCurve[*it]=-1;
02567   mDesc1->split3DCurveWithPlane(origin,vec,eps,cut3DCurve);
02568   std::vector< std::pair<int,int> > cut3DSurf(mDesc2->getNumberOfCells());
02569   AssemblyForSplitFrom3DCurve(cut3DCurve,nodes,mDesc2->getNodalConnectivity()->getConstPointer(),mDesc2->getNodalConnectivityIndex()->getConstPointer(),
02570                               mDesc1->getNodalConnectivity()->getConstPointer(),mDesc1->getNodalConnectivityIndex()->getConstPointer(),
02571                               desc1->getConstPointer(),descIndx1->getConstPointer(),cut3DSurf);
02572   std::vector<int> conn,connI,cellIds2;
02573   connI.push_back(0);
02574   subMesh->assemblyForSplitFrom3DSurf(cut3DSurf,desc2->getConstPointer(),descIndx2->getConstPointer(),conn,connI,cellIds2);
02575   if(cellIds2.empty())
02576     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D : No 3D cells in this intercepts the specified plane !");
02577   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret=MEDCouplingUMesh::New("Slice3D",2);
02578   ret->setCoords(mDesc1->getCoords());
02579   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c=DataArrayInt::New();
02580   c->alloc((int)conn.size(),1); std::copy(conn.begin(),conn.end(),c->getPointer());
02581   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cI=DataArrayInt::New();
02582   cI->alloc((int)connI.size(),1); std::copy(connI.begin(),connI.end(),cI->getPointer());
02583   ret->setConnectivity(c,cI,true);
02584   cellIds=candidates->selectByTupleId(&cellIds2[0],&cellIds2[0]+cellIds2.size());
02585   ret->incrRef();
02586   return ret;
02587 }
02588 
02603 MEDCouplingUMesh *MEDCouplingUMesh::buildSlice3DSurf(const double *origin, const double *vec, double eps, DataArrayInt *&cellIds) const throw(INTERP_KERNEL::Exception)
02604 {
02605   checkFullyDefined();
02606   if(getMeshDimension()!=2 || getSpaceDimension()!=3)
02607     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D works on umeshes with meshdim equal to 2 and spaceDim equal to 3 !");
02608   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> candidates=getCellIdsCrossingPlane(origin,vec,eps);
02609   if(candidates->empty())
02610     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D : No 3D cells in this intercepts the specified plane considering bounding boxes !");
02611   std::vector<int> nodes;
02612   std::vector<int> cellIds1D;
02613   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> subMesh=static_cast<MEDCouplingUMesh*>(buildPartOfMySelf(candidates->begin(),candidates->end(),false));
02614   subMesh->findNodesOnPlane(origin,vec,eps,nodes);
02615   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> desc1=DataArrayInt::New();
02616   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> descIndx1=DataArrayInt::New();
02617   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> revDesc1=DataArrayInt::New();
02618   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> revDescIndx1=DataArrayInt::New();
02619   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> mDesc1=subMesh->buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1);//meshDim==1 spaceDim==3
02620   mDesc1->fillCellIdsToKeepFromNodeIds(&nodes[0],&nodes[0]+nodes.size(),true,cellIds1D);
02621   //
02622   std::vector<int> cut3DCurve(mDesc1->getNumberOfCells(),-2);
02623   for(std::vector<int>::const_iterator it=cellIds1D.begin();it!=cellIds1D.end();it++)
02624     cut3DCurve[*it]=-1;
02625   mDesc1->split3DCurveWithPlane(origin,vec,eps,cut3DCurve);
02626   int ncellsSub=subMesh->getNumberOfCells();
02627   std::vector< std::pair<int,int> > cut3DSurf(ncellsSub);
02628   AssemblyForSplitFrom3DCurve(cut3DCurve,nodes,subMesh->getNodalConnectivity()->getConstPointer(),subMesh->getNodalConnectivityIndex()->getConstPointer(),
02629                               mDesc1->getNodalConnectivity()->getConstPointer(),mDesc1->getNodalConnectivityIndex()->getConstPointer(),
02630                               desc1->getConstPointer(),descIndx1->getConstPointer(),cut3DSurf);
02631   std::vector<int> conn,connI,cellIds2; connI.push_back(0);
02632   const int *nodal=subMesh->getNodalConnectivity()->getConstPointer();
02633   const int *nodalI=subMesh->getNodalConnectivityIndex()->getConstPointer();
02634   for(int i=0;i<ncellsSub;i++)
02635     {
02636       if(cut3DSurf[i].first!=-1 && cut3DSurf[i].second!=-1)
02637         {
02638           if(cut3DSurf[i].first!=-2)
02639             {
02640               conn.push_back((int)INTERP_KERNEL::NORM_SEG2); conn.push_back(cut3DSurf[i].first); conn.push_back(cut3DSurf[i].second);
02641               connI.push_back((int)conn.size());
02642               cellIds2.push_back(i);
02643             }
02644           else
02645             {
02646               int cellId3DSurf=cut3DSurf[i].second;
02647               int offset=nodalI[cellId3DSurf]+1;
02648               int nbOfEdges=nodalI[cellId3DSurf+1]-offset;
02649               for(int j=0;j<nbOfEdges;j++)
02650                 {
02651                   conn.push_back((int)INTERP_KERNEL::NORM_SEG2); conn.push_back(nodal[offset+j]); conn.push_back(nodal[offset+(j+1)%nbOfEdges]);
02652                   connI.push_back((int)conn.size());
02653                   cellIds2.push_back(cellId3DSurf);
02654                 }
02655             }
02656         }
02657     }
02658   if(cellIds2.empty())
02659     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3DSurf : No 3DSurf cells in this intercepts the specified plane !");
02660   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret=MEDCouplingUMesh::New("Slice3DSurf",1);
02661   ret->setCoords(mDesc1->getCoords());
02662   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c=DataArrayInt::New();
02663   c->alloc((int)conn.size(),1); std::copy(conn.begin(),conn.end(),c->getPointer());
02664   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cI=DataArrayInt::New();
02665   cI->alloc((int)connI.size(),1); std::copy(connI.begin(),connI.end(),cI->getPointer());
02666   ret->setConnectivity(c,cI,true);
02667   cellIds=candidates->selectByTupleId(&cellIds2[0],&cellIds2[0]+cellIds2.size());
02668   ret->incrRef();
02669   return ret;
02670 }
02671 
02678 DataArrayInt *MEDCouplingUMesh::getCellIdsCrossingPlane(const double *origin, const double *vec, double eps) const throw(INTERP_KERNEL::Exception)
02679 {
02680   checkFullyDefined();
02681   if(getSpaceDimension()!=3)
02682     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D works on umeshes with spaceDim equal to 3 !");
02683   double normm=sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
02684   if(normm<1e-6)
02685     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getCellIdsCrossingPlane : parameter 'vec' should have a norm2 greater than 1e-6 !");
02686   double vec2[3];
02687   vec2[0]=vec[1]; vec2[1]=-vec[0]; vec2[2]=0.;//vec2 is the result of cross product of vec with (0,0,1)
02688   double angle=acos(vec[2]/normm);
02689   std::vector<int> cellIds;
02690   double bbox[6];
02691   if(angle>eps)
02692     {
02693       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo=_coords->deepCpy();
02694       MEDCouplingPointSet::Rotate3DAlg(origin,vec2,angle,coo->getNumberOfTuples(),coo->getPointer());
02695       MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> mw=clone(false);//false -> shallow copy
02696       mw->setCoords(coo);
02697       mw->getBoundingBox(bbox);
02698       bbox[4]=origin[2]-eps; bbox[5]=origin[2]+eps;
02699       mw->getCellsInBoundingBox(bbox,eps,cellIds);
02700     }
02701   else
02702     {
02703       getBoundingBox(bbox);
02704       bbox[4]=origin[2]-eps; bbox[5]=origin[2]+eps;
02705       getCellsInBoundingBox(bbox,eps,cellIds);
02706     }
02707   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
02708   ret->alloc((int)cellIds.size(),1);
02709   std::copy(cellIds.begin(),cellIds.end(),ret->getPointer());
02710   ret->incrRef();
02711   return ret;
02712 }
02713 
02721 bool MEDCouplingUMesh::isContiguous1D() const throw(INTERP_KERNEL::Exception)
02722 {
02723   if(getMeshDimension()!=1)
02724     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::isContiguous1D : this method has a sense only for 1D mesh !");
02725   int nbCells=getNumberOfCells();
02726   if(nbCells<1)
02727     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::isContiguous1D : this method has a sense for non empty mesh !");
02728   const int *connI=_nodal_connec_index->getConstPointer();
02729   const int *conn=_nodal_connec->getConstPointer();
02730   int ref=conn[connI[0]+2];
02731   for(int i=1;i<nbCells;i++)
02732     {
02733       if(conn[connI[i]+1]!=ref)
02734         return false;
02735       ref=conn[connI[i]+2];
02736     }
02737   return true;
02738 }
02739 
02748 void MEDCouplingUMesh::project1D(const double *pt, const double *v, double eps, double *res) const
02749 {
02750   if(getMeshDimension()!=1)
02751     throw INTERP_KERNEL::Exception("Expected a umesh with meshDim == 1 for project1D !");
02752    if(_types.size()!=1 || *(_types.begin())!=INTERP_KERNEL::NORM_SEG2)
02753      throw INTERP_KERNEL::Exception("Expected a umesh with only NORM_SEG2 type of elements for project1D !");
02754    if(getSpaceDimension()!=3)
02755      throw INTERP_KERNEL::Exception("Expected a umesh with spaceDim==3 for project1D !");
02756    MEDCouplingFieldDouble *f=buildDirectionVectorField();
02757    const double *fPtr=f->getArray()->getConstPointer();
02758    double tmp[3];
02759    for(int i=0;i<getNumberOfCells();i++)
02760      {
02761        const double *tmp1=fPtr+3*i;
02762        tmp[0]=tmp1[1]*v[2]-tmp1[2]*v[1];
02763        tmp[1]=tmp1[2]*v[0]-tmp1[0]*v[2];
02764        tmp[2]=tmp1[0]*v[1]-tmp1[1]*v[0];
02765        double n1=INTERP_KERNEL::norm<3>(tmp);
02766        n1/=INTERP_KERNEL::norm<3>(tmp1);
02767        if(n1>eps)
02768          {
02769            f->decrRef();
02770            throw INTERP_KERNEL::Exception("UMesh::Projection 1D failed !");
02771          }
02772      }
02773    const double *coo=getCoords()->getConstPointer();
02774    for(int i=0;i<getNumberOfNodes();i++)
02775      {
02776        std::transform(coo+i*3,coo+i*3+3,pt,tmp,std::minus<double>());
02777        std::transform(tmp,tmp+3,v,tmp,std::multiplies<double>());
02778        res[i]=std::accumulate(tmp,tmp+3,0.);
02779      }
02780    f->decrRef();
02781 }
02782 
02789 int MEDCouplingUMesh::getCellContainingPoint(const double *pos, double eps) const
02790 {
02791   std::vector<int> elts;
02792   getCellsContainingPoint(pos,eps,elts);
02793   if(elts.empty())
02794     return -1;
02795   return elts.front();
02796 }
02797 
02803 void MEDCouplingUMesh::getCellsContainingPoint(const double *pos, double eps, std::vector<int>& elts) const
02804 {
02805   std::vector<int> eltsIndex;
02806   getCellsContainingPoints(pos,1,eps,elts,eltsIndex);
02807 }
02808 
02810 
02811 namespace ParaMEDMEM
02812 {
02813   template<const int SPACEDIMM>
02814   class DummyClsMCUG
02815   {
02816   public:
02817     static const int MY_SPACEDIM=SPACEDIMM;
02818     static const int MY_MESHDIM=8;
02819     typedef int MyConnType;
02820     static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
02821     // begin
02822     // useless, but for windows compilation ...
02823     const double* getCoordinatesPtr() const { return 0; }
02824     const int* getConnectivityPtr() const { return 0; }
02825     const int* getConnectivityIndexPtr() const { return 0; }
02826     INTERP_KERNEL::NormalizedCellType getTypeOfElement(int) const { return (INTERP_KERNEL::NormalizedCellType)0; }
02827     // end
02828   };
02829 
02830   INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge(INTERP_KERNEL::NormalizedCellType typ, std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >& mapp2, const int *bg)
02831   {
02832     INTERP_KERNEL::Edge *ret=0;
02833     switch(typ)
02834       {
02835       case INTERP_KERNEL::NORM_SEG2:
02836         {
02837           ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[1]].first);
02838           break;
02839         }
02840       case INTERP_KERNEL::NORM_SEG3:
02841         {
02842           INTERP_KERNEL::EdgeLin *e1=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[2]].first);
02843           INTERP_KERNEL::EdgeLin *e2=new INTERP_KERNEL::EdgeLin(mapp2[bg[2]].first,mapp2[bg[1]].first);
02844           INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
02845           bool colinearity=inters.areColinears();
02846           delete e1; delete e2;
02847           if(colinearity)
02848             ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]].first,mapp2[bg[1]].first);
02849           else
02850             ret=new INTERP_KERNEL::EdgeArcCircle(mapp2[bg[0]].first,mapp2[bg[2]].first,mapp2[bg[1]].first);
02851           mapp2[bg[2]].second=false;
02852           break;
02853         }
02854       default:
02855         throw INTERP_KERNEL::Exception("MEDCouplingUMeshBuildQPFromEdge : Expecting a mesh with spaceDim==2 and meshDim==1 !");
02856       }
02857     return ret;
02858   }
02859 
02865   INTERP_KERNEL::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMesh(const MEDCouplingUMesh *mDesc, const std::vector<int>& candidates, std::map<INTERP_KERNEL::Node *,int>& mapp) throw(INTERP_KERNEL::Exception)
02866   {
02867     mapp.clear();
02868     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.
02869     const double *coo=mDesc->getCoords()->getConstPointer();
02870     const int *c=mDesc->getNodalConnectivity()->getConstPointer();
02871     const int *cI=mDesc->getNodalConnectivityIndex()->getConstPointer();
02872     std::set<int> s;
02873     for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
02874       s.insert(c+cI[*it]+1,c+cI[(*it)+1]);
02875     for(std::set<int>::const_iterator it2=s.begin();it2!=s.end();it2++)
02876       {
02877         INTERP_KERNEL::Node *n=new INTERP_KERNEL::Node(coo[2*(*it2)],coo[2*(*it2)+1]);
02878         mapp2[*it2]=std::pair<INTERP_KERNEL::Node *,bool>(n,true);
02879       }
02880     INTERP_KERNEL::QuadraticPolygon *ret=new INTERP_KERNEL::QuadraticPolygon;
02881     for(std::vector<int>::const_iterator it=candidates.begin();it!=candidates.end();it++)
02882       {
02883         INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)c[cI[*it]];
02884         ret->pushBack(MEDCouplingUMeshBuildQPFromEdge(typ,mapp2,c+cI[*it]+1));
02885       }
02886     for(std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >::const_iterator it2=mapp2.begin();it2!=mapp2.end();it2++)
02887       {
02888         if((*it2).second.second)
02889           mapp[(*it2).second.first]=(*it2).first;
02890         ((*it2).second.first)->decrRef();
02891       }
02892     return ret;
02893   }
02894 
02895   INTERP_KERNEL::Node *MEDCouplingUMeshBuildQPNode(int nodeId, const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo)
02896   {
02897     if(nodeId>=offset2)
02898       {
02899         int locId=nodeId-offset2;
02900         return new INTERP_KERNEL::Node(addCoo[2*locId],addCoo[2*locId+1]);
02901       }
02902     if(nodeId>=offset1)
02903       {
02904         int locId=nodeId-offset1;
02905         return new INTERP_KERNEL::Node(coo2[2*locId],coo2[2*locId+1]);
02906       }
02907     return new INTERP_KERNEL::Node(coo1[2*nodeId],coo1[2*nodeId+1]);
02908   }
02909 
02910   void MEDCouplingUMeshBuildQPFromMesh3(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo,
02911                                         const int *desc1Bg, const int *desc1End, const std::vector<std::vector<int> >& intesctEdges1,
02912                                         /*output*/std::map<INTERP_KERNEL::Node *,int>& mapp, std::map<int,INTERP_KERNEL::Node *>& mappRev)
02913   {
02914     for(const int *desc1=desc1Bg;desc1!=desc1End;desc1++)
02915       {
02916         int eltId1=abs(*desc1)-1;
02917         for(std::vector<int>::const_iterator it1=intesctEdges1[eltId1].begin();it1!=intesctEdges1[eltId1].end();it1++)
02918           {
02919             std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.find(*it1);
02920             if(it==mappRev.end())
02921               {
02922                 INTERP_KERNEL::Node *node=MEDCouplingUMeshBuildQPNode(*it1,coo1,offset1,coo2,offset2,addCoo);
02923                 mapp[node]=*it1;
02924                 mappRev[*it1]=node;
02925               }
02926           }
02927       }
02928   }
02929 }
02930 
02932 
02933 template<int SPACEDIM>
02934 void MEDCouplingUMesh::getCellsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints,
02935                                                    double eps, std::vector<int>& elts, std::vector<int>& eltsIndex) const
02936 {
02937   std::vector<double> bbox;
02938   eltsIndex.resize(nbOfPoints+1);
02939   eltsIndex[0]=0;
02940   elts.clear();
02941   getBoundingBoxForBBTree(bbox);
02942   int nbOfCells=getNumberOfCells();
02943   const int *conn=_nodal_connec->getConstPointer();
02944   const int *connI=_nodal_connec_index->getConstPointer();
02945   double bb[2*SPACEDIM];
02946   BBTree<SPACEDIM,int> myTree(&bbox[0],0,0,nbOfCells,-eps);
02947   for(int i=0;i<nbOfPoints;i++)
02948     {
02949       eltsIndex[i+1]=eltsIndex[i];
02950       for(int j=0;j<SPACEDIM;j++)
02951         {
02952           bb[2*j]=pos[SPACEDIM*i+j];
02953           bb[2*j+1]=pos[SPACEDIM*i+j];
02954         }
02955       std::vector<int> candidates;
02956       myTree.getIntersectingElems(bb,candidates);
02957       for(std::vector<int>::const_iterator iter=candidates.begin();iter!=candidates.end();iter++)
02958         {
02959           int sz=connI[(*iter)+1]-connI[*iter]-1;
02960           if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCUG<SPACEDIM> >::isElementContainsPoint(pos+i*SPACEDIM,
02961                                                                                                (INTERP_KERNEL::NormalizedCellType)conn[connI[*iter]],
02962                                                                                                coords,conn+connI[*iter]+1,sz,eps))
02963             {
02964               eltsIndex[i+1]++;
02965               elts.push_back(*iter);
02966             }
02967         }
02968     }
02969 }
02970 
02980 void MEDCouplingUMesh::getCellsContainingPoints(const double *pos, int nbOfPoints, double eps,
02981                                                 std::vector<int>& elts, std::vector<int>& eltsIndex) const
02982 {
02983   int spaceDim=getSpaceDimension();
02984   int mDim=getMeshDimension();
02985   if(spaceDim==3)
02986     {
02987       if(mDim==3)
02988         {
02989           const double *coords=_coords->getConstPointer();
02990           getCellsContainingPointsAlg<3>(coords,pos,nbOfPoints,eps,elts,eltsIndex);
02991         }
02992       /*else if(mDim==2)
02993         {
02994           
02995         }*/
02996       else
02997         throw INTERP_KERNEL::Exception("For spaceDim==3 only meshDim==3 implemented for getelementscontainingpoints !");
02998     }
02999   else if(spaceDim==2)
03000     {
03001       if(mDim==2)
03002         {
03003           const double *coords=_coords->getConstPointer();
03004           getCellsContainingPointsAlg<2>(coords,pos,nbOfPoints,eps,elts,eltsIndex);
03005         }
03006       else
03007         throw INTERP_KERNEL::Exception("For spaceDim==2 only meshDim==2 implemented for getelementscontainingpoints !");
03008     }
03009   else if(spaceDim==1)
03010     {
03011       if(mDim==1)
03012         {
03013           const double *coords=_coords->getConstPointer();
03014           getCellsContainingPointsAlg<1>(coords,pos,nbOfPoints,eps,elts,eltsIndex);
03015         }
03016       else
03017         throw INTERP_KERNEL::Exception("For spaceDim==1 only meshDim==1 implemented for getelementscontainingpoints !");
03018     }
03019   else
03020     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getCellsContainingPoints : not managed for mdim not in [1,2,3] !");
03021 }
03022 
03029 void MEDCouplingUMesh::checkButterflyCells(std::vector<int>& cells, double eps) const
03030 {
03031   const char msg[]="Butterfly detection work only for 2D cells with spaceDim==2 or 3!";
03032   if(getMeshDimension()!=2)
03033     throw INTERP_KERNEL::Exception(msg);
03034   int spaceDim=getSpaceDimension();
03035   if(spaceDim!=2 && spaceDim!=3)
03036     throw INTERP_KERNEL::Exception(msg);
03037   const int *conn=_nodal_connec->getConstPointer();
03038   const int *connI=_nodal_connec_index->getConstPointer();
03039   int nbOfCells=getNumberOfCells();
03040   std::vector<double> cell2DinS2;
03041   for(int i=0;i<nbOfCells;i++)
03042     {
03043       int offset=connI[i];
03044       int nbOfNodesForCell=connI[i+1]-offset-1;
03045       if(nbOfNodesForCell<=3)
03046         continue;
03047       bool isQuad=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[offset]).isQuadratic();
03048       project2DCellOnXY(conn+offset+1,conn+connI[i+1],cell2DinS2);
03049       if(isButterfly2DCell(cell2DinS2,isQuad,eps))
03050         cells.push_back(i);
03051       cell2DinS2.clear();
03052     }
03053 }
03054 
03061 void MEDCouplingUMesh::findAndCorrectBadOriented3DExtrudedCells(std::vector<int>& cells) throw(INTERP_KERNEL::Exception)
03062 {
03063   const char msg[]="check3DCellsWellOriented detection works only for 3D cells !";
03064   if(getMeshDimension()!=3)
03065     throw INTERP_KERNEL::Exception(msg);
03066   int spaceDim=getSpaceDimension();
03067   if(spaceDim!=3)
03068     throw INTERP_KERNEL::Exception(msg);
03069   //
03070   int nbOfCells=getNumberOfCells();
03071   int *conn=_nodal_connec->getPointer();
03072   const int *connI=_nodal_connec_index->getConstPointer();
03073   const double *coo=getCoords()->getConstPointer();
03074   double vec0[3],vec1[3];
03075   for(int i=0;i<nbOfCells;i++)
03076     {
03077       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[connI[i]]);
03078       if(cm.isExtruded() && !cm.isDynamic() && !cm.isQuadratic())
03079         {
03080           INTERP_KERNEL::AutoPtr<int> tmp=new int[connI[i+1]-connI[i]-1];
03081           int nbOfNodes=cm.fillSonCellNodalConnectivity(0,conn+connI[i]+1,tmp);
03082           INTERP_KERNEL::areaVectorOfPolygon<int,INTERP_KERNEL::ALL_C_MODE>(tmp,nbOfNodes,coo,vec0);
03083           const double *pt0=coo+3*conn[connI[i]+1];
03084           const double *pt1=coo+3*conn[connI[i]+nbOfNodes+1];
03085           vec1[0]=pt0[0]-pt1[0]; vec1[1]=pt0[1]-pt1[1]; vec1[2]=pt0[2]-pt1[2];
03086           double dot=vec0[0]*vec1[0]+vec0[1]*vec1[1]+vec0[2]*vec1[2];
03087           if(dot<0)
03088             {
03089               cells.push_back(i);
03090               std::copy(conn+connI[i]+1,conn+connI[i+1],(int *)tmp);
03091               for(int j=1;j<nbOfNodes;j++)
03092                 {
03093                   conn[connI[i]+1+j]=tmp[nbOfNodes-j];
03094                   conn[connI[i]+1+j+nbOfNodes]=tmp[nbOfNodes+nbOfNodes-j];
03095                 }
03096             }
03097         }
03098     }
03099 }
03100 
03109 MEDCouplingUMesh *MEDCouplingUMesh::buildExtrudedMesh(const MEDCouplingUMesh *mesh1D, int policy)
03110 {
03111   checkFullyDefined();
03112   mesh1D->checkFullyDefined();
03113   if(!mesh1D->isContiguous1D())
03114     throw INTERP_KERNEL::Exception("buildExtrudedMesh : 1D mesh passed in parameter is not contiguous !");
03115   if(getSpaceDimension()!=mesh1D->getSpaceDimension())
03116     throw INTERP_KERNEL::Exception("Invalid call to buildExtrudedMesh this and mesh1D must have same dimension !");
03117   if((getMeshDimension()!=2 || getSpaceDimension()!=3) && (getMeshDimension()!=1 || getSpaceDimension()!=2))
03118     throw INTERP_KERNEL::Exception("Invalid 'this' for buildExtrudedMesh method : must be (meshDim==2 and spaceDim==3) or (meshDim==1 and spaceDim==2) !");
03119   if(mesh1D->getMeshDimension()!=1)
03120     throw INTERP_KERNEL::Exception("Invalid 'mesh1D' for buildExtrudedMesh method : must be meshDim==1 !");
03121   bool isQuad=false;
03122   if(isPresenceOfQuadratic())
03123     {
03124       if(mesh1D->isFullyQuadratic())
03125         isQuad=true;
03126       else
03127         throw INTERP_KERNEL::Exception("Invalid 2D mesh and 1D mesh because 2D mesh has quadratic cells and 1D is not fully quadratic !");
03128     }
03129   zipCoords();
03130   int oldNbOfNodes=getNumberOfNodes();
03131   DataArrayDouble *newCoords=0;
03132   switch(policy)
03133     {
03134     case 0:
03135       {
03136         newCoords=fillExtCoordsUsingTranslation(mesh1D,isQuad);
03137         break;
03138       }
03139     case 1:
03140       {
03141         newCoords=fillExtCoordsUsingTranslAndAutoRotation(mesh1D,isQuad);
03142         break;
03143       }
03144     default:
03145       throw INTERP_KERNEL::Exception("Not implemented extrusion policy : must be in (0) !");
03146     }
03147   setCoords(newCoords);
03148   newCoords->decrRef();
03149   MEDCouplingUMesh *ret=buildExtrudedMeshFromThisLowLev(oldNbOfNodes,isQuad);
03150   updateTime();
03151   return ret;
03152 }
03153 
03165 void MEDCouplingUMesh::split3DCurveWithPlane(const double *origin, const double *vec, double eps, std::vector<int>& cut3DCurve) throw(INTERP_KERNEL::Exception)
03166 {
03167   checkFullyDefined();
03168   if(getMeshDimension()!=1 || getSpaceDimension()!=3)
03169     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane works on umeshes with meshdim equal to 1 and spaceDim equal to 3 !");
03170   int ncells=getNumberOfCells();
03171   int nnodes=getNumberOfNodes();
03172   double vec2[3],vec3[3],vec4[3];
03173   double normm=sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]);
03174   if(normm<1e-6)
03175     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane : parameter 'vec' should have a norm2 greater than 1e-6 !");
03176   vec2[0]=vec[0]/normm; vec2[1]=vec[1]/normm; vec2[2]=vec[2]/normm;
03177   const int *conn=_nodal_connec->getConstPointer();
03178   const int *connI=_nodal_connec_index->getConstPointer();
03179   const double *coo=_coords->getConstPointer();
03180   std::vector<double> addCoo;
03181   for(int i=0;i<ncells;i++)
03182     {
03183       if(conn[connI[i]]==(int)INTERP_KERNEL::NORM_SEG2)
03184         {
03185           if(cut3DCurve[i]==-2)
03186             {
03187               int st=conn[connI[i]+1],endd=conn[connI[i]+2];
03188               vec3[0]=coo[3*endd]-coo[3*st]; vec3[1]=coo[3*endd+1]-coo[3*st+1]; vec3[2]=coo[3*endd+2]-coo[3*st+2];
03189               double normm2=sqrt(vec3[0]*vec3[0]+vec3[1]*vec3[1]+vec3[2]*vec3[2]);
03190               double colin=std::abs((vec3[0]*vec2[0]+vec3[1]*vec2[1]+vec3[2]*vec2[2])/normm2);
03191               if(colin>eps)//if colin<=eps -> current SEG2 is colinear to the input plane
03192                 {
03193                   const double *st2=coo+3*st;
03194                   vec4[0]=st2[0]-origin[0]; vec4[1]=st2[1]-origin[1]; vec4[2]=st2[2]-origin[2];
03195                   double pos=-(vec4[0]*vec2[0]+vec4[1]*vec2[1]+vec4[2]*vec2[2])/((vec3[0]*vec2[0]+vec3[1]*vec2[1]+vec3[2]*vec2[2]));
03196                   if(pos>eps && pos<1-eps)
03197                     {
03198                       int nNode=((int)addCoo.size())/3;
03199                       vec4[0]=st2[0]+pos*vec3[0]; vec4[1]=st2[1]+pos*vec3[1]; vec4[2]=st2[2]+pos*vec3[2];
03200                       addCoo.insert(addCoo.end(),vec4,vec4+3);
03201                       cut3DCurve[i]=nnodes+nNode;
03202                     }
03203                 }
03204             }
03205         }
03206       else
03207         throw INTERP_KERNEL::Exception("MEDCouplingUMesh::split3DCurveWithPlane : this method is only available for linear cell (NORM_SEG2) !");
03208     }
03209   if(!addCoo.empty())
03210     {
03211       int newNbOfNodes=nnodes+((int)addCoo.size())/3;
03212       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo2=DataArrayDouble::New();
03213       coo2->alloc(newNbOfNodes,3);
03214       double *tmp=coo2->getPointer();
03215       tmp=std::copy(_coords->begin(),_coords->end(),tmp);
03216       std::copy(addCoo.begin(),addCoo.end(),tmp);
03217       DataArrayDouble::SetArrayIn(coo2,_coords);
03218     }
03219 }
03220 
03226 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslation(const MEDCouplingUMesh *mesh1D, bool isQuad) const
03227 {
03228   int oldNbOfNodes=getNumberOfNodes();
03229   int nbOf1DCells=mesh1D->getNumberOfCells();
03230   int spaceDim=getSpaceDimension();
03231   DataArrayDouble *ret=DataArrayDouble::New();
03232   std::vector<bool> isQuads;
03233   int nbOfLevsInVec=isQuad?2*nbOf1DCells+1:nbOf1DCells+1;
03234   ret->alloc(oldNbOfNodes*nbOfLevsInVec,spaceDim);
03235   double *retPtr=ret->getPointer();
03236   const double *coords=getCoords()->getConstPointer();
03237   double *work=std::copy(coords,coords+spaceDim*oldNbOfNodes,retPtr);
03238   std::vector<int> v;
03239   std::vector<double> c;
03240   double vec[3];
03241   v.reserve(3);
03242   c.reserve(6);
03243   for(int i=0;i<nbOf1DCells;i++)
03244     {
03245       v.resize(0);
03246       mesh1D->getNodeIdsOfCell(i,v);
03247       c.resize(0);
03248       mesh1D->getCoordinatesOfNode(v[isQuad?2:1],c);
03249       mesh1D->getCoordinatesOfNode(v[0],c);
03250       std::transform(c.begin(),c.begin()+spaceDim,c.begin()+spaceDim,vec,std::minus<double>());
03251       for(int j=0;j<oldNbOfNodes;j++)
03252         work=std::transform(vec,vec+spaceDim,retPtr+spaceDim*(i*oldNbOfNodes+j),work,std::plus<double>());
03253       if(isQuad)
03254         {
03255           c.resize(0);
03256           mesh1D->getCoordinatesOfNode(v[1],c);
03257           mesh1D->getCoordinatesOfNode(v[0],c);
03258           std::transform(c.begin(),c.begin()+spaceDim,c.begin()+spaceDim,vec,std::minus<double>());
03259           for(int j=0;j<oldNbOfNodes;j++)
03260             work=std::transform(vec,vec+spaceDim,retPtr+spaceDim*(i*oldNbOfNodes+j),work,std::plus<double>());
03261         }
03262     }
03263   ret->copyStringInfoFrom(*getCoords());
03264   return ret;
03265 }
03266 
03272 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation(const MEDCouplingUMesh *mesh1D, bool isQuad) const throw(INTERP_KERNEL::Exception)
03273 {
03274   if(mesh1D->getSpaceDimension()==2)
03275     return fillExtCoordsUsingTranslAndAutoRotation2D(mesh1D,isQuad);
03276   if(mesh1D->getSpaceDimension()==3)
03277     return fillExtCoordsUsingTranslAndAutoRotation3D(mesh1D,isQuad);
03278   throw INTERP_KERNEL::Exception("Not implemented rotation and translation alg. for spacedim other than 2 and 3 !");
03279 }
03280 
03286 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D(const MEDCouplingUMesh *mesh1D, bool isQuad) const throw(INTERP_KERNEL::Exception)
03287 {
03288   if(isQuad)
03289     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D : not implemented for quadratic cells !");
03290   int oldNbOfNodes=getNumberOfNodes();
03291   int nbOf1DCells=mesh1D->getNumberOfCells();
03292   if(nbOf1DCells<2)
03293     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D : impossible to detect any angle of rotation ! Change extrusion policy 1->0 !");
03294   DataArrayDouble *ret=DataArrayDouble::New();
03295   int nbOfLevsInVec=nbOf1DCells+1;
03296   ret->alloc(oldNbOfNodes*nbOfLevsInVec,2);
03297   double *retPtr=ret->getPointer();
03298   retPtr=std::copy(getCoords()->getConstPointer(),getCoords()->getConstPointer()+getCoords()->getNbOfElems(),retPtr);
03299   MEDCouplingUMesh *tmp=MEDCouplingUMesh::New();
03300   DataArrayDouble *tmp2=getCoords()->deepCpy();
03301   tmp->setCoords(tmp2);
03302   tmp2->decrRef();
03303   const double *coo1D=mesh1D->getCoords()->getConstPointer();
03304   const int *conn1D=mesh1D->getNodalConnectivity()->getConstPointer();
03305   const int *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer();
03306   for(int i=1;i<nbOfLevsInVec;i++)
03307     {
03308       const double *begin=coo1D+2*conn1D[connI1D[i-1]+1];
03309       const double *end=coo1D+2*conn1D[connI1D[i-1]+2];
03310       const double *third=i+1<nbOfLevsInVec?coo1D+2*conn1D[connI1D[i]+2]:coo1D+2*conn1D[connI1D[i-2]+1];
03311       const double vec[2]={end[0]-begin[0],end[1]-begin[1]};
03312       tmp->translate(vec);
03313       double tmp3[2],radius,alpha,alpha0;
03314       const double *p0=i+1<nbOfLevsInVec?begin:third;
03315       const double *p1=i+1<nbOfLevsInVec?end:begin;
03316       const double *p2=i+1<nbOfLevsInVec?third:end;
03317       INTERP_KERNEL::EdgeArcCircle::GetArcOfCirclePassingThru(p0,p1,p2,tmp3,radius,alpha,alpha0);
03318       double cosangle=i+1<nbOfLevsInVec?(p0[0]-tmp3[0])*(p1[0]-tmp3[0])+(p0[1]-tmp3[1])*(p1[1]-tmp3[1]):(p2[0]-tmp3[0])*(p1[0]-tmp3[0])+(p2[1]-tmp3[1])*(p1[1]-tmp3[1]);
03319       double angle=acos(cosangle/(radius*radius));
03320       tmp->rotate(end,0,angle);
03321       retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr);
03322     }
03323   tmp->decrRef();
03324   return ret;
03325 }
03326 
03332 DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D(const MEDCouplingUMesh *mesh1D, bool isQuad) const throw(INTERP_KERNEL::Exception)
03333 {
03334   if(isQuad)
03335     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D : not implemented for quadratic cells !");
03336   int oldNbOfNodes=getNumberOfNodes();
03337   int nbOf1DCells=mesh1D->getNumberOfCells();
03338   if(nbOf1DCells<2)
03339     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D : impossible to detect any angle of rotation ! Change extrusion policy 1->0 !");
03340   DataArrayDouble *ret=DataArrayDouble::New();
03341   int nbOfLevsInVec=nbOf1DCells+1;
03342   ret->alloc(oldNbOfNodes*nbOfLevsInVec,3);
03343   double *retPtr=ret->getPointer();
03344   retPtr=std::copy(getCoords()->getConstPointer(),getCoords()->getConstPointer()+getCoords()->getNbOfElems(),retPtr);
03345   MEDCouplingUMesh *tmp=MEDCouplingUMesh::New();
03346   DataArrayDouble *tmp2=getCoords()->deepCpy();
03347   tmp->setCoords(tmp2);
03348   tmp2->decrRef();
03349   const double *coo1D=mesh1D->getCoords()->getConstPointer();
03350   const int *conn1D=mesh1D->getNodalConnectivity()->getConstPointer();
03351   const int *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer();
03352   for(int i=1;i<nbOfLevsInVec;i++)
03353     {
03354       const double *begin=coo1D+3*conn1D[connI1D[i-1]+1];
03355       const double *end=coo1D+3*conn1D[connI1D[i-1]+2];
03356       const double *third=i+1<nbOfLevsInVec?coo1D+3*conn1D[connI1D[i]+2]:coo1D+3*conn1D[connI1D[i-2]+1];
03357       const double vec[3]={end[0]-begin[0],end[1]-begin[1],end[2]-begin[2]};
03358       tmp->translate(vec);
03359       double tmp3[2],radius,alpha,alpha0;
03360       const double *p0=i+1<nbOfLevsInVec?begin:third;
03361       const double *p1=i+1<nbOfLevsInVec?end:begin;
03362       const double *p2=i+1<nbOfLevsInVec?third:end;
03363       double vecPlane[3]={
03364         (p1[1]-p0[1])*(p2[2]-p1[2])-(p1[2]-p0[2])*(p2[1]-p1[1]),
03365         (p1[2]-p0[2])*(p2[0]-p1[0])-(p1[0]-p0[0])*(p2[2]-p1[2]),
03366         (p1[0]-p0[0])*(p2[1]-p1[1])-(p1[1]-p0[1])*(p2[0]-p1[0]),
03367       };
03368       double norm=sqrt(vecPlane[0]*vecPlane[0]+vecPlane[1]*vecPlane[1]+vecPlane[2]*vecPlane[2]);
03369       if(norm>1.e-7)
03370         {
03371           vecPlane[0]/=norm; vecPlane[1]/=norm; vecPlane[2]/=norm;
03372           double norm2=sqrt(vecPlane[0]*vecPlane[0]+vecPlane[1]*vecPlane[1]);
03373           double vec2[2]={vecPlane[1]/norm2,-vecPlane[0]/norm2};
03374           double s2=norm2;
03375           double c2=cos(asin(s2));
03376           double m[3][3]={
03377             {vec2[0]*vec2[0]*(1-c2)+c2, vec2[0]*vec2[1]*(1-c2), vec2[1]*s2},
03378             {vec2[0]*vec2[1]*(1-c2), vec2[1]*vec2[1]*(1-c2)+c2, -vec2[0]*s2},
03379             {-vec2[1]*s2, vec2[0]*s2, c2}
03380           };
03381           double p0r[3]={m[0][0]*p0[0]+m[0][1]*p0[1]+m[0][2]*p0[2], m[1][0]*p0[0]+m[1][1]*p0[1]+m[1][2]*p0[2], m[2][0]*p0[0]+m[2][1]*p0[1]+m[2][2]*p0[2]};
03382           double p1r[3]={m[0][0]*p1[0]+m[0][1]*p1[1]+m[0][2]*p1[2], m[1][0]*p1[0]+m[1][1]*p1[1]+m[1][2]*p1[2], m[2][0]*p1[0]+m[2][1]*p1[1]+m[2][2]*p1[2]};
03383           double p2r[3]={m[0][0]*p2[0]+m[0][1]*p2[1]+m[0][2]*p2[2], m[1][0]*p2[0]+m[1][1]*p2[1]+m[1][2]*p2[2], m[2][0]*p2[0]+m[2][1]*p2[1]+m[2][2]*p2[2]};
03384           INTERP_KERNEL::EdgeArcCircle::GetArcOfCirclePassingThru(p0r,p1r,p2r,tmp3,radius,alpha,alpha0);
03385           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]);
03386           double angle=acos(cosangle/(radius*radius));
03387           tmp->rotate(end,vecPlane,angle);
03388           
03389         }
03390       retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr);
03391     }
03392   tmp->decrRef();
03393   return ret;
03394 }
03395 
03402 MEDCouplingUMesh *MEDCouplingUMesh::buildExtrudedMeshFromThisLowLev(int nbOfNodesOf1Lev, bool isQuad) const
03403 {
03404   int nbOf1DCells=getNumberOfNodes()/nbOfNodesOf1Lev-1;
03405   int nbOf2DCells=getNumberOfCells();
03406   int nbOf3DCells=nbOf2DCells*nbOf1DCells;
03407   MEDCouplingUMesh *ret=MEDCouplingUMesh::New("Extruded",getMeshDimension()+1);
03408   const int *conn=_nodal_connec->getConstPointer();
03409   const int *connI=_nodal_connec_index->getConstPointer();
03410   DataArrayInt *newConn=DataArrayInt::New();
03411   DataArrayInt *newConnI=DataArrayInt::New();
03412   newConnI->alloc(nbOf3DCells+1,1);
03413   int *newConnIPtr=newConnI->getPointer();
03414   *newConnIPtr++=0;
03415   std::vector<int> newc;
03416   for(int j=0;j<nbOf2DCells;j++)
03417     {
03418       AppendExtrudedCell(conn+connI[j],conn+connI[j+1],nbOfNodesOf1Lev,isQuad,newc);
03419       *newConnIPtr++=(int)newc.size();
03420     }
03421   newConn->alloc((int)(newc.size())*nbOf1DCells,1);
03422   int *newConnPtr=newConn->getPointer();
03423   int deltaPerLev=isQuad?2*nbOfNodesOf1Lev:nbOfNodesOf1Lev;
03424   newConnIPtr=newConnI->getPointer();
03425   for(int iz=0;iz<nbOf1DCells;iz++)
03426     {
03427       if(iz!=0)
03428         std::transform(newConnIPtr+1,newConnIPtr+1+nbOf2DCells,newConnIPtr+1+iz*nbOf2DCells,std::bind2nd(std::plus<int>(),newConnIPtr[iz*nbOf2DCells]));
03429       for(std::vector<int>::const_iterator iter=newc.begin();iter!=newc.end();iter++,newConnPtr++)
03430         {
03431           int icell=(int)(iter-newc.begin());
03432           if(std::find(newConnIPtr,newConnIPtr+nbOf2DCells,icell)==newConnIPtr+nbOf2DCells)
03433             {
03434               if(*iter!=-1)
03435                 *newConnPtr=(*iter)+iz*deltaPerLev;
03436               else
03437                 *newConnPtr=-1;
03438             }
03439           else
03440             *newConnPtr=(*iter);
03441         }
03442     }
03443   ret->setConnectivity(newConn,newConnI,true);
03444   newConn->decrRef();
03445   newConnI->decrRef();
03446   ret->setCoords(getCoords());
03447   return ret;
03448 }
03449 
03453 bool MEDCouplingUMesh::isFullyQuadratic() const
03454 {
03455   checkFullyDefined();
03456   bool ret=true;
03457   int nbOfCells=getNumberOfCells();
03458   for(int i=0;i<nbOfCells && ret;i++)
03459     {
03460       INTERP_KERNEL::NormalizedCellType type=getTypeOfCell(i);
03461       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
03462       ret=cm.isQuadratic();
03463     }
03464   return ret;
03465 }
03466 
03470 bool MEDCouplingUMesh::isPresenceOfQuadratic() const
03471 {
03472   checkFullyDefined();
03473   bool ret=false;
03474   int nbOfCells=getNumberOfCells();
03475   for(int i=0;i<nbOfCells && !ret;i++)
03476     {
03477       INTERP_KERNEL::NormalizedCellType type=getTypeOfCell(i);
03478       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
03479       ret=cm.isQuadratic();
03480     }
03481   return ret;
03482 }
03483 
03488 void MEDCouplingUMesh::convertQuadraticCellsToLinear() throw(INTERP_KERNEL::Exception)
03489 {
03490   checkFullyDefined();
03491   int nbOfCells=getNumberOfCells();
03492   int delta=0;
03493   for(int i=0;i<nbOfCells;i++)
03494     {
03495       INTERP_KERNEL::NormalizedCellType type=getTypeOfCell(i);
03496       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
03497       if(cm.isQuadratic())
03498         {
03499           INTERP_KERNEL::NormalizedCellType typel=cm.getLinearType();
03500           const INTERP_KERNEL::CellModel& cml=INTERP_KERNEL::CellModel::GetCellModel(typel);
03501           delta+=cm.getNumberOfNodes()-cml.getNumberOfNodes();
03502         }
03503     }
03504   if(delta==0)
03505     return ;
03506   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConn=DataArrayInt::New();
03507   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConnI=DataArrayInt::New();
03508   newConn->alloc(getMeshLength()-delta,1);
03509   newConnI->alloc(nbOfCells+1,1);
03510   const int *icptr=_nodal_connec->getConstPointer();
03511   const int *iciptr=_nodal_connec_index->getConstPointer();
03512   int *ocptr=newConn->getPointer();
03513   int *ociptr=newConnI->getPointer();
03514   *ociptr=0;
03515   _types.clear();
03516   for(int i=0;i<nbOfCells;i++,ociptr++)
03517     {
03518       INTERP_KERNEL::NormalizedCellType type=getTypeOfCell(i);
03519       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
03520       if(!cm.isQuadratic())
03521         {
03522           _types.insert(type);
03523           ocptr=std::copy(icptr+iciptr[i],icptr+iciptr[i+1],ocptr);
03524           ociptr[1]=ociptr[0]+iciptr[i+1]-iciptr[i];
03525         }
03526       else
03527         {
03528           INTERP_KERNEL::NormalizedCellType typel=cm.getLinearType();
03529           _types.insert(typel);
03530           const INTERP_KERNEL::CellModel& cml=INTERP_KERNEL::CellModel::GetCellModel(typel);
03531           int newNbOfNodes=cml.getNumberOfNodes();
03532           *ocptr++=(int)typel;
03533           ocptr=std::copy(icptr+iciptr[i]+1,icptr+iciptr[i]+newNbOfNodes+1,ocptr);
03534           ociptr[1]=ociptr[0]+newNbOfNodes+1;
03535         }
03536     }
03537   setConnectivity(newConn,newConnI,false);
03538 }
03539 
03548 void MEDCouplingUMesh::tessellate2D(double eps) throw(INTERP_KERNEL::Exception)
03549 {
03550   checkFullyDefined();
03551   if(getMeshDimension()!=2 || getSpaceDimension()!=2)  
03552     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2D works on umeshes with meshdim equal to 2 and spaceDim equal to 2 too!");
03553   double epsa=fabs(eps);
03554   if(epsa<std::numeric_limits<double>::min())
03555     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2DCurve : epsilon is null ! Please specify a higher epsilon. If too tiny it can lead to a huge amount of nodes and memory !");
03556   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> desc1=DataArrayInt::New();
03557   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> descIndx1=DataArrayInt::New();
03558   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> revDesc1=DataArrayInt::New();
03559   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> revDescIndx1=DataArrayInt::New();
03560   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> mDesc=buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
03561   revDesc1=0; revDescIndx1=0;
03562   mDesc->tessellate2DCurve(eps);
03563   subDivide2DMesh(mDesc->_nodal_connec->getConstPointer(),mDesc->_nodal_connec_index->getConstPointer(),desc1->getConstPointer(),descIndx1->getConstPointer());
03564   setCoords(mDesc->getCoords());
03565 }
03566 
03575 void MEDCouplingUMesh::tessellate2DCurve(double eps) throw(INTERP_KERNEL::Exception)
03576 {
03577   checkFullyDefined();
03578   if(getMeshDimension()!=1 || getSpaceDimension()!=2)
03579     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2DCurve works on umeshes with meshdim equal to 1 and spaceDim equal to 2 too!");
03580   double epsa=fabs(eps);
03581   if(epsa<std::numeric_limits<double>::min())
03582     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2DCurve : epsilon is null ! Please specify a higher epsilon. If too tiny it can lead to a huge amount of nodes and memory !");
03583   INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=1.e-10;
03584   int nbCells=getNumberOfCells();
03585   int nbNodes=getNumberOfNodes();
03586   const int *conn=_nodal_connec->getConstPointer();
03587   const int *connI=_nodal_connec_index->getConstPointer();
03588   const double *coords=_coords->getConstPointer();
03589   std::vector<double> addCoo;
03590   std::vector<int> newConn;
03591   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConnI=DataArrayInt::New();
03592   newConnI->alloc(nbCells+1,1);
03593   int *newConnIPtr=newConnI->getPointer();
03594   *newConnIPtr=0;
03595   int tmp1[3];
03596   INTERP_KERNEL::Node *tmp2[3];
03597   std::set<INTERP_KERNEL::NormalizedCellType> types;
03598   for(int i=0;i<nbCells;i++,newConnIPtr++)
03599     {
03600       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[connI[i]]);
03601       if(cm.isQuadratic())
03602         {//assert(connI[i+1]-connI[i]-1==3)
03603           tmp1[0]=conn[connI[i]+1+0]; tmp1[1]=conn[connI[i]+1+1]; tmp1[2]=conn[connI[i]+1+2];
03604           tmp2[0]=new INTERP_KERNEL::Node(coords[2*tmp1[0]],coords[2*tmp1[0]+1]);
03605           tmp2[1]=new INTERP_KERNEL::Node(coords[2*tmp1[1]],coords[2*tmp1[1]+1]);
03606           tmp2[2]=new INTERP_KERNEL::Node(coords[2*tmp1[2]],coords[2*tmp1[2]+1]);
03607           INTERP_KERNEL::EdgeArcCircle *eac=INTERP_KERNEL::EdgeArcCircle::BuildFromNodes(tmp2[0],tmp2[2],tmp2[1]);
03608           if(eac)
03609             {
03610               eac->tesselate(tmp1,nbNodes,epsa,newConn,addCoo);
03611               types.insert((INTERP_KERNEL::NormalizedCellType)newConn[newConnIPtr[0]]);
03612               delete eac;
03613               newConnIPtr[1]=(int)newConn.size();
03614             }
03615           else
03616             {
03617               types.insert(INTERP_KERNEL::NORM_SEG2);
03618               newConn.push_back(INTERP_KERNEL::NORM_SEG2);
03619               newConn.insert(newConn.end(),conn+connI[i]+1,conn+connI[i]+3);
03620               newConnIPtr[1]=newConnIPtr[0]+3;
03621             }
03622         }
03623       else
03624         {
03625           types.insert((INTERP_KERNEL::NormalizedCellType)conn[connI[i]]);
03626           newConn.insert(newConn.end(),conn+connI[i],conn+connI[i+1]);
03627           newConnIPtr[1]=newConnIPtr[0]+3;
03628         }
03629     }
03630   if(addCoo.empty() && ((int)newConn.size())==_nodal_connec->getNumberOfTuples())//nothing happens during tasselation : no update needed
03631     return ;
03632   _types=types;
03633   DataArrayInt::SetArrayIn(newConnI,_nodal_connec_index);
03634   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConnArr=DataArrayInt::New();
03635   newConnArr->alloc((int)newConn.size(),1);
03636   std::copy(newConn.begin(),newConn.end(),newConnArr->getPointer());
03637   DataArrayInt::SetArrayIn(newConnArr,_nodal_connec);
03638   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=DataArrayDouble::New();
03639   newCoords->alloc(nbNodes+((int)addCoo.size())/2,2);
03640   double *work=std::copy(_coords->begin(),_coords->end(),newCoords->getPointer());
03641   std::copy(addCoo.begin(),addCoo.end(),work);
03642   DataArrayDouble::SetArrayIn(newCoords,_coords);
03643   updateTime();
03644 }
03645 
03655 DataArrayInt *MEDCouplingUMesh::simplexize(int policy) throw(INTERP_KERNEL::Exception)
03656 {
03657   switch(policy)
03658     {
03659     case 0:
03660       return simplexizePol0();
03661     case 1:
03662       return simplexizePol1();
03663     default:
03664       throw INTERP_KERNEL::Exception("MEDCouplingUMesh::simplexize : unrecognized policy ! Must be 0 or 1 !");
03665     }
03666 }
03667 
03668 bool MEDCouplingUMesh::areOnlySimplexCells() const throw(INTERP_KERNEL::Exception)
03669 {
03670   checkFullyDefined();
03671   if(getMeshDimension()<1)
03672     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::areOnlySimplexCells : only available with meshes having a meshdim >= 1 !");
03673   int nbCells=getNumberOfCells();
03674   const int *conn=_nodal_connec->getConstPointer();
03675   const int *connI=_nodal_connec_index->getConstPointer();
03676   for(int i=0;i<nbCells;i++)
03677     {
03678       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[connI[i]]);
03679       if(!cm.isSimplex())
03680         return false;
03681     }
03682   return true;
03683 }
03684 
03688 DataArrayInt *MEDCouplingUMesh::simplexizePol0() throw(INTERP_KERNEL::Exception)
03689 {
03690   if(getMeshDimension()!=2)
03691     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::simplexizePol0 : this policy is only available for mesh with meshdim == 2 !");
03692   int nbOfCells=getNumberOfCells();
03693   DataArrayInt *ret=DataArrayInt::New();
03694   int nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_QUAD4);
03695   ret->alloc(nbOfCells+nbOfCutCells,1);
03696   if(nbOfCutCells==0)
03697     {
03698       ret->iota(0);
03699       return ret;
03700     }
03701   int *retPt=ret->getPointer();
03702   DataArrayInt *newConn=DataArrayInt::New();
03703   DataArrayInt *newConnI=DataArrayInt::New();
03704   newConnI->alloc(nbOfCells+nbOfCutCells+1,1);
03705   newConn->alloc(getMeshLength()+3*nbOfCutCells,1);
03706   int *pt=newConn->getPointer();
03707   int *ptI=newConnI->getPointer();
03708   ptI[0]=0;
03709   const int *oldc=_nodal_connec->getConstPointer();
03710   const int *ci=_nodal_connec_index->getConstPointer();
03711   for(int i=0;i<nbOfCells;i++,ci++)
03712     {
03713       if((INTERP_KERNEL::NormalizedCellType)oldc[ci[0]]==INTERP_KERNEL::NORM_QUAD4)
03714         {
03715           const int tmp[8]={(int)INTERP_KERNEL::NORM_TRI3,oldc[ci[0]+1],oldc[ci[0]+2],oldc[ci[0]+3],
03716                             (int)INTERP_KERNEL::NORM_TRI3,oldc[ci[0]+1],oldc[ci[0]+3],oldc[ci[0]+4]};
03717           pt=std::copy(tmp,tmp+8,pt);
03718           ptI[1]=ptI[0]+4;
03719           ptI[2]=ptI[0]+8;
03720           *retPt++=i;
03721           *retPt++=i;
03722           ptI+=2;
03723         }
03724       else
03725         {
03726           pt=std::copy(oldc+ci[0],oldc+ci[1],pt);
03727           ptI[1]=ptI[0]+ci[1]-ci[0];
03728           ptI++;
03729           *retPt++=i;
03730         }
03731     }
03732   _nodal_connec->decrRef();
03733   _nodal_connec=newConn;
03734   _nodal_connec_index->decrRef();
03735   _nodal_connec_index=newConnI;
03736   computeTypes();
03737   updateTime();
03738   return ret;
03739 }
03740 
03744 DataArrayInt *MEDCouplingUMesh::simplexizePol1() throw(INTERP_KERNEL::Exception)
03745 {
03746   if(getMeshDimension()!=2)
03747     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::simplexizePol0 : this policy is only available for mesh with meshdim == 2 !");
03748   int nbOfCells=getNumberOfCells();
03749   DataArrayInt *ret=DataArrayInt::New();
03750   int nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_QUAD4);
03751   ret->alloc(nbOfCells+nbOfCutCells,1);
03752   if(nbOfCutCells==0)
03753     {
03754       ret->iota(0);
03755       return ret;
03756     }
03757   int *retPt=ret->getPointer();
03758   DataArrayInt *newConn=DataArrayInt::New();
03759   DataArrayInt *newConnI=DataArrayInt::New();
03760   newConnI->alloc(nbOfCells+nbOfCutCells+1,1);
03761   newConn->alloc(getMeshLength()+3*nbOfCutCells,1);
03762   int *pt=newConn->getPointer();
03763   int *ptI=newConnI->getPointer();
03764   ptI[0]=0;
03765   const int *oldc=_nodal_connec->getConstPointer();
03766   const int *ci=_nodal_connec_index->getConstPointer();
03767   for(int i=0;i<nbOfCells;i++,ci++)
03768     {
03769       if((INTERP_KERNEL::NormalizedCellType)oldc[ci[0]]==INTERP_KERNEL::NORM_QUAD4)
03770         {
03771           const int tmp[8]={(int)INTERP_KERNEL::NORM_TRI3,oldc[ci[0]+1],oldc[ci[0]+2],oldc[ci[0]+4],
03772                             (int)INTERP_KERNEL::NORM_TRI3,oldc[ci[0]+2],oldc[ci[0]+3],oldc[ci[0]+4]};
03773           pt=std::copy(tmp,tmp+8,pt);
03774           ptI[1]=ptI[0]+4;
03775           ptI[2]=ptI[0]+8;
03776           *retPt++=i;
03777           *retPt++=i;
03778           ptI+=2;
03779         }
03780       else
03781         {
03782           pt=std::copy(oldc+ci[0],oldc+ci[1],pt);
03783           ptI[1]=ptI[0]+ci[1]-ci[0];
03784           ptI++;
03785           *retPt++=i;
03786         }
03787     }
03788   _nodal_connec->decrRef();
03789   _nodal_connec=newConn;
03790   _nodal_connec_index->decrRef();
03791   _nodal_connec_index=newConnI;
03792   computeTypes();
03793   updateTime();
03794   return ret;
03795 }
03796 
03805 void MEDCouplingUMesh::subDivide2DMesh(const int *nodeSubdived, const int *nodeIndxSubdived, const int *desc, const int *descIndex) throw(INTERP_KERNEL::Exception)
03806 {
03807   checkFullyDefined();
03808   if(getMeshDimension()!=2)
03809     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::subDivide2DMesh : works only on umesh with meshdim==2 !");
03810   int nbOfCells=getNumberOfCells();
03811   int *connI=_nodal_connec_index->getPointer();
03812   int newConnLgth=0;
03813   for(int i=0;i<nbOfCells;i++,connI++)
03814     {
03815       int offset=descIndex[i];
03816       int nbOfEdges=descIndex[i+1]-offset;
03817       //
03818       bool ddirect=desc[offset+nbOfEdges-1]>0;
03819       int eedgeId=std::abs(desc[offset+nbOfEdges-1])-1;
03820       int ref=ddirect?nodeSubdived[nodeIndxSubdived[eedgeId+1]-1]:nodeSubdived[nodeIndxSubdived[eedgeId]+1];
03821       for(int j=0;j<nbOfEdges;j++)
03822         {
03823           bool direct=desc[offset+j]>0;
03824           int edgeId=std::abs(desc[offset+j])-1;
03825           if(!INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodeSubdived[nodeIndxSubdived[edgeId]]).isQuadratic())
03826             {
03827               int id1=nodeSubdived[nodeIndxSubdived[edgeId]+1];
03828               int id2=nodeSubdived[nodeIndxSubdived[edgeId+1]-1];
03829               int ref2=direct?id1:id2;
03830               if(ref==ref2)
03831                 {
03832                   int nbOfSubNodes=nodeIndxSubdived[edgeId+1]-nodeIndxSubdived[edgeId]-1;
03833                   newConnLgth+=nbOfSubNodes-1;
03834                   ref=direct?id2:id1;
03835                 }
03836               else
03837                 {
03838                   std::ostringstream oss; oss << "MEDCouplingUMesh::subDivide2DMesh : On polygon #" << i << " edgeid #" << j << " subedges mismatch : end subedge k!=start subedge k+1 !";
03839                   throw INTERP_KERNEL::Exception(oss.str().c_str());
03840                 }
03841             }
03842           else
03843             {
03844               throw INTERP_KERNEL::Exception("MEDCouplingUMesh::subDivide2DMesh : this method only subdivides into linear edges !");
03845             }
03846         }
03847       newConnLgth++;//+1 is for cell type
03848       connI[1]=newConnLgth;
03849     }
03850   //
03851   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConn=DataArrayInt::New();
03852   newConn->alloc(newConnLgth,1);
03853   int *work=newConn->getPointer();
03854   for(int i=0;i<nbOfCells;i++)
03855     {
03856       *work++=INTERP_KERNEL::NORM_POLYGON;
03857       int offset=descIndex[i];
03858       int nbOfEdges=descIndex[i+1]-offset;
03859       for(int j=0;j<nbOfEdges;j++)
03860         {
03861           bool direct=desc[offset+j]>0;
03862           int edgeId=std::abs(desc[offset+j])-1;
03863           if(direct)
03864             work=std::copy(nodeSubdived+nodeIndxSubdived[edgeId]+1,nodeSubdived+nodeIndxSubdived[edgeId+1]-1,work);
03865           else
03866             {
03867               int nbOfSubNodes=nodeIndxSubdived[edgeId+1]-nodeIndxSubdived[edgeId]-1;
03868               std::reverse_iterator<const int *> it(nodeSubdived+nodeIndxSubdived[edgeId+1]);
03869               work=std::copy(it,it+nbOfSubNodes-1,work);
03870             }
03871         }
03872     }
03873   DataArrayInt::SetArrayIn(newConn,_nodal_connec);
03874   _types.clear();
03875   if(nbOfCells>0)
03876     _types.insert(INTERP_KERNEL::NORM_POLYGON);
03877 }
03878 
03888 void MEDCouplingUMesh::convertDegeneratedCells() throw(INTERP_KERNEL::Exception)
03889 {
03890   checkFullyDefined();
03891   if(getMeshDimension()<=1)
03892     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::convertDegeneratedCells works on umeshes with meshdim equals to 2 or 3 !");
03893   int nbOfCells=getNumberOfCells();
03894   if(nbOfCells<1)
03895     return ;
03896   int initMeshLgth=getMeshLength();
03897   int *conn=_nodal_connec->getPointer();
03898   int *index=_nodal_connec_index->getPointer();
03899   int posOfCurCell=0;
03900   int newPos=0;
03901   int lgthOfCurCell;
03902   for(int i=0;i<nbOfCells;i++)
03903     {
03904       lgthOfCurCell=index[i+1]-posOfCurCell;
03905       INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)conn[posOfCurCell];
03906       int newLgth;
03907       INTERP_KERNEL::NormalizedCellType newType=INTERP_KERNEL::CellSimplify::simplifyDegeneratedCell(type,conn+posOfCurCell+1,lgthOfCurCell-1,
03908                                                                                                      conn+newPos+1,newLgth);
03909       conn[newPos]=newType;
03910       newPos+=newLgth+1;
03911       posOfCurCell=index[i+1];
03912       index[i+1]=newPos;
03913     }
03914   if(newPos!=initMeshLgth)
03915     _nodal_connec->reAlloc(newPos);
03916   computeTypes();
03917 }
03918 
03925 void MEDCouplingUMesh::are2DCellsNotCorrectlyOriented(const double *vec, bool polyOnly, std::vector<int>& cells) const throw(INTERP_KERNEL::Exception)
03926 {
03927   if(getMeshDimension()!=2 || getSpaceDimension()!=3)
03928     throw INTERP_KERNEL::Exception("Invalid mesh to apply are2DCellsNotCorrectlyOriented on it : must be meshDim==2 and spaceDim==3 !");
03929   int nbOfCells=getNumberOfCells();
03930   const int *conn=_nodal_connec->getConstPointer();
03931   const int *connI=_nodal_connec_index->getConstPointer();
03932   const double *coordsPtr=_coords->getConstPointer();
03933   for(int i=0;i<nbOfCells;i++)
03934     {
03935       INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)conn[connI[i]];
03936       if(!polyOnly || (type==INTERP_KERNEL::NORM_POLYGON || type==INTERP_KERNEL::NORM_QPOLYG))
03937         {
03938           bool isQuadratic=INTERP_KERNEL::CellModel::GetCellModel(type).isQuadratic();
03939           if(!IsPolygonWellOriented(isQuadratic,vec,conn+connI[i]+1,conn+connI[i+1],coordsPtr))
03940             cells.push_back(i);
03941         }
03942     }
03943 }
03944 
03950 void MEDCouplingUMesh::orientCorrectly2DCells(const double *vec, bool polyOnly) throw(INTERP_KERNEL::Exception)
03951 {
03952   if(getMeshDimension()!=2 || getSpaceDimension()!=3)
03953     throw INTERP_KERNEL::Exception("Invalid mesh to apply orientCorrectly2DCells on it : must be meshDim==2 and spaceDim==3 !");
03954   int nbOfCells=getNumberOfCells();
03955   int *conn=_nodal_connec->getPointer();
03956   const int *connI=_nodal_connec_index->getConstPointer();
03957   const double *coordsPtr=_coords->getConstPointer();
03958   bool isModified=false;
03959   for(int i=0;i<nbOfCells;i++)
03960     {
03961       INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)conn[connI[i]];
03962       if(!polyOnly || (type==INTERP_KERNEL::NORM_POLYGON || type==INTERP_KERNEL::NORM_QPOLYG))
03963         {
03964           bool isQuadratic=INTERP_KERNEL::CellModel::GetCellModel(type).isQuadratic();
03965           if(!IsPolygonWellOriented(isQuadratic,vec,conn+connI[i]+1,conn+connI[i+1],coordsPtr))
03966             {
03967               isModified=true;
03968               std::vector<int> tmp(connI[i+1]-connI[i]-2);
03969               std::copy(conn+connI[i]+2,conn+connI[i+1],tmp.rbegin());
03970               std::copy(tmp.begin(),tmp.end(),conn+connI[i]+2);
03971             }
03972         }
03973     }
03974   if(isModified)
03975     _nodal_connec->declareAsNew();
03976   updateTime();
03977 }
03978 
03984 void MEDCouplingUMesh::arePolyhedronsNotCorrectlyOriented(std::vector<int>& cells) const throw(INTERP_KERNEL::Exception)
03985 {
03986   if(getMeshDimension()!=3 || getSpaceDimension()!=3)
03987     throw INTERP_KERNEL::Exception("Invalid mesh to apply arePolyhedronsNotCorrectlyOriented on it : must be meshDim==3 and spaceDim==3 !");
03988   int nbOfCells=getNumberOfCells();
03989   const int *conn=_nodal_connec->getConstPointer();
03990   const int *connI=_nodal_connec_index->getConstPointer();
03991   const double *coordsPtr=_coords->getConstPointer();
03992   for(int i=0;i<nbOfCells;i++)
03993     {
03994       INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)conn[connI[i]];
03995       if(type==INTERP_KERNEL::NORM_POLYHED)
03996         {
03997           if(!IsPolyhedronWellOriented(conn+connI[i]+1,conn+connI[i+1],coordsPtr))
03998             cells.push_back(i);
03999         }
04000     }
04001 }
04002 
04007 void MEDCouplingUMesh::orientCorrectlyPolyhedrons() throw(INTERP_KERNEL::Exception)
04008 {
04009   if(getMeshDimension()!=3 || getSpaceDimension()!=3)
04010     throw INTERP_KERNEL::Exception("Invalid mesh to apply orientCorrectlyPolyhedrons on it : must be meshDim==3 and spaceDim==3 !");
04011   int nbOfCells=getNumberOfCells();
04012   int *conn=_nodal_connec->getPointer();
04013   const int *connI=_nodal_connec_index->getConstPointer();
04014   const double *coordsPtr=_coords->getConstPointer();
04015   bool isModified=false;
04016   for(int i=0;i<nbOfCells;i++)
04017     {
04018       INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)conn[connI[i]];
04019       if(type==INTERP_KERNEL::NORM_POLYHED)
04020         if(!IsPolyhedronWellOriented(conn+connI[i]+1,conn+connI[i+1],coordsPtr))
04021           {
04022             TryToCorrectPolyhedronOrientation(conn+connI[i]+1,conn+connI[i+1],coordsPtr);
04023             isModified=true;
04024           }
04025     }
04026   if(isModified)
04027     _nodal_connec->declareAsNew();
04028   updateTime();
04029 }
04030 
04038 void MEDCouplingUMesh::getFastAveragePlaneOfThis(double *vec, double *pos) const throw(INTERP_KERNEL::Exception)
04039 {
04040   if(getMeshDimension()!=2 || getSpaceDimension()!=3)
04041     throw INTERP_KERNEL::Exception("Invalid mesh to apply getFastAveragePlaneOfThis on it : must be meshDim==2 and spaceDim==3 !");
04042   const int *conn=_nodal_connec->getConstPointer();
04043   const int *connI=_nodal_connec_index->getConstPointer();
04044   const double *coordsPtr=_coords->getConstPointer();
04045   INTERP_KERNEL::areaVectorOfPolygon<int,INTERP_KERNEL::ALL_C_MODE>(conn+1,connI[1]-connI[0]-1,coordsPtr,vec);
04046   std::copy(coordsPtr+3*conn[1],coordsPtr+3*conn[1]+3,pos);
04047 }
04048 
04055 MEDCouplingFieldDouble *MEDCouplingUMesh::getEdgeRatioField() const throw(INTERP_KERNEL::Exception)
04056 {
04057   checkCoherency();
04058   int spaceDim=getSpaceDimension();
04059   int meshDim=getMeshDimension();
04060   if(spaceDim!=2 && spaceDim!=3)
04061     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getEdgeRatioField : SpaceDimension must be equal to 2 or 3 !");
04062   if(meshDim!=2 && meshDim!=3)
04063     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getEdgeRatioField : MeshDimension must be equal to 2 or 3 !");
04064   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
04065   ret->setMesh(this);
04066   int nbOfCells=getNumberOfCells();
04067   DataArrayDouble *arr=DataArrayDouble::New();
04068   arr->alloc(nbOfCells,1);
04069   double *pt=arr->getPointer();
04070   ret->setArray(arr);//In case of throw to avoid mem leaks arr will be used after decrRef.
04071   arr->decrRef();
04072   const int *conn=_nodal_connec->getConstPointer();
04073   const int *connI=_nodal_connec_index->getConstPointer();
04074   const double *coo=_coords->getConstPointer();
04075   double tmp[12];
04076   for(int i=0;i<nbOfCells;i++,pt++)
04077     {
04078       INTERP_KERNEL::NormalizedCellType t=(INTERP_KERNEL::NormalizedCellType)*conn;
04079       switch(t)
04080         {
04081           case INTERP_KERNEL::NORM_TRI3:
04082             {
04083               FillInCompact3DMode(spaceDim,3,conn+1,coo,tmp);
04084               *pt=INTERP_KERNEL::triEdgeRatio(tmp);
04085               break;
04086             }
04087           case INTERP_KERNEL::NORM_QUAD4:
04088             {
04089               FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
04090               *pt=INTERP_KERNEL::quadEdgeRatio(tmp);
04091               break;
04092             }
04093           case INTERP_KERNEL::NORM_TETRA4:
04094             {
04095               FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
04096               *pt=INTERP_KERNEL::tetraEdgeRatio(tmp);
04097               break;
04098             }
04099         default:
04100           throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getEdgeRatioField : A cell with not manged type (NORM_TRI3, NORM_QUAD4 and NORM_TETRA4) has been detected !");
04101         }
04102       conn+=connI[i+1]-connI[i];
04103     }
04104   ret->setName("EdgeRatio");
04105   ret->incrRef();
04106   return ret;
04107 }
04108 
04115 MEDCouplingFieldDouble *MEDCouplingUMesh::getAspectRatioField() const throw(INTERP_KERNEL::Exception)
04116 {
04117   checkCoherency();
04118   int spaceDim=getSpaceDimension();
04119   int meshDim=getMeshDimension();
04120   if(spaceDim!=2 && spaceDim!=3)
04121     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getAspectRatioField : SpaceDimension must be equal to 2 or 3 !");
04122   if(meshDim!=2 && meshDim!=3)
04123     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getAspectRatioField : MeshDimension must be equal to 2 or 3 !");
04124   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
04125   ret->setMesh(this);
04126   int nbOfCells=getNumberOfCells();
04127   DataArrayDouble *arr=DataArrayDouble::New();
04128   arr->alloc(nbOfCells,1);
04129   double *pt=arr->getPointer();
04130   ret->setArray(arr);//In case of throw to avoid mem leaks arr will be used after decrRef.
04131   arr->decrRef();
04132   const int *conn=_nodal_connec->getConstPointer();
04133   const int *connI=_nodal_connec_index->getConstPointer();
04134   const double *coo=_coords->getConstPointer();
04135   double tmp[12];
04136   for(int i=0;i<nbOfCells;i++,pt++)
04137     {
04138       INTERP_KERNEL::NormalizedCellType t=(INTERP_KERNEL::NormalizedCellType)*conn;
04139       switch(t)
04140         {
04141           case INTERP_KERNEL::NORM_TRI3:
04142             {
04143               FillInCompact3DMode(spaceDim,3,conn+1,coo,tmp);
04144               *pt=INTERP_KERNEL::triAspectRatio(tmp);
04145               break;
04146             }
04147           case INTERP_KERNEL::NORM_QUAD4:
04148             {
04149               FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
04150               *pt=INTERP_KERNEL::quadAspectRatio(tmp);
04151               break;
04152             }
04153           case INTERP_KERNEL::NORM_TETRA4:
04154             {
04155               FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
04156               *pt=INTERP_KERNEL::tetraAspectRatio(tmp);
04157               break;
04158             }
04159         default:
04160           throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getAspectRatioField : A cell with not manged type (NORM_TRI3, NORM_QUAD4 and NORM_TETRA4) has been detected !");
04161         }
04162       conn+=connI[i+1]-connI[i];
04163     }
04164   ret->setName("AspectRatio");
04165   ret->incrRef();
04166   return ret;
04167 }
04168 
04175 MEDCouplingFieldDouble *MEDCouplingUMesh::getWarpField() const throw(INTERP_KERNEL::Exception)
04176 {
04177   checkCoherency();
04178   int spaceDim=getSpaceDimension();
04179   int meshDim=getMeshDimension();
04180   if(spaceDim!=3)
04181     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getWarpField : SpaceDimension must be equal to 3 !");
04182   if(meshDim!=2)
04183     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getWarpField : MeshDimension must be equal to 2 !");
04184   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
04185   ret->setMesh(this);
04186   int nbOfCells=getNumberOfCells();
04187   DataArrayDouble *arr=DataArrayDouble::New();
04188   arr->alloc(nbOfCells,1);
04189   double *pt=arr->getPointer();
04190   ret->setArray(arr);//In case of throw to avoid mem leaks arr will be used after decrRef.
04191   arr->decrRef();
04192   const int *conn=_nodal_connec->getConstPointer();
04193   const int *connI=_nodal_connec_index->getConstPointer();
04194   const double *coo=_coords->getConstPointer();
04195   double tmp[12];
04196   for(int i=0;i<nbOfCells;i++,pt++)
04197     {
04198       INTERP_KERNEL::NormalizedCellType t=(INTERP_KERNEL::NormalizedCellType)*conn;
04199       switch(t)
04200         {
04201           case INTERP_KERNEL::NORM_QUAD4:
04202             {
04203               FillInCompact3DMode(3,4,conn+1,coo,tmp);
04204               *pt=INTERP_KERNEL::quadWarp(tmp);
04205               break;
04206             }
04207         default:
04208           throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getWarpField : A cell with not manged type (NORM_QUAD4) has been detected !");
04209         }
04210       conn+=connI[i+1]-connI[i];
04211     }
04212   ret->setName("Warp");
04213   ret->incrRef();
04214   return ret;
04215 }
04216 
04223 MEDCouplingFieldDouble *MEDCouplingUMesh::getSkewField() const throw(INTERP_KERNEL::Exception)
04224 {
04225   checkCoherency();
04226   int spaceDim=getSpaceDimension();
04227   int meshDim=getMeshDimension();
04228   if(spaceDim!=3)
04229     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getSkewField : SpaceDimension must be equal to 3 !");
04230   if(meshDim!=2)
04231     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getSkewField : MeshDimension must be equal to 2 !");
04232   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
04233   ret->setMesh(this);
04234   int nbOfCells=getNumberOfCells();
04235   DataArrayDouble *arr=DataArrayDouble::New();
04236   arr->alloc(nbOfCells,1);
04237   double *pt=arr->getPointer();
04238   ret->setArray(arr);//In case of throw to avoid mem leaks arr will be used after decrRef.
04239   arr->decrRef();
04240   const int *conn=_nodal_connec->getConstPointer();
04241   const int *connI=_nodal_connec_index->getConstPointer();
04242   const double *coo=_coords->getConstPointer();
04243   double tmp[12];
04244   for(int i=0;i<nbOfCells;i++,pt++)
04245     {
04246       INTERP_KERNEL::NormalizedCellType t=(INTERP_KERNEL::NormalizedCellType)*conn;
04247       switch(t)
04248         {
04249           case INTERP_KERNEL::NORM_QUAD4:
04250             {
04251               FillInCompact3DMode(3,4,conn+1,coo,tmp);
04252               *pt=INTERP_KERNEL::quadSkew(tmp);
04253               break;
04254             }
04255         default:
04256           throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getSkewField : A cell with not manged type (NORM_QUAD4) has been detected !");
04257         }
04258       conn+=connI[i+1]-connI[i];
04259     }
04260   ret->setName("Skew");
04261   ret->incrRef();
04262   return ret;
04263 }
04264 
04269 void MEDCouplingUMesh::getBoundingBoxForBBTree(std::vector<double>& bbox) const
04270 {
04271   int spaceDim=getSpaceDimension();
04272   int nbOfCells=getNumberOfCells();
04273   bbox.resize(2*nbOfCells*spaceDim);
04274   for(int i=0;i<nbOfCells*spaceDim;i++)
04275     {
04276       bbox[2*i]=std::numeric_limits<double>::max();
04277       bbox[2*i+1]=-std::numeric_limits<double>::max();
04278     }
04279   const double *coordsPtr=_coords->getConstPointer();
04280   const int *conn=_nodal_connec->getConstPointer();
04281   const int *connI=_nodal_connec_index->getConstPointer();
04282   for(int i=0;i<nbOfCells;i++)
04283     {
04284       int offset=connI[i]+1;
04285       int nbOfNodesForCell=connI[i+1]-offset;
04286       for(int j=0;j<nbOfNodesForCell;j++)
04287         {
04288           int nodeId=conn[offset+j];
04289           if(nodeId>=0)
04290             for(int k=0;k<spaceDim;k++)
04291               {
04292                 bbox[2*spaceDim*i+2*k]=std::min(bbox[2*spaceDim*i+2*k],coordsPtr[spaceDim*nodeId+k]);
04293                 bbox[2*spaceDim*i+2*k+1]=std::max(bbox[2*spaceDim*i+2*k+1],coordsPtr[spaceDim*nodeId+k]);
04294               }
04295         }
04296     }
04297 }
04298 
04300 
04301 namespace ParaMEDMEMImpl
04302 {
04303   class ConnReader
04304   {
04305   public:
04306     ConnReader(const int *c, int val):_conn(c),_val(val) { }
04307     bool operator() (const int& pos) { return _conn[pos]!=_val; }
04308   private:
04309     const int *_conn;
04310     int _val;
04311   };
04312 
04313   class ConnReader2
04314   {
04315   public:
04316     ConnReader2(const int *c, int val):_conn(c),_val(val) { }
04317     bool operator() (const int& pos) { return _conn[pos]==_val; }
04318   private:
04319     const int *_conn;
04320     int _val;
04321   };
04322 }
04323 
04325 
04334 std::vector<int> MEDCouplingUMesh::getDistributionOfTypes() const throw(INTERP_KERNEL::Exception)
04335 {
04336   checkConnectivityFullyDefined();
04337   const int *conn=_nodal_connec->getConstPointer();
04338   const int *connI=_nodal_connec_index->getConstPointer();
04339   const int *work=connI;
04340   int nbOfCells=getNumberOfCells();
04341   std::size_t n=getAllTypes().size();
04342   std::vector<int> ret(3*n,0); //ret[3*k+2]==0 because it has no sense here
04343   std::set<INTERP_KERNEL::NormalizedCellType> types;
04344   for(std::size_t i=0;work!=connI+nbOfCells;i++)
04345     {
04346       INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn[*work];
04347       if(types.find(typ)!=types.end())
04348         {
04349           std::ostringstream oss; oss << "MEDCouplingUMesh::getDistributionOfTypes : Type " << INTERP_KERNEL::CellModel::GetCellModel(typ).getRepr();
04350           oss << " is not contiguous !";
04351           throw INTERP_KERNEL::Exception(oss.str().c_str());
04352         }
04353       types.insert(typ);
04354       ret[3*i]=typ;
04355       const int *work2=std::find_if(work+1,connI+nbOfCells,ParaMEDMEMImpl::ConnReader(conn,typ));
04356       ret[3*i+1]=(int)std::distance(work,work2);
04357       work=work2;
04358     }
04359   return ret;
04360 }
04361 
04379 DataArrayInt *MEDCouplingUMesh::checkTypeConsistencyAndContig(const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
04380 {
04381   if(code.empty())
04382     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::checkTypeConsistencyAndContig : code is empty, should not !");
04383   std::size_t sz=code.size();
04384   std::size_t n=sz/3;
04385   if(sz%3!=0)
04386     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::checkTypeConsistencyAndContig : code size is NOT %3 !");
04387   std::vector<INTERP_KERNEL::NormalizedCellType> types;
04388   int nb=0;
04389   for(std::size_t i=0;i<n;i++)
04390     if(std::find(types.begin(),types.end(),(INTERP_KERNEL::NormalizedCellType)code[3*i])==types.end())
04391       {
04392         types.push_back((INTERP_KERNEL::NormalizedCellType)code[3*i]);
04393         nb+=code[3*i+1];
04394         if(_types.find((INTERP_KERNEL::NormalizedCellType)code[3*i])==_types.end())
04395           throw INTERP_KERNEL::Exception("MEDCouplingUMesh::checkTypeConsistencyAndContig : expected geo types not in this !");
04396       }
04397   if(types.size()!=n)
04398     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::checkTypeConsistencyAndContig : code contains duplication of types in unstructured mesh !");
04399   if(idsPerType.empty())
04400     {
04401       if(!checkConsecutiveCellTypesAndOrder(&types[0],&types[0]+types.size()))
04402         throw INTERP_KERNEL::Exception("MEDCouplingUMesh::checkTypeConsistencyAndContig : non contiguous type !");
04403       if(types.size()==_types.size())
04404         return 0;
04405     }
04406   DataArrayInt *ret=DataArrayInt::New();
04407   ret->alloc(nb,1);
04408   int *retPtr=ret->getPointer();
04409   const int *connI=_nodal_connec_index->getConstPointer();
04410   const int *conn=_nodal_connec->getConstPointer();
04411   int nbOfCells=getNumberOfCells();
04412   const int *i=connI;
04413   int kk=0;
04414   for(std::vector<INTERP_KERNEL::NormalizedCellType>::const_iterator it=types.begin();it!=types.end();it++,kk++)
04415     {
04416       i=std::find_if(i,connI+nbOfCells,ParaMEDMEMImpl::ConnReader2(conn,(int)(*it)));
04417       int offset=(int)std::distance(connI,i);
04418       if(code[3*kk+2]==-1)
04419         {
04420           const int *j=std::find_if(i+1,connI+nbOfCells,ParaMEDMEMImpl::ConnReader(conn,(int)(*it)));
04421           std::size_t pos2=std::distance(i,j);
04422           for(std::size_t k=0;k<pos2;k++)
04423             *retPtr++=(int)k+offset;
04424           i=j;
04425         }
04426       else
04427         {
04428           retPtr=std::transform(idsPerType[code[3*kk+2]]->getConstPointer(),idsPerType[code[3*kk+2]]->getConstPointer()+idsPerType[code[3*kk+2]]->getNbOfElems(),
04429                                 retPtr,std::bind2nd(std::plus<int>(),offset));
04430         }
04431     }
04432   return ret;
04433 }
04434 
04441 void MEDCouplingUMesh::splitProfilePerType(const DataArrayInt *profile, std::vector<int>& code, std::vector<DataArrayInt *>& idsInPflPerType, std::vector<DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
04442 {
04443   if(profile->getNumberOfComponents()!=1)
04444     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::splitProfilePerType : input profile should have exactly one component !");
04445   checkConnectivityFullyDefined();
04446   const int *conn=_nodal_connec->getConstPointer();
04447   const int *connI=_nodal_connec_index->getConstPointer();
04448   int nbOfCells=getNumberOfCells();
04449   std::vector<INTERP_KERNEL::NormalizedCellType> types;
04450   std::vector<int> typeRangeVals(1);
04451   for(const int *i=connI;i!=connI+nbOfCells;)
04452     {
04453       INTERP_KERNEL::NormalizedCellType curType=(INTERP_KERNEL::NormalizedCellType)conn[*i];
04454       if(std::find(types.begin(),types.end(),curType)!=types.end())
04455         {
04456           throw INTERP_KERNEL::Exception("MEDCouplingUMesh::splitProfilePerType : current mesh is not sorted by type !");
04457         }
04458       types.push_back(curType);
04459       i=std::find_if(i+1,connI+nbOfCells,ParaMEDMEMImpl::ConnReader(conn,(int)curType));
04460       typeRangeVals.push_back((int)std::distance(connI,i));
04461     }
04462   //
04463   DataArrayInt *castArr=0,*rankInsideCast=0,*castsPresent=0;
04464   profile->splitByValueRange(&typeRangeVals[0],&typeRangeVals[0]+typeRangeVals.size(),castArr,rankInsideCast,castsPresent);
04465   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp0=castArr;
04466   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp1=rankInsideCast;
04467   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp2=castsPresent;
04468   //
04469   int nbOfCastsFinal=castsPresent->getNumberOfTuples();
04470   code.resize(3*nbOfCastsFinal);
04471   std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayInt> > idsInPflPerType2;
04472   std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayInt> > idsPerType2;
04473   for(int i=0;i<nbOfCastsFinal;i++)
04474     {
04475       int castId=castsPresent->getIJ(i,0);
04476       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp3=castArr->getIdsEqual(castId);
04477       idsInPflPerType2.push_back(tmp3);
04478       code[3*i]=(int)types[castId];
04479       code[3*i+1]=tmp3->getNumberOfTuples();
04480       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp4=rankInsideCast->selectByTupleId(tmp3->getConstPointer(),tmp3->getConstPointer()+tmp3->getNumberOfTuples());
04481       if(tmp4->getNumberOfTuples()!=typeRangeVals[castId+1]-typeRangeVals[castId] || !tmp4->isIdentity())
04482         {
04483           tmp4->copyStringInfoFrom(*profile);
04484           idsPerType2.push_back(tmp4);
04485           code[3*i+2]=(int)idsPerType2.size()-1;
04486         }
04487       else
04488         {
04489           code[3*i+2]=-1;
04490         }
04491     }
04492   std::size_t sz2=idsInPflPerType2.size();
04493   idsInPflPerType.resize(sz2);
04494   for(std::size_t i=0;i<sz2;i++)
04495     {
04496       DataArrayInt *locDa=idsInPflPerType2[i];
04497       locDa->incrRef();
04498       idsInPflPerType[i]=locDa;
04499     }
04500   std::size_t sz=idsPerType2.size();
04501   idsPerType.resize(sz);
04502   for(std::size_t i=0;i<sz;i++)
04503     {
04504       DataArrayInt *locDa=idsPerType2[i];
04505       locDa->incrRef();
04506       idsPerType[i]=locDa;
04507     }
04508 }
04509 
04516 MEDCouplingUMesh *MEDCouplingUMesh::emulateMEDMEMBDC(const MEDCouplingUMesh *nM1LevMesh, DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *&revDesc, DataArrayInt *&revDescIndx, DataArrayInt *& nM1LevMeshIds, DataArrayInt *&meshnM1Old2New) const throw(INTERP_KERNEL::Exception)
04517 {
04518   checkFullyDefined();
04519   nM1LevMesh->checkFullyDefined();
04520   if(getMeshDimension()-1!=nM1LevMesh->getMeshDimension())
04521     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::emulateMEDMEMBDC : The mesh passed as first argument should have a meshDim equal to this->getMeshDimension()-1 !" );
04522   if(_coords!=nM1LevMesh->getCoords())
04523     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::emulateMEDMEMBDC : 'this' and mesh in first argument should share the same coords : Use tryToShareSameCoords method !");
04524   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp0=DataArrayInt::New();
04525   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp1=DataArrayInt::New();
04526   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret1=buildDescendingConnectivity(desc,descIndx,tmp0,tmp1);
04527   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret0=ret1->sortCellsInMEDFileFrmt();
04528   desc->transformWithIndArr(ret0->getConstPointer(),ret0->getConstPointer()+ret0->getNbOfElems());
04529   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> tmp=MEDCouplingUMesh::New();
04530   tmp->setConnectivity(tmp0,tmp1);
04531   tmp->renumberCells(ret0->getConstPointer(),false);
04532   revDesc=tmp->getNodalConnectivity();
04533   revDescIndx=tmp->getNodalConnectivityIndex();
04534   DataArrayInt *ret=0;
04535   if(!ret1->areCellsIncludedIn(nM1LevMesh,2,ret))
04536     {
04537       int tmp2;
04538       ret->getMaxValue(tmp2);
04539       ret->decrRef();
04540       std::ostringstream oss; oss << "MEDCouplingUMesh::emulateMEDMEMBDC : input N-1 mesh present a cell not in descending mesh ... Id of cell is " << tmp2 << " !";
04541       throw INTERP_KERNEL::Exception(oss.str().c_str());
04542     }
04543   nM1LevMeshIds=ret;
04544   //
04545   revDesc->incrRef();
04546   revDescIndx->incrRef();
04547   ret1->incrRef();
04548   ret0->incrRef();
04549   meshnM1Old2New=ret0;
04550   return ret1;
04551 }
04552 
04559 DataArrayInt *MEDCouplingUMesh::sortCellsInMEDFileFrmt() throw(INTERP_KERNEL::Exception)
04560 {
04561   static const int N=19;
04562   static const INTERP_KERNEL::NormalizedCellType MEDMEM_ORDER[N] = { INTERP_KERNEL::NORM_POINT1, INTERP_KERNEL::NORM_SEG2, INTERP_KERNEL::NORM_SEG3, INTERP_KERNEL::NORM_TRI3, INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6, INTERP_KERNEL::NORM_QUAD8, INTERP_KERNEL::NORM_TETRA4, INTERP_KERNEL::NORM_PYRA5, INTERP_KERNEL::NORM_PENTA6, INTERP_KERNEL::NORM_HEXA8, INTERP_KERNEL::NORM_HEXGP12, INTERP_KERNEL::NORM_TETRA10, INTERP_KERNEL::NORM_PYRA13, INTERP_KERNEL::NORM_PENTA15, INTERP_KERNEL::NORM_HEXA20, INTERP_KERNEL::NORM_POLYGON, INTERP_KERNEL::NORM_QPOLYG, INTERP_KERNEL::NORM_POLYHED };
04563   checkConnectivityFullyDefined();
04564   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=getRenumArrForConsecutiveCellTypesSpec(MEDMEM_ORDER,MEDMEM_ORDER+N);
04565   renumberCells(ret->getConstPointer(),false);
04566   ret->incrRef();
04567   return ret;
04568 }
04569 
04574 bool MEDCouplingUMesh::checkConsecutiveCellTypes() const
04575 {
04576   checkFullyDefined();
04577   const int *conn=_nodal_connec->getConstPointer();
04578   const int *connI=_nodal_connec_index->getConstPointer();
04579   int nbOfCells=getNumberOfCells();
04580   std::set<INTERP_KERNEL::NormalizedCellType> types;
04581   for(const int *i=connI;i!=connI+nbOfCells;)
04582     {
04583       INTERP_KERNEL::NormalizedCellType curType=(INTERP_KERNEL::NormalizedCellType)conn[*i];
04584       if(types.find(curType)!=types.end())
04585         return false;
04586       types.insert(curType);
04587       i=std::find_if(i+1,connI+nbOfCells,ParaMEDMEMImpl::ConnReader(conn,(int)curType));
04588     }
04589   return true;
04590 }
04591 
04596 bool MEDCouplingUMesh::checkConsecutiveCellTypesAndOrder(const INTERP_KERNEL::NormalizedCellType *orderBg, const INTERP_KERNEL::NormalizedCellType *orderEnd) const
04597 {
04598   checkFullyDefined();
04599   const int *conn=_nodal_connec->getConstPointer();
04600   const int *connI=_nodal_connec_index->getConstPointer();
04601   int nbOfCells=getNumberOfCells();
04602   int lastPos=-1;
04603   for(const int *i=connI;i!=connI+nbOfCells;)
04604     {
04605       INTERP_KERNEL::NormalizedCellType curType=(INTERP_KERNEL::NormalizedCellType)conn[*i];
04606       int pos=(int)std::distance(orderBg,std::find(orderBg,orderEnd,curType));
04607       if(pos<=lastPos)
04608         return false;
04609       lastPos=pos;
04610       i=std::find_if(i+1,connI+nbOfCells,ParaMEDMEMImpl::ConnReader(conn,(int)curType));
04611     }
04612   return true;
04613 }
04614 
04620 DataArrayInt *MEDCouplingUMesh::getLevArrPerCellTypes(const INTERP_KERNEL::NormalizedCellType *orderBg, const INTERP_KERNEL::NormalizedCellType *orderEnd, DataArrayInt *&nbPerType) const throw(INTERP_KERNEL::Exception)
04621 {
04622   checkConnectivityFullyDefined();
04623   int nbOfCells=getNumberOfCells();
04624   const int *conn=_nodal_connec->getConstPointer();
04625   const int *connI=_nodal_connec_index->getConstPointer();
04626   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmpa=DataArrayInt::New();
04627   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmpb=DataArrayInt::New();
04628   tmpa->alloc(nbOfCells,1);
04629   tmpb->alloc((int)std::distance(orderBg,orderEnd),1);
04630   tmpb->fillWithZero();
04631   int *tmp=tmpa->getPointer();
04632   int *tmp2=tmpb->getPointer();
04633   for(const int *i=connI;i!=connI+nbOfCells;i++)
04634     {
04635       const INTERP_KERNEL::NormalizedCellType *where=std::find(orderBg,orderEnd,(INTERP_KERNEL::NormalizedCellType)conn[*i]);
04636       if(where!=orderEnd)
04637         {
04638           int pos=(int)std::distance(orderBg,where);
04639           tmp2[pos]++;
04640           tmp[std::distance(connI,i)]=pos;
04641         }
04642       else
04643         {
04644           const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[*i]);
04645           std::ostringstream oss; oss << "MEDCouplingUMesh::getLevArrPerCellTypes : Cell #" << std::distance(connI,i);
04646           oss << " has a type " << cm.getRepr() << " not in input array of type !";
04647           throw INTERP_KERNEL::Exception(oss.str().c_str());
04648         }
04649     }
04650   nbPerType=tmpb;
04651   tmpa->incrRef();
04652   tmpb->incrRef();
04653   return tmpa;
04654 }
04655 
04662 DataArrayInt *MEDCouplingUMesh::getRenumArrForConsecutiveCellTypesSpec(const INTERP_KERNEL::NormalizedCellType *orderBg, const INTERP_KERNEL::NormalizedCellType *orderEnd) const throw(INTERP_KERNEL::Exception)
04663 {
04664   DataArrayInt *nbPerType=0;
04665   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmpa=getLevArrPerCellTypes(orderBg,orderEnd,nbPerType);
04666   nbPerType->decrRef();
04667   return tmpa->buildPermArrPerLevel();
04668 }
04669 
04678 DataArrayInt *MEDCouplingUMesh::rearrange2ConsecutiveCellTypes()
04679 {
04680   checkFullyDefined();
04681   computeTypes();
04682   const int *conn=_nodal_connec->getConstPointer();
04683   const int *connI=_nodal_connec_index->getConstPointer();
04684   int nbOfCells=getNumberOfCells();
04685   std::vector<INTERP_KERNEL::NormalizedCellType> types;
04686   for(const int *i=connI;i!=connI+nbOfCells && (types.size()!=_types.size());)
04687     if(std::find(types.begin(),types.end(),(INTERP_KERNEL::NormalizedCellType)conn[*i])==types.end())
04688       {
04689         INTERP_KERNEL::NormalizedCellType curType=(INTERP_KERNEL::NormalizedCellType)conn[*i];
04690         types.push_back(curType);
04691         for(i++;i!=connI+nbOfCells && (INTERP_KERNEL::NormalizedCellType)conn[*i]==curType;i++);
04692       }
04693   DataArrayInt *ret=DataArrayInt::New();
04694   ret->alloc(nbOfCells,1);
04695   int *retPtr=ret->getPointer();
04696   std::fill(retPtr,retPtr+nbOfCells,-1);
04697   int newCellId=0;
04698   for(std::vector<INTERP_KERNEL::NormalizedCellType>::const_iterator iter=types.begin();iter!=types.end();iter++)
04699     {
04700       for(const int *i=connI;i!=connI+nbOfCells;i++)
04701         if((INTERP_KERNEL::NormalizedCellType)conn[*i]==(*iter))
04702           retPtr[std::distance(connI,i)]=newCellId++;
04703     }
04704   renumberCells(retPtr,false);
04705   return ret;
04706 }
04707 
04713 std::vector<MEDCouplingUMesh *> MEDCouplingUMesh::splitByType() const
04714 {
04715   checkFullyDefined();
04716   const int *conn=_nodal_connec->getConstPointer();
04717   const int *connI=_nodal_connec_index->getConstPointer();
04718   int nbOfCells=getNumberOfCells();
04719   std::vector<MEDCouplingUMesh *> ret;
04720   for(const int *i=connI;i!=connI+nbOfCells;)
04721     {
04722       INTERP_KERNEL::NormalizedCellType curType=(INTERP_KERNEL::NormalizedCellType)conn[*i];
04723       int beginCellId=(int)std::distance(connI,i);
04724       i=std::find_if(i+1,connI+nbOfCells,ParaMEDMEMImpl::ConnReader(conn,(int)curType));
04725       int endCellId=(int)std::distance(connI,i);
04726       int sz=endCellId-beginCellId;
04727       int *cells=new int[sz];
04728       for(int j=0;j<sz;j++)
04729         cells[j]=beginCellId+j;
04730       MEDCouplingUMesh *m=(MEDCouplingUMesh *)buildPartOfMySelf(cells,cells+sz,true);
04731       delete [] cells;
04732       ret.push_back(m);
04733     }
04734   return ret;
04735 }
04736 
04741 DataArrayInt *MEDCouplingUMesh::keepCellIdsByType(INTERP_KERNEL::NormalizedCellType type, const int *begin, const int *end) const throw(INTERP_KERNEL::Exception)
04742 {
04743   checkFullyDefined();
04744   std::vector<int> r;
04745   const int *conn=_nodal_connec->getConstPointer();
04746   const int *connIndex=_nodal_connec_index->getConstPointer();
04747   for(const int *w=begin;w!=end;w++)
04748     if((INTERP_KERNEL::NormalizedCellType)conn[connIndex[*w]]==type)
04749       r.push_back(*w);
04750   DataArrayInt *ret=DataArrayInt::New();
04751   ret->alloc((int)r.size(),1);
04752   std::copy(r.begin(),r.end(),ret->getPointer());
04753   return ret;
04754 }
04755 
04760 DataArrayInt *MEDCouplingUMesh::convertCellArrayPerGeoType(const DataArrayInt *da) const throw(INTERP_KERNEL::Exception)
04761 {
04762   checkFullyDefined();
04763   const int *conn=_nodal_connec->getConstPointer();
04764   const int *connI=_nodal_connec_index->getConstPointer();
04765   int nbOfCells=getNumberOfCells();
04766   std::set<INTERP_KERNEL::NormalizedCellType> types=getAllTypes();
04767   int *tmp=new int[nbOfCells];
04768   for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator iter=types.begin();iter!=types.end();iter++)
04769     {
04770       int j=0;
04771       for(const int *i=connI;i!=connI+nbOfCells;i++)
04772         if((INTERP_KERNEL::NormalizedCellType)conn[*i]==(*iter))
04773           tmp[std::distance(connI,i)]=j++;
04774     }
04775   DataArrayInt *ret=DataArrayInt::New();
04776   ret->alloc(da->getNumberOfTuples(),da->getNumberOfComponents());
04777   ret->copyStringInfoFrom(*da);
04778   int *retPtr=ret->getPointer();
04779   const int *daPtr=da->getConstPointer();
04780   int nbOfElems=da->getNbOfElems();
04781   for(int k=0;k<nbOfElems;k++)
04782     retPtr[k]=tmp[daPtr[k]];
04783   delete [] tmp;
04784   return ret;
04785 }
04786 
04792 MEDCouplingUMesh *MEDCouplingUMesh::keepSpecifiedCells(INTERP_KERNEL::NormalizedCellType type, const int *idsPerGeoTypeBg, const int *idsPerGeoTypeEnd) const
04793 {
04794   std::vector<int> idsTokeep;
04795   int nbOfCells=getNumberOfCells();
04796   int j=0;
04797   for(int i=0;i<nbOfCells;i++)
04798     if(getTypeOfCell(i)!=type)
04799       idsTokeep.push_back(i);
04800     else
04801       {
04802         if(std::find(idsPerGeoTypeBg,idsPerGeoTypeEnd,j)!=idsPerGeoTypeEnd)
04803           idsTokeep.push_back(i);
04804         j++;
04805       }
04806   MEDCouplingPointSet *ret=buildPartOfMySelf(&idsTokeep[0],&idsTokeep[0]+idsTokeep.size(),true);
04807   MEDCouplingUMesh *ret2=dynamic_cast<MEDCouplingUMesh *>(ret);
04808   if(!ret2)
04809     {
04810       ret->decrRef();
04811       return 0;
04812     }
04813   ret2->copyTinyInfoFrom(this);
04814   return ret2;
04815 }
04816 
04821 std::vector<bool> MEDCouplingUMesh::getQuadraticStatus() const throw(INTERP_KERNEL::Exception)
04822 {
04823   int ncell=getNumberOfCells();
04824   std::vector<bool> ret(ncell);
04825   const int *cI=getNodalConnectivityIndex()->getConstPointer();
04826   const int *c=getNodalConnectivity()->getConstPointer();
04827   for(int i=0;i<ncell;i++)
04828     {
04829       INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)c[cI[i]];
04830       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
04831       ret[i]=cm.isQuadratic();
04832     }
04833   return ret;
04834 }
04835 
04839 MEDCouplingMesh *MEDCouplingUMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
04840 {
04841   if(other->getType()!=UNSTRUCTURED)
04842     throw INTERP_KERNEL::Exception("Merge of umesh only available with umesh each other !");
04843   const MEDCouplingUMesh *otherC=static_cast<const MEDCouplingUMesh *>(other);
04844   return MergeUMeshes(this,otherC);
04845 }
04846 
04852 DataArrayDouble *MEDCouplingUMesh::getBarycenterAndOwner() const
04853 {
04854   DataArrayDouble *ret=DataArrayDouble::New();
04855   int spaceDim=getSpaceDimension();
04856   int nbOfCells=getNumberOfCells();
04857   ret->alloc(nbOfCells,spaceDim);
04858   ret->copyStringInfoFrom(*getCoords());
04859   double *ptToFill=ret->getPointer();
04860   double *tmp=new double[spaceDim];
04861   const int *nodal=_nodal_connec->getConstPointer();
04862   const int *nodalI=_nodal_connec_index->getConstPointer();
04863   const double *coor=_coords->getConstPointer();
04864   for(int i=0;i<nbOfCells;i++)
04865     {
04866       INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)nodal[nodalI[i]];
04867       INTERP_KERNEL::computeBarycenter2<int,INTERP_KERNEL::ALL_C_MODE>(type,nodal+nodalI[i]+1,nodalI[i+1]-nodalI[i]-1,coor,spaceDim,ptToFill);
04868       ptToFill+=spaceDim;
04869     }
04870   delete [] tmp;
04871   return ret;
04872 }
04873 
04879 DataArrayDouble *MEDCouplingUMesh::getPartBarycenterAndOwner(const int *begin, const int *end) const
04880 {
04881   DataArrayDouble *ret=DataArrayDouble::New();
04882   int spaceDim=getSpaceDimension();
04883   int nbOfTuple=(int)std::distance(begin,end);
04884   ret->alloc(nbOfTuple,spaceDim);
04885   double *ptToFill=ret->getPointer();
04886   double *tmp=new double[spaceDim];
04887   const int *nodal=_nodal_connec->getConstPointer();
04888   const int *nodalI=_nodal_connec_index->getConstPointer();
04889   const double *coor=_coords->getConstPointer();
04890   for(const int *w=begin;w!=end;w++)
04891     {
04892       INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)nodal[nodalI[*w]];
04893       INTERP_KERNEL::computeBarycenter2<int,INTERP_KERNEL::ALL_C_MODE>(type,nodal+nodalI[*w]+1,nodalI[*w+1]-nodalI[*w]-1,coor,spaceDim,ptToFill);
04894       ptToFill+=spaceDim;
04895     }
04896   delete [] tmp;
04897   return ret;
04898 }
04899 
04904 MEDCouplingUMesh *MEDCouplingUMesh::Build0DMeshFromCoords(DataArrayDouble *da) throw(INTERP_KERNEL::Exception)
04905 {
04906   if(!da)
04907     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Build0DMeshFromCoords : instance of DataArrayDouble must be not null !");
04908   da->checkAllocated();
04909   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret=MEDCouplingUMesh::New(da->getName().c_str(),0);
04910   ret->setCoords(da);
04911   int nbOfTuples=da->getNumberOfTuples();
04912   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c=DataArrayInt::New();
04913   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cI=DataArrayInt::New();
04914   c->alloc(2*nbOfTuples,1);
04915   cI->alloc(nbOfTuples+1,1);
04916   int *cp=c->getPointer();
04917   int *cip=cI->getPointer();
04918   *cip++=0;
04919   for(int i=0;i<nbOfTuples;i++)
04920     {
04921       *cp++=INTERP_KERNEL::NORM_POINT1;
04922       *cp++=i;
04923       *cip++=2*(i+1);
04924     }
04925   ret->setConnectivity(c,cI,true);
04926   ret->incrRef();
04927   return ret;
04928 }
04929 
04934 MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshes(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception)
04935 {
04936   std::vector<const MEDCouplingUMesh *> tmp(2);
04937   tmp[0]=const_cast<MEDCouplingUMesh *>(mesh1); tmp[1]=const_cast<MEDCouplingUMesh *>(mesh2);
04938   return MergeUMeshes(tmp);
04939 }
04940 
04947 MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshes(std::vector<const MEDCouplingUMesh *>& a) throw(INTERP_KERNEL::Exception)
04948 {
04949   std::size_t sz=a.size();
04950   if(sz==0)
04951     return MergeUMeshesLL(a);
04952   for(std::size_t ii=0;ii<sz;ii++)
04953     if(!a[ii])
04954       {
04955         std::ostringstream oss; oss << "MEDCouplingUMesh::MergeUMeshes : item #" << ii << " in input array of size "<< sz << " is empty !";
04956         throw INTERP_KERNEL::Exception(oss.str().c_str());
04957       }
04958   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > bb(sz);
04959   std::vector< const MEDCouplingUMesh * > aa(sz);
04960   int spaceDim=-3;
04961   for(std::size_t i=0;i<sz && spaceDim==-3;i++)
04962     {
04963       const MEDCouplingUMesh *cur=a[i];
04964       const DataArrayDouble *coo=cur->getCoords();
04965       if(coo)
04966         spaceDim=coo->getNumberOfComponents();
04967     }
04968   if(spaceDim==-3)
04969     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::MergeUMeshes : no spaceDim specified ! unable to perform merge !");
04970   for(std::size_t i=0;i<sz;i++)
04971     {
04972       bb[i]=a[i]->buildSetInstanceFromThis(spaceDim);
04973       aa[i]=bb[i];
04974     }
04975   return MergeUMeshesLL(aa);
04976 }
04977 
04979 
04980 MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshesLL(std::vector<const MEDCouplingUMesh *>& a) throw(INTERP_KERNEL::Exception)
04981 {
04982   if(a.empty())
04983     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::MergeUMeshes : input array must be NON EMPTY !");
04984   std::vector<const MEDCouplingUMesh *>::const_iterator it=a.begin();
04985   int meshDim=(*it)->getMeshDimension();
04986   int nbOfCells=(*it)->getNumberOfCells();
04987   int meshLgth=(*it++)->getMeshLength();
04988   for(;it!=a.end();it++)
04989     {
04990       if(meshDim!=(*it)->getMeshDimension())
04991         throw INTERP_KERNEL::Exception("Mesh dimensions mismatches, MergeUMeshes impossible !");
04992       nbOfCells+=(*it)->getNumberOfCells();
04993       meshLgth+=(*it)->getMeshLength();
04994     }
04995   std::vector<const MEDCouplingPointSet *> aps(a.size());
04996   std::copy(a.begin(),a.end(),aps.begin());
04997   DataArrayDouble *pts=MergeNodesArray(aps);
04998   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret=MEDCouplingUMesh::New("merge",meshDim);
04999   ret->setCoords(pts);
05000   pts->decrRef();
05001   DataArrayInt *c=DataArrayInt::New();
05002   c->alloc(meshLgth,1);
05003   int *cPtr=c->getPointer();
05004   DataArrayInt *cI=DataArrayInt::New();
05005   cI->alloc(nbOfCells+1,1);
05006   int *cIPtr=cI->getPointer();
05007   *cIPtr++=0;
05008   int offset=0;
05009   int offset2=0;
05010   for(it=a.begin();it!=a.end();it++)
05011     {
05012       int curNbOfCell=(*it)->getNumberOfCells();
05013       const int *curCI=(*it)->_nodal_connec_index->getConstPointer();
05014       const int *curC=(*it)->_nodal_connec->getConstPointer();
05015       cIPtr=std::transform(curCI+1,curCI+curNbOfCell+1,cIPtr,std::bind2nd(std::plus<int>(),offset));
05016       for(int j=0;j<curNbOfCell;j++)
05017         {
05018           const int *src=curC+curCI[j];
05019           *cPtr++=*src++;
05020           for(;src!=curC+curCI[j+1];src++,cPtr++)
05021             {
05022               if(*src!=-1)
05023                 *cPtr=*src+offset2;
05024               else
05025                 *cPtr=-1;
05026             }
05027         }
05028       offset+=curCI[curNbOfCell];
05029       offset2+=(*it)->getNumberOfNodes();
05030     }
05031   //
05032   ret->setConnectivity(c,cI,true);
05033   c->decrRef();
05034   cI->decrRef();
05035   ret->incrRef();
05036   return ret;
05037 }
05038 
05040 
05045 MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshesOnSameCoords(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception)
05046 {
05047   std::vector<const MEDCouplingUMesh *> tmp(2);
05048   tmp[0]=mesh1; tmp[1]=mesh2;
05049   return MergeUMeshesOnSameCoords(tmp);
05050 }
05051 
05056 MEDCouplingUMesh *MEDCouplingUMesh::MergeUMeshesOnSameCoords(const std::vector<const MEDCouplingUMesh *>& meshes)
05057 {
05058   if(meshes.empty())
05059     throw INTERP_KERNEL::Exception("meshes input parameter is expected to be non empty.");
05060   for(std::size_t ii=0;ii<meshes.size();ii++)
05061     if(!meshes[ii])
05062       {
05063         std::ostringstream oss; oss << "MEDCouplingUMesh::MergeUMeshesOnSameCoords : item #" << ii << " in input array of size "<< meshes.size() << " is empty !";;
05064         throw INTERP_KERNEL::Exception(oss.str().c_str());
05065       }
05066   const DataArrayDouble *coords=meshes.front()->getCoords();
05067   int meshDim=meshes.front()->getMeshDimension();
05068   std::vector<const MEDCouplingUMesh *>::const_iterator iter=meshes.begin();
05069   int meshLgth=0;
05070   int meshIndexLgth=0;
05071   for(;iter!=meshes.end();iter++)
05072     {
05073       if(coords!=(*iter)->getCoords())
05074         throw INTERP_KERNEL::Exception("meshes does not share the same coords ! Try using tryToShareSameCoords method !");
05075       if(meshDim!=(*iter)->getMeshDimension())
05076         throw INTERP_KERNEL::Exception("Mesh dimensions mismatches, FuseUMeshesOnSameCoords impossible !");
05077       meshLgth+=(*iter)->getMeshLength();
05078       meshIndexLgth+=(*iter)->getNumberOfCells();
05079     }
05080   DataArrayInt *nodal=DataArrayInt::New();
05081   nodal->alloc(meshLgth,1);
05082   int *nodalPtr=nodal->getPointer();
05083   DataArrayInt *nodalIndex=DataArrayInt::New();
05084   nodalIndex->alloc(meshIndexLgth+1,1);
05085   int *nodalIndexPtr=nodalIndex->getPointer();
05086   int offset=0;
05087   for(iter=meshes.begin();iter!=meshes.end();iter++)
05088     {
05089       const int *nod=(*iter)->getNodalConnectivity()->getConstPointer();
05090       const int *index=(*iter)->getNodalConnectivityIndex()->getConstPointer();
05091       int nbOfCells=(*iter)->getNumberOfCells();
05092       int meshLgth2=(*iter)->getMeshLength();
05093       nodalPtr=std::copy(nod,nod+meshLgth2,nodalPtr);
05094       if(iter!=meshes.begin())
05095         nodalIndexPtr=std::transform(index+1,index+nbOfCells+1,nodalIndexPtr,std::bind2nd(std::plus<int>(),offset));
05096       else
05097         nodalIndexPtr=std::copy(index,index+nbOfCells+1,nodalIndexPtr);
05098       offset+=meshLgth2;
05099     }
05100   MEDCouplingUMesh *ret=MEDCouplingUMesh::New();
05101   ret->setName("merge");
05102   ret->setMeshDimension(meshDim);
05103   ret->setConnectivity(nodal,nodalIndex,true);
05104   ret->setCoords(coords);
05105   nodalIndex->decrRef();
05106   nodal->decrRef();
05107   return ret;
05108 }
05109 
05122 MEDCouplingUMesh *MEDCouplingUMesh::FuseUMeshesOnSameCoords(const std::vector<const MEDCouplingUMesh *>& meshes, int compType, std::vector<DataArrayInt *>& corr)
05123 {
05124   //All checks are delegated to MergeUMeshesOnSameCoords
05125   MEDCouplingUMesh *ret=MergeUMeshesOnSameCoords(meshes);
05126   DataArrayInt *o2n=ret->zipConnectivityTraducer(compType);
05127   corr.resize(meshes.size());
05128   std::size_t nbOfMeshes=meshes.size();
05129   int offset=0;
05130   const int *o2nPtr=o2n->getConstPointer();
05131   for(std::size_t i=0;i<nbOfMeshes;i++)
05132     {
05133       DataArrayInt *tmp=DataArrayInt::New();
05134       int curNbOfCells=meshes[i]->getNumberOfCells();
05135       tmp->alloc(curNbOfCells,1);
05136       std::copy(o2nPtr+offset,o2nPtr+offset+curNbOfCells,tmp->getPointer());
05137       offset+=curNbOfCells;
05138       tmp->setName(meshes[i]->getName());
05139       corr[i]=tmp;
05140     }
05141   o2n->decrRef();
05142   return ret;
05143 }
05144 
05157 void MEDCouplingUMesh::PutUMeshesOnSameAggregatedCoords(const std::vector<MEDCouplingUMesh *>& meshes) throw(INTERP_KERNEL::Exception)
05158 {
05159   std::size_t sz=meshes.size();
05160   if(sz==0 || sz==1)
05161     return;
05162   std::vector< const DataArrayDouble * > coords(meshes.size());
05163   std::vector< const DataArrayDouble * >::iterator it2=coords.begin();
05164   for(std::vector<MEDCouplingUMesh *>::const_iterator it=meshes.begin();it!=meshes.end();it++,it2++)
05165     {
05166       if((*it))
05167         {
05168           (*it)->checkConnectivityFullyDefined();
05169           const DataArrayDouble *coo=(*it)->getCoords();
05170           if(coo)
05171             *it2=coo;
05172           else
05173             {
05174               std::ostringstream oss; oss << " MEDCouplingUMesh::PutUMeshesOnSameAggregatedCoords : Item #" << std::distance(meshes.begin(),it) << " inside the vector of length " << meshes.size();
05175               oss << " has no coordinate array defined !";
05176               throw INTERP_KERNEL::Exception(oss.str().c_str());
05177             }
05178         }
05179       else
05180         {
05181           std::ostringstream oss; oss << " MEDCouplingUMesh::PutUMeshesOnSameAggregatedCoords : Item #" << std::distance(meshes.begin(),it) << " inside the vector of length " << meshes.size();
05182           oss << " is null !";
05183           throw INTERP_KERNEL::Exception(oss.str().c_str());
05184         }
05185     }
05186   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> res=DataArrayDouble::Aggregate(coords);
05187   std::vector<MEDCouplingUMesh *>::const_iterator it=meshes.begin();
05188   int offset=(*it)->getNumberOfNodes();
05189   (*it++)->setCoords(res);
05190   for(;it!=meshes.end();it++)
05191     {
05192       int oldNumberOfNodes=(*it)->getNumberOfNodes();
05193       (*it)->setCoords(res);
05194       (*it)->shiftNodeNumbersInConn(offset);
05195       offset+=oldNumberOfNodes;
05196     }
05197 }
05198 
05213 void MEDCouplingUMesh::MergeNodesOnUMeshesSharingSameCoords(const std::vector<MEDCouplingUMesh *>& meshes, double eps) throw(INTERP_KERNEL::Exception)
05214 {
05215   if(meshes.empty())
05216     return ;
05217   std::set<const DataArrayDouble *> s;
05218   for(std::vector<MEDCouplingUMesh *>::const_iterator it=meshes.begin();it!=meshes.end();it++)
05219     {
05220       if(*it)
05221         s.insert((*it)->getCoords());
05222       else
05223         {
05224           std::ostringstream oss; oss << "MEDCouplingUMesh::MergeNodesOnUMeshesSharingSameCoords : In input vector of unstructured meshes of size " << meshes.size() << " the element #" << std::distance(meshes.begin(),it) << " is null !";
05225           throw INTERP_KERNEL::Exception(oss.str().c_str());
05226         }
05227     }
05228   if(s.size()!=1)
05229     {
05230       std::ostringstream oss; oss << "MEDCouplingUMesh::MergeNodesOnUMeshesSharingSameCoords : In input vector of unstructured meshes of size " << meshes.size() << ", it appears that they do not share the same instance of DataArrayDouble for coordiantes ! tryToShareSameCoordsPermute method can help to reach that !";
05231       throw INTERP_KERNEL::Exception(oss.str().c_str());
05232     }
05233   const DataArrayDouble *coo=*(s.begin());
05234   if(!coo)
05235     return;
05236   //
05237   DataArrayInt *comm,*commI;
05238   coo->findCommonTuples(eps,-1,comm,commI);
05239   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp1(comm),tmp2(commI);
05240   int oldNbOfNodes=coo->getNumberOfTuples();
05241   int newNbOfNodes;
05242   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> o2n=DataArrayInt::BuildOld2NewArrayFromSurjectiveFormat2(oldNbOfNodes,comm,commI,newNbOfNodes);
05243   if(oldNbOfNodes==newNbOfNodes)
05244     return ;
05245   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> newCoords=coo->renumberAndReduce(o2n->getConstPointer(),newNbOfNodes);
05246   for(std::vector<MEDCouplingUMesh *>::const_iterator it=meshes.begin();it!=meshes.end();it++)
05247     {
05248       (*it)->renumberNodesInConn(o2n->getConstPointer());
05249       (*it)->setCoords(newCoords);
05250     } 
05251 }
05252 
05259 void MEDCouplingUMesh::AppendExtrudedCell(const int *connBg, const int *connEnd, int nbOfNodesPerLev, bool isQuad, std::vector<int>& ret)
05260 {
05261   INTERP_KERNEL::NormalizedCellType flatType=(INTERP_KERNEL::NormalizedCellType)connBg[0];
05262   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(flatType);
05263   ret.push_back(cm.getExtrudedType());
05264   int deltaz=isQuad?2*nbOfNodesPerLev:nbOfNodesPerLev;
05265   switch(flatType)
05266     {
05267     case INTERP_KERNEL::NORM_POINT1:
05268       {
05269         ret.push_back(connBg[1]);
05270         ret.push_back(connBg[1]+nbOfNodesPerLev);
05271         break;
05272       }
05273     case INTERP_KERNEL::NORM_SEG2:
05274       {
05275         int conn[4]={connBg[1],connBg[2],connBg[2]+deltaz,connBg[1]+deltaz};
05276         ret.insert(ret.end(),conn,conn+4);
05277         break;
05278       }
05279     case INTERP_KERNEL::NORM_SEG3:
05280       {
05281         int conn[8]={connBg[1],connBg[3],connBg[3]+deltaz,connBg[1]+deltaz,connBg[2],connBg[3]+nbOfNodesPerLev,connBg[2]+deltaz,connBg[1]+nbOfNodesPerLev};
05282         ret.insert(ret.end(),conn,conn+8);
05283         break;
05284       }
05285     case INTERP_KERNEL::NORM_QUAD4:
05286       {
05287         int conn[8]={connBg[1],connBg[2],connBg[3],connBg[4],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4]+deltaz};
05288         ret.insert(ret.end(),conn,conn+8);
05289         break;
05290       }
05291     case INTERP_KERNEL::NORM_TRI3:
05292       {
05293         int conn[6]={connBg[1],connBg[2],connBg[3],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz};
05294         ret.insert(ret.end(),conn,conn+6);
05295         break;
05296       }
05297     case INTERP_KERNEL::NORM_TRI6:
05298       {
05299         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,
05300                       connBg[1]+nbOfNodesPerLev,connBg[2]+nbOfNodesPerLev,connBg[3]+nbOfNodesPerLev};
05301         ret.insert(ret.end(),conn,conn+15);
05302         break;
05303       }
05304     case INTERP_KERNEL::NORM_QUAD8:
05305       {
05306         int conn[20]={
05307           connBg[1],connBg[2],connBg[3],connBg[4],connBg[1]+deltaz,connBg[2]+deltaz,connBg[3]+deltaz,connBg[4]+deltaz,
05308           connBg[5],connBg[6],connBg[7],connBg[8],connBg[5]+deltaz,connBg[6]+deltaz,connBg[7]+deltaz,connBg[8]+deltaz,
05309           connBg[1]+nbOfNodesPerLev,connBg[2]+nbOfNodesPerLev,connBg[3]+nbOfNodesPerLev,connBg[4]+nbOfNodesPerLev
05310         };
05311         ret.insert(ret.end(),conn,conn+20);
05312         break;
05313       }
05314     case INTERP_KERNEL::NORM_POLYGON:
05315       {
05316         std::back_insert_iterator< std::vector<int> > ii(ret);
05317         std::copy(connBg+1,connEnd,ii);
05318         *ii++=-1;
05319         std::reverse_iterator<const int *> rConnBg(connEnd);
05320         std::reverse_iterator<const int *> rConnEnd(connBg+1);
05321         std::transform(rConnBg,rConnEnd,ii,std::bind2nd(std::plus<int>(),deltaz));
05322         std::size_t nbOfRadFaces=std::distance(connBg+1,connEnd);
05323         for(std::size_t i=0;i<nbOfRadFaces;i++)
05324           {
05325             *ii++=-1;
05326             int conn[4]={connBg[(i+1)%nbOfRadFaces+1],connBg[i+1],connBg[i+1]+deltaz,connBg[(i+1)%nbOfRadFaces+1]+deltaz};
05327             std::copy(conn,conn+4,ii);
05328           }
05329         break;
05330       }
05331     default:
05332       throw INTERP_KERNEL::Exception("A flat type has been detected that has not its extruded representation !");
05333     }
05334 }
05335 
05339 bool MEDCouplingUMesh::IsPolygonWellOriented(bool isQuadratic, const double *vec, const int *begin, const int *end, const double *coords)
05340 {
05341   double v[3]={0.,0.,0.};
05342   std::size_t sz=std::distance(begin,end);
05343   if(isQuadratic)
05344     sz/=2;
05345   for(std::size_t i=0;i<sz;i++)
05346     {
05347       v[0]+=coords[3*begin[i]+1]*coords[3*begin[(i+1)%sz]+2]-coords[3*begin[i]+2]*coords[3*begin[(i+1)%sz]+1];
05348       v[1]+=coords[3*begin[i]+2]*coords[3*begin[(i+1)%sz]]-coords[3*begin[i]]*coords[3*begin[(i+1)%sz]+2];
05349       v[2]+=coords[3*begin[i]]*coords[3*begin[(i+1)%sz]+1]-coords[3*begin[i]+1]*coords[3*begin[(i+1)%sz]];
05350     }
05351   return vec[0]*v[0]+vec[1]*v[1]+vec[2]*v[2]>0.;
05352 }
05353 
05357 bool MEDCouplingUMesh::IsPolyhedronWellOriented(const int *begin, const int *end, const double *coords)
05358 {
05359   std::vector<std::pair<int,int> > edges;
05360   std::size_t nbOfFaces=std::count(begin,end,-1)+1;
05361   const int *bgFace=begin;
05362   for(std::size_t i=0;i<nbOfFaces;i++)
05363     {
05364       const int *endFace=std::find(bgFace+1,end,-1);
05365       std::size_t nbOfEdgesInFace=std::distance(bgFace,endFace);
05366       for(std::size_t j=0;j<nbOfEdgesInFace;j++)
05367         {
05368           std::pair<int,int> p1(bgFace[j],bgFace[(j+1)%nbOfEdgesInFace]);
05369           if(std::find(edges.begin(),edges.end(),p1)!=edges.end())
05370             return false;
05371           edges.push_back(p1);
05372         }
05373       bgFace=endFace+1;
05374     }
05375   return INTERP_KERNEL::calculateVolumeForPolyh2<int,INTERP_KERNEL::ALL_C_MODE>(begin,(int)std::distance(begin,end),coords)>-EPS_FOR_POLYH_ORIENTATION;
05376 }
05377 
05382 void MEDCouplingUMesh::TryToCorrectPolyhedronOrientation(int *begin, int *end, const double *coords) throw(INTERP_KERNEL::Exception)
05383 {
05384   std::vector<std::pair<int,int> > edges;
05385   std::size_t nbOfFaces=std::count(begin,end,-1)+1;
05386   int *bgFace=begin;
05387   std::vector<bool> isPerm(nbOfFaces);
05388   for(std::size_t i=0;i<nbOfFaces;i++)
05389     {
05390       int *endFace=std::find(bgFace+1,end,-1);
05391       std::size_t nbOfEdgesInFace=std::distance(bgFace,endFace);
05392       for(std::size_t l=0;l<nbOfEdgesInFace;l++)
05393         {
05394           std::pair<int,int> p1(bgFace[l],bgFace[(l+1)%nbOfEdgesInFace]);
05395           edges.push_back(p1);
05396         }
05397       int *bgFace2=endFace+1;
05398       for(std::size_t k=i+1;k<nbOfFaces;k++)
05399         {
05400           int *endFace2=std::find(bgFace2+1,end,-1);
05401           std::size_t nbOfEdgesInFace2=std::distance(bgFace2,endFace2);
05402           for(std::size_t j=0;j<nbOfEdgesInFace2;j++)
05403             {
05404               std::pair<int,int> p2(bgFace2[j],bgFace2[(j+1)%nbOfEdgesInFace2]);
05405               if(std::find(edges.begin(),edges.end(),p2)!=edges.end())
05406                 {
05407                   if(isPerm[k])
05408                     throw INTERP_KERNEL::Exception("Fail to repare polyhedron ! Polyedron looks bad !");
05409                   std::vector<int> tmp(nbOfEdgesInFace2-1);
05410                   std::copy(bgFace2+1,endFace2,tmp.rbegin());
05411                   std::copy(tmp.begin(),tmp.end(),bgFace2+1);
05412                   isPerm[k]=true;
05413                   continue;
05414                 }
05415             }
05416           bgFace2=endFace2+1;
05417         }
05418       bgFace=endFace+1;
05419     }
05420   if(INTERP_KERNEL::calculateVolumeForPolyh2<int,INTERP_KERNEL::ALL_C_MODE>(begin,(int)std::distance(begin,end),coords)<-EPS_FOR_POLYH_ORIENTATION)
05421     {//not lucky ! The first face was not correctly oriented : reorient all faces...
05422       bgFace=begin;
05423       for(std::size_t i=0;i<nbOfFaces;i++)
05424         {
05425           int *endFace=std::find(bgFace+1,end,-1);
05426           std::size_t nbOfEdgesInFace=std::distance(bgFace,endFace);
05427           std::vector<int> tmp(nbOfEdgesInFace-1);
05428           std::copy(bgFace+1,endFace,tmp.rbegin());
05429           std::copy(tmp.begin(),tmp.end(),bgFace+1);
05430           bgFace=endFace+1;
05431         }
05432     }
05433 }
05434 
05439 void MEDCouplingUMesh::FillInCompact3DMode(int spaceDim, int nbOfNodesInCell, const int *conn, const double *coo, double *zipFrmt) throw(INTERP_KERNEL::Exception)
05440 {
05441   double *w=zipFrmt;
05442   if(spaceDim==3)
05443     for(int i=0;i<nbOfNodesInCell;i++)
05444       w=std::copy(coo+3*conn[i],coo+3*conn[i]+3,w);
05445   else if(spaceDim==2)
05446     {
05447       for(int i=0;i<nbOfNodesInCell;i++)
05448         {
05449           w=std::copy(coo+2*conn[i],coo+2*conn[i]+2,w);
05450           *w++=0.;
05451         }
05452     }
05453   else
05454     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::FillInCompact3DMode : Invalid spaceDim specified : must be 2 or 3 !");
05455 }
05456 
05457 void MEDCouplingUMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData) const throw(INTERP_KERNEL::Exception)
05458 {
05459   if(getNumberOfCells()<=0)
05460     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::writeVTK : the unstructured mesh has no cells !");
05461   static const int PARAMEDMEM2VTKTYPETRADUCER[INTERP_KERNEL::NORM_MAXTYPE+1]={1,3,21,5,9,7,22,-1,23,-1,-1,-1,-1,-1,10,14,13,-1,12,-1,24,-1,16,27,-1,26,-1,-1,-1,-1,25,42,-1,4};
05462   ofs << "  <" << getVTKDataSetType() << ">\n";
05463   ofs << "    <Piece NumberOfPoints=\"" << getNumberOfNodes() << "\" NumberOfCells=\"" << getNumberOfCells() << "\">\n";
05464   ofs << "      <PointData>\n" << pointData << std::endl;
05465   ofs << "      </PointData>\n";
05466   ofs << "      <CellData>\n" << cellData << std::endl;
05467   ofs << "      </CellData>\n";
05468   ofs << "      <Points>\n";
05469   if(getSpaceDimension()==3)
05470     _coords->writeVTK(ofs,8,"Points");
05471   else
05472     {
05473       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo=_coords->changeNbOfComponents(3,0.);
05474       coo->writeVTK(ofs,8,"Points");
05475     }
05476   ofs << "      </Points>\n";
05477   ofs << "      <Cells>\n";
05478   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c0=_nodal_connec_index->buildComplement(_nodal_connec->getNumberOfTuples()+1);
05479   c0=_nodal_connec->selectByTupleId(c0->begin(),c0->end());
05480   c0->writeVTK(ofs,8,"Int64","connectivity");
05481   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c1=_nodal_connec_index->deltaShiftIndex();
05482   c1->applyLin(1,-1);
05483   c1->computeOffsets2();
05484   c1=c1->selectByTupleId2(1,c1->getNumberOfTuples(),1);
05485   c1->writeVTK(ofs,8,"Int64","offsets");
05486   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c2=_nodal_connec->selectByTupleId(_nodal_connec_index->getConstPointer(),_nodal_connec_index->getConstPointer()+getNumberOfCells());
05487   c2->transformWithIndArr(PARAMEDMEM2VTKTYPETRADUCER,PARAMEDMEM2VTKTYPETRADUCER+INTERP_KERNEL::NORM_MAXTYPE);
05488   c2->writeVTK(ofs,8,"UInt8","types");
05489   ofs << "      </Cells>\n";
05490   ofs << "    </Piece>\n";
05491   ofs << "  </" << getVTKDataSetType() << ">\n";
05492 }
05493 
05494 std::string MEDCouplingUMesh::getVTKDataSetType() const throw(INTERP_KERNEL::Exception)
05495 {
05496   return std::string("UnstructuredGrid");
05497 }
05498 
05500 
05501 MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2) throw(INTERP_KERNEL::Exception)
05502 {
05503   m1->checkFullyDefined();
05504   m2->checkFullyDefined();
05505   if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2)
05506     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2  with meshdim equal to 2 and spaceDim equal to 2 too!");
05507   std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
05508   MEDCouplingUMesh *m1Desc=0,*m2Desc=0;
05509   DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0;
05510   std::vector<double> addCoo,addCoordsQuadratic;
05511   INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
05512   INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
05513   IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,colinear2, subDiv2,m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
05514                               m2Desc,desc2,descIndx2,revDesc2,revDescIndx2,addCoo);
05515   revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef();
05516   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(desc2),dd4(descIndx2);
05517   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> dd5(m1Desc),dd6(m2Desc);
05518   std::vector< std::vector<int> > intersectEdge2;
05519   BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2);
05520   subDiv2.clear(); dd5=0; dd6=0;
05521   std::vector<int> cr,crI;
05522   std::vector<int> cNb1,cNb2;
05523   BuildIntersecting2DCellsFromEdges(eps,m1,desc1->getConstPointer(),descIndx1->getConstPointer(),intersectEdge1,colinear2,m2,desc2->getConstPointer(),descIndx2->getConstPointer(),intersectEdge2,addCoo,
05524                                     /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2);
05525   //
05526   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> addCooDa=DataArrayDouble::New();
05527   addCooDa->alloc((int)(addCoo.size())/2,2);
05528   std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
05529   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> addCoordsQuadraticDa=DataArrayDouble::New();
05530   addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2);
05531   std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
05532   std::vector<const DataArrayDouble *> coordss(4);
05533   coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
05534   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo=DataArrayDouble::Aggregate(coordss);
05535   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret=MEDCouplingUMesh::New("Intersect2D",2);
05536   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn=DataArrayInt::New(); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
05537   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> connI=DataArrayInt::New(); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
05538   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c1=DataArrayInt::New(); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer()); cellNb1=c1;
05539   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c2=DataArrayInt::New(); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer()); cellNb2=c2;
05540   ret->setConnectivity(conn,connI,true);
05541   ret->setCoords(coo);
05542   ret->incrRef(); c1->incrRef(); c2->incrRef();
05543   return ret;
05544 }
05545 
05547 
05548 void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1,
05549                                                          const std::vector<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
05550                                                          const MEDCouplingUMesh *m2, const int *desc2, const int *descIndx2, const std::vector<std::vector<int> >& intesctEdges2,
05551                                                          const std::vector<double>& addCoords,
05552                                                          std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI, std::vector<int>& cNb1, std::vector<int>& cNb2)
05553 {
05554   static const int SPACEDIM=2;
05555   std::vector<double> bbox1,bbox2;
05556   const double *coo1=m1->getCoords()->getConstPointer();
05557   const int *conn1=m1->getNodalConnectivity()->getConstPointer();
05558   const int *connI1=m1->getNodalConnectivityIndex()->getConstPointer();
05559   int offset1=m1->getNumberOfNodes();
05560   const double *coo2=m2->getCoords()->getConstPointer();
05561   const int *conn2=m2->getNodalConnectivity()->getConstPointer();
05562   const int *connI2=m2->getNodalConnectivityIndex()->getConstPointer();
05563   int offset2=offset1+m2->getNumberOfNodes();
05564   int offset3=offset2+((int)addCoords.size())/2;
05565   m1->getBoundingBoxForBBTree(bbox1);
05566   m2->getBoundingBoxForBBTree(bbox2);
05567   BBTree<SPACEDIM,int> myTree(&bbox2[0],0,0,m2->getNumberOfCells(),eps);
05568   int ncell1=m1->getNumberOfCells();
05569   crI.push_back(0);
05570   for(int i=0;i<ncell1;i++)
05571     {
05572       std::vector<int> candidates2;
05573       myTree.getIntersectingElems(&bbox1[i*2*SPACEDIM],candidates2);
05574       std::map<INTERP_KERNEL::Node *,int> mapp;
05575       std::map<int,INTERP_KERNEL::Node *> mappRev;
05576       INTERP_KERNEL::QuadraticPolygon pol1;
05577       INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]];
05578       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
05579       MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
05580       pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
05581                                    desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
05582       std::vector<int> crTmp,crITmp;
05583       crITmp.push_back(crI.back());
05584       for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++)
05585         {
05586           INTERP_KERNEL::QuadraticPolygon pol2;
05587           pol1.initLocations();
05588           MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
05589           INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
05590           const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
05591           pol2.buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
05592                                         pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
05593           //MEDCouplingUMeshAssignOnLoc(pol1,pol2,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,colinear2);
05594           pol1.buildPartitionsAbs(pol2,mapp,i,*it2,offset3,addCoordsQuadratic,cr,crI,cNb1,cNb2);
05595         }
05596       if(!crTmp.empty())
05597         {
05598           cr.insert(cr.end(),crTmp.begin(),crTmp.end());
05599           crI.insert(crI.end(),crITmp.begin()+1,crITmp.end());
05600         }
05601       for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
05602         (*it).second->decrRef();
05603     }
05604 }
05605 
05610 void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
05611                                                    std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
05612                                                    MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
05613                                                    MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2,
05614                                                    std::vector<double>& addCoo) throw(INTERP_KERNEL::Exception)
05615 {
05616   static const int SPACEDIM=2;
05617   desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New();
05618   desc2=DataArrayInt::New();
05619   descIndx2=DataArrayInt::New();
05620   revDesc2=DataArrayInt::New();
05621   revDescIndx2=DataArrayInt::New();
05622   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
05623   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd5(desc2),dd6(descIndx2),dd7(revDesc2),dd8(revDescIndx2);
05624   m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
05625   m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2);
05626   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> dd9(m1Desc),dd10(m2Desc);
05627   const int *c1=m1Desc->getNodalConnectivity()->getConstPointer();
05628   const int *ci1=m1Desc->getNodalConnectivityIndex()->getConstPointer();
05629   std::vector<double> bbox1,bbox2;
05630   m1Desc->getBoundingBoxForBBTree(bbox1);
05631   m2Desc->getBoundingBoxForBBTree(bbox2);
05632   int ncell1=m1Desc->getNumberOfCells();
05633   int ncell2=m2Desc->getNumberOfCells();
05634   intersectEdge1.resize(ncell1);
05635   colinear2.resize(ncell2);
05636   subDiv2.resize(ncell2);
05637   BBTree<SPACEDIM,int> myTree(&bbox2[0],0,0,m2Desc->getNumberOfCells(),-eps);
05638   std::vector<int> candidates1(1);
05639   int offset1=m1->getNumberOfNodes();
05640   int offset2=offset1+m2->getNumberOfNodes();
05641   for(int i=0;i<ncell1;i++)
05642     {
05643       std::vector<int> candidates2;
05644       myTree.getIntersectingElems(&bbox1[i*2*SPACEDIM],candidates2);
05645       if(!candidates2.empty())
05646         {
05647           std::map<INTERP_KERNEL::Node *,int> map1,map2;
05648           INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2);
05649           candidates1[0]=i;
05650           INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1);
05651           pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo);
05652           delete pol2;
05653           delete pol1;
05654         }
05655       else
05656         intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i+1]);
05657     }
05658   m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef();
05659   m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef();
05660 }
05661 
05673 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) throw(INTERP_KERNEL::Exception)
05674 {
05675   int offset1=m1->getNumberOfNodes();
05676   int ncell=m2->getNumberOfCells();
05677   const int *c=m2->getNodalConnectivity()->getConstPointer();
05678   const int *cI=m2->getNodalConnectivityIndex()->getConstPointer();
05679   const double *coo=m2->getCoords()->getConstPointer();
05680   const double *cooBis=m1->getCoords()->getConstPointer();
05681   int offset2=offset1+m2->getNumberOfNodes();
05682   intersectEdge.resize(ncell);
05683   for(int i=0;i<ncell;i++,cI++)
05684     {
05685       const std::vector<int>& divs=subDiv[i];
05686       int nnode=cI[1]-cI[0]-1;
05687       std::map<int, std::pair<INTERP_KERNEL::Node *,bool> > mapp2;
05688       std::map<INTERP_KERNEL::Node *, int> mapp22;
05689       for(int j=0;j<nnode;j++)
05690         {
05691           INTERP_KERNEL::Node *nn=new INTERP_KERNEL::Node(coo[2*c[(*cI)+j+1]],coo[2*c[(*cI)+j+1]+1]);
05692           int nnid=c[(*cI)+j+1];
05693           mapp2[nnid]=std::pair<INTERP_KERNEL::Node *,bool>(nn,true);
05694           mapp22[nn]=nnid+offset1;
05695         }
05696       INTERP_KERNEL::Edge *e=MEDCouplingUMeshBuildQPFromEdge((INTERP_KERNEL::NormalizedCellType)c[*cI],mapp2,c+(*cI)+1);
05697       for(std::map<int, std::pair<INTERP_KERNEL::Node *,bool> >::const_iterator it=mapp2.begin();it!=mapp2.end();it++)
05698         ((*it).second.first)->decrRef();
05699       std::vector<INTERP_KERNEL::Node *> addNodes(divs.size());
05700       std::map<INTERP_KERNEL::Node *,int> mapp3;
05701       for(std::size_t j=0;j<divs.size();j++)
05702         {
05703           int id=divs[j];
05704           INTERP_KERNEL::Node *tmp=0;
05705           if(id<offset1)
05706             tmp=new INTERP_KERNEL::Node(cooBis[2*id],cooBis[2*id+1]);
05707           else if(id<offset2)
05708             tmp=new INTERP_KERNEL::Node(coo[2*(id-offset1)],coo[2*(id-offset1)+1]);//if it happens, bad news mesh 'm2' is non conform.
05709           else
05710             tmp=new INTERP_KERNEL::Node(addCoo[2*(id-offset2)],addCoo[2*(id-offset2)+1]);
05711           addNodes[j]=tmp;
05712           mapp3[tmp]=id;
05713         }
05714       e->sortIdsAbs(addNodes,mapp22,mapp3,intersectEdge[i]);
05715       for(std::vector<INTERP_KERNEL::Node *>::const_iterator it=addNodes.begin();it!=addNodes.end();it++)
05716         (*it)->decrRef();
05717       e->decrRef();
05718     }
05719 }
05720 
05735 void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector<int>& cut3DCurve, std::vector<int>& nodesOnPlane, const int *nodal3DSurf, const int *nodalIndx3DSurf,
05736                                                    const int *nodal3DCurve, const int *nodalIndx3DCurve,
05737                                                    const int *desc, const int *descIndx, 
05738                                                    std::vector< std::pair<int,int> >& cut3DSurf) throw(INTERP_KERNEL::Exception)
05739 {
05740   std::set<int> nodesOnP(nodesOnPlane.begin(),nodesOnPlane.end());
05741   int nbOf3DSurfCell=(int)cut3DSurf.size();
05742   for(int i=0;i<nbOf3DSurfCell;i++)
05743     {
05744       std::vector<int> res;
05745       int offset=descIndx[i];
05746       int nbOfSeg=descIndx[i+1]-offset;
05747       for(int j=0;j<nbOfSeg;j++)
05748         {
05749           int edgeId=desc[offset+j];
05750           int status=cut3DCurve[edgeId];
05751           if(status!=-2)
05752             {
05753               if(status>-1)
05754                 res.push_back(status);
05755               else
05756                 {
05757                   res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+1]);
05758                   res.push_back(nodal3DCurve[nodalIndx3DCurve[edgeId]+2]);
05759                 }
05760             }
05761         }
05762       switch(res.size())
05763         {
05764         case 2:
05765           {
05766             cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
05767             break;
05768           }
05769         case 1:
05770         case 0:
05771           {
05772             std::set<int> s1(nodal3DSurf+nodalIndx3DSurf[i]+1,nodal3DSurf+nodalIndx3DSurf[i+1]);
05773             std::set_intersection(nodesOnP.begin(),nodesOnP.end(),s1.begin(),s1.end(),std::back_insert_iterator< std::vector<int> >(res));
05774             if(res.size()==2)
05775               {
05776                 cut3DSurf[i].first=res[0]; cut3DSurf[i].second=res[1];
05777               }
05778             else
05779               {
05780                 cut3DSurf[i].first=-1; cut3DSurf[i].second=-1;
05781               }
05782             break;
05783           }
05784         default:
05785           {// case when plane is on a multi colinear edge of a polyhedron
05786             if((int)res.size()==2*nbOfSeg)
05787               {
05788                 cut3DSurf[i].first=-2; cut3DSurf[i].second=i;
05789               }
05790             else
05791               throw INTERP_KERNEL::Exception("MEDCouplingUMesh::AssemblyPointsFrom3DCurve : unexpected situation !");
05792           }
05793         }
05794     }
05795 }
05796 
05806 void MEDCouplingUMesh::assemblyForSplitFrom3DSurf(const std::vector< std::pair<int,int> >& cut3DSurf,
05807                                                   const int *desc, const int *descIndx,
05808                                                   std::vector<int>& nodalRes, std::vector<int>& nodalResIndx, std::vector<int>& cellIds) const throw(INTERP_KERNEL::Exception)
05809 {
05810   checkFullyDefined();
05811   if(getMeshDimension()!=3 || getSpaceDimension()!=3)
05812     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::assemblyForSplitFrom3DSurf works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!");
05813   const int *nodal3D=_nodal_connec->getConstPointer();
05814   const int *nodalIndx3D=_nodal_connec_index->getConstPointer();
05815   int nbOfCells=getNumberOfCells();
05816   for(int i=0;i<nbOfCells;i++)
05817     {
05818       std::map<int, std::set<int> > m;
05819       int offset=descIndx[i];
05820       int nbOfFaces=descIndx[i+1]-offset;
05821       int start=-1;
05822       int end=-1;
05823       for(int j=0;j<nbOfFaces;j++)
05824         {
05825           const std::pair<int,int>& p=cut3DSurf[desc[offset+j]];
05826           if(p.first!=-1 && p.second!=-1)
05827             {
05828               if(p.first!=-2)
05829                 {
05830                   start=p.first; end=p.second;
05831                   m[p.first].insert(p.second);
05832                   m[p.second].insert(p.first);
05833                 }
05834               else
05835                 {
05836                   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodal3D[nodalIndx3D[i]]);
05837                   int sz=nodalIndx3D[i+1]-nodalIndx3D[i]-1;
05838                   INTERP_KERNEL::AutoPtr<int> tmp=new int[sz];
05839                   INTERP_KERNEL::NormalizedCellType cmsId;
05840                   unsigned nbOfNodesSon=cm.fillSonCellNodalConnectivity2(j,nodal3D+nodalIndx3D[i]+1,sz,tmp,cmsId);
05841                   start=tmp[0]; end=tmp[nbOfNodesSon-1];
05842                   for(unsigned k=0;k<nbOfNodesSon;k++)
05843                     {
05844                       m[tmp[k]].insert(tmp[(k+1)%nbOfNodesSon]);
05845                       m[tmp[(k+1)%nbOfNodesSon]].insert(tmp[k]);
05846                     }
05847                 }
05848             }
05849         }
05850       if(m.empty())
05851         continue;
05852       std::vector<int> conn(1,(int)INTERP_KERNEL::NORM_POLYGON);
05853       int prev=end;
05854       while(end!=start)
05855         {
05856           std::map<int, std::set<int> >::const_iterator it=m.find(start);
05857           const std::set<int>& s=(*it).second;
05858           std::set<int> s2; s2.insert(prev);
05859           std::set<int> s3;
05860           std::set_difference(s.begin(),s.end(),s2.begin(),s2.end(),inserter(s3,s3.begin()));
05861           if(s3.size()==1)
05862             {
05863               int val=*s3.begin();
05864               conn.push_back(start);
05865               prev=start;
05866               start=val;
05867             }
05868           else
05869             start=end;
05870         }
05871       conn.push_back(end);
05872       if(conn.size()>3)
05873         {
05874           nodalRes.insert(nodalRes.end(),conn.begin(),conn.end());
05875           nodalResIndx.push_back((int)nodalRes.size());
05876           cellIds.push_back(i);
05877         }
05878     }
05879 }
05880 
05884 void MEDCouplingUMesh::BuildUnionOf2DMesh(const std::vector<int>& conn2D, const std::vector<int>& connI2D, std::vector<int>& polyUnion)
05885 {
05886 }
05887 
05888 MEDCouplingUMeshCellIterator::MEDCouplingUMeshCellIterator(MEDCouplingUMesh *mesh):_mesh(mesh),_cell(new MEDCouplingUMeshCell(mesh)),
05889                                                                                    _own_cell(true),_cell_id(-1),_nb_cell(0)
05890 {
05891   if(mesh)
05892     {
05893       mesh->incrRef();
05894       _nb_cell=mesh->getNumberOfCells();
05895     }
05896 }
05897 
05898 MEDCouplingUMeshCellIterator::~MEDCouplingUMeshCellIterator()
05899 {
05900   if(_mesh)
05901     _mesh->decrRef();
05902   if(_own_cell)
05903     delete _cell;
05904 }
05905 
05906 MEDCouplingUMeshCellIterator::MEDCouplingUMeshCellIterator(MEDCouplingUMesh *mesh, MEDCouplingUMeshCell *itc, int bg, int end):_mesh(mesh),_cell(itc),
05907                                                                                                                                _own_cell(false),_cell_id(bg-1),
05908                                                                                                                                _nb_cell(end)
05909 {
05910   if(mesh)
05911     mesh->incrRef();
05912 }
05913 
05914 MEDCouplingUMeshCell *MEDCouplingUMeshCellIterator::nextt()
05915 {
05916   _cell_id++;
05917   if(_cell_id<_nb_cell)
05918     {
05919       _cell->next();
05920       return _cell;
05921     }
05922   else
05923     return 0;
05924 }
05925 
05926 MEDCouplingUMeshCellByTypeEntry::MEDCouplingUMeshCellByTypeEntry(MEDCouplingUMesh *mesh):_mesh(mesh)
05927 {
05928   if(_mesh)
05929     _mesh->incrRef();
05930 }
05931 
05932 MEDCouplingUMeshCellByTypeIterator *MEDCouplingUMeshCellByTypeEntry::iterator()
05933 {
05934   return new MEDCouplingUMeshCellByTypeIterator(_mesh);
05935 }
05936 
05937 MEDCouplingUMeshCellByTypeEntry::~MEDCouplingUMeshCellByTypeEntry()
05938 {
05939   if(_mesh)
05940     _mesh->decrRef();
05941 }
05942 
05943 MEDCouplingUMeshCellEntry::MEDCouplingUMeshCellEntry(MEDCouplingUMesh *mesh,  INTERP_KERNEL::NormalizedCellType type, MEDCouplingUMeshCell *itc, int bg, int end):_mesh(mesh),_type(type),
05944                                                                                                                                                                   _itc(itc),
05945                                                                                                                                                                   _bg(bg),_end(end)
05946 {
05947   if(_mesh)
05948     _mesh->incrRef();
05949 }
05950 
05951 MEDCouplingUMeshCellEntry::~MEDCouplingUMeshCellEntry()
05952 {
05953   if(_mesh)
05954     _mesh->decrRef();
05955 }
05956 
05957 INTERP_KERNEL::NormalizedCellType MEDCouplingUMeshCellEntry::getType() const
05958 {
05959   return _type;
05960 }
05961 
05962 int MEDCouplingUMeshCellEntry::getNumberOfElems() const
05963 {
05964   return _end-_bg;
05965 }
05966 
05967 MEDCouplingUMeshCellIterator *MEDCouplingUMeshCellEntry::iterator()
05968 {
05969   return new MEDCouplingUMeshCellIterator(_mesh,_itc,_bg,_end);
05970 }
05971 
05972 MEDCouplingUMeshCellByTypeIterator::MEDCouplingUMeshCellByTypeIterator(MEDCouplingUMesh *mesh):_mesh(mesh),_cell(new MEDCouplingUMeshCell(mesh)),_cell_id(0),_nb_cell(0)
05973 {
05974   if(mesh)
05975     {
05976       mesh->incrRef();
05977       _nb_cell=mesh->getNumberOfCells();
05978     }
05979 }
05980 
05981 MEDCouplingUMeshCellByTypeIterator::~MEDCouplingUMeshCellByTypeIterator()
05982 {
05983   if(_mesh)
05984     _mesh->decrRef();
05985   delete _cell;
05986 }
05987 
05988 MEDCouplingUMeshCellEntry *MEDCouplingUMeshCellByTypeIterator::nextt()
05989 {
05990   const int *c=_mesh->getNodalConnectivity()->getConstPointer();
05991   const int *ci=_mesh->getNodalConnectivityIndex()->getConstPointer();
05992   if(_cell_id<_nb_cell)
05993     {
05994       INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)c[ci[_cell_id]];
05995       int nbOfElems=(int)std::distance(ci+_cell_id,std::find_if(ci+_cell_id,ci+_nb_cell,ParaMEDMEMImpl::ConnReader(c,type)));
05996       int startId=_cell_id;
05997       _cell_id+=nbOfElems;
05998       return new MEDCouplingUMeshCellEntry(_mesh,type,_cell,startId,_cell_id);
05999     }
06000   else
06001     return 0;
06002 }
06003 
06004 MEDCouplingUMeshCell::MEDCouplingUMeshCell(MEDCouplingUMesh *mesh):_conn(0),_conn_indx(0),_conn_lgth(NOTICABLE_FIRST_VAL)
06005 {
06006   if(mesh)
06007     {
06008       _conn=mesh->getNodalConnectivity()->getPointer();
06009       _conn_indx=mesh->getNodalConnectivityIndex()->getPointer();
06010     }
06011 }
06012 
06013 void MEDCouplingUMeshCell::next()
06014 {
06015   if(_conn_lgth!=NOTICABLE_FIRST_VAL)
06016     {
06017       _conn+=_conn_lgth;
06018       _conn_indx++;
06019     }
06020   _conn_lgth=_conn_indx[1]-_conn_indx[0];
06021 }
06022 
06023 std::string MEDCouplingUMeshCell::repr() const
06024 {
06025   if(_conn_lgth!=NOTICABLE_FIRST_VAL)
06026     {
06027       std::ostringstream oss; oss << "Cell Type " << INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)_conn[0]).getRepr();
06028       oss << " : ";
06029       std::copy(_conn+1,_conn+_conn_lgth,std::ostream_iterator<int>(oss," "));
06030       return oss.str();
06031     }
06032   else
06033     return std::string("MEDCouplingUMeshCell::repr : Invalid pos");
06034 }
06035 
06036 INTERP_KERNEL::NormalizedCellType MEDCouplingUMeshCell::getType() const
06037 {
06038   if(_conn_lgth!=NOTICABLE_FIRST_VAL)
06039     return (INTERP_KERNEL::NormalizedCellType)_conn[0];
06040   else
06041     return INTERP_KERNEL::NORM_ERROR;
06042 }
06043 
06044 const int *MEDCouplingUMeshCell::getAllConn(int& lgth) const
06045 {
06046   lgth=_conn_lgth;
06047   if(_conn_lgth!=NOTICABLE_FIRST_VAL)
06048     return _conn;
06049   else
06050     return 0;
06051 }