[EDF23738] : multi geometric type in a single mesh having structure elements on each...

[tools/medcoupling.git] / src / MEDLoader / MEDFileBlowStrEltUp.cxx

diff --git a/src/MEDLoader/MEDFileBlowStrEltUp.cxx b/src/MEDLoader/MEDFileBlowStrEltUp.cxx
index c2403ff3e6f9494b16c1e1df49d3d96541740ac4..a94a7d7609d1a3c33fec4de8b333762bcb05da70 100644 (file)
--- a/src/MEDLoader/MEDFileBlowStrEltUp.cxx
+++ b/src/MEDLoader/MEDFileBlowStrEltUp.cxx
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2017  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2021  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
@@ -23,6 +23,7 @@
 #include "MEDFileFieldVisitor.hxx"
 #include "MEDCouplingPartDefinition.hxx"
 #include "MCAuto.txx"
+#include <numeric>
 
 using namespace MEDCoupling;
 
@@ -94,36 +95,35 @@ MCAuto<MEDFileEltStruct4Mesh> MEDFileBlowStrEltUp::dealWithMEDBALLInMesh(const M
   const DataArrayDouble *coo(mesh->getCoords());
   if(!coo)
     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLInMesh : null coords !");
-  MCAuto<DataArrayInt> conn(zeStr->getConn());
+  MCAuto<DataArrayIdType> conn(zeStr->getConn());
   if(conn.isNull())
     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLInMesh : null connectivity !");
   conn->checkAllocated();
   if(conn->getNumberOfComponents()!=1)
     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLInMesh : excepted to be single compo !");
-  int nbCells(conn->getNumberOfTuples());
   MCAuto<DataArrayDouble> connOut(coo->selectByTupleIdSafe(conn->begin(),conn->end()));
   MCAuto<MEDCouplingUMesh> mcOut(MEDCouplingUMesh::Build0DMeshFromCoords(connOut));
   mcOut->setName(BuildNewMeshName(mesh->getName(),MED_BALL_STR));
   mOut->setMeshAtLevel(0,mcOut);
-  const DataArrayInt *ff1(mesh->getFamilyFieldAtLevel(1));
+  const DataArrayIdType *ff1(mesh->getFamilyFieldAtLevel(1));
   if(ff1)
     {
-      MCAuto<DataArrayInt> ff1o(ff1->selectByTupleIdSafe(conn->begin(),conn->end()));
+      MCAuto<DataArrayIdType> ff1o(ff1->selectByTupleIdSafe(conn->begin(),conn->end()));
       mOut->setFamilyFieldArr(1,ff1o);
     }
-  const DataArrayInt *nf1(mesh->getNumberFieldAtLevel(1));
+  const DataArrayIdType *nf1(mesh->getNumberFieldAtLevel(1));
   if(nf1)
     {
-      MCAuto<DataArrayInt> nf1o(nf1->selectByTupleIdSafe(conn->begin(),conn->end()));
+      MCAuto<DataArrayIdType> nf1o(nf1->selectByTupleIdSafe(conn->begin(),conn->end()));
       mOut->setRenumFieldArr(1,nf1o);
     }
   MCAuto<MEDFileUMeshPerTypeCommon> md(zeStr->getMeshDef());
-  const DataArrayInt *ff0(md->getFam());
+  const DataArrayIdType *ff0(md->getFam());
   if(ff0)
-    mOut->setFamilyFieldArr(0,const_cast<DataArrayInt *>(ff0));
-  const DataArrayInt *nf0(md->getNum());
+    mOut->setFamilyFieldArr(0,const_cast<DataArrayIdType *>(ff0));
+  const DataArrayIdType *nf0(md->getNum());
   if(nf0)
-    mOut->setRenumFieldArr(0,const_cast<DataArrayInt *>(nf0));
+    mOut->setRenumFieldArr(0,const_cast<DataArrayIdType *>(nf0));
   mOut->copyFamGrpMapsFrom(*mesh);
   const std::vector< MCAuto<DataArray> >& vars(zeStr->getVars());
   for(std::vector< MCAuto<DataArray> >::const_iterator it=vars.begin();it!=vars.end();it++)
@@ -222,7 +222,7 @@ void MEDFileBlowStrEltUp::dealWithMEDBALLSInFields(const MEDFileFields *fs, cons
           std::vector<std::string> pfls(zeGuideForPfl->getPflsReallyUsed());
           if(pfls.size()>=2)
             throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLSInFields : drink less coffee");
-          MCAuto<DataArrayInt> pflMyLove;
+          MCAuto<DataArrayIdType> pflMyLove;
           if(pfls.size()==1)
             pflMyLove.takeRef(zeGuideForPfl->getProfile(pfls[0]));
           // Yeah we have pfls
@@ -336,7 +336,7 @@ public:
 private:
   void checkUniqueLoc(const std::string& loc) const;
   static MCAuto<DataArrayDouble> BuildMeshFromAngleVrille(INTERP_KERNEL::NormalizedCellType gt, const DataArrayDouble *angleDeVrille, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs);
-  static MCAuto<DataArrayDouble> BuildMeshFromEpaisseur(INTERP_KERNEL::NormalizedCellType gt, const DataArrayDouble *thikness, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs);
+  static MCAuto<DataArrayDouble> BuildMeshFromEpaisseur(INTERP_KERNEL::NormalizedCellType gt, const DataArrayDouble *thickness, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs);
   static MCAuto<DataArrayDouble> BuildMeshPipeSEG3(const DataArrayDouble *angle, const DataArrayDouble *scale, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs);
   static MCAuto<MEDCouplingUMesh> BuildMeshCommon(INTERP_KERNEL::NormalizedCellType gt, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs, MCAuto<DataArrayDouble>& ptsForLoc);
   static MCAuto<DataArrayDouble> BuildMeshFromStructure(INTERP_KERNEL::NormalizedCellType gt, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs);
@@ -385,7 +385,7 @@ void LocInfo::checkUniqueLoc(const std::string& loc) const
 
 MCAuto<MEDCouplingUMesh> LocInfo::BuildMeshCommon(INTERP_KERNEL::NormalizedCellType gt, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs, MCAuto<DataArrayDouble>& ptsForLoc)
 {
-  MCAuto<DataArrayInt> conn(zeStr->getConn());
+  MCAuto<DataArrayIdType> conn(zeStr->getConn());
   conn=conn->deepCopy(); conn->rearrange(1);
   MCAuto<MEDCouplingUMesh> geoMesh;
   {
@@ -397,7 +397,7 @@ MCAuto<MEDCouplingUMesh> LocInfo::BuildMeshCommon(INTERP_KERNEL::NormalizedCellT
   //
   if(!pfl.empty())
     {
-      const DataArrayInt *pflArr(globs->getProfile(pfl));
+      const DataArrayIdType *pflArr(globs->getProfile(pfl));
       geoMesh=geoMesh->buildPartOfMySelf(pflArr->begin(),pflArr->end(),true);
     }
   //
@@ -417,7 +417,7 @@ MCAuto<DataArrayDouble> LocInfo::BuildMeshFromAngleVrille(INTERP_KERNEL::Normali
   MCConstAuto<DataArrayDouble> angleVrille;
   if(!pfl.empty())
     {
-      const DataArrayInt *pflArr(globs->getProfile(pfl));
+      const DataArrayIdType *pflArr(globs->getProfile(pfl));
       angleVrille=angleDeVrille->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
     }
   else
@@ -425,12 +425,12 @@ MCAuto<DataArrayDouble> LocInfo::BuildMeshFromAngleVrille(INTERP_KERNEL::Normali
   //
   MCAuto<MEDCouplingFieldDouble> dir(geoMesh->buildDirectionVectorField());
   MCAuto<DataArrayDouble> rot(dir->getArray()->fromCartToSpher());
-  int nbCompo(ptsForLoc->getNumberOfComponents());
+  std::size_t nbCompo(ptsForLoc->getNumberOfComponents());
   MCAuto<DataArrayDouble> secPts(section->getCoords()->changeNbOfComponents(nbCompo,0.));
-  int nbSecPts(secPts->getNumberOfTuples()),nbCells(geoMesh->getNumberOfCells()),nbg(loc.getGaussWeights().size());
+  mcIdType nbSecPts(secPts->getNumberOfTuples()),nbCells(geoMesh->getNumberOfCells()),nbg(ToIdType(loc.getGaussWeights().size()));
   {
     const int TAB[3]={2,0,1};
-    std::vector<int> v(TAB,TAB+3);
+    std::vector<std::size_t> v(TAB,TAB+3);
     secPts=secPts->keepSelectedComponents(v);
   }
   const double CENTER[3]={0.,0.,0.},AX0[3]={0.,0.,1.};
@@ -448,7 +448,7 @@ MCAuto<DataArrayDouble> LocInfo::BuildMeshFromAngleVrille(INTERP_KERNEL::Normali
       for(int l=0;l<nbg;l++)
         {
           MCAuto<DataArrayDouble> p2(p->deepCopy());
-          for(int k=0;k<nbCompo;k++)
+          for(std::size_t k=0;k<nbCompo;k++)
             p2->applyLin(1.,ptsForLoc->getIJ(j*nbg+l,k),k);
           arrs[j*nbg+l]=p2;
         }
@@ -458,32 +458,31 @@ MCAuto<DataArrayDouble> LocInfo::BuildMeshFromAngleVrille(INTERP_KERNEL::Normali
   return resu;
 }
 
-MCAuto<DataArrayDouble> LocInfo::BuildMeshFromEpaisseur(INTERP_KERNEL::NormalizedCellType gt, const DataArrayDouble *thikness, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs)
+MCAuto<DataArrayDouble> LocInfo::BuildMeshFromEpaisseur(INTERP_KERNEL::NormalizedCellType gt, const DataArrayDouble *thickness, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs)
 {
-#if __cplusplus >= 201103L
   MCAuto<DataArrayDouble> ptsForLoc;
   MCAuto<MEDCouplingUMesh> geoMesh(BuildMeshCommon(gt,pfl,loc,zeStr,mesh,section,globs,ptsForLoc));
-  int nbSecPts(section->getNumberOfNodes()),nbCells(geoMesh->getNumberOfCells()),nbg(loc.getGaussWeights().size());
-  MCConstAuto<DataArrayDouble> zeThikness;
+  mcIdType nbCells(geoMesh->getNumberOfCells()),nbg(ToIdType(loc.getGaussWeights().size()));
+  MCConstAuto<DataArrayDouble> zeThickness;
   if(!pfl.empty())
     {
-      const DataArrayInt *pflArr(globs->getProfile(pfl));
-      zeThikness=thikness->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
+      const DataArrayIdType *pflArr(globs->getProfile(pfl));
+      zeThickness=thickness->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
     }
   else
-    zeThikness.takeRef(thikness);
+    zeThickness.takeRef(thickness);
   MCAuto<DataArrayDouble> orthoArr;
   {
     MCAuto<MEDCouplingFieldDouble> ortho(geoMesh->buildOrthogonalField());
     orthoArr.takeRef(ortho->getArray());
   }
-  int nbCompo(orthoArr->getNumberOfComponents());
+  mcIdType nbCompo(ToIdType(orthoArr->getNumberOfComponents()));
   MCAuto<DataArrayDouble> secPts(section->getCoords()->duplicateEachTupleNTimes(nbCompo));
   secPts->rearrange(nbCompo);
   std::vector< MCAuto<DataArrayDouble> > arrs(nbCells*nbg);
-  for(int j=0;j<nbCells;j++)
+  for(mcIdType j=0;j<nbCells;j++)
     {
-      double thck(zeThikness->getIJ(j,0)),eccentricity(zeThikness->getIJ(j,1));
+      double thck(zeThickness->getIJ(j,0)),eccentricity(zeThickness->getIJ(j,1));
       MCAuto<DataArrayDouble> fact(DataArrayDouble::New()),fact2(DataArrayDouble::New()); fact->alloc(1,nbCompo); fact2->alloc(1,nbCompo);
       std::copy(orthoArr->begin()+j*nbCompo,orthoArr->begin()+(j+1)*nbCompo,fact->getPointer());
       std::copy(orthoArr->begin()+j*nbCompo,orthoArr->begin()+(j+1)*nbCompo,fact2->getPointer());
@@ -502,9 +501,6 @@ MCAuto<DataArrayDouble> LocInfo::BuildMeshFromEpaisseur(INTERP_KERNEL::Normalize
   std::vector<const DataArrayDouble *> arrs2(VecAutoToVecOfCstPt(arrs));
   MCAuto<DataArrayDouble> resu(DataArrayDouble::Aggregate(arrs2));
   return resu;
-#else
-  throw INTERP_KERNEL::Exception("Broken news : 10% off for C++11 compiler :)");
-#endif
 }
 
 MCAuto<DataArrayDouble> LocInfo::BuildMeshPipeSEG3(const DataArrayDouble *angle, const DataArrayDouble *scale, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs)
@@ -512,11 +508,11 @@ MCAuto<DataArrayDouble> LocInfo::BuildMeshPipeSEG3(const DataArrayDouble *angle,
   static const char MSG1[]="BuildMeshPipeSEG3 : not recognized pattern ! Send mail to anthony.geay@edf.fr with corresponding MED file !";
   MCAuto<DataArrayDouble> ptsForLoc;
   MCAuto<MEDCouplingUMesh> geoMesh(BuildMeshCommon(INTERP_KERNEL::NORM_SEG3,pfl,loc,zeStr,mesh,section,globs,ptsForLoc));
-  int nbSecPts(section->getNumberOfNodes()),nbCells(geoMesh->getNumberOfCells()),nbg(loc.getGaussWeights().size());
+  mcIdType nbSecPts(section->getNumberOfNodes()),nbCells(geoMesh->getNumberOfCells()),nbg(ToIdType(loc.getGaussWeights().size()));
   MCConstAuto<DataArrayDouble> zeAngle,zeScale;
   if(!pfl.empty())
     {
-      const DataArrayInt *pflArr(globs->getProfile(pfl));
+      const DataArrayIdType *pflArr(globs->getProfile(pfl));
       zeAngle=angle->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
       zeScale=scale->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
     }
@@ -534,20 +530,32 @@ MCAuto<DataArrayDouble> LocInfo::BuildMeshPipeSEG3(const DataArrayDouble *angle,
     dir=geoMesh2->buildDirectionVectorField();
   }
   MCAuto<DataArrayDouble> rot(dir->getArray()->fromCartToSpher());
-  int nbCompo(ptsForLoc->getNumberOfComponents());
+  std::size_t nbCompo(ptsForLoc->getNumberOfComponents());
   MCAuto<DataArrayDouble> secPts(section->getCoords()->changeNbOfComponents(nbCompo,0.));
   {
     const int TAB[3]={2,0,1};
-    std::vector<int> v(TAB,TAB+3);
+    std::vector<std::size_t> v(TAB,TAB+3);
     secPts=secPts->keepSelectedComponents(v);
   }
   const double CENTER[3]={0.,0.,0.},AX0[3]={0.,0.,1.};
   double AX1[3]; AX1[2]=0.;
   std::vector< MCAuto<DataArrayDouble> > arrs(nbCells*nbg);
-  for(int j=0;j<nbCells;j++)
+  for(mcIdType j=0;j<nbCells;j++)
     {
+      constexpr int DIM=3;
       MCAuto<DataArrayDouble> p(secPts->deepCopy());
       double ang0(rot->getIJ(j,2));
+      double rmin(zeScale->getIJ(j,0)),rmax(zeScale->getIJ(j,1));
+      {
+        auto pt(p->rwBegin());
+        for(int i=0;i<nbSecPts;i++)
+          {
+            auto nrm(sqrt(std::accumulate(pt,pt+DIM,0.,[](double sum, double v) { return sum+v*v; } )));
+            auto sca((rmin+2.*(nrm-0.5)*(rmax-rmin))/nrm);
+            std::for_each(pt,pt+3,[sca](double& val) { val*=sca; } );
+            std::advance(pt,DIM);
+          }
+      }
       DataArrayDouble::Rotate3DAlg(CENTER,AX0,ang0,nbSecPts,p->begin(),p->getPointer());
       AX1[0]=-sin(ang0); AX1[1]=cos(ang0);// rot Oy around OZ
       double ang1(M_PI/2.-rot->getIJ(j,1));
@@ -557,7 +565,7 @@ MCAuto<DataArrayDouble> LocInfo::BuildMeshPipeSEG3(const DataArrayDouble *angle,
           MCAuto<DataArrayDouble> p3(p->deepCopy());
           DataArrayDouble::Rotate3DAlg(CENTER,dir->getArray()->begin()+j*3,zeAngle->getIJ(j,l),nbSecPts,p3->begin(),p3->getPointer());
           MCAuto<DataArrayDouble> p2(p3->deepCopy());
-          for(int k=0;k<nbCompo;k++)
+          for(std::size_t k=0;k<nbCompo;k++)
             p2->applyLin(1.,ptsForLoc->getIJ(j*nbg+l,k),k);
           arrs[j*nbg+l]=p2;
         }
@@ -639,11 +647,11 @@ MCAuto<MEDFileUMesh> LocInfo::generateNonClassicalData(int zePos, const MEDFileU
         if(um->getNumberOfCells()!=1)
           throw INTERP_KERNEL::Exception(MSG1);
         gt=um->getTypeOfCell(0);
-        std::vector<int> v;
+        std::vector<mcIdType> v;
         um->getNodeIdsOfCell(0,v);
         std::size_t sz2(v.size());
         for(std::size_t j=0;j<sz2;j++)
-          if(v[j]!=j)
+          if(v[j]!=ToIdType(j))
             throw INTERP_KERNEL::Exception(MSG1);
       }
       const std::vector< MCAuto<MEDFileEltStruct4Mesh> >& strs(mesh->getAccessOfUndergroundEltStrs());