src/MEDLoader/MEDFileBlowStrEltUp.cxx

   1 // Copyright (C) 2007-2023  CEA, EDF
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (EDF R&D)
  20
  21 #include "MEDFileBlowStrEltUp.hxx"
  22 #include "MEDCouplingFieldDouble.hxx"
  23 #include "MEDFileFieldVisitor.hxx"
  24 #include "MEDCouplingPartDefinition.hxx"
  25 #include "MCAuto.txx"
  26 #include <numeric>
  27
  28 using namespace MEDCoupling;
  29
  30 const char MEDFileBlowStrEltUp::MED_BALL_STR[]="MED_BALL";
  31
  32 MEDFileBlowStrEltUp::MEDFileBlowStrEltUp(const MEDFileFields *fsOnlyOnSE, const MEDFileMeshes *ms, const MEDFileStructureElements *ses)
  33 {
  34   if(!fsOnlyOnSE || !ms || !ses)
  35     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp constructor : NULL input pointer !");
  36   _ms.takeRef(ms); _ses.takeRef(ses);
  37   std::vector< std::pair<std::string,std::string> > ps;
  38   fsOnlyOnSE->getMeshSENames(ps);
  39   std::size_t sz(ps.size());
  40   _elts.resize(sz);
  41   for(std::size_t i=0;i<sz;i++)
  42     {
  43       const std::pair<std::string,std::string>& p(ps[i]);
  44       MCAuto<MEDFileFields> f(fsOnlyOnSE->partOfThisLyingOnSpecifiedMeshSEName(p.first,p.second));
  45       _elts[i]=f;
  46     }
  47   for(std::size_t i=0;i<sz;i++)
  48     {
  49       const std::pair<std::string,std::string>& p(ps[i]);
  50       MEDFileMesh *mesh(_ms->getMeshWithName(p.first));
  51       if(!mesh)
  52         throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp : NULL mesh !");
  53       MEDFileUMesh *umesh(dynamic_cast<MEDFileUMesh *>(mesh));
  54       if(!umesh)
  55         throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp : Blow up of Stru Elt not managed yet for unstructured meshes !");
  56     }
  57 }
  58
  59 /*!
  60  * \param [in] mesh - The mesh containing structure element called \a seName. After the call of this method the Structure elements parts will be removed.
  61  * \param [out] mOut - the physical mesh of the structure element \a seName in mesh \a mesh
  62  * \param [out] fsOut - the list of var attribute of structure element \a seName - \b WARNING no time steps here
  63  */
  64 MCAuto<MEDFileEltStruct4Mesh> MEDFileBlowStrEltUp::dealWithSEInMesh(const std::string& seName, MEDFileUMesh *mesh, MCAuto<MEDFileUMesh>& mOut, MCAuto<MEDFileFields>& fsOut) const
  65 {
  66   if(!mesh)
  67     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithSEInMesh : null pointer !");
  68   if(seName==MED_BALL_STR)
  69     {
  70       MCAuto<MEDFileEltStruct4Mesh> ret(dealWithMEDBALLInMesh(mesh,mOut,fsOut));
  71       mesh->killStructureElements();
  72       return ret;
  73     }
  74   throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithSEInMesh : only MED_BALL is managed for the moment, but if you are interested please send spec to anthony.geay@edf.fr !");
  75 }
  76
  77 MCAuto<MEDFileEltStruct4Mesh> MEDFileBlowStrEltUp::dealWithMEDBALLInMesh(const MEDFileUMesh *mesh, MCAuto<MEDFileUMesh>& mOut, MCAuto<MEDFileFields>& fsOut) const
  78 {
  79   mOut=MEDFileUMesh::New(); fsOut=MEDFileFields::New();
  80   const std::vector< MCAuto<MEDFileEltStruct4Mesh> >& strs(mesh->getAccessOfUndergroundEltStrs());
  81   MCAuto<MEDFileEltStruct4Mesh> zeStr;
  82   for(std::vector< MCAuto<MEDFileEltStruct4Mesh> >::const_iterator it=strs.begin();it!=strs.end();it++)
  83     {
  84       if((*it)->getGeoTypeName()==MED_BALL_STR)
  85         {
  86           zeStr=*it;
  87           break;
  88         }
  89     }
  90   if(zeStr.isNull())
  91     {
  92       std::ostringstream oss; oss << "MEDFileBlowStrEltUp::dealWithMEDBALLInMesh : no geo type with name " <<  MED_BALL_STR << " in " << mesh->getName() << " !";
  93       throw INTERP_KERNEL::Exception(oss.str());
  94     }
  95   const DataArrayDouble *coo(mesh->getCoords());
  96   if(!coo)
  97     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLInMesh : null coords !");
  98   MCAuto<DataArrayIdType> conn(zeStr->getConn());
  99   if(conn.isNull())
 100     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLInMesh : null connectivity !");
 101   conn->checkAllocated();
 102   if(conn->getNumberOfComponents()!=1)
 103     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLInMesh : excepted to be single compo !");
 104   MCAuto<DataArrayDouble> connOut(coo->selectByTupleIdSafe(conn->begin(),conn->end()));
 105   MCAuto<MEDCouplingUMesh> mcOut(MEDCouplingUMesh::Build0DMeshFromCoords(connOut));
 106   mcOut->setName(BuildNewMeshName(mesh->getName(),MED_BALL_STR));
 107   mOut->setMeshAtLevel(0,mcOut);
 108   const DataArrayIdType *ff1(mesh->getFamilyFieldAtLevel(1));
 109   if(ff1)
 110     {
 111       MCAuto<DataArrayIdType> ff1o(ff1->selectByTupleIdSafe(conn->begin(),conn->end()));
 112       mOut->setFamilyFieldArr(1,ff1o);
 113     }
 114   const DataArrayIdType *nf1(mesh->getNumberFieldAtLevel(1));
 115   if(nf1)
 116     {
 117       MCAuto<DataArrayIdType> nf1o(nf1->selectByTupleIdSafe(conn->begin(),conn->end()));
 118       mOut->setRenumFieldArr(1,nf1o);
 119     }
 120   MCAuto<MEDFileUMeshPerTypeCommon> md(zeStr->getMeshDef());
 121   const DataArrayIdType *ff0(md->getFam());
 122   if(ff0)
 123     mOut->setFamilyFieldArr(0,const_cast<DataArrayIdType *>(ff0));
 124   const DataArrayIdType *nf0(md->getNum());
 125   if(nf0)
 126     mOut->setRenumFieldArr(0,const_cast<DataArrayIdType *>(nf0));
 127   mOut->copyFamGrpMapsFrom(*mesh);
 128   const std::vector< MCAuto<DataArray> >& vars(zeStr->getVars());
 129   for(std::vector< MCAuto<DataArray> >::const_iterator it=vars.begin();it!=vars.end();it++)
 130     {
 131       const DataArray *elt(*it);
 132       if(!elt)
 133         continue;
 134       {
 135         const DataArrayDouble *eltC(dynamic_cast<const DataArrayDouble *>(elt));
 136         if(eltC)
 137           {
 138             MCAuto<MEDFileFieldMultiTS> fmts(MEDFileFieldMultiTS::New());
 139             MCAuto<MEDFileField1TS> f1ts(MEDFileField1TS::New());
 140             MCAuto<MEDCouplingFieldDouble> f(MEDCouplingFieldDouble::New(ON_NODES));
 141             f->setMesh(mcOut);
 142             f->setArray(const_cast<DataArrayDouble *>(eltC));
 143             f->setName(eltC->getName());
 144             f1ts->setFieldNoProfileSBT(f);
 145             fmts->pushBackTimeStep(f1ts);
 146             fsOut->pushField(fmts);
 147           }
 148       }
 149     }
 150   return zeStr;
 151 }
 152
 153 /*!
 154  * \param [in] fs - fields lying all on same mesh and on same structure element
 155  * \param [in] zeStr - ze structure of current structure element
 156  * \param [in] varAtt - fields containing var att of current structure element. WARNING at this stage the number of iteration are equal to one for each field in \a varAtt
 157  * \param [out] zeOutputs - ze fields that are the concatenation of fields in \a fs transformed and those in \a varAtt normalized in time space
 158  */
 159 void MEDFileBlowStrEltUp::dealWithSEInFields(const std::string& seName, const MEDFileFields *fs, const MEDFileEltStruct4Mesh *zeStr, const MEDFileFields *varAtt, MEDFileFields *zeOutputs) const
 160 {
 161   if(!fs)
 162     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithSEInFields : null pointer !");
 163   if(seName==MED_BALL_STR)
 164     {
 165       dealWithMEDBALLSInFields(fs,zeStr,varAtt,zeOutputs);
 166       return ;
 167     }
 168   throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithSEInFields : only MED_BALL is managed for the moment, but if you are interested please send spec to anthony.geay@edf.fr !");
 169 }
 170
 171 void MEDFileBlowStrEltUp::dealWithMEDBALLSInFields(const MEDFileFields *fs, const MEDFileEltStruct4Mesh *zeStr, const MEDFileFields *varAtt, MEDFileFields *zeOutputs) const
 172 {
 173   int nbf(fs->getNumberOfFields());
 174   std::vector< MCAuto<MEDFileAnyTypeFieldMultiTS> > elts0;
 175   std::vector< MEDFileAnyTypeFieldMultiTS * > elts1;
 176   std::string zeMeshName;
 177   for(int i=0;i<nbf;i++)
 178     {
 179       MCAuto<MEDFileAnyTypeFieldMultiTS> elt(fs->getFieldAtPos(i));
 180       MCAuto<MEDFileAnyTypeFieldMultiTS> eltOut(elt->buildNewEmpty());
 181       int nbTS(elt->getNumberOfTS());
 182       for(int j=0;j<nbTS;j++)
 183         {
 184           MCAuto<MEDFileAnyTypeField1TS> eltt(elt->getTimeStepAtPos(j));
 185           MCAuto<MEDFileAnyTypeField1TS> elttOut(eltt->deepCopy());
 186           std::string meshName(eltt->getMeshName());
 187           zeMeshName=BuildNewMeshName(meshName,MED_BALL_STR);
 188           elttOut->setMeshName(zeMeshName);
 189           elttOut->convertMedBallIntoClassic();
 190           eltOut->pushBackTimeStep(elttOut);
 191         }
 192       elts0.push_back(eltOut); elts1.push_back(eltOut);
 193     }
 194   //
 195   const MEDFileMesh *zeCurrentMesh(_ms->getMeshWithName(zeMeshName));
 196   //
 197   std::size_t ii(0);
 198   std::vector< std::vector<MEDFileAnyTypeFieldMultiTS *> > sp(MEDFileAnyTypeFieldMultiTS::SplitIntoCommonTimeSeries(elts1));
 199   for(std::vector< std::vector<MEDFileAnyTypeFieldMultiTS *> >::const_iterator it0=sp.begin();it0!=sp.end();it0++,ii++)
 200     {
 201       std::vector< MCAuto<MEDFileFastCellSupportComparator> > fsc;
 202       std::vector< std::vector<MEDFileAnyTypeFieldMultiTS *> > sp2(MEDFileAnyTypeFieldMultiTS::SplitPerCommonSupport(*it0,zeCurrentMesh,fsc));
 203       std::size_t jj(0);
 204       for(std::vector< std::vector<MEDFileAnyTypeFieldMultiTS *> >::const_iterator it1=sp2.begin();it1!=sp2.end();it1++,jj++)
 205         {
 206           for(std::vector<MEDFileAnyTypeFieldMultiTS *>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 207             zeOutputs->pushField(*it2);
 208           // The most exciting part. Users that put profiles on struct elements part of fields. Reduce var att.
 209           if((*it1).size()<1)
 210             throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLSInFields : take a deep breath !");
 211           MCAuto<MEDFileAnyTypeField1TS> zeGuideForPfl;// This var is the reference for pfl management.
 212           {
 213             if(!(*it1)[0])
 214               throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLSInFields : take a deep breath 2 !");
 215             int pdm((*it1)[0]->getNumberOfTS());
 216             if(pdm<1)
 217               throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLSInFields : take a deep breath 3 !");
 218             zeGuideForPfl=(*it1)[0]->getTimeStepAtPos(0);
 219           }
 220           if(zeGuideForPfl.isNull())
 221             throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLSInFields : take a deep breath 4 !");
 222           std::vector<std::string> pfls(zeGuideForPfl->getPflsReallyUsed());
 223           if(pfls.size()>=2)
 224             throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLSInFields : drink less coffee");
 225           MCAuto<DataArrayIdType> pflMyLove;
 226           if(pfls.size()==1)
 227             pflMyLove.takeRef(zeGuideForPfl->getProfile(pfls[0]));
 228           // Yeah we have pfls
 229           std::vector<double> t2s;
 230           std::vector< std::pair<int,int> > t1s((*it1)[0]->getTimeSteps(t2s));
 231           std::size_t nbTS3(t2s.size());
 232           int nbf2(varAtt->getNumberOfFields());
 233           for(int i=0;i<nbf2;i++)
 234             {
 235               MCAuto<MEDFileAnyTypeFieldMultiTS> elt(varAtt->getFieldAtPos(i));
 236               int nbTS2(elt->getNumberOfTS());
 237               if(nbTS2!=1)
 238                 throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLSInFields : internal error ! The dealWithMEDBALLInMesh is expected to return a single TS !");
 239               MCAuto<MEDFileAnyTypeField1TS> elt2(elt->getTimeStepAtPos(0));
 240               MCAuto<MEDFileAnyTypeFieldMultiTS> elt4(elt->buildNewEmpty());
 241               for(std::size_t j=0;j<nbTS3;j++)
 242                 {
 243                   MCAuto<MEDFileAnyTypeField1TS> elt3(elt2->deepCopy());
 244                   elt3->setTime(t1s[j].first,t1s[j].second,t2s[j]);
 245                   elt3->setName(BuildVarAttName(ii,sp.size(),jj,sp2.size(),elt3->getName()));
 246                   if(pflMyLove.isNotNull())
 247                     elt3->makeReduction(INTERP_KERNEL::NORM_ERROR,ON_NODES,pflMyLove);
 248                   elt4->pushBackTimeStep(elt3);
 249                 }
 250               zeOutputs->pushField(elt4);
 251             }
 252         }
 253     }
 254 }
 255
 256 void MEDFileBlowStrEltUp::generate(MEDFileMeshes *msOut, MEDFileFields *allZeOutFields)
 257 {
 258   for(std::vector< MCAuto<MEDFileFields> >::iterator elt=_elts.begin();elt!=_elts.end();elt++)
 259     {
 260       std::vector< std::pair<std::string,std::string> > ps;
 261       (*elt)->getMeshSENames(ps);
 262       if(ps.size()!=1)
 263         throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::generateMeshes : internal error !");
 264       MEDFileMesh *mesh(_ms->getMeshWithName(ps[0].first));
 265       if(!mesh)
 266         throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::generateMeshes : NULL mesh !");
 267       MEDFileUMesh *umesh(dynamic_cast<MEDFileUMesh *>(mesh));
 268       if(!umesh)
 269         throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::generateMeshes : Blow up of Stru Elt not managed yet for unstructured meshes !");
 270       //
 271       MCAuto<MEDFileFields> classicalSEFields(splitFieldsPerLoc(*elt,umesh,msOut,allZeOutFields));
 272       if(classicalSEFields.isNotNull())
 273         {
 274           MCAuto<MEDFileUMesh> mOut;
 275           MCAuto<MEDFileFields> fsOut1;
 276           MCAuto<MEDFileEltStruct4Mesh> zeStr(dealWithSEInMesh(ps[0].second,umesh,mOut,fsOut1));
 277           msOut->pushMesh(mOut);
 278           dealWithSEInFields(ps[0].second,classicalSEFields,zeStr,fsOut1,allZeOutFields);
 279         }
 280     }
 281 }
 282
 283 std::string MEDFileBlowStrEltUp::BuildNewMeshName(const std::string& meshName, const std::string& seName)
 284 {
 285   std::ostringstream mNameOut;
 286   mNameOut << meshName << "_" << seName;
 287   return mNameOut.str();
 288 }
 289
 290 std::string MEDFileBlowStrEltUp::BuildVarAttName(std::size_t iPart, std::size_t totINbParts, std::size_t jPart, std::size_t totJNbParts, const std::string& name)
 291 {
 292   if(totINbParts==1 && totJNbParts==1)
 293     return name;
 294   std::ostringstream oss;
 295   oss << name << "@" << iPart << "@" << jPart;
 296   return oss.str();
 297 }
 298
 299 void MEDFileBlowStrEltUp::DealWithSE(MEDFileFields *fs, MEDFileMeshes *ms, const MEDFileStructureElements *ses)
 300 {
 301   MCAuto<MEDFileFields> fsSEOnly(fs->partOfThisOnStructureElements());
 302   fs->killStructureElements();
 303   MEDFileBlowStrEltUp bu(fsSEOnly,ms,ses);
 304   bu.generate(ms,fs);
 305   fs->killStructureElementsInGlobs();
 306 }
 307
 308 //
 309
 310 class FieldWalker2
 311 {
 312 public:
 313   FieldWalker2(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd);
 314   std::string getLoc() const { return _loc; }
 315   std::string getPfl() const { return _pfl; }
 316   INTERP_KERNEL::NormalizedCellType getGeoType() const { return _ct; }
 317   bool isClassic() const { return _is_classic; }
 318   bool operator!=(const FieldWalker2& other) const;
 319   bool operator==(const FieldWalker2& other) const;
 320   const SlicePartDefinition *getPartDef() const { return _pd; }
 321 private:
 322   std::string _loc;
 323   std::string _pfl;
 324   INTERP_KERNEL::NormalizedCellType _ct;
 325   bool _is_classic;
 326   MCAuto<SlicePartDefinition> _pd;
 327 };
 328
 329 class LocInfo
 330 {
 331 public:
 332   LocInfo() { }
 333   LocInfo(const std::vector<FieldWalker2>& fw);
 334   bool operator==(const LocInfo& other) const { return _locs==other._locs && _pfl==other._pfl; }
 335   void push(const std::string& loc, const std::string& pfl) { checkUniqueLoc(loc); _locs.push_back(loc); _pfl.push_back(pfl); }
 336   MCAuto<MEDFileUMesh> generateNonClassicalData(int zePos, const MEDFileUMesh *mesh, const MEDFileFieldGlobsReal *globs) const;
 337   const PartDefinition *getPartDef() const { return _pd; }
 338 private:
 339   void checkUniqueLoc(const std::string& loc) const;
 340   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);
 341   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);
 342   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);
 343   static MCAuto<DataArrayDouble> BuildMeshPipeSEG4(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);
 344   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);
 345   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);
 346 public:
 347   static const char ANGLE_DE_VRILLE[];
 348   static const char ANGLE[];
 349   static const char SCALE[];
 350   static const char EPAISSEUR[];
 351 private:
 352   std::vector<std::string> _locs;
 353   std::vector<std::string> _pfl;
 354   std::vector<INTERP_KERNEL::NormalizedCellType> _cts;
 355   MCAuto<PartDefinition> _pd;
 356 };
 357
 358 const char LocInfo::ANGLE_DE_VRILLE[]="ANGLE DE VRILLE";
 359
 360 const char LocInfo::ANGLE[]="ANGLE";
 361
 362 const char LocInfo::SCALE[]="SCALE";
 363
 364 const char LocInfo::EPAISSEUR[]="EPAISSEUR";
 365
 366 LocInfo::LocInfo(const std::vector<FieldWalker2>& fw)
 367 {
 368   std::size_t sz(fw.size());
 369   _locs.resize(sz); _pfl.resize(sz); _cts.resize(sz);
 370   if(sz>0)
 371     _pd=fw[0].getPartDef()->deepCopy();
 372   for(std::size_t i=0;i<sz;i++)
 373     {
 374       _locs[i]=fw[i].getLoc();
 375       _pfl[i]=fw[i].getPfl();
 376       _cts[i]=fw[i].getGeoType();
 377       if(i>0)
 378         _pd=(*_pd)+(*(fw[i].getPartDef()));
 379     }
 380 }
 381
 382 void LocInfo::checkUniqueLoc(const std::string& loc) const
 383 {
 384   if(std::find(_locs.begin(),_locs.end(),loc)!=_locs.end())
 385     {
 386       std::ostringstream oss; oss << "LocInfo::checkUniqueLoc : loc \"" << loc << "\" already exists !";
 387       throw INTERP_KERNEL::Exception(oss.str());
 388     }
 389 }
 390
 391 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)
 392 {
 393   MCAuto<DataArrayIdType> conn(zeStr->getConn());
 394   conn=conn->deepCopy(); conn->rearrange(1);
 395   MCAuto<MEDCouplingUMesh> geoMesh;
 396   {
 397     MCAuto<MEDCoupling1SGTUMesh> umesh(MEDCoupling1SGTUMesh::New("",gt));
 398     umesh->setCoords(mesh->getCoords());
 399     umesh->setNodalConnectivity(conn);
 400     geoMesh=umesh->buildUnstructured();
 401   }
 402   //
 403   if(!pfl.empty())
 404     {
 405       const DataArrayIdType *pflArr(globs->getProfile(pfl));
 406       geoMesh=geoMesh->buildPartOfMySelf(pflArr->begin(),pflArr->end(),true);
 407     }
 408   //
 409   MCAuto<MEDCouplingFieldDouble> fakeF(MEDCouplingFieldDouble::New(ON_GAUSS_PT));
 410   fakeF->setMesh(geoMesh);
 411   fakeF->setGaussLocalizationOnType(gt,loc.getRefCoords(),loc.getGaussCoords(),loc.getGaussWeights());
 412   ptsForLoc=fakeF->getLocalizationOfDiscr();
 413   //
 414   return geoMesh;
 415 }
 416
 417 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)
 418 {
 419   MCAuto<DataArrayDouble> ptsForLoc;
 420   MCAuto<MEDCouplingUMesh> geoMesh(BuildMeshCommon(gt,pfl,loc,zeStr,mesh,section,globs,ptsForLoc));
 421   //
 422   MCConstAuto<DataArrayDouble> angleVrille;
 423   if(!pfl.empty())
 424     {
 425       const DataArrayIdType *pflArr(globs->getProfile(pfl));
 426       angleVrille=angleDeVrille->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
 427     }
 428   else
 429     angleVrille.takeRef(angleDeVrille);
 430   //
 431   MCAuto<MEDCouplingFieldDouble> dir(geoMesh->buildDirectionVectorField());
 432   MCAuto<DataArrayDouble> rot(dir->getArray()->fromCartToSpher());
 433   std::size_t nbCompo(ptsForLoc->getNumberOfComponents());
 434   MCAuto<DataArrayDouble> secPts(section->getCoords()->changeNbOfComponents(nbCompo,0.));
 435   mcIdType nbSecPts(secPts->getNumberOfTuples()),nbCells(geoMesh->getNumberOfCells()),nbg(ToIdType(loc.getGaussWeights().size()));
 436   {
 437     const int TAB[3]={2,0,1};
 438     std::vector<std::size_t> v(TAB,TAB+3);
 439     secPts=secPts->keepSelectedComponents(v);
 440   }
 441   const double CENTER[3]={0.,0.,0.},AX0[3]={0.,0.,1.};
 442   double AX1[3]; AX1[2]=0.;
 443   std::vector< MCAuto<DataArrayDouble> > arrs(nbCells*nbg);
 444   for(int j=0;j<nbCells;j++)
 445     {
 446       MCAuto<DataArrayDouble> p(secPts->deepCopy());
 447       double ang0(rot->getIJ(j,2));
 448       DataArrayDouble::Rotate3DAlg(CENTER,AX0,ang0,nbSecPts,p->begin(),p->getPointer());
 449       AX1[0]=-sin(ang0); AX1[1]=cos(ang0);// rot Oy around OZ
 450       double ang1(M_PI/2.-rot->getIJ(j,1));
 451       DataArrayDouble::Rotate3DAlg(CENTER,AX1,-ang1,nbSecPts,p->begin(),p->getPointer());
 452       DataArrayDouble::Rotate3DAlg(CENTER,dir->getArray()->begin()+j*3,angleVrille->getIJ(j,0),nbSecPts,p->begin(),p->getPointer());
 453       for(int l=0;l<nbg;l++)
 454         {
 455           MCAuto<DataArrayDouble> p2(p->deepCopy());
 456           for(std::size_t k=0;k<nbCompo;k++)
 457             p2->applyLin(1.,ptsForLoc->getIJ(j*nbg+l,k),k);
 458           arrs[j*nbg+l]=p2;
 459         }
 460     }
 461   std::vector<const DataArrayDouble *> arrs2(VecAutoToVecOfCstPt(arrs));
 462   MCAuto<DataArrayDouble> resu(DataArrayDouble::Aggregate(arrs2));
 463   return resu;
 464 }
 465
 466 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)
 467 {
 468   MCAuto<DataArrayDouble> ptsForLoc;
 469   MCAuto<MEDCouplingUMesh> geoMesh(BuildMeshCommon(gt,pfl,loc,zeStr,mesh,section,globs,ptsForLoc));
 470   mcIdType nbCells(geoMesh->getNumberOfCells()),nbg(ToIdType(loc.getGaussWeights().size()));
 471   MCConstAuto<DataArrayDouble> zeThickness;
 472   if(!pfl.empty())
 473     {
 474       const DataArrayIdType *pflArr(globs->getProfile(pfl));
 475       zeThickness=thickness->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
 476     }
 477   else
 478     zeThickness.takeRef(thickness);
 479   MCAuto<DataArrayDouble> orthoArr;
 480   {
 481     MCAuto<MEDCouplingFieldDouble> ortho(geoMesh->buildOrthogonalField());
 482     orthoArr.takeRef(ortho->getArray());
 483   }
 484   mcIdType nbCompo(ToIdType(orthoArr->getNumberOfComponents()));
 485   MCAuto<DataArrayDouble> secPts(section->getCoords()->duplicateEachTupleNTimes(nbCompo));
 486   secPts->rearrange(nbCompo);
 487   std::vector< MCAuto<DataArrayDouble> > arrs(nbCells*nbg);
 488   for(mcIdType j=0;j<nbCells;j++)
 489     {
 490       double thck(zeThickness->getIJ(j,0)),eccentricity(zeThickness->getIJ(j,1));
 491       MCAuto<DataArrayDouble> fact(DataArrayDouble::New()),fact2(DataArrayDouble::New()); fact->alloc(1,nbCompo); fact2->alloc(1,nbCompo);
 492       std::copy(orthoArr->begin()+j*nbCompo,orthoArr->begin()+(j+1)*nbCompo,fact->getPointer());
 493       std::copy(orthoArr->begin()+j*nbCompo,orthoArr->begin()+(j+1)*nbCompo,fact2->getPointer());
 494       std::for_each(fact2->rwBegin(),fact2->rwEnd(),[eccentricity](double& val){ val*=eccentricity; });
 495       std::transform(fact->begin(),fact->end(),fact->getPointer(),[thck](const double& val){ return val*thck/2.; });
 496       MCAuto<DataArrayDouble> p(DataArrayDouble::Multiply(secPts,fact));
 497       p=DataArrayDouble::Add(p,fact2);
 498       for(int l=0;l<nbg;l++)
 499         {
 500           MCAuto<DataArrayDouble> p2(p->deepCopy());
 501           for(int k=0;k<nbCompo;k++)
 502             p2->applyLin(1.,ptsForLoc->getIJ(j*nbg+l,k),k);
 503           arrs[j*nbg+l]=p2;
 504         }
 505     }
 506   std::vector<const DataArrayDouble *> arrs2(VecAutoToVecOfCstPt(arrs));
 507   MCAuto<DataArrayDouble> resu(DataArrayDouble::Aggregate(arrs2));
 508   return resu;
 509 }
 510
 511 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 {
 513   static const char MSG1[]="BuildMeshPipeSEG3 : not recognized pattern ! Send mail to anthony.geay@edf.fr with corresponding MED file !";
 514   MCAuto<DataArrayDouble> ptsForLoc;
 515   MCAuto<MEDCouplingUMesh> geoMesh(BuildMeshCommon(INTERP_KERNEL::NORM_SEG3,pfl,loc,zeStr,mesh,section,globs,ptsForLoc));
 516   mcIdType nbSecPts(section->getNumberOfNodes()),nbCells(geoMesh->getNumberOfCells()),nbg(ToIdType(loc.getGaussWeights().size()));
 517   MCConstAuto<DataArrayDouble> zeAngle,zeScale;
 518   if(!pfl.empty())
 519     {
 520       const DataArrayIdType *pflArr(globs->getProfile(pfl));
 521       zeAngle=angle->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
 522       zeScale=scale->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
 523     }
 524   else
 525     {
 526       zeAngle.takeRef(angle);
 527       zeScale.takeRef(scale);
 528     }
 529   if(zeAngle->getNumberOfComponents()!=3 || zeScale->getNumberOfComponents()!=2 || nbg!=3)
 530     throw INTERP_KERNEL::Exception(MSG1);
 531   MCAuto<MEDCouplingFieldDouble> dir;
 532   {
 533     MCAuto<MEDCouplingUMesh> geoMesh2(geoMesh->deepCopy());
 534     geoMesh2->convertQuadraticCellsToLinear();
 535     dir=geoMesh2->buildDirectionVectorField();
 536   }
 537   MCAuto<DataArrayDouble> rot(dir->getArray()->fromCartToSpher());
 538   std::size_t nbCompo(ptsForLoc->getNumberOfComponents());
 539   MCAuto<DataArrayDouble> secPts(section->getCoords()->changeNbOfComponents(nbCompo,0.));
 540   {
 541     const int TAB[3]={2,0,1};
 542     std::vector<std::size_t> v(TAB,TAB+3);
 543     secPts=secPts->keepSelectedComponents(v);
 544   }
 545   const double CENTER[3]={0.,0.,0.},AX0[3]={0.,0.,1.};
 546   double AX1[3]; AX1[2]=0.;
 547   std::vector< MCAuto<DataArrayDouble> > arrs(nbCells*nbg);
 548   for(mcIdType j=0;j<nbCells;j++)
 549     {
 550       constexpr int DIM=3;
 551       MCAuto<DataArrayDouble> p(secPts->deepCopy());
 552       double ang0(rot->getIJ(j,2));
 553       double rmin(zeScale->getIJ(j,0)),rmax(zeScale->getIJ(j,1));
 554       {
 555         auto pt(p->rwBegin());
 556         for(int i=0;i<nbSecPts;i++)
 557           {
 558             auto nrm(sqrt(std::accumulate(pt,pt+DIM,0.,[](double sum, double v) { return sum+v*v; } )));
 559             auto sca((rmin+2.*(nrm-0.5)*(rmax-rmin))/nrm);
 560             std::for_each(pt,pt+3,[sca](double& val) { val*=sca; } );
 561             std::advance(pt,DIM);
 562           }
 563       }
 564       DataArrayDouble::Rotate3DAlg(CENTER,AX0,ang0,nbSecPts,p->begin(),p->getPointer());
 565       AX1[0]=-sin(ang0); AX1[1]=cos(ang0);// rot Oy around OZ
 566       double ang1(M_PI/2.-rot->getIJ(j,1));
 567       DataArrayDouble::Rotate3DAlg(CENTER,AX1,-ang1,nbSecPts,p->begin(),p->getPointer());
 568       for(int l=0;l<3;l++)
 569         {
 570           MCAuto<DataArrayDouble> p3(p->deepCopy());
 571           DataArrayDouble::Rotate3DAlg(CENTER,dir->getArray()->begin()+j*3,zeAngle->getIJ(j,l),nbSecPts,p3->begin(),p3->getPointer());
 572           MCAuto<DataArrayDouble> p2(p3->deepCopy());
 573           for(std::size_t k=0;k<nbCompo;k++)
 574             p2->applyLin(1.,ptsForLoc->getIJ(j*nbg+l,k),k);
 575           arrs[j*nbg+l]=p2;
 576         }
 577     }
 578   std::vector<const DataArrayDouble *> arrs2(VecAutoToVecOfCstPt(arrs));
 579   MCAuto<DataArrayDouble> resu(DataArrayDouble::Aggregate(arrs2));
 580   return resu;
 581 }
 582
 583 MCAuto<DataArrayDouble> LocInfo::BuildMeshPipeSEG4(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)
 584 {
 585   static const char MSG1[]="BuildMeshPipeSEG4 : not recognized pattern ! Send mail to anthony.geay@edf.fr with corresponding MED file !";
 586   MCAuto<DataArrayDouble> ptsForLoc;
 587   MCAuto<MEDCouplingUMesh> geoMesh(BuildMeshCommon(INTERP_KERNEL::NORM_SEG4,pfl,loc,zeStr,mesh,section,globs,ptsForLoc));
 588   mcIdType nbSecPts(section->getNumberOfNodes()),nbCells(geoMesh->getNumberOfCells()),nbg(ToIdType(loc.getGaussWeights().size()));
 589   MCConstAuto<DataArrayDouble> zeAngle,zeScale;
 590   if(!pfl.empty())
 591     {
 592       const DataArrayIdType *pflArr(globs->getProfile(pfl));
 593       zeAngle=angle->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
 594       zeScale=scale->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
 595     }
 596   else
 597     {
 598       zeAngle.takeRef(angle);
 599       zeScale.takeRef(scale);
 600     }
 601   if(zeAngle->getNumberOfComponents()!=4 || zeScale->getNumberOfComponents()!=2 || nbg!=3)
 602     throw INTERP_KERNEL::Exception(MSG1);
 603   MCAuto<MEDCouplingFieldDouble> dir;
 604   {
 605     MCAuto<MEDCouplingUMesh> geoMesh2(geoMesh->deepCopy());
 606     geoMesh2->convertQuadraticCellsToLinear();
 607     dir=geoMesh2->buildDirectionVectorField();
 608   }
 609   MCAuto<DataArrayDouble> rot(dir->getArray()->fromCartToSpher());
 610   std::size_t nbCompo(ptsForLoc->getNumberOfComponents());
 611   MCAuto<DataArrayDouble> secPts(section->getCoords()->changeNbOfComponents(nbCompo,0.));
 612   {
 613     const int TAB[3]={2,0,1};
 614     std::vector<std::size_t> v(TAB,TAB+3);
 615     secPts=secPts->keepSelectedComponents(v);
 616   }
 617   const double CENTER[3]={0.,0.,0.},AX0[3]={0.,0.,1.};
 618   double AX1[3]; AX1[2]=0.;
 619   std::vector< MCAuto<DataArrayDouble> > arrs(nbCells*nbg);
 620   for(mcIdType j=0;j<nbCells;j++)
 621     {
 622       constexpr int DIM=3;
 623       MCAuto<DataArrayDouble> p(secPts->deepCopy());
 624       double ang0(rot->getIJ(j,2));
 625       double rmin(zeScale->getIJ(j,0)),rmax(zeScale->getIJ(j,1));
 626       {
 627         auto pt(p->rwBegin());
 628         for(int i=0;i<nbSecPts;i++)
 629           {
 630             auto nrm(sqrt(std::accumulate(pt,pt+DIM,0.,[](double sum, double v) { return sum+v*v; } )));
 631             auto sca((rmin+2.*(nrm-0.5)*(rmax-rmin))/nrm);
 632             std::for_each(pt,pt+3,[sca](double& val) { val*=sca; } );
 633             std::advance(pt,DIM);
 634           }
 635       }
 636       DataArrayDouble::Rotate3DAlg(CENTER,AX0,ang0,nbSecPts,p->begin(),p->getPointer());
 637       AX1[0]=-sin(ang0); AX1[1]=cos(ang0);// rot Oy around OZ
 638       double ang1(M_PI/2.-rot->getIJ(j,1));
 639       DataArrayDouble::Rotate3DAlg(CENTER,AX1,-ang1,nbSecPts,p->begin(),p->getPointer());
 640       for(int l=0;l<3;l++)
 641         {
 642           MCAuto<DataArrayDouble> p3(p->deepCopy());
 643           DataArrayDouble::Rotate3DAlg(CENTER,dir->getArray()->begin()+j*3,zeAngle->getIJ(j,l),nbSecPts,p3->begin(),p3->getPointer());
 644           MCAuto<DataArrayDouble> p2(p3->deepCopy());
 645           for(std::size_t k=0;k<nbCompo;k++)
 646             p2->applyLin(1.,ptsForLoc->getIJ(j*nbg+l,k),k);
 647           arrs[j*nbg+l]=p2;
 648         }
 649     }
 650   std::vector<const DataArrayDouble *> arrs2(VecAutoToVecOfCstPt(arrs));
 651   MCAuto<DataArrayDouble> resu(DataArrayDouble::Aggregate(arrs2));
 652   return resu;
 653 }
 654
 655 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)
 656 {
 657   static const char MSG1[]="BuildMeshFromStructure : not recognized pattern ! Send mail to anthony.geay@edf.fr with corresponding MED file !";
 658   const std::vector< MCAuto<DataArray> >& vars(zeStr->getVars());
 659   if(vars.size()==1)
 660     {
 661       MCAuto<DataArray> zeArr(vars[0]);
 662       if(zeArr.isNull())
 663         throw INTERP_KERNEL::Exception(MSG1);
 664       MCAuto<DataArrayDouble> zeArr2(DynamicCast<DataArray,DataArrayDouble>(zeArr));
 665       if(zeArr2.isNull())
 666         throw INTERP_KERNEL::Exception(MSG1);
 667       if(zeArr2->getName()==ANGLE_DE_VRILLE || zeArr2->getName()==ANGLE)
 668         return BuildMeshFromAngleVrille(gt,zeArr2,pfl,loc,zeStr,mesh,section,globs);
 669       if(zeArr2->getName()==EPAISSEUR || zeArr2->getName()==SCALE)
 670         return BuildMeshFromEpaisseur(gt,zeArr2,pfl,loc,zeStr,mesh,section,globs);
 671     }
 672   if(vars.size()==2)
 673     {
 674       MCAuto<DataArray> zeArr0(vars[0]),zeArr1(vars[1]);
 675       if(zeArr0.isNull() || zeArr1.isNull())
 676         throw INTERP_KERNEL::Exception(MSG1);
 677       MCAuto<DataArrayDouble> zeArr00(DynamicCastSafe<DataArray,DataArrayDouble>(zeArr0)),zeArr11(DynamicCastSafe<DataArray,DataArrayDouble>(zeArr1));
 678       MCAuto<DataArrayDouble> angle,scale;
 679       if(zeArr00->getName()==ANGLE)
 680         angle=zeArr00;
 681       if(zeArr00->getName()==SCALE)
 682         scale=zeArr00;
 683       if(zeArr11->getName()==ANGLE)
 684         angle=zeArr11;
 685       if(zeArr11->getName()==SCALE)
 686         scale=zeArr11;
 687       if(angle.isNull() || scale.isNull())
 688         throw INTERP_KERNEL::Exception(MSG1);
 689       switch(gt)
 690       {
 691         case INTERP_KERNEL::NORM_SEG3:
 692           {
 693             return BuildMeshPipeSEG3(angle,scale,pfl,loc,zeStr,mesh,section,globs);
 694           }
 695         case INTERP_KERNEL::NORM_SEG4:
 696           {
 697             return BuildMeshPipeSEG4(angle,scale,pfl,loc,zeStr,mesh,section,globs);
 698           }
 699         default:
 700           throw INTERP_KERNEL::Exception(MSG1);
 701       }
 702     }
 703   throw INTERP_KERNEL::Exception(MSG1);
 704 }
 705
 706 MCAuto<MEDFileUMesh> LocInfo::generateNonClassicalData(int zePos, const MEDFileUMesh *mesh, const MEDFileFieldGlobsReal *globs) const
 707 {
 708   static const char MSG1[]="LocInfo::generateNonClassicalData : no spec for GAUSS on StructureElement with more than one cell !";
 709   std::size_t sz(_locs.size());
 710   std::vector< MCAuto<DataArrayDouble> > arrs(sz);
 711   for(std::size_t i=0;i<sz;i++)
 712     {
 713       const MEDFileFieldLoc& loc(globs->getLocalization(_locs[i]));
 714       const MEDFileGTKeeper *gtk(loc.getUndergroundGTKeeper());
 715       const MEDFileGTKeeperDyn *gtk2(dynamic_cast<const MEDFileGTKeeperDyn *>(gtk));
 716       if(!gtk2)
 717         throw INTERP_KERNEL::Exception("LocInfo::generateNonClassicalData : internal error !");
 718       const MEDFileUMesh *meshLoc(gtk2->getMesh()),*section(gtk2->getSection());
 719       const MEDFileStructureElement *se(gtk2->getSE());
 720       MCAuto<MEDCouplingUMesh> um(meshLoc->getMeshAtLevel(0));
 721       INTERP_KERNEL::NormalizedCellType gt(_cts[i]);
 722       {
 723         std::vector<int> nel(meshLoc->getNonEmptyLevels());
 724         if(nel.size()!=1)
 725           throw INTERP_KERNEL::Exception(MSG1);
 726         if(nel[0]!=0)
 727           throw INTERP_KERNEL::Exception(MSG1);
 728         mcIdType zePos(-1);
 729         for(mcIdType icell = 0 ; icell < um->getNumberOfCells() ; ++icell)
 730           if( gt == um->getTypeOfCell(icell) )
 731             zePos = icell;
 732         if(zePos == -1)
 733           throw INTERP_KERNEL::Exception(MSG1);
 734         std::vector<mcIdType> v;
 735         um->getNodeIdsOfCell(zePos,v);
 736         std::size_t sz2(v.size());
 737         for(std::size_t j=0;j<sz2;j++)
 738           if(v[j]!=ToIdType(j))
 739             throw INTERP_KERNEL::Exception(MSG1);
 740       }
 741       const std::vector< MCAuto<MEDFileEltStruct4Mesh> >& strs(mesh->getAccessOfUndergroundEltStrs());
 742       MCAuto<MEDFileEltStruct4Mesh> zeStr;
 743       for(std::vector< MCAuto<MEDFileEltStruct4Mesh> >::const_iterator it=strs.begin();it!=strs.end();it++)
 744         {
 745           if((*it)->getGeoTypeName()==se->getName())
 746             {
 747               zeStr=*it;
 748               break;
 749             }
 750         }
 751       if(zeStr.isNull())
 752         {
 753           std::ostringstream oss; oss << "LocInfo::generateNonClassicalData :  : no geo type with name " <<  se->getName() << " in " << mesh->getName() << " !";
 754           throw INTERP_KERNEL::Exception(oss.str());
 755         }
 756       arrs[i]=BuildMeshFromStructure(gt,_pfl[i],loc,zeStr,mesh,section,globs);
 757     }
 758   std::vector<const DataArrayDouble *> arrs2(VecAutoToVecOfCstPt(arrs));
 759   MCAuto<DataArrayDouble> resu(DataArrayDouble::Aggregate(arrs2));
 760   MCAuto<MEDFileUMesh> ret(MEDFileUMesh::New());
 761   ret->setCoords(resu);
 762   std::ostringstream meshName; meshName << mesh->getName() << "_on_" << sz << "_sections" << "_" << zePos;
 763   ret->setName(meshName.str());
 764   return ret;
 765 }
 766
 767 FieldWalker2::FieldWalker2(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd)
 768 {
 769   _loc=pmptpd->getLocalization();
 770   _pfl=pmptpd->getProfile();
 771   _ct=pmptpd->getGeoTypeStatic();
 772   _is_classic=pmptpd->getType()!=ON_GAUSS_PT;
 773   _pd=SlicePartDefinition::New(pmptpd->getStart(),pmptpd->getEnd(),1);
 774 }
 775
 776 bool FieldWalker2::operator!=(const FieldWalker2& other) const
 777 {
 778   return !((*this)==other);
 779 }
 780
 781 bool FieldWalker2::operator==(const FieldWalker2& other) const
 782 {
 783   bool ret2(false);
 784   {
 785     std::string tmp;
 786     ret2=_pd->isEqual(other._pd,tmp);
 787   }
 788   bool ret(_loc==other._loc && _pfl==other._pfl && _is_classic==other._is_classic && ret2);
 789   return ret;
 790 }
 791
 792 class FieldWalker1
 793 {
 794 public:
 795   FieldWalker1(const MEDFileAnyTypeField1TSWithoutSDA *ts):_ts(ts),_pm_pt(0),_nb_mesh(0) { }
 796   void newMeshEntry(const MEDFileFieldPerMesh *fpm);
 797   void endMeshEntry(const MEDFileFieldPerMesh *fpm) { }
 798   void newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt);
 799   void endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt);
 800   void newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd);
 801   void checkOK(const FieldWalker1& other) const;
 802   bool isClassical() const;
 803   std::vector<FieldWalker2> getNonClassicalData() const { return _fw; }
 804 private:
 805   const MEDFileAnyTypeField1TSWithoutSDA *_ts;
 806   const MEDFileFieldPerMeshPerTypeDyn *_pm_pt;
 807   std::vector<FieldWalker2> _fw;
 808   int _nb_mesh;
 809 };
 810
 811 void FieldWalker1::newMeshEntry(const MEDFileFieldPerMesh *fpm)
 812 {
 813   if(_nb_mesh++==1)
 814     throw INTERP_KERNEL::Exception("FieldWalker1::newMeshEntry : multi mesh not supported !");
 815 }
 816
 817 void FieldWalker1::newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt)
 818 {
 819   if(_pm_pt)
 820     throw INTERP_KERNEL::Exception("FieldWalker1::newPerMeshPerTypeEntry : multi SE loc not managed yet !");
 821   const MEDFileFieldPerMeshPerTypeDyn *pmpt2(dynamic_cast<const MEDFileFieldPerMeshPerTypeDyn *>(pmpt));
 822   if(!pmpt2)
 823     throw INTERP_KERNEL::Exception("newPerMeshPerTypeEntry : internal error !");
 824   _pm_pt=pmpt2;
 825 }
 826
 827 void FieldWalker1::endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *)
 828 {
 829   isClassical();
 830 }
 831
 832 void FieldWalker1::newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd)
 833 {
 834   _fw.push_back(FieldWalker2(pmptpd));
 835 }
 836
 837 void FieldWalker1::checkOK(const FieldWalker1& other) const
 838 {
 839   std::size_t sz(_fw.size());
 840   if(other._fw.size()!=sz)
 841     throw INTERP_KERNEL::Exception("checkOK : not OK because size are not the same !");
 842   for(std::size_t i=0;i<sz;i++)
 843     if(_fw[i]!=other._fw[i])
 844       throw INTERP_KERNEL::Exception("checkOK : not OK because an element mismatches !");
 845 }
 846
 847 bool FieldWalker1::isClassical() const
 848 {
 849   if(_fw.empty())
 850     throw INTERP_KERNEL::Exception("FieldWalker1::endPerMeshPerTypeEntry : internal error !");
 851   std::size_t ic(0),inc(0);
 852   for(std::vector<FieldWalker2>::const_iterator it=_fw.begin();it!=_fw.end();it++)
 853     {
 854       if((*it).isClassic())
 855         ic++;
 856       else
 857         inc++;
 858     }
 859   if(ic!=0 && inc!=0)
 860     throw INTERP_KERNEL::Exception("FieldWalker1::endPerMeshPerTypeEntry : mix is not allowed yet !");
 861   return inc==0;
 862 }
 863
 864 class FieldWalker
 865 {
 866 public:
 867   FieldWalker(const MEDFileAnyTypeFieldMultiTSWithoutSDA *f):_f(f) { }
 868   void newTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts);
 869   void endTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts);
 870   void newMeshEntry(const MEDFileFieldPerMesh *fpm);
 871   void endMeshEntry(const MEDFileFieldPerMesh *fpm);
 872   void newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt);
 873   void endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt);
 874   void newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd);
 875 public:
 876   bool isEmpty() const;
 877   bool isClassical() const;
 878   const MEDFileAnyTypeFieldMultiTSWithoutSDA *field() const { return _f; }
 879   std::vector<FieldWalker2> getNonClassicalData() const { return _fw_prev->getNonClassicalData(); }
 880 private:
 881   const MEDFileAnyTypeFieldMultiTSWithoutSDA *_f;
 882   mutable INTERP_KERNEL::AutoCppPtr<FieldWalker1> _fw;
 883   mutable INTERP_KERNEL::AutoCppPtr<FieldWalker1> _fw_prev;
 884 };
 885
 886 bool FieldWalker::isEmpty() const
 887 {
 888   return _fw_prev.isNull();
 889 }
 890
 891 bool FieldWalker::isClassical() const
 892 {
 893   return _fw_prev->isClassical();
 894 }
 895
 896 void FieldWalker::newTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts)
 897 {
 898   _fw=new FieldWalker1(ts);
 899 }
 900
 901 void FieldWalker::endTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts)
 902 {
 903   if(_fw_prev.isNull())
 904     _fw_prev=new FieldWalker1(*_fw);
 905   else
 906     _fw_prev->checkOK(*_fw);
 907   _fw=0;
 908 }
 909
 910 void FieldWalker::newMeshEntry(const MEDFileFieldPerMesh *fpm)
 911 {
 912   _fw->newMeshEntry(fpm);
 913 }
 914
 915 void FieldWalker::endMeshEntry(const MEDFileFieldPerMesh *fpm)
 916 {
 917   _fw->endMeshEntry(fpm);
 918 }
 919
 920 void FieldWalker::newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt)
 921 {
 922   _fw->newPerMeshPerTypeEntry(pmpt);
 923 }
 924
 925 void FieldWalker::endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt)
 926 {
 927   _fw->endPerMeshPerTypeEntry(pmpt);
 928 }
 929
 930 void FieldWalker::newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd)
 931 {
 932   _fw->newPerMeshPerTypePerDisc(pmptpd);
 933 }
 934
 935 // this class splits fields into same
 936 class LocSpliter : public MEDFileFieldVisitor
 937 {
 938 public:
 939   LocSpliter(const MEDFileFieldGlobsReal *globs):_globs(globs),_fw(0) { }
 940   MCAuto<MEDFileFields> getClassical() const { return _classical; }
 941   void generateNonClassicalData(const MEDFileUMesh *mesh, std::vector< MCAuto<MEDFileFields> >& outFields, std::vector< MCAuto<MEDFileUMesh> >& outMeshes) const;
 942 private:
 943   void newFieldEntry(const MEDFileAnyTypeFieldMultiTSWithoutSDA *field);
 944   void endFieldEntry(const MEDFileAnyTypeFieldMultiTSWithoutSDA *field);
 945   //
 946   void newTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts);
 947   void endTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts);
 948   //
 949   void newMeshEntry(const MEDFileFieldPerMesh *fpm);
 950   void endMeshEntry(const MEDFileFieldPerMesh *fpm);
 951   //
 952   void newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt);
 953   void endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt);
 954   //
 955   void newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd);
 956 private:
 957   const MEDFileFieldGlobsReal *_globs;
 958   std::vector< LocInfo > _locs;
 959   std::vector< MCAuto<MEDFileFields> > _fields_on_locs;//size of _locs== size of _fields_on_locs
 960   MCAuto<MEDFileFields> _classical;
 961 private:
 962   mutable INTERP_KERNEL::AutoCppPtr<FieldWalker> _fw;
 963 };
 964
 965 void LocSpliter::newFieldEntry(const MEDFileAnyTypeFieldMultiTSWithoutSDA *field)
 966 {
 967   _fw=new FieldWalker(field);
 968 }
 969
 970 void LocSpliter::endFieldEntry(const MEDFileAnyTypeFieldMultiTSWithoutSDA *field)
 971 {
 972   if(_fw->isEmpty())
 973     return ;
 974   MCAuto<MEDFileAnyTypeFieldMultiTS> f(MEDFileAnyTypeFieldMultiTS::BuildNewInstanceFromContent(const_cast<MEDFileAnyTypeFieldMultiTSWithoutSDA *>(field)));
 975   if(_fw->isClassical())
 976     {
 977       if(_classical.isNull())
 978         {
 979           _classical=MEDFileFields::New();
 980           _classical->shallowCpyGlobs(*_globs);
 981         }
 982       _classical->pushField(f);
 983     }
 984   else
 985     {
 986       std::vector<FieldWalker2> fw2(_fw->getNonClassicalData());
 987       LocInfo elt(fw2);
 988       std::vector< LocInfo >::iterator it(std::find(_locs.begin(),_locs.end(),elt));
 989       if(it==_locs.end())
 990         {
 991           _locs.push_back(elt);
 992           MCAuto<MEDFileFields> zeF(MEDFileFields::New());
 993           zeF->shallowCpyGlobs(*_globs);
 994           zeF->pushField(f);
 995           _fields_on_locs.push_back(zeF);
 996         }
 997       else
 998         {
 999           MCAuto<MEDFileFields> zeF(_fields_on_locs[std::distance(_locs.begin(),it)]);
1000           zeF->pushField(f);
1001         }
1002     }
1003 }
1004
1005 void LocSpliter::generateNonClassicalData(const MEDFileUMesh *mesh, std::vector< MCAuto<MEDFileFields> >& outFields, std::vector< MCAuto<MEDFileUMesh> >& outMeshes) const
1006 {
1007   int i(0);
1008   for(std::vector<LocInfo>::const_iterator it=_locs.begin();it!=_locs.end();it++,i++)
1009     {
1010       MCAuto<MEDFileUMesh> m((*it).generateNonClassicalData(i,mesh,_globs));
1011       outMeshes.push_back(m);
1012       MCAuto<MEDCouplingUMesh> mcm(MEDCouplingUMesh::Build0DMeshFromCoords(m->getCoords()));
1013       mcm->setName(m->getName());
1014       MCAuto<MEDFileFields> fs(_fields_on_locs[i]);
1015       MCAuto<MEDFileFields> outFs(MEDFileFields::New());
1016       for(int j=0;j<fs->getNumberOfFields();j++)
1017         {
1018           MCAuto<MEDFileAnyTypeFieldMultiTS> fmtsNC(fs->getFieldAtPos(j));
1019           MCAuto<MEDFileFieldMultiTS> fmts(DynamicCastSafe<MEDFileAnyTypeFieldMultiTS,MEDFileFieldMultiTS>(fmtsNC));
1020           MCAuto<MEDFileFieldMultiTS> outFmts(MEDFileFieldMultiTS::New());
1021           for(int k=0;k<fmts->getNumberOfTS();k++)
1022             {
1023               MCAuto<MEDFileField1TS> outF1t(MEDFileField1TS::New());
1024               MCAuto<MEDFileField1TS> f1ts(fmts->getTimeStepAtPos(k));
1025               int t2,t3;
1026               double t1(f1ts->getTime(t2,t3));
1027               MCAuto<MEDCouplingFieldDouble> mcf(MEDCouplingFieldDouble::New(ON_NODES));
1028               MCAuto<DataArrayDouble> arr,arr2;
1029               arr.takeRef(f1ts->getUndergroundDataArray());
1030               arr2=arr->selectPartDef((*it).getPartDef());
1031               mcf->setArray(arr2);
1032               mcf->setTime(t1,t2,t3);
1033               mcf->setName(f1ts->getName());
1034               mcf->setMesh(mcm);
1035               outF1t->setFieldNoProfileSBT(mcf);
1036               outFmts->pushBackTimeStep(outF1t);
1037             }
1038           outFs->pushField(outFmts);
1039         }
1040       outFields.push_back(outFs);
1041     }
1042 }
1043
1044 void LocSpliter::newTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts)
1045 {
1046   _fw->newTimeStepEntry(ts);
1047 }
1048
1049 void LocSpliter::endTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts)
1050 {
1051   _fw->endTimeStepEntry(ts);
1052 }
1053
1054 void LocSpliter::newMeshEntry(const MEDFileFieldPerMesh *fpm)
1055 {
1056   _fw->newMeshEntry(fpm);
1057 }
1058
1059 void LocSpliter::endMeshEntry(const MEDFileFieldPerMesh *fpm)
1060 {
1061   _fw->endMeshEntry(fpm);
1062 }
1063
1064 void LocSpliter::newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt)
1065 {
1066   _fw->newPerMeshPerTypeEntry(pmpt);
1067 }
1068
1069 void LocSpliter::endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt)
1070 {
1071   _fw->endPerMeshPerTypeEntry(pmpt);
1072 }
1073
1074 void LocSpliter::newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd)
1075 {
1076   _fw->newPerMeshPerTypePerDisc(pmptpd);
1077 }
1078
1079 void MEDFileBlowStrEltUp::DealWithConflictNames(MEDFileAnyTypeFieldMultiTS *fmtsToAdd, const MEDFileFields *fs)
1080 {
1081   std::vector<std::string> fnames(fs->getFieldsNames());
1082   for(int i=0;i<1000;i++)
1083     {
1084       std::ostringstream oss; oss << fmtsToAdd->getName();
1085       if(i>=1)
1086         oss << "_" << i-1;
1087       if(std::find(fnames.begin(),fnames.end(),oss.str())==fnames.end())
1088         {
1089           fmtsToAdd->setName(oss.str());
1090           return ;
1091         }
1092     }
1093   throw INTERP_KERNEL::Exception("DealWithConflictNames : Eh eh interesting !");
1094 }
1095
1096 MCAuto<MEDFileFields> MEDFileBlowStrEltUp::splitFieldsPerLoc(const MEDFileFields *fields, const MEDFileUMesh *mesh, MEDFileMeshes *msOut, MEDFileFields *allZeOutFields)
1097 {
1098   LocSpliter ls(fields);
1099   fields->accept(ls);
1100   std::vector< MCAuto<MEDFileFields> > outFields;
1101   std::vector< MCAuto<MEDFileUMesh> > outMeshes;
1102   ls.generateNonClassicalData(mesh,outFields,outMeshes);
1103   for(std::vector< MCAuto<MEDFileFields> >::iterator it=outFields.begin();it!=outFields.end();it++)
1104     {
1105       for(int j=0;j<(*it)->getNumberOfFields();j++)
1106         {
1107           MCAuto<MEDFileAnyTypeFieldMultiTS> fmts((*it)->getFieldAtPos(j));
1108           //DealWithConflictNames(fmts,allZeOutFields);// uncomment to have a writable data structure
1109           allZeOutFields->pushField(fmts);
1110         }
1111     }
1112   for(std::vector< MCAuto<MEDFileUMesh> >::iterator it=outMeshes.begin();it!=outMeshes.end();it++)
1113     msOut->pushMesh(*it);
1114   return ls.getClassical();
1115 }