X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FMEDLoader%2FMEDFileBlowStrEltUp.cxx;h=6390279e3ad614dc32391aaa539863bd9ebb3dce;hb=5ea9370443568d9dd6c42c970cbc24c1b6ae6c28;hp=9e6ee4f45ce6b4cd0d7eb7c6993ac4d3a22f98bd;hpb=14f30db13a9749dd4a47f5c14944c7331e7390d1;p=tools%2Fmedcoupling.git

diff --git a/src/MEDLoader/MEDFileBlowStrEltUp.cxx b/src/MEDLoader/MEDFileBlowStrEltUp.cxx
index 9e6ee4f45..6390279e3 100644
--- a/src/MEDLoader/MEDFileBlowStrEltUp.cxx
+++ b/src/MEDLoader/MEDFileBlowStrEltUp.cxx
@@ -21,7 +21,9 @@
 #include "MEDFileBlowStrEltUp.hxx"
 #include "MEDCouplingFieldDouble.hxx"
 #include "MEDFileFieldVisitor.hxx"
+#include "MEDCouplingPartDefinition.hxx"
 #include "MCAuto.txx"
+#include <numeric>
 
 using namespace MEDCoupling;
 
@@ -315,10 +317,12 @@ public:
   bool isClassic() const { return _is_classic; }
   bool operator!=(const FieldWalker2& other) const;
   bool operator==(const FieldWalker2& other) const;
+  const SlicePartDefinition *getPartDef() const { return _pd; }
 private:
   std::string _loc;
   std::string _pfl;
   bool _is_classic;
+  MCAuto<SlicePartDefinition> _pd;
 };
 
 class LocInfo
@@ -329,27 +333,45 @@ public:
   bool operator==(const LocInfo& other) const { return _locs==other._locs && _pfl==other._pfl; }
   void push(const std::string& loc, const std::string& pfl) { checkUniqueLoc(loc); _locs.push_back(loc); _pfl.push_back(pfl); }
   MCAuto<MEDFileUMesh> generateNonClassicalData(int zePos, const MEDFileUMesh *mesh, const MEDFileFieldGlobsReal *globs) const;
+  const PartDefinition *getPartDef() const { return _pd; }
 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 *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);
 public:
   static const char ANGLE_DE_VRILLE[];
+  static const char ANGLE[];
+  static const char SCALE[];
+  static const char EPAISSEUR[];
 private:
   std::vector<std::string> _locs;
   std::vector<std::string> _pfl;
+  MCAuto<PartDefinition> _pd;
 };
 
 const char LocInfo::ANGLE_DE_VRILLE[]="ANGLE DE VRILLE";
 
+const char LocInfo::ANGLE[]="ANGLE";
+
+const char LocInfo::SCALE[]="SCALE";
+
+const char LocInfo::EPAISSEUR[]="EPAISSEUR";
+
 LocInfo::LocInfo(const std::vector<FieldWalker2>& fw)
 {
   std::size_t sz(fw.size());
   _locs.resize(sz); _pfl.resize(sz);
+  if(sz>0)
+    _pd=fw[0].getPartDef()->deepCopy();
   for(std::size_t i=0;i<sz;i++)
     {
       _locs[i]=fw[i].getLoc();
       _pfl[i]=fw[i].getPfl();
+      if(i>0)
+        _pd=(*_pd)+(*(fw[i].getPartDef()));
     }
 }
 
@@ -362,7 +384,7 @@ void LocInfo::checkUniqueLoc(const std::string& loc) const
     }
 }
 
-MCAuto<DataArrayDouble> LocInfo::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)
+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());
   conn=conn->deepCopy(); conn->rearrange(1);
@@ -374,27 +396,39 @@ MCAuto<DataArrayDouble> LocInfo::BuildMeshFromAngleVrille(INTERP_KERNEL::Normali
     geoMesh=umesh->buildUnstructured();
   }
   //
-  MCConstAuto<DataArrayDouble> angleVrille;
   if(!pfl.empty())
     {
       const DataArrayInt *pflArr(globs->getProfile(pfl));
       geoMesh=geoMesh->buildPartOfMySelf(pflArr->begin(),pflArr->end(),true);
-      angleVrille=angleDeVrille->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
     }
-  else
-    angleVrille.takeRef(angleDeVrille);
   //
   MCAuto<MEDCouplingFieldDouble> fakeF(MEDCouplingFieldDouble::New(ON_GAUSS_PT));
   fakeF->setMesh(geoMesh);
-  int nbg(loc.getGaussWeights().size());
   fakeF->setGaussLocalizationOnType(gt,loc.getRefCoords(),loc.getGaussCoords(),loc.getGaussWeights());
-  MCAuto<DataArrayDouble> ptsForLoc(fakeF->getLocalizationOfDiscr());
+  ptsForLoc=fakeF->getLocalizationOfDiscr();
+  //
+  return geoMesh;
+}
+
+MCAuto<DataArrayDouble> LocInfo::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)
+{
+  MCAuto<DataArrayDouble> ptsForLoc;
+  MCAuto<MEDCouplingUMesh> geoMesh(BuildMeshCommon(gt,pfl,loc,zeStr,mesh,section,globs,ptsForLoc));
+  //
+  MCConstAuto<DataArrayDouble> angleVrille;
+  if(!pfl.empty())
+    {
+      const DataArrayInt *pflArr(globs->getProfile(pfl));
+      angleVrille=angleDeVrille->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
+    }
+  else
+    angleVrille.takeRef(angleDeVrille);
   //
   MCAuto<MEDCouplingFieldDouble> dir(geoMesh->buildDirectionVectorField());
   MCAuto<DataArrayDouble> rot(dir->getArray()->fromCartToSpher());
-  int nbCells(geoMesh->getNumberOfCells()),nbCompo(ptsForLoc->getNumberOfComponents());
+  int nbCompo(ptsForLoc->getNumberOfComponents());
   MCAuto<DataArrayDouble> secPts(section->getCoords()->changeNbOfComponents(nbCompo,0.));
-  int nbSecPts(secPts->getNumberOfTuples());
+  int nbSecPts(secPts->getNumberOfTuples()),nbCells(geoMesh->getNumberOfCells()),nbg(loc.getGaussWeights().size());
   {
     const int TAB[3]={2,0,1};
     std::vector<int> v(TAB,TAB+3);
@@ -425,21 +459,168 @@ MCAuto<DataArrayDouble> LocInfo::BuildMeshFromAngleVrille(INTERP_KERNEL::Normali
   return resu;
 }
 
+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)
+{
+  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> zeThickness;
+  if(!pfl.empty())
+    {
+      const DataArrayInt *pflArr(globs->getProfile(pfl));
+      zeThickness=thickness->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
+    }
+  else
+    zeThickness.takeRef(thickness);
+  MCAuto<DataArrayDouble> orthoArr;
+  {
+    MCAuto<MEDCouplingFieldDouble> ortho(geoMesh->buildOrthogonalField());
+    orthoArr.takeRef(ortho->getArray());
+  }
+  int nbCompo(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++)
+    {
+      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());
+      std::for_each(fact2->rwBegin(),fact2->rwEnd(),[eccentricity](double& val){ val*=eccentricity; });
+      std::transform(fact->begin(),fact->end(),fact->getPointer(),[thck](const double& val){ return val*thck/2.; });
+      MCAuto<DataArrayDouble> p(DataArrayDouble::Multiply(secPts,fact));
+      p=DataArrayDouble::Add(p,fact2);
+      for(int l=0;l<nbg;l++)
+        {
+          MCAuto<DataArrayDouble> p2(p->deepCopy());
+          for(int k=0;k<nbCompo;k++)
+            p2->applyLin(1.,ptsForLoc->getIJ(j*nbg+l,k),k);
+          arrs[j*nbg+l]=p2;
+        }
+    }
+  std::vector<const DataArrayDouble *> arrs2(VecAutoToVecOfCstPt(arrs));
+  MCAuto<DataArrayDouble> resu(DataArrayDouble::Aggregate(arrs2));
+  return resu;
+}
+
+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)
+{
+  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());
+  MCConstAuto<DataArrayDouble> zeAngle,zeScale;
+  if(!pfl.empty())
+    {
+      const DataArrayInt *pflArr(globs->getProfile(pfl));
+      zeAngle=angle->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
+      zeScale=scale->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
+    }
+  else
+    {
+      zeAngle.takeRef(angle);
+      zeScale.takeRef(scale);
+    }
+  if(zeAngle->getNumberOfComponents()!=3 || zeScale->getNumberOfComponents()!=2 || nbg!=3)
+    throw INTERP_KERNEL::Exception(MSG1);
+  MCAuto<MEDCouplingFieldDouble> dir;
+  {
+    MCAuto<MEDCouplingUMesh> geoMesh2(geoMesh->deepCopy());
+    geoMesh2->convertQuadraticCellsToLinear();
+    dir=geoMesh2->buildDirectionVectorField();
+  }
+  MCAuto<DataArrayDouble> rot(dir->getArray()->fromCartToSpher());
+  int 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);
+    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++)
+    {
+      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));
+      DataArrayDouble::Rotate3DAlg(CENTER,AX1,-ang1,nbSecPts,p->begin(),p->getPointer());
+      for(int l=0;l<3;l++)
+        {
+          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++)
+            p2->applyLin(1.,ptsForLoc->getIJ(j*nbg+l,k),k);
+          arrs[j*nbg+l]=p2;
+        }
+    }
+  std::vector<const DataArrayDouble *> arrs2(VecAutoToVecOfCstPt(arrs));
+  MCAuto<DataArrayDouble> resu(DataArrayDouble::Aggregate(arrs2));
+  return resu;
+}
+
 MCAuto<DataArrayDouble> LocInfo::BuildMeshFromStructure(INTERP_KERNEL::NormalizedCellType gt, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs)
 {
   static const char MSG1[]="BuildMeshFromStructure : not recognized pattern ! Send mail to anthony.geay@edf.fr with corresponding MED file !";
   const std::vector< MCAuto<DataArray> >& vars(zeStr->getVars());
-  if(vars.size()!=1)
-    throw INTERP_KERNEL::Exception(MSG1);
-  MCAuto<DataArray> zeArr(vars[0]);
-  if(zeArr.isNull())
-    throw INTERP_KERNEL::Exception(MSG1);
-  MCAuto<DataArrayDouble> zeArr2(DynamicCast<DataArray,DataArrayDouble>(zeArr));
-  if(zeArr2.isNull())
-    throw INTERP_KERNEL::Exception(MSG1);
-  if(zeArr2->getName()!=ANGLE_DE_VRILLE)
-    throw INTERP_KERNEL::Exception(MSG1);
-  return BuildMeshFromAngleVrille(gt,zeArr2,pfl,loc,zeStr,mesh,section,globs);
+  if(vars.size()==1)
+    {
+      MCAuto<DataArray> zeArr(vars[0]);
+      if(zeArr.isNull())
+        throw INTERP_KERNEL::Exception(MSG1);
+      MCAuto<DataArrayDouble> zeArr2(DynamicCast<DataArray,DataArrayDouble>(zeArr));
+      if(zeArr2.isNull())
+        throw INTERP_KERNEL::Exception(MSG1);
+      if(zeArr2->getName()==ANGLE_DE_VRILLE || zeArr2->getName()==ANGLE)
+        return BuildMeshFromAngleVrille(gt,zeArr2,pfl,loc,zeStr,mesh,section,globs);
+      if(zeArr2->getName()==EPAISSEUR || zeArr2->getName()==SCALE)
+        return BuildMeshFromEpaisseur(gt,zeArr2,pfl,loc,zeStr,mesh,section,globs);
+    }
+  if(vars.size()==2)
+    {
+      MCAuto<DataArray> zeArr0(vars[0]),zeArr1(vars[1]);
+      if(zeArr0.isNull() || zeArr1.isNull())
+        throw INTERP_KERNEL::Exception(MSG1);
+      MCAuto<DataArrayDouble> zeArr00(DynamicCastSafe<DataArray,DataArrayDouble>(zeArr0)),zeArr11(DynamicCastSafe<DataArray,DataArrayDouble>(zeArr1));
+      switch(gt)
+      {
+        case INTERP_KERNEL::NORM_SEG3:
+          {
+            MCAuto<DataArrayDouble> angle,scale;
+            if(zeArr00->getName()==ANGLE)
+              angle=zeArr00;
+            if(zeArr00->getName()==SCALE)
+              scale=zeArr00;
+            if(zeArr11->getName()==ANGLE)
+              angle=zeArr11;
+            if(zeArr11->getName()==SCALE)
+              scale=zeArr11;
+            if(angle.isNull() || scale.isNull())
+              throw INTERP_KERNEL::Exception(MSG1);
+            return BuildMeshPipeSEG3(angle,scale,pfl,loc,zeStr,mesh,section,globs);
+          }
+        default:
+          throw INTERP_KERNEL::Exception(MSG1);
+      }
+    }
+  throw INTERP_KERNEL::Exception(MSG1);
 }
 
 MCAuto<MEDFileUMesh> LocInfo::generateNonClassicalData(int zePos, const MEDFileUMesh *mesh, const MEDFileFieldGlobsReal *globs) const
@@ -505,17 +686,22 @@ FieldWalker2::FieldWalker2(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd)
   _loc=pmptpd->getLocalization();
   _pfl=pmptpd->getProfile();
   _is_classic=pmptpd->getType()!=ON_GAUSS_PT;
+  _pd=SlicePartDefinition::New(pmptpd->getStart(),pmptpd->getEnd(),1);
 }
 
 bool FieldWalker2::operator!=(const FieldWalker2& other) const
 {
-  bool ret(_loc==other._loc && _pfl==other._pfl && _is_classic==other._is_classic);
-  return !ret;
+  return !((*this)==other);
 }
 
 bool FieldWalker2::operator==(const FieldWalker2& other) const
 {
-  bool ret(_loc==other._loc && _pfl==other._pfl && _is_classic==other._is_classic);
+  bool ret2(false);
+  {
+    std::string tmp;
+    ret2=_pd->isEqual(other._pd,tmp);
+  }
+  bool ret(_loc==other._loc && _pfl==other._pfl && _is_classic==other._is_classic && ret2);
   return ret;
 }
 
@@ -755,7 +941,10 @@ void LocSpliter::generateNonClassicalData(const MEDFileUMesh *mesh, std::vector<
               int t2,t3;
               double t1(f1ts->getTime(t2,t3));
               MCAuto<MEDCouplingFieldDouble> mcf(MEDCouplingFieldDouble::New(ON_NODES));
-              mcf->setArray(f1ts->getUndergroundDataArray());
+              MCAuto<DataArrayDouble> arr,arr2;
+              arr.takeRef(f1ts->getUndergroundDataArray());
+              arr2=arr->selectPartDef((*it).getPartDef());
+              mcf->setArray(arr2);
               mcf->setTime(t1,t2,t3);
               mcf->setName(f1ts->getName());
               mcf->setMesh(mcm);
@@ -832,7 +1021,7 @@ MCAuto<MEDFileFields> MEDFileBlowStrEltUp::splitFieldsPerLoc(const MEDFileFields
       for(int j=0;j<(*it)->getNumberOfFields();j++)
         {
           MCAuto<MEDFileAnyTypeFieldMultiTS> fmts((*it)->getFieldAtPos(j));
-          DealWithConflictNames(fmts,allZeOutFields);
+          //DealWithConflictNames(fmts,allZeOutFields);// uncomment to have a writable data structure
           allZeOutFields->pushField(fmts);
         }
     }