Salome HOME
Some useful tools for MEDCouplingMappedExtrudedMesh users that needs to store it...
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingStructuredMesh.cxx
index 1e2186855f2412843fee22ca0679d26bb8e6fba9..aa83f7f1591ed91a6fe754afc9cd62c144cb5398 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2014  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2016  CEA/DEN, EDF R&D
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 #include "MEDCouplingMemArray.hxx"
 #include "MEDCoupling1GTUMesh.hxx"
 #include "MEDCouplingUMesh.hxx"
+#include "MEDCouplingIMesh.hxx"//tony to throw when optimization will be performed in AssignPartOfFieldOfDoubleUsing
 
 #include <numeric>
 
-using namespace ParaMEDMEM;
+using namespace MEDCoupling;
 
 MEDCouplingStructuredMesh::MEDCouplingStructuredMesh()
 {
@@ -96,7 +97,7 @@ int MEDCouplingStructuredMesh::getNumberOfCellsWithType(INTERP_KERNEL::Normalize
 
 DataArrayInt *MEDCouplingStructuredMesh::giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
 {
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
+  MCAuto<DataArrayInt> ret=DataArrayInt::New();
   if(getTypeOfCell(0)==type)
     {
       ret->alloc(getNumberOfCells(),1);
@@ -110,7 +111,7 @@ DataArrayInt *MEDCouplingStructuredMesh::giveCellsWithType(INTERP_KERNEL::Normal
 DataArrayInt *MEDCouplingStructuredMesh::computeNbOfNodesPerCell() const
 {
   int nbCells=getNumberOfCells();
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
+  MCAuto<DataArrayInt> ret=DataArrayInt::New();
   ret->alloc(nbCells,1);
   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(getTypeOfCell(0));
   ret->fillWithValue((int)cm.getNumberOfNodes());
@@ -120,7 +121,7 @@ DataArrayInt *MEDCouplingStructuredMesh::computeNbOfNodesPerCell() const
 DataArrayInt *MEDCouplingStructuredMesh::computeNbOfFacesPerCell() const
 {
   int nbCells=getNumberOfCells();
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
+  MCAuto<DataArrayInt> ret=DataArrayInt::New();
   ret->alloc(nbCells,1);
   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(getTypeOfCell(0));
   ret->fillWithValue((int)cm.getNumberOfSons());
@@ -315,15 +316,15 @@ void MEDCouplingStructuredMesh::splitProfilePerType(const DataArrayInt *profile,
     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::splitProfilePerType : input profile is NULL or not allocated !");
   if(profile->getNumberOfComponents()!=1)
     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::splitProfilePerType : input profile should have exactly one component !");
-  int nbTuples=profile->getNumberOfTuples();
+  int nbTuples(profile->getNumberOfTuples());
   int nbOfCells=getNumberOfCells();
   code.resize(3); idsInPflPerType.resize(1);
   code[0]=(int)getTypeOfCell(0); code[1]=nbOfCells;
   idsInPflPerType.resize(1);
-  if(profile->isIdentity() && nbTuples==nbOfCells)
+  if(profile->isIota(nbOfCells))
     {
       code[2]=-1;
-      idsInPflPerType[0]=0;
+      idsInPflPerType[0]=profile->deepCopy();
       idsPerType.clear();
       return ;
     }
@@ -331,7 +332,7 @@ void MEDCouplingStructuredMesh::splitProfilePerType(const DataArrayInt *profile,
   code[2]=0;
   profile->checkAllIdsInRange(0,nbOfCells);
   idsPerType.resize(1);
-  idsPerType[0]=profile->deepCpy();
+  idsPerType[0]=profile->deepCopy();
   idsInPflPerType[0]=DataArrayInt::Range(0,nbTuples,1);
 }
 
@@ -346,11 +347,11 @@ MEDCoupling1SGTUMesh *MEDCouplingStructuredMesh::build1SGTUnstructured() const
   int meshDim(getMeshDimension()),spaceDim(getSpaceDimensionOnNodeStruct());
   if((meshDim<0 || meshDim>3) || (spaceDim<0 || spaceDim>3))
     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::build1SGTUnstructured : meshdim and spacedim must be in [1,2,3] !");
-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords(getCoordinatesAndOwner());
+  MCAuto<DataArrayDouble> coords(getCoordinatesAndOwner());
   int ns[3];
   getNodeGridStructure(ns);
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(Build1GTNodalConnectivity(ns,ns+spaceDim));
-  MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> ret(MEDCoupling1SGTUMesh::New(getName(),GetGeoTypeGivenMeshDimension(meshDim)));
+  MCAuto<DataArrayInt> conn(Build1GTNodalConnectivity(ns,ns+spaceDim));
+  MCAuto<MEDCoupling1SGTUMesh> ret(MEDCoupling1SGTUMesh::New(getName(),GetGeoTypeGivenMeshDimension(meshDim)));
   ret->setNodalConnectivity(conn); ret->setCoords(coords);
   try
     { ret->copyTinyInfoFrom(this); }
@@ -369,11 +370,11 @@ MEDCoupling1SGTUMesh *MEDCouplingStructuredMesh::build1SGTSubLevelMesh() const
   int meshDim(getMeshDimension());
   if(meshDim<1 || meshDim>3)
     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::build1SGTSubLevelMesh : meshdim must be in [2,3] !");
-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords(getCoordinatesAndOwner());
+  MCAuto<DataArrayDouble> coords(getCoordinatesAndOwner());
   int ns[3];
   getNodeGridStructure(ns);
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(Build1GTNodalConnectivityOfSubLevelMesh(ns,ns+meshDim));
-  MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> ret(MEDCoupling1SGTUMesh::New(getName(),GetGeoTypeGivenMeshDimension(meshDim-1)));
+  MCAuto<DataArrayInt> conn(Build1GTNodalConnectivityOfSubLevelMesh(ns,ns+meshDim));
+  MCAuto<MEDCoupling1SGTUMesh> ret(MEDCoupling1SGTUMesh::New(getName(),GetGeoTypeGivenMeshDimension(meshDim-1)));
   ret->setNodalConnectivity(conn); ret->setCoords(coords);
   return ret.retn();
 }
@@ -386,7 +387,7 @@ MEDCoupling1SGTUMesh *MEDCouplingStructuredMesh::build1SGTSubLevelMesh() const
  */
 MEDCouplingUMesh *MEDCouplingStructuredMesh::buildUnstructured() const
 {
-  MEDCouplingAutoRefCountObjectPtr<MEDCoupling1SGTUMesh> ret0(build1SGTUnstructured());
+  MCAuto<MEDCoupling1SGTUMesh> ret0(build1SGTUnstructured());
   return ret0->buildUnstructured();
 }
 
@@ -414,14 +415,14 @@ MEDCouplingMesh *MEDCouplingStructuredMesh::buildPartAndReduceNodes(const int *s
   std::vector< std::pair<int,int> > cellPartFormat,nodePartFormat;
   if(IsPartStructured(start,end,cgs,cellPartFormat))
     {
-      MEDCouplingAutoRefCountObjectPtr<MEDCouplingStructuredMesh> ret(buildStructuredSubPart(cellPartFormat));
+      MCAuto<MEDCouplingStructuredMesh> ret(buildStructuredSubPart(cellPartFormat));
       nodePartFormat=cellPartFormat;
       for(std::vector< std::pair<int,int> >::iterator it=nodePartFormat.begin();it!=nodePartFormat.end();it++)
         (*it).second++;
-      MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp1(BuildExplicitIdsFrom(getNodeGridStructure(),nodePartFormat));
-      MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp2(DataArrayInt::New()); tmp2->alloc(getNumberOfNodes(),1);
+      MCAuto<DataArrayInt> tmp1(BuildExplicitIdsFrom(getNodeGridStructure(),nodePartFormat));
+      MCAuto<DataArrayInt> tmp2(DataArrayInt::New()); tmp2->alloc(getNumberOfNodes(),1);
       tmp2->fillWithValue(-1);
-      MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp3(DataArrayInt::New()); tmp3->alloc(tmp1->getNumberOfTuples(),1); tmp3->iota(0);
+      MCAuto<DataArrayInt> tmp3(DataArrayInt::New()); tmp3->alloc(tmp1->getNumberOfTuples(),1); tmp3->iota(0);
       tmp2->setPartOfValues3(tmp3,tmp1->begin(),tmp1->end(),0,1,1);
       arr=tmp2.retn();
       return ret.retn();
@@ -631,7 +632,7 @@ DataArrayInt *MEDCouplingStructuredMesh::Build1GTNodalConnectivity(const int *no
   {
     case 0:
       {
-        MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New());
+        MCAuto<DataArrayInt> conn(DataArrayInt::New());
         conn->alloc(1,1); conn->setIJ(0,0,0);
         return conn.retn();
       }
@@ -673,7 +674,7 @@ DataArrayInt *MEDCouplingStructuredMesh::ComputeCornersGhost(const std::vector<i
   if(ghostLev<0)
     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ComputeCornersGhost : ghost lev must be >= 0 !");
   std::size_t dim(st.size());
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New());
+  MCAuto<DataArrayInt> ret(DataArrayInt::New());
   switch(dim)
   {
     case 1:
@@ -962,7 +963,7 @@ std::vector< std::vector<int> > MEDCouplingStructuredMesh::ComputeSignaturePerAx
 DataArrayInt *MEDCouplingStructuredMesh::Build1GTNodalConnectivity1D(const int *nodeStBg)
 {
   int nbOfCells(*nodeStBg-1);
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New());
+  MCAuto<DataArrayInt> conn(DataArrayInt::New());
   conn->alloc(2*nbOfCells,1);
   int *cp=conn->getPointer();
   for(int i=0;i<nbOfCells;i++)
@@ -977,7 +978,7 @@ DataArrayInt *MEDCouplingStructuredMesh::Build1GTNodalConnectivity2D(const int *
 {
   int n1=nodeStBg[0]-1;
   int n2=nodeStBg[1]-1;
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New());
+  MCAuto<DataArrayInt> conn(DataArrayInt::New());
   conn->alloc(4*n1*n2,1);
   int *cp=conn->getPointer();
   int pos=0;
@@ -997,7 +998,7 @@ DataArrayInt *MEDCouplingStructuredMesh::Build1GTNodalConnectivity3D(const int *
   int n1=nodeStBg[0]-1;
   int n2=nodeStBg[1]-1;
   int n3=nodeStBg[2]-1;
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New());
+  MCAuto<DataArrayInt> conn(DataArrayInt::New());
   conn->alloc(8*n1*n2*n3,1);
   int *cp=conn->getPointer();
   int pos=0;
@@ -1023,7 +1024,7 @@ DataArrayInt *MEDCouplingStructuredMesh::Build1GTNodalConnectivityOfSubLevelMesh
   std::vector<int> ngs(3);
   int n0(nodeStBg[0]-1),n1(nodeStBg[1]-1),n2(nodeStBg[2]-1); ngs[0]=n0; ngs[1]=n1; ngs[2]=n2;
   int off0(nodeStBg[0]),off1(nodeStBg[0]*nodeStBg[1]);
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New());
+  MCAuto<DataArrayInt> conn(DataArrayInt::New());
   conn->alloc(4*GetNumberOfCellsOfSubLevelMesh(ngs,3));
   int *cp(conn->getPointer());
   //X
@@ -1188,7 +1189,7 @@ DataArrayInt *MEDCouplingStructuredMesh::Build1GTNodalConnectivityOfSubLevelMesh
   std::vector<int> ngs(2);
   int n0(nodeStBg[0]-1),n1(nodeStBg[1]-1); ngs[0]=n0; ngs[1]=n1;
   int off0(nodeStBg[0]);
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn(DataArrayInt::New());
+  MCAuto<DataArrayInt> conn(DataArrayInt::New());
   conn->alloc(2*GetNumberOfCellsOfSubLevelMesh(ngs,2));
   int *cp(conn->getPointer());
   //X
@@ -1270,9 +1271,55 @@ int MEDCouplingStructuredMesh::getNumberOfNodes() const
   return ret;
 }
 
-void MEDCouplingStructuredMesh::GetPosFromId(int nodeId, int meshDim, const int *split, int *res)
+/*!
+ * This method returns for a cell which id is \a cellId the location (locX,locY,locZ) of this cell in \a this.
+ * 
+ * \param [in] cellId
+ * \return - A vector of size this->getMeshDimension()
+ * \throw if \a cellId not in [ 0, this->getNumberOfCells() )
+ */
+std::vector<int> MEDCouplingStructuredMesh::getLocationFromCellId(int cellId) const
 {
-  int work=nodeId;
+  int meshDim(getMeshDimension());
+  std::vector<int> ret(meshDim);
+  std::vector<int> struc(getCellGridStructure());
+  int nbCells(std::accumulate(struc.begin(),struc.end(),1,std::multiplies<int>()));
+  if(cellId<0 || cellId>=nbCells)
+    {
+      std::ostringstream oss; oss << "MEDCouplingStructuredMesh::getLocationFromCellId : Input cell id (" << cellId << ") is invalid ! Should be in [0," << nbCells << ") !";
+      throw INTERP_KERNEL::Exception(oss.str().c_str());
+    }
+  std::vector<int> spt(GetSplitVectFromStruct(struc));
+  GetPosFromId(cellId,meshDim,&spt[0],&ret[0]);
+  return ret;
+}
+
+/*!
+ * This method returns for a node which id is \a nodeId the location (locX,locY,locZ) of this node in \a this.
+ * 
+ * \param [in] nodeId
+ * \return - A vector of size this->getSpaceDimension()
+ * \throw if \a cellId not in [ 0, this->getNumberOfNodes() )
+ */
+std::vector<int> MEDCouplingStructuredMesh::getLocationFromNodeId(int nodeId) const
+{
+  int spaceDim(getSpaceDimension());
+  std::vector<int> ret(spaceDim);
+  std::vector<int> struc(getNodeGridStructure());
+  int nbNodes(std::accumulate(struc.begin(),struc.end(),1,std::multiplies<int>()));
+  if(nodeId<0 || nodeId>=nbNodes)
+    {
+      std::ostringstream oss; oss << "MEDCouplingStructuredMesh::getLocationFromNodeId : Input node id (" << nodeId << ") is invalid ! Should be in [0," << nbNodes << ") !";
+      throw INTERP_KERNEL::Exception(oss.str().c_str());
+    }
+  std::vector<int> spt(GetSplitVectFromStruct(struc));
+  GetPosFromId(nodeId,spaceDim,&spt[0],&ret[0]);
+  return ret;
+}
+
+void MEDCouplingStructuredMesh::GetPosFromId(int eltId, int meshDim, const int *split, int *res)
+{
+  int work(eltId);
   for(int i=meshDim-1;i>=0;i--)
     {
       int pos=work/split[i];
@@ -1469,6 +1516,8 @@ std::vector< std::pair<int,int> > MEDCouplingStructuredMesh::GetCompactFrmtFromD
 /*!
  * This method returns the intersection zone of two ranges (in compact format) \a r1 and \a r2.
  * This method will throw exception if on one axis the intersection is empty.
+ *
+ * \sa AreRangesIntersect
  */
 std::vector< std::pair<int,int> > MEDCouplingStructuredMesh::IntersectRanges(const std::vector< std::pair<int,int> >& r1, const std::vector< std::pair<int,int> >& r2)
 {
@@ -1499,6 +1548,36 @@ std::vector< std::pair<int,int> > MEDCouplingStructuredMesh::IntersectRanges(con
   return ret;
 }
 
+/*!
+ * This method states if \a r1 and \a r2 do overlap of not. If yes you can call IntersectRanges to know the intersection area.
+ *
+ * \sa IntersectRanges
+ */
+bool MEDCouplingStructuredMesh::AreRangesIntersect(const std::vector< std::pair<int,int> >& r1, const std::vector< std::pair<int,int> >& r2)
+{
+  std::size_t sz(r1.size());
+  if(sz!=r2.size())
+    throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::AreRangesIntersect : the two ranges must have the same dimension !");
+  for(std::size_t i=0;i<sz;i++)
+    {
+      if(r1[i].first>r1[i].second)
+        {
+          std::ostringstream oss; oss << "MEDCouplingStructuredMesh::AreRangesIntersect : On axis " << i << " of range r1, end is before start !";
+          throw INTERP_KERNEL::Exception(oss.str().c_str());
+        }
+      if(r2[i].first>r2[i].second)
+        {
+          std::ostringstream oss; oss << "MEDCouplingStructuredMesh::AreRangesIntersect : On axis " << i << " of range r2, end is before start !";
+          throw INTERP_KERNEL::Exception(oss.str().c_str());
+        }
+      if(r1[i].second<=r2[i].first)
+        return false;
+      if(r1[i].first>=r2[i].second)
+        return false;
+    }
+  return true;
+}
+
 /*!
  * This method is close to BuildExplicitIdsFrom except that instead of returning a DataArrayInt instance containing explicit ids it
  * enable elems in the vector of booleans (for performance reasons). As it is method for performance, this method is \b not
@@ -1634,7 +1713,7 @@ DataArrayDouble *MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom(const std::
     throw INTERP_KERNEL::Exception("MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom : invalid size of input array of double regarding the structure !");
   std::vector<int> dims(GetDimensionsFromCompactFrmt(partCompactFormat));
   int nbOfTuplesOfOutField(DeduceNumberOfGivenStructure(dims)),nbComp(fieldOfDbl->getNumberOfComponents());
-  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfTuplesOfOutField,nbComp);
+  MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfTuplesOfOutField,nbComp);
   ret->copyStringInfoFrom(*fieldOfDbl);
   double *ptRet(ret->getPointer());
   const double *fieldOfDblPtr(fieldOfDbl->begin());
@@ -1676,6 +1755,20 @@ DataArrayDouble *MEDCouplingStructuredMesh::ExtractFieldOfDoubleFrom(const std::
   return ret.retn();
 }
 
+/*!
+ * This method assign a part of values in \a fieldOfDbl using entirely values of \b other.
+ *
+ * \param [in] st - the structure of \a fieldOfDbl.
+ * \param [in,out] fieldOfDbl - the array that will be partially filled using \a other.
+ * \param [in] partCompactFormat - the specification of the part.
+ * \param [in] other - the array that will be used to fill \a fieldOfDbl.
+ */
+void MEDCouplingStructuredMesh::AssignPartOfFieldOfDoubleUsing(const std::vector<int>& st, DataArrayDouble *fieldOfDbl, const std::vector< std::pair<int,int> >& partCompactFormat, const DataArrayDouble *other)
+{//to be optimized
+  std::vector<int> facts(st.size(),1.);
+  MEDCouplingIMesh::CondenseFineToCoarse(st,other,partCompactFormat,facts,fieldOfDbl);
+}
+
 /*!
  * This method changes the reference of a part of structured mesh \a partOfBigInAbs define in absolute reference to a new reference \a bigInAbs.
  * So this method only performs a translation by doing \a partOfBigRelativeToBig = \a partOfBigInAbs - \a bigInAbs
@@ -1816,7 +1909,7 @@ DataArrayInt *MEDCouplingStructuredMesh::BuildExplicitIdsFrom(const std::vector<
       dims[i]=partCompactFormat[i].second-partCompactFormat[i].first;
       nbOfItems*=dims[i];
     }
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret(DataArrayInt::New());
+  MCAuto<DataArrayInt> ret(DataArrayInt::New());
   ret->alloc(nbOfItems,1);
   int *pt(ret->getPointer());
   switch(st.size())