src/MEDLoader/MEDFileBlowStrEltUp.cxx

   1 // Copyright (C) 2007-2017  CEA/DEN, EDF R&D
   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 "MCAuto.txx"
  25
  26 using namespace MEDCoupling;
  27
  28 const char MEDFileBlowStrEltUp::MED_BALL_STR[]="MED_BALL";
  29
  30 MEDFileBlowStrEltUp::MEDFileBlowStrEltUp(const MEDFileFields *fsOnlyOnSE, const MEDFileMeshes *ms, const MEDFileStructureElements *ses)
  31 {
  32   if(!fsOnlyOnSE || !ms || !ses)
  33     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp constructor : NULL input pointer !");
  34   _ms.takeRef(ms); _ses.takeRef(ses);
  35   std::vector< std::pair<std::string,std::string> > ps;
  36   fsOnlyOnSE->getMeshSENames(ps);
  37   std::size_t sz(ps.size());
  38   _elts.resize(sz);
  39   for(std::size_t i=0;i<sz;i++)
  40     {
  41       const std::pair<std::string,std::string>& p(ps[i]);
  42       MCAuto<MEDFileFields> f(fsOnlyOnSE->partOfThisLyingOnSpecifiedMeshSEName(p.first,p.second));
  43       _elts[i]=f;
  44     }
  45   for(std::size_t i=0;i<sz;i++)
  46     {
  47       const std::pair<std::string,std::string>& p(ps[i]);
  48       MEDFileMesh *mesh(_ms->getMeshWithName(p.first));
  49       if(!mesh)
  50         throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp : NULL mesh !");
  51       MEDFileUMesh *umesh(dynamic_cast<MEDFileUMesh *>(mesh));
  52       if(!umesh)
  53         throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp : Blow up of Stru Elt not managed yet for unstructured meshes !");
  54     }
  55 }
  56
  57 /*!
  58  * \param [in] mesh - The mesh containing structure element called \a seName. After the call of this method the Structure elements parts will be removed.
  59  * \param [out] mOut - the physical mesh of the structure element \a seName in mesh \a mesh
  60  * \param [out] fsOut - the list of var attribute of structure element \a seName - \b WARNING no time steps here
  61  */
  62 MCAuto<MEDFileEltStruct4Mesh> MEDFileBlowStrEltUp::dealWithSEInMesh(const std::string& seName, MEDFileUMesh *mesh, MCAuto<MEDFileUMesh>& mOut, MCAuto<MEDFileFields>& fsOut) const
  63 {
  64   if(!mesh)
  65     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithSEInMesh : null pointer !");
  66   if(seName==MED_BALL_STR)
  67     {
  68       MCAuto<MEDFileEltStruct4Mesh> ret(dealWithMEDBALLInMesh(mesh,mOut,fsOut));
  69       mesh->killStructureElements();
  70       return ret;
  71     }
  72   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 !");
  73 }
  74
  75 MCAuto<MEDFileEltStruct4Mesh> MEDFileBlowStrEltUp::dealWithMEDBALLInMesh(const MEDFileUMesh *mesh, MCAuto<MEDFileUMesh>& mOut, MCAuto<MEDFileFields>& fsOut) const
  76 {
  77   mOut=MEDFileUMesh::New(); fsOut=MEDFileFields::New();
  78   const std::vector< MCAuto<MEDFileEltStruct4Mesh> >& strs(mesh->getAccessOfUndergroundEltStrs());
  79   MCAuto<MEDFileEltStruct4Mesh> zeStr;
  80   for(std::vector< MCAuto<MEDFileEltStruct4Mesh> >::const_iterator it=strs.begin();it!=strs.end();it++)
  81     {
  82       if((*it)->getGeoTypeName()==MED_BALL_STR)
  83         {
  84           zeStr=*it;
  85           break;
  86         }
  87     }
  88   if(zeStr.isNull())
  89     {
  90       std::ostringstream oss; oss << "MEDFileBlowStrEltUp::dealWithMEDBALLInMesh : no geo type with name " <<  MED_BALL_STR << " in " << mesh->getName() << " !";
  91       throw INTERP_KERNEL::Exception(oss.str());
  92     }
  93   const DataArrayDouble *coo(mesh->getCoords());
  94   if(!coo)
  95     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLInMesh : null coords !");
  96   MCAuto<DataArrayInt> conn(zeStr->getConn());
  97   if(conn.isNull())
  98     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLInMesh : null connectivity !");
  99   conn->checkAllocated();
 100   if(conn->getNumberOfComponents()!=1)
 101     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLInMesh : excepted to be single compo !");
 102   int nbCells(conn->getNumberOfTuples());
 103   MCAuto<DataArrayDouble> connOut(coo->selectByTupleIdSafe(conn->begin(),conn->end()));
 104   MCAuto<MEDCouplingUMesh> mcOut(MEDCouplingUMesh::Build0DMeshFromCoords(connOut));
 105   mcOut->setName(BuildNewMeshName(mesh->getName(),MED_BALL_STR));
 106   mOut->setMeshAtLevel(0,mcOut);
 107   const DataArrayInt *ff1(mesh->getFamilyFieldAtLevel(1));
 108   if(ff1)
 109     {
 110       MCAuto<DataArrayInt> ff1o(ff1->selectByTupleIdSafe(conn->begin(),conn->end()));
 111       mOut->setFamilyFieldArr(1,ff1o);
 112     }
 113   const DataArrayInt *nf1(mesh->getNumberFieldAtLevel(1));
 114   if(nf1)
 115     {
 116       MCAuto<DataArrayInt> nf1o(nf1->selectByTupleIdSafe(conn->begin(),conn->end()));
 117       mOut->setRenumFieldArr(1,nf1o);
 118     }
 119   MCAuto<MEDFileUMeshPerTypeCommon> md(zeStr->getMeshDef());
 120   const DataArrayInt *ff0(md->getFam());
 121   if(ff0)
 122     mOut->setFamilyFieldArr(0,const_cast<DataArrayInt *>(ff0));
 123   const DataArrayInt *nf0(md->getNum());
 124   if(nf0)
 125     mOut->setRenumFieldArr(0,const_cast<DataArrayInt *>(nf0));
 126   mOut->copyFamGrpMapsFrom(*mesh);
 127   const std::vector< MCAuto<DataArray> >& vars(zeStr->getVars());
 128   for(std::vector< MCAuto<DataArray> >::const_iterator it=vars.begin();it!=vars.end();it++)
 129     {
 130       const DataArray *elt(*it);
 131       if(!elt)
 132         continue;
 133       {
 134         const DataArrayDouble *eltC(dynamic_cast<const DataArrayDouble *>(elt));
 135         if(eltC)
 136           {
 137             MCAuto<MEDFileFieldMultiTS> fmts(MEDFileFieldMultiTS::New());
 138             MCAuto<MEDFileField1TS> f1ts(MEDFileField1TS::New());
 139             MCAuto<MEDCouplingFieldDouble> f(MEDCouplingFieldDouble::New(ON_NODES));
 140             f->setMesh(mcOut);
 141             f->setArray(const_cast<DataArrayDouble *>(eltC));
 142             f->setName(eltC->getName());
 143             f1ts->setFieldNoProfileSBT(f);
 144             fmts->pushBackTimeStep(f1ts);
 145             fsOut->pushField(fmts);
 146           }
 147       }
 148     }
 149   return zeStr;
 150 }
 151
 152 /*!
 153  * \param [in] fs - fields lying all on same mesh and on same structure element
 154  * \param [in] zeStr - ze structure of current structure element
 155  * \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
 156  * \param [out] zeOutputs - ze fields that are the concatenation of fields in \a fs transformed and those in \a varAtt normalized in time space
 157  */
 158 void MEDFileBlowStrEltUp::dealWithSEInFields(const std::string& seName, const MEDFileFields *fs, const MEDFileEltStruct4Mesh *zeStr, const MEDFileFields *varAtt, MEDFileFields *zeOutputs) const
 159 {
 160   if(!fs)
 161     throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithSEInFields : null pointer !");
 162   if(seName==MED_BALL_STR)
 163     {
 164       dealWithMEDBALLSInFields(fs,zeStr,varAtt,zeOutputs);
 165       return ;
 166     }
 167   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 !");
 168 }
 169
 170 void MEDFileBlowStrEltUp::dealWithMEDBALLSInFields(const MEDFileFields *fs, const MEDFileEltStruct4Mesh *zeStr, const MEDFileFields *varAtt, MEDFileFields *zeOutputs) const
 171 {
 172   int nbf(fs->getNumberOfFields());
 173   std::vector< MCAuto<MEDFileAnyTypeFieldMultiTS> > elts0;
 174   std::vector< MEDFileAnyTypeFieldMultiTS * > elts1;
 175   std::string zeMeshName;
 176   for(int i=0;i<nbf;i++)
 177     {
 178       MCAuto<MEDFileAnyTypeFieldMultiTS> elt(fs->getFieldAtPos(i));
 179       MCAuto<MEDFileAnyTypeFieldMultiTS> eltOut(elt->buildNewEmpty());
 180       int nbTS(elt->getNumberOfTS());
 181       for(int j=0;j<nbTS;j++)
 182         {
 183           MCAuto<MEDFileAnyTypeField1TS> eltt(elt->getTimeStepAtPos(j));
 184           MCAuto<MEDFileAnyTypeField1TS> elttOut(eltt->deepCopy());
 185           std::string meshName(eltt->getMeshName());
 186           zeMeshName=BuildNewMeshName(meshName,MED_BALL_STR);
 187           elttOut->setMeshName(zeMeshName);
 188           elttOut->convertMedBallIntoClassic();
 189           eltOut->pushBackTimeStep(elttOut);
 190         }
 191       elts0.push_back(eltOut); elts1.push_back(eltOut);
 192     }
 193   //
 194   const MEDFileMesh *zeCurrentMesh(_ms->getMeshWithName(zeMeshName));
 195   //
 196   std::size_t ii(0);
 197   std::vector< std::vector<MEDFileAnyTypeFieldMultiTS *> > sp(MEDFileAnyTypeFieldMultiTS::SplitIntoCommonTimeSeries(elts1));
 198   for(std::vector< std::vector<MEDFileAnyTypeFieldMultiTS *> >::const_iterator it0=sp.begin();it0!=sp.end();it0++,ii++)
 199     {
 200       std::vector< MCAuto<MEDFileFastCellSupportComparator> > fsc;
 201       std::vector< std::vector<MEDFileAnyTypeFieldMultiTS *> > sp2(MEDFileAnyTypeFieldMultiTS::SplitPerCommonSupport(*it0,zeCurrentMesh,fsc));
 202       std::size_t jj(0);
 203       for(std::vector< std::vector<MEDFileAnyTypeFieldMultiTS *> >::const_iterator it1=sp2.begin();it1!=sp2.end();it1++,jj++)
 204         {
 205           for(std::vector<MEDFileAnyTypeFieldMultiTS *>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 206             zeOutputs->pushField(*it2);
 207           // The most exciting part. Users that put profiles on struct elements part of fields. Reduce var att.
 208           if((*it1).size()<1)
 209             throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLSInFields : take a deep breath !");
 210           MCAuto<MEDFileAnyTypeField1TS> zeGuideForPfl;// This var is the reference for pfl management.
 211           {
 212             if(!(*it1)[0])
 213               throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLSInFields : take a deep breath 2 !");
 214             int pdm((*it1)[0]->getNumberOfTS());
 215             if(pdm<1)
 216               throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLSInFields : take a deep breath 3 !");
 217             zeGuideForPfl=(*it1)[0]->getTimeStepAtPos(0);
 218           }
 219           if(zeGuideForPfl.isNull())
 220             throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLSInFields : take a deep breath 4 !");
 221           std::vector<std::string> pfls(zeGuideForPfl->getPflsReallyUsed());
 222           if(pfls.size()>=2)
 223             throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLSInFields : drink less coffee");
 224           MCAuto<DataArrayInt> pflMyLove;
 225           if(pfls.size()==1)
 226             pflMyLove.takeRef(zeGuideForPfl->getProfile(pfls[0]));
 227           // Yeah we have pfls
 228           std::vector<double> t2s;
 229           std::vector< std::pair<int,int> > t1s((*it1)[0]->getTimeSteps(t2s));
 230           std::size_t nbTS3(t2s.size());
 231           int nbf2(varAtt->getNumberOfFields());
 232           for(int i=0;i<nbf2;i++)
 233             {
 234               MCAuto<MEDFileAnyTypeFieldMultiTS> elt(varAtt->getFieldAtPos(i));
 235               int nbTS2(elt->getNumberOfTS());
 236               if(nbTS2!=1)
 237                 throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::dealWithMEDBALLSInFields : internal error ! The dealWithMEDBALLInMesh is expected to return a single TS !");
 238               MCAuto<MEDFileAnyTypeField1TS> elt2(elt->getTimeStepAtPos(0));
 239               MCAuto<MEDFileAnyTypeFieldMultiTS> elt4(elt->buildNewEmpty());
 240               for(std::size_t j=0;j<nbTS3;j++)
 241                 {
 242                   MCAuto<MEDFileAnyTypeField1TS> elt3(elt2->deepCopy());
 243                   elt3->setTime(t1s[j].first,t1s[j].second,t2s[j]);
 244                   elt3->setName(BuildVarAttName(ii,sp.size(),jj,sp2.size(),elt3->getName()));
 245                   if(pflMyLove.isNotNull())
 246                     elt3->makeReduction(INTERP_KERNEL::NORM_ERROR,ON_NODES,pflMyLove);
 247                   elt4->pushBackTimeStep(elt3);
 248                 }
 249               zeOutputs->pushField(elt4);
 250             }
 251         }
 252     }
 253 }
 254
 255 void MEDFileBlowStrEltUp::generate(MEDFileMeshes *msOut, MEDFileFields *allZeOutFields)
 256 {
 257   for(std::vector< MCAuto<MEDFileFields> >::iterator elt=_elts.begin();elt!=_elts.end();elt++)
 258     {
 259       std::vector< std::pair<std::string,std::string> > ps;
 260       (*elt)->getMeshSENames(ps);
 261       if(ps.size()!=1)
 262         throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::generateMeshes : internal error !");
 263       MEDFileMesh *mesh(_ms->getMeshWithName(ps[0].first));
 264       if(!mesh)
 265         throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::generateMeshes : NULL mesh !");
 266       MEDFileUMesh *umesh(dynamic_cast<MEDFileUMesh *>(mesh));
 267       if(!umesh)
 268         throw INTERP_KERNEL::Exception("MEDFileBlowStrEltUp::generateMeshes : Blow up of Stru Elt not managed yet for unstructured meshes !");
 269       //
 270       MCAuto<MEDFileFields> classicalSEFields(splitFieldsPerLoc(*elt,umesh,msOut,allZeOutFields));
 271       if(classicalSEFields.isNotNull())
 272         {
 273           MCAuto<MEDFileUMesh> mOut;
 274           MCAuto<MEDFileFields> fsOut1;
 275           MCAuto<MEDFileEltStruct4Mesh> zeStr(dealWithSEInMesh(ps[0].second,umesh,mOut,fsOut1));
 276           msOut->pushMesh(mOut);
 277           dealWithSEInFields(ps[0].second,classicalSEFields,zeStr,fsOut1,allZeOutFields);
 278         }
 279     }
 280 }
 281
 282 std::string MEDFileBlowStrEltUp::BuildNewMeshName(const std::string& meshName, const std::string& seName)
 283 {
 284   std::ostringstream mNameOut;
 285   mNameOut << meshName << "_" << seName;
 286   return mNameOut.str();
 287 }
 288
 289 std::string MEDFileBlowStrEltUp::BuildVarAttName(std::size_t iPart, std::size_t totINbParts, std::size_t jPart, std::size_t totJNbParts, const std::string& name)
 290 {
 291   if(totINbParts==1 && totJNbParts==1)
 292     return name;
 293   std::ostringstream oss;
 294   oss << name << "@" << iPart << "@" << jPart;
 295   return oss.str();
 296 }
 297
 298 void MEDFileBlowStrEltUp::DealWithSE(MEDFileFields *fs, MEDFileMeshes *ms, const MEDFileStructureElements *ses)
 299 {
 300   MCAuto<MEDFileFields> fsSEOnly(fs->partOfThisOnStructureElements());
 301   fs->killStructureElements();
 302   MEDFileBlowStrEltUp bu(fsSEOnly,ms,ses);
 303   bu.generate(ms,fs);
 304   fs->killStructureElementsInGlobs();
 305 }
 306
 307 //
 308
 309 class FieldWalker2
 310 {
 311 public:
 312   FieldWalker2(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd);
 313   std::string getLoc() const { return _loc; }
 314   std::string getPfl() const { return _pfl; }
 315   bool isClassic() const { return _is_classic; }
 316   bool operator!=(const FieldWalker2& other) const;
 317   bool operator==(const FieldWalker2& other) const;
 318 private:
 319   std::string _loc;
 320   std::string _pfl;
 321   bool _is_classic;
 322 };
 323
 324 class LocInfo
 325 {
 326 public:
 327   LocInfo() { }
 328   LocInfo(const std::vector<FieldWalker2>& fw);
 329   bool operator==(const LocInfo& other) const { return _locs==other._locs && _pfl==other._pfl; }
 330   void push(const std::string& loc, const std::string& pfl) { checkUniqueLoc(loc); _locs.push_back(loc); _pfl.push_back(pfl); }
 331   MCAuto<MEDFileUMesh> generateNonClassicalData(int zePos, const MEDFileUMesh *mesh, const MEDFileFieldGlobsReal *globs) const;
 332 private:
 333   void checkUniqueLoc(const std::string& loc) const;
 334   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);
 335   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);
 336 public:
 337   static const char ANGLE_DE_VRILLE[];
 338 private:
 339   std::vector<std::string> _locs;
 340   std::vector<std::string> _pfl;
 341 };
 342
 343 const char LocInfo::ANGLE_DE_VRILLE[]="ANGLE DE VRILLE";
 344
 345 LocInfo::LocInfo(const std::vector<FieldWalker2>& fw)
 346 {
 347   std::size_t sz(fw.size());
 348   _locs.resize(sz); _pfl.resize(sz);
 349   for(std::size_t i=0;i<sz;i++)
 350     {
 351       _locs[i]=fw[i].getLoc();
 352       _pfl[i]=fw[i].getPfl();
 353     }
 354 }
 355
 356 void LocInfo::checkUniqueLoc(const std::string& loc) const
 357 {
 358   if(std::find(_locs.begin(),_locs.end(),loc)!=_locs.end())
 359     {
 360       std::ostringstream oss; oss << "LocInfo::checkUniqueLoc : loc \"" << loc << "\" already exists !";
 361       throw INTERP_KERNEL::Exception(oss.str());
 362     }
 363 }
 364
 365 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)
 366 {
 367   MCAuto<DataArrayInt> conn(zeStr->getConn());
 368   conn=conn->deepCopy(); conn->rearrange(1);
 369   MCAuto<MEDCouplingUMesh> geoMesh;
 370   {
 371     MCAuto<MEDCoupling1SGTUMesh> umesh(MEDCoupling1SGTUMesh::New("",gt));
 372     umesh->setCoords(mesh->getCoords());
 373     umesh->setNodalConnectivity(conn);
 374     geoMesh=umesh->buildUnstructured();
 375   }
 376   //
 377   MCConstAuto<DataArrayDouble> angleVrille;
 378   if(!pfl.empty())
 379     {
 380       const DataArrayInt *pflArr(globs->getProfile(pfl));
 381       geoMesh=geoMesh->buildPartOfMySelf(pflArr->begin(),pflArr->end(),true);
 382       angleVrille=angleDeVrille->selectByTupleIdSafe(pflArr->begin(),pflArr->end());
 383     }
 384   else
 385     angleVrille.takeRef(angleDeVrille);
 386   //
 387   MCAuto<MEDCouplingFieldDouble> fakeF(MEDCouplingFieldDouble::New(ON_GAUSS_PT));
 388   fakeF->setMesh(geoMesh);
 389   int nbg(loc.getGaussWeights().size());
 390   fakeF->setGaussLocalizationOnType(gt,loc.getRefCoords(),loc.getGaussCoords(),loc.getGaussWeights());
 391   MCAuto<DataArrayDouble> ptsForLoc(fakeF->getLocalizationOfDiscr());
 392   //
 393   MCAuto<MEDCouplingFieldDouble> dir(geoMesh->buildDirectionVectorField());
 394   MCAuto<DataArrayDouble> rot(dir->getArray()->fromCartToSpher());
 395   int nbCells(geoMesh->getNumberOfCells()),nbCompo(ptsForLoc->getNumberOfComponents());
 396   MCAuto<DataArrayDouble> secPts(section->getCoords()->changeNbOfComponents(nbCompo,0.));
 397   int nbSecPts(secPts->getNumberOfTuples());
 398   {
 399     const int TAB[3]={2,0,1};
 400     std::vector<int> v(TAB,TAB+3);
 401     secPts=secPts->keepSelectedComponents(v);
 402   }
 403   const double CENTER[3]={0.,0.,0.},AX0[3]={0.,0.,1.};
 404   double AX1[3]; AX1[2]=0.;
 405   std::vector< MCAuto<DataArrayDouble> > arrs(nbCells*nbg);
 406   for(int j=0;j<nbCells;j++)
 407     {
 408       MCAuto<DataArrayDouble> p(secPts->deepCopy());
 409       double ang0(rot->getIJ(j,2));
 410       DataArrayDouble::Rotate3DAlg(CENTER,AX0,ang0,nbSecPts,p->begin(),p->getPointer());
 411       AX1[0]=-sin(ang0); AX1[1]=cos(ang0);// rot Oy around OZ
 412       double ang1(M_PI/2.-rot->getIJ(j,1));
 413       DataArrayDouble::Rotate3DAlg(CENTER,AX1,-ang1,nbSecPts,p->begin(),p->getPointer());
 414       DataArrayDouble::Rotate3DAlg(CENTER,dir->getArray()->begin()+j*3,angleVrille->getIJ(j,0),nbSecPts,p->begin(),p->getPointer());
 415       for(int l=0;l<nbg;l++)
 416         {
 417           MCAuto<DataArrayDouble> p2(p->deepCopy());
 418           for(int k=0;k<nbCompo;k++)
 419             p2->applyLin(1.,ptsForLoc->getIJ(j*nbg+l,k),k);
 420           arrs[j*nbg+l]=p2;
 421         }
 422     }
 423   std::vector<const DataArrayDouble *> arrs2(VecAutoToVecOfCstPt(arrs));
 424   MCAuto<DataArrayDouble> resu(DataArrayDouble::Aggregate(arrs2));
 425   return resu;
 426 }
 427
 428 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)
 429 {
 430   static const char MSG1[]="BuildMeshFromStructure : not recognized pattern ! Send mail to anthony.geay@edf.fr with corresponding MED file !";
 431   const std::vector< MCAuto<DataArray> >& vars(zeStr->getVars());
 432   if(vars.size()!=1)
 433     throw INTERP_KERNEL::Exception(MSG1);
 434   MCAuto<DataArray> zeArr(vars[0]);
 435   if(zeArr.isNull())
 436     throw INTERP_KERNEL::Exception(MSG1);
 437   MCAuto<DataArrayDouble> zeArr2(DynamicCast<DataArray,DataArrayDouble>(zeArr));
 438   if(zeArr2.isNull())
 439     throw INTERP_KERNEL::Exception(MSG1);
 440   if(zeArr2->getName()!=ANGLE_DE_VRILLE)
 441     throw INTERP_KERNEL::Exception(MSG1);
 442   return BuildMeshFromAngleVrille(gt,zeArr2,pfl,loc,zeStr,mesh,section,globs);
 443 }
 444
 445 MCAuto<MEDFileUMesh> LocInfo::generateNonClassicalData(int zePos, const MEDFileUMesh *mesh, const MEDFileFieldGlobsReal *globs) const
 446 {
 447   static const char MSG1[]="LocInfo::generateNonClassicalData : no spec for GAUSS on StructureElement with more than one cell !";
 448   std::size_t sz(_locs.size());
 449   std::vector< MCAuto<DataArrayDouble> > arrs(sz);
 450   for(std::size_t i=0;i<sz;i++)
 451     {
 452       const MEDFileFieldLoc& loc(globs->getLocalization(_locs[i]));
 453       const MEDFileGTKeeper *gtk(loc.getUndergroundGTKeeper());
 454       const MEDFileGTKeeperDyn *gtk2(dynamic_cast<const MEDFileGTKeeperDyn *>(gtk));
 455       if(!gtk2)
 456         throw INTERP_KERNEL::Exception("LocInfo::generateNonClassicalData : internal error !");
 457       const MEDFileUMesh *meshLoc(gtk2->getMesh()),*section(gtk2->getSection());
 458       const MEDFileStructureElement *se(gtk2->getSE());
 459       INTERP_KERNEL::NormalizedCellType gt;
 460       {
 461         std::vector<int> nel(meshLoc->getNonEmptyLevels());
 462         if(nel.size()!=1)
 463           throw INTERP_KERNEL::Exception(MSG1);
 464         if(nel[0]!=0)
 465           throw INTERP_KERNEL::Exception(MSG1);
 466         MCAuto<MEDCouplingUMesh> um(meshLoc->getMeshAtLevel(0));
 467         if(um->getNumberOfCells()!=1)
 468           throw INTERP_KERNEL::Exception(MSG1);
 469         gt=um->getTypeOfCell(0);
 470         std::vector<int> v;
 471         um->getNodeIdsOfCell(0,v);
 472         std::size_t sz2(v.size());
 473         for(std::size_t j=0;j<sz2;j++)
 474           if(v[j]!=j)
 475             throw INTERP_KERNEL::Exception(MSG1);
 476       }
 477       const std::vector< MCAuto<MEDFileEltStruct4Mesh> >& strs(mesh->getAccessOfUndergroundEltStrs());
 478       MCAuto<MEDFileEltStruct4Mesh> zeStr;
 479       for(std::vector< MCAuto<MEDFileEltStruct4Mesh> >::const_iterator it=strs.begin();it!=strs.end();it++)
 480         {
 481           if((*it)->getGeoTypeName()==se->getName())
 482             {
 483               zeStr=*it;
 484               break;
 485             }
 486         }
 487       if(zeStr.isNull())
 488         {
 489           std::ostringstream oss; oss << "LocInfo::generateNonClassicalData :  : no geo type with name " <<  se->getName() << " in " << mesh->getName() << " !";
 490           throw INTERP_KERNEL::Exception(oss.str());
 491         }
 492       arrs[i]=BuildMeshFromStructure(gt,_pfl[i],loc,zeStr,mesh,section,globs);
 493     }
 494   std::vector<const DataArrayDouble *> arrs2(VecAutoToVecOfCstPt(arrs));
 495   MCAuto<DataArrayDouble> resu(DataArrayDouble::Aggregate(arrs2));
 496   MCAuto<MEDFileUMesh> ret(MEDFileUMesh::New());
 497   ret->setCoords(resu);
 498   std::ostringstream meshName; meshName << mesh->getName() << "_on_" << sz << "_sections" << "_" << zePos;
 499   ret->setName(meshName.str());
 500   return ret;
 501 }
 502
 503 FieldWalker2::FieldWalker2(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd)
 504 {
 505   _loc=pmptpd->getLocalization();
 506   _pfl=pmptpd->getProfile();
 507   _is_classic=pmptpd->getType()!=ON_GAUSS_PT;
 508 }
 509
 510 bool FieldWalker2::operator!=(const FieldWalker2& other) const
 511 {
 512   bool ret(_loc==other._loc && _pfl==other._pfl && _is_classic==other._is_classic);
 513   return !ret;
 514 }
 515
 516 bool FieldWalker2::operator==(const FieldWalker2& other) const
 517 {
 518   bool ret(_loc==other._loc && _pfl==other._pfl && _is_classic==other._is_classic);
 519   return ret;
 520 }
 521
 522 class FieldWalker1
 523 {
 524 public:
 525   FieldWalker1(const MEDFileAnyTypeField1TSWithoutSDA *ts):_ts(ts),_pm_pt(0),_nb_mesh(0) { }
 526   void newMeshEntry(const MEDFileFieldPerMesh *fpm);
 527   void endMeshEntry(const MEDFileFieldPerMesh *fpm) { }
 528   void newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt);
 529   void endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt);
 530   void newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd);
 531   void checkOK(const FieldWalker1& other) const;
 532   bool isClassical() const;
 533   std::vector<FieldWalker2> getNonClassicalData() const { return _fw; }
 534 private:
 535   const MEDFileAnyTypeField1TSWithoutSDA *_ts;
 536   const MEDFileFieldPerMeshPerTypeDyn *_pm_pt;
 537   std::vector<FieldWalker2> _fw;
 538   int _nb_mesh;
 539 };
 540
 541 void FieldWalker1::newMeshEntry(const MEDFileFieldPerMesh *fpm)
 542 {
 543   if(_nb_mesh++==1)
 544     throw INTERP_KERNEL::Exception("FieldWalker1::newMeshEntry : multi mesh not supported !");
 545 }
 546
 547 void FieldWalker1::newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt)
 548 {
 549   if(_pm_pt)
 550     throw INTERP_KERNEL::Exception("FieldWalker1::newPerMeshPerTypeEntry : multi SE loc not managed yet !");
 551   const MEDFileFieldPerMeshPerTypeDyn *pmpt2(dynamic_cast<const MEDFileFieldPerMeshPerTypeDyn *>(pmpt));
 552   if(!pmpt2)
 553     throw INTERP_KERNEL::Exception("newPerMeshPerTypeEntry : internal error !");
 554   _pm_pt=pmpt2;
 555 }
 556
 557 void FieldWalker1::endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *)
 558 {
 559   isClassical();
 560 }
 561
 562 void FieldWalker1::newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd)
 563 {
 564   _fw.push_back(FieldWalker2(pmptpd));
 565 }
 566
 567 void FieldWalker1::checkOK(const FieldWalker1& other) const
 568 {
 569   std::size_t sz(_fw.size());
 570   if(other._fw.size()!=sz)
 571     throw INTERP_KERNEL::Exception("checkOK : not OK because size are not the same !");
 572   for(std::size_t i=0;i<sz;i++)
 573     if(_fw[i]!=other._fw[i])
 574       throw INTERP_KERNEL::Exception("checkOK : not OK because an element mismatches !");
 575 }
 576
 577 bool FieldWalker1::isClassical() const
 578 {
 579   if(_fw.empty())
 580     throw INTERP_KERNEL::Exception("FieldWalker1::endPerMeshPerTypeEntry : internal error !");
 581   std::size_t ic(0),inc(0);
 582   for(std::vector<FieldWalker2>::const_iterator it=_fw.begin();it!=_fw.end();it++)
 583     {
 584       if((*it).isClassic())
 585         ic++;
 586       else
 587         inc++;
 588     }
 589   if(ic!=0 && inc!=0)
 590     throw INTERP_KERNEL::Exception("FieldWalker1::endPerMeshPerTypeEntry : mix is not allowed yet !");
 591   return inc==0;
 592 }
 593
 594 class FieldWalker
 595 {
 596 public:
 597   FieldWalker(const MEDFileAnyTypeFieldMultiTSWithoutSDA *f):_f(f) { }
 598   void newTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts);
 599   void endTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts);
 600   void newMeshEntry(const MEDFileFieldPerMesh *fpm);
 601   void endMeshEntry(const MEDFileFieldPerMesh *fpm);
 602   void newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt);
 603   void endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt);
 604   void newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd);
 605 public:
 606   bool isEmpty() const;
 607   bool isClassical() const;
 608   const MEDFileAnyTypeFieldMultiTSWithoutSDA *field() const { return _f; }
 609   std::vector<FieldWalker2> getNonClassicalData() const { return _fw_prev->getNonClassicalData(); }
 610 private:
 611   const MEDFileAnyTypeFieldMultiTSWithoutSDA *_f;
 612   mutable INTERP_KERNEL::AutoCppPtr<FieldWalker1> _fw;
 613   mutable INTERP_KERNEL::AutoCppPtr<FieldWalker1> _fw_prev;
 614 };
 615
 616 bool FieldWalker::isEmpty() const
 617 {
 618   return _fw_prev.isNull();
 619 }
 620
 621 bool FieldWalker::isClassical() const
 622 {
 623   return _fw_prev->isClassical();
 624 }
 625
 626 void FieldWalker::newTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts)
 627 {
 628   _fw=new FieldWalker1(ts);
 629 }
 630
 631 void FieldWalker::endTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts)
 632 {
 633   if(_fw_prev.isNull())
 634     _fw_prev=new FieldWalker1(*_fw);
 635   else
 636     _fw_prev->checkOK(*_fw);
 637   _fw=0;
 638 }
 639
 640 void FieldWalker::newMeshEntry(const MEDFileFieldPerMesh *fpm)
 641 {
 642   _fw->newMeshEntry(fpm);
 643 }
 644
 645 void FieldWalker::endMeshEntry(const MEDFileFieldPerMesh *fpm)
 646 {
 647   _fw->endMeshEntry(fpm);
 648 }
 649
 650 void FieldWalker::newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt)
 651 {
 652   _fw->newPerMeshPerTypeEntry(pmpt);
 653 }
 654
 655 void FieldWalker::endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt)
 656 {
 657   _fw->endPerMeshPerTypeEntry(pmpt);
 658 }
 659
 660 void FieldWalker::newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd)
 661 {
 662   _fw->newPerMeshPerTypePerDisc(pmptpd);
 663 }
 664
 665 // this class splits fields into same
 666 class LocSpliter : public MEDFileFieldVisitor
 667 {
 668 public:
 669   LocSpliter(const MEDFileFieldGlobsReal *globs):_globs(globs),_fw(0) { }
 670   MCAuto<MEDFileFields> getClassical() const { return _classical; }
 671   void generateNonClassicalData(const MEDFileUMesh *mesh, std::vector< MCAuto<MEDFileFields> >& outFields, std::vector< MCAuto<MEDFileUMesh> >& outMeshes) const;
 672 private:
 673   void newFieldEntry(const MEDFileAnyTypeFieldMultiTSWithoutSDA *field);
 674   void endFieldEntry(const MEDFileAnyTypeFieldMultiTSWithoutSDA *field);
 675   //
 676   void newTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts);
 677   void endTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts);
 678   //
 679   void newMeshEntry(const MEDFileFieldPerMesh *fpm);
 680   void endMeshEntry(const MEDFileFieldPerMesh *fpm);
 681   //
 682   void newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt);
 683   void endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt);
 684   //
 685   void newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd);
 686 private:
 687   const MEDFileFieldGlobsReal *_globs;
 688   std::vector< LocInfo > _locs;
 689   std::vector< MCAuto<MEDFileFields> > _fields_on_locs;//size of _locs== size of _fields_on_locs
 690   MCAuto<MEDFileFields> _classical;
 691 private:
 692   mutable INTERP_KERNEL::AutoCppPtr<FieldWalker> _fw;
 693 };
 694
 695 void LocSpliter::newFieldEntry(const MEDFileAnyTypeFieldMultiTSWithoutSDA *field)
 696 {
 697   _fw=new FieldWalker(field);
 698 }
 699
 700 void LocSpliter::endFieldEntry(const MEDFileAnyTypeFieldMultiTSWithoutSDA *field)
 701 {
 702   if(_fw->isEmpty())
 703     return ;
 704   MCAuto<MEDFileAnyTypeFieldMultiTS> f(MEDFileAnyTypeFieldMultiTS::BuildNewInstanceFromContent(const_cast<MEDFileAnyTypeFieldMultiTSWithoutSDA *>(field)));
 705   if(_fw->isClassical())
 706     {
 707       if(_classical.isNull())
 708         {
 709           _classical=MEDFileFields::New();
 710           _classical->shallowCpyGlobs(*_globs);
 711         }
 712       _classical->pushField(f);
 713     }
 714   else
 715     {
 716       std::vector<FieldWalker2> fw2(_fw->getNonClassicalData());
 717       LocInfo elt(fw2);
 718       std::vector< LocInfo >::iterator it(std::find(_locs.begin(),_locs.end(),elt));
 719       if(it==_locs.end())
 720         {
 721           _locs.push_back(elt);
 722           MCAuto<MEDFileFields> zeF(MEDFileFields::New());
 723           zeF->shallowCpyGlobs(*_globs);
 724           zeF->pushField(f);
 725           _fields_on_locs.push_back(zeF);
 726         }
 727       else
 728         {
 729           MCAuto<MEDFileFields> zeF(_fields_on_locs[std::distance(_locs.begin(),it)]);
 730           zeF->pushField(f);
 731         }
 732     }
 733 }
 734
 735 void LocSpliter::generateNonClassicalData(const MEDFileUMesh *mesh, std::vector< MCAuto<MEDFileFields> >& outFields, std::vector< MCAuto<MEDFileUMesh> >& outMeshes) const
 736 {
 737   int i(0);
 738   for(std::vector<LocInfo>::const_iterator it=_locs.begin();it!=_locs.end();it++,i++)
 739     {
 740       MCAuto<MEDFileUMesh> m((*it).generateNonClassicalData(i,mesh,_globs));
 741       outMeshes.push_back(m);
 742       MCAuto<MEDCouplingUMesh> mcm(MEDCouplingUMesh::Build0DMeshFromCoords(m->getCoords()));
 743       mcm->setName(m->getName());
 744       MCAuto<MEDFileFields> fs(_fields_on_locs[i]);
 745       MCAuto<MEDFileFields> outFs(MEDFileFields::New());
 746       for(int j=0;j<fs->getNumberOfFields();j++)
 747         {
 748           MCAuto<MEDFileAnyTypeFieldMultiTS> fmtsNC(fs->getFieldAtPos(j));
 749           MCAuto<MEDFileFieldMultiTS> fmts(DynamicCastSafe<MEDFileAnyTypeFieldMultiTS,MEDFileFieldMultiTS>(fmtsNC));
 750           MCAuto<MEDFileFieldMultiTS> outFmts(MEDFileFieldMultiTS::New());
 751           for(int k=0;k<fmts->getNumberOfTS();k++)
 752             {
 753               MCAuto<MEDFileField1TS> outF1t(MEDFileField1TS::New());
 754               MCAuto<MEDFileField1TS> f1ts(fmts->getTimeStepAtPos(k));
 755               int t2,t3;
 756               double t1(f1ts->getTime(t2,t3));
 757               MCAuto<MEDCouplingFieldDouble> mcf(MEDCouplingFieldDouble::New(ON_NODES));
 758               mcf->setArray(f1ts->getUndergroundDataArray());
 759               mcf->setTime(t1,t2,t3);
 760               mcf->setName(f1ts->getName());
 761               mcf->setMesh(mcm);
 762               outF1t->setFieldNoProfileSBT(mcf);
 763               outFmts->pushBackTimeStep(outF1t);
 764             }
 765           outFs->pushField(outFmts);
 766         }
 767       outFields.push_back(outFs);
 768     }
 769 }
 770
 771 void LocSpliter::newTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts)
 772 {
 773   _fw->newTimeStepEntry(ts);
 774 }
 775
 776 void LocSpliter::endTimeStepEntry(const MEDFileAnyTypeField1TSWithoutSDA *ts)
 777 {
 778   _fw->endTimeStepEntry(ts);
 779 }
 780
 781 void LocSpliter::newMeshEntry(const MEDFileFieldPerMesh *fpm)
 782 {
 783   _fw->newMeshEntry(fpm);
 784 }
 785
 786 void LocSpliter::endMeshEntry(const MEDFileFieldPerMesh *fpm)
 787 {
 788   _fw->endMeshEntry(fpm);
 789 }
 790
 791 void LocSpliter::newPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt)
 792 {
 793   _fw->newPerMeshPerTypeEntry(pmpt);
 794 }
 795
 796 void LocSpliter::endPerMeshPerTypeEntry(const MEDFileFieldPerMeshPerTypeCommon *pmpt)
 797 {
 798   _fw->endPerMeshPerTypeEntry(pmpt);
 799 }
 800
 801 void LocSpliter::newPerMeshPerTypePerDisc(const MEDFileFieldPerMeshPerTypePerDisc *pmptpd)
 802 {
 803   _fw->newPerMeshPerTypePerDisc(pmptpd);
 804 }
 805
 806 void MEDFileBlowStrEltUp::DealWithConflictNames(MEDFileAnyTypeFieldMultiTS *fmtsToAdd, const MEDFileFields *fs)
 807 {
 808   std::vector<std::string> fnames(fs->getFieldsNames());
 809   for(int i=0;i<1000;i++)
 810     {
 811       std::ostringstream oss; oss << fmtsToAdd->getName();
 812       if(i>=1)
 813         oss << "_" << i-1;
 814       if(std::find(fnames.begin(),fnames.end(),oss.str())==fnames.end())
 815         {
 816           fmtsToAdd->setName(oss.str());
 817           return ;
 818         }
 819     }
 820   throw INTERP_KERNEL::Exception("DealWithConflictNames : Eh eh interesting !");
 821 }
 822
 823 MCAuto<MEDFileFields> MEDFileBlowStrEltUp::splitFieldsPerLoc(const MEDFileFields *fields, const MEDFileUMesh *mesh, MEDFileMeshes *msOut, MEDFileFields *allZeOutFields)
 824 {
 825   LocSpliter ls(fields);
 826   fields->accept(ls);
 827   std::vector< MCAuto<MEDFileFields> > outFields;
 828   std::vector< MCAuto<MEDFileUMesh> > outMeshes;
 829   ls.generateNonClassicalData(mesh,outFields,outMeshes);
 830   for(std::vector< MCAuto<MEDFileFields> >::iterator it=outFields.begin();it!=outFields.end();it++)
 831     {
 832       for(int j=0;j<(*it)->getNumberOfFields();j++)
 833         {
 834           MCAuto<MEDFileAnyTypeFieldMultiTS> fmts((*it)->getFieldAtPos(j));
 835           DealWithConflictNames(fmts,allZeOutFields);
 836           allZeOutFields->pushField(fmts);
 837         }
 838     }
 839   for(std::vector< MCAuto<MEDFileUMesh> >::iterator it=outMeshes.begin();it!=outMeshes.end();it++)
 840     msOut->pushMesh(*it);
 841   return ls.getClassical();
 842 }