src/MEDWrapper/Base/MED_Algorithm.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 //
  23 #include "MED_Algorithm.hxx"
  24 #include "MED_Wrapper.hxx"
  25
  26 #include "MED_Utilities.hxx"
  27
  28 #ifdef _DEBUG_
  29 static int MYDEBUG = 0;
  30 static int MYVALUEDEBUG = 0;
  31 #else
  32 // static int MYDEBUG = 0;
  33 // static int MYVALUEDEBUG = 0;
  34 #endif
  35
  36 namespace MED
  37 {
  38   //---------------------------------------------------------------
  39   TEntity2TGeom2ElemInfo
  40   GetEntity2TGeom2ElemInfo(const PWrapper& theWrapper,
  41                            const PMeshInfo& theMeshInfo,
  42                            const MED::TEntityInfo& theEntityInfo)
  43   {
  44     MSG(MYDEBUG,"GetElemsByEntity(...)");
  45     TEntity2TGeom2ElemInfo anEntity2TGeom2ElemInfo;
  46     MED::TEntityInfo::const_iterator anIter = theEntityInfo.begin();
  47     PElemInfo anElemInfo;
  48     TErr anErr;
  49     for(; anIter != theEntityInfo.end(); anIter++){
  50       const EEntiteMaillage& anEntity = anIter->first;
  51       const TGeom2Size& aGeom2Size = anIter->second;
  52       TGeom2ElemInfo& aGeom2ElemInfo = anEntity2TGeom2ElemInfo[anEntity];
  53
  54       if(anEntity == eNOEUD){
  55         aGeom2ElemInfo[ePOINT1] = theWrapper->GetPElemInfo(theMeshInfo);
  56         continue;
  57       }
  58
  59       TGeom2Size::const_iterator anIter2 = aGeom2Size.begin();
  60       for(; anIter2 != aGeom2Size.end(); anIter2++){
  61         const EGeometrieElement& aGeom = anIter2->first;
  62         aGeom2ElemInfo[aGeom] = theWrapper->GetPElemInfo(theMeshInfo,anEntity,aGeom,MED::eNOD,&anErr);
  63       }
  64     }
  65     ADDMSG(MYDEBUG,"\n");
  66     return anEntity2TGeom2ElemInfo;
  67   }
  68
  69
  70   //---------------------------------------------------------------
  71   TFamilyInfoSet
  72   GetFamilyInfoSet(const PWrapper& theWrapper,
  73                    const PMeshInfo& theMeshInfo)
  74   {
  75     MSG(MYDEBUG,"GetFamilies(...)");
  76     TErr anErr;
  77     TFamilyInfoSet aFamilyInfoSet;
  78     TInt aNbFam = theWrapper->GetNbFamilies(*theMeshInfo);
  79     INITMSG(MYDEBUG,"GetNbFamilies() = "<<aNbFam<<"\n");
  80     for(TInt iFam = 1; iFam <= aNbFam; iFam++){
  81       PFamilyInfo aFamilyInfo = theWrapper->GetPFamilyInfo(theMeshInfo,iFam,&anErr);
  82       if(anErr >= 0)
  83         aFamilyInfoSet.insert(aFamilyInfo);
  84     }
  85     ADDMSG(MYDEBUG,"\n");
  86     return aFamilyInfoSet;
  87   }
  88
  89
  90   //---------------------------------------------------------------
  91   TGroupInfo
  92   GetGroupInfo(const TFamilyInfoSet& theFamilyInfoSet)
  93   {
  94     MSG(MYDEBUG,"GetFamiliesByGroup(...)");
  95     TGroupInfo aGroup;
  96     TFamilyInfoSet::const_iterator anIter = theFamilyInfoSet.begin();
  97     for(; anIter != theFamilyInfoSet.end(); anIter++){
  98       const PFamilyInfo& aFamilyInfo = *anIter;
  99       TInt aNbGroup = aFamilyInfo->GetNbGroup();
 100       for(TInt iGroup = 0; iGroup < aNbGroup; iGroup++){
 101         aGroup[aFamilyInfo->GetGroupName(iGroup)].insert(aFamilyInfo);
 102       }
 103     }
 104
 105 #ifdef _DEBUG_
 106     if(MYDEBUG){
 107       TGroupInfo::const_iterator anIter = aGroup.begin();
 108       for(; anIter != aGroup.end(); anIter++){
 109         const std::string& aName = anIter->first;
 110         INITMSG(MYDEBUG,"aGroupName = '"<<aName<<"'\n");
 111         const TFamilyInfoSet& aFamilyInfoSet = anIter->second;
 112         TFamilyInfoSet::const_iterator anFamIter = aFamilyInfoSet.begin();
 113         for(; anFamIter != aFamilyInfoSet.end(); anFamIter++){
 114           const PFamilyInfo& aFamilyInfo = *anFamIter;
 115           INITMSG(MYDEBUG,"aFamilyName = '"<<aFamilyInfo->GetName()<<"'\n");
 116         }
 117       }
 118       ADDMSG(MYDEBUG,"\n");
 119     }
 120 #endif
 121
 122     return aGroup;
 123   }
 124
 125
 126   //---------------------------------------------------------------
 127   TFieldInfo2TimeStampInfoSet
 128   GetFieldInfo2TimeStampInfoSet(const PWrapper& theWrapper,
 129                                 const PMeshInfo& theMeshInfo,
 130                                 const MED::TEntityInfo& theEntityInfo)
 131   {
 132     MSG(MYDEBUG,"GetFieldsByEntity(...)");
 133     TFieldInfo2TimeStampInfoSet aFieldInfo2TimeStampInfoSet;
 134     TInt aNbFields = theWrapper->GetNbFields();
 135     INITMSG(MYDEBUG,"GetNbFields() = "<<aNbFields<<"\n");
 136     for(TInt iField = 1; iField <= aNbFields; iField++){
 137       PFieldInfo aFieldInfo = theWrapper->GetPFieldInfo(theMeshInfo,iField);
 138       INITMSG(MYDEBUG,"aFieldName = '"<<aFieldInfo->GetName()<<
 139               "'; aNbComp = "<<aFieldInfo->GetNbComp()<<"; ");
 140       TGeom2Size aGeom2Size;
 141       EEntiteMaillage anEntity = EEntiteMaillage(-1);
 142       TInt aNbTimeStamps = theWrapper->GetNbTimeStamps(aFieldInfo,theEntityInfo,anEntity,aGeom2Size);
 143       ADDMSG(MYDEBUG,"anEntity = "<<anEntity<<"; GetNbTimeStamps = "<<aNbTimeStamps<<"\n");
 144       for(TInt iTimeStamp = 1; iTimeStamp <= aNbTimeStamps; iTimeStamp++){
 145         PTimeStampInfo aTimeStamp =
 146           theWrapper->GetPTimeStampInfo(aFieldInfo,anEntity,aGeom2Size,iTimeStamp);
 147         aFieldInfo2TimeStampInfoSet[aFieldInfo].insert(aTimeStamp);
 148         INITMSG(MYDEBUG,
 149                 "aDt = "<<aTimeStamp->GetDt()<<
 150                 ", Unit = \'"<<aTimeStamp->GetUnitDt()<<"\n");
 151       }
 152     }
 153     ADDMSG(MYDEBUG,"\n");
 154     return aFieldInfo2TimeStampInfoSet;
 155   }
 156
 157
 158   //---------------------------------------------------------------
 159   TEntite2TFieldInfo2TimeStampInfoSet
 160   GetEntite2TFieldInfo2TimeStampInfoSet(const TFieldInfo2TimeStampInfoSet& theFieldInfo2TimeStampInfoSet)
 161   {
 162     TEntite2TFieldInfo2TimeStampInfoSet anEntite2TFieldInfo2TimeStampInfoSet;
 163     TFieldInfo2TimeStampInfoSet::const_iterator anIter = theFieldInfo2TimeStampInfoSet.begin();
 164     for(; anIter != theFieldInfo2TimeStampInfoSet.end(); anIter++){
 165       const TTimeStampInfoSet& aTimeStampInfoSet = anIter->second;
 166       //const PFieldInfo& aFieldInfo = anIter->first;
 167       if(aTimeStampInfoSet.empty())
 168         continue;
 169       const PTimeStampInfo& aTimeStampInfo = *aTimeStampInfoSet.begin();
 170       anEntite2TFieldInfo2TimeStampInfoSet[ConvertEntity(aTimeStampInfo->GetEntity())].insert(*anIter);
 171     }
 172     return anEntite2TFieldInfo2TimeStampInfoSet;
 173   }
 174
 175
 176   //---------------------------------------------------------------
 177   bool
 178   operator<(const TFamilyTSize& theLeft, const TFamilyTSize& theRight)
 179   {
 180     const MED::PFamilyInfo& aLeftInfo = boost::get<0>(theLeft);
 181     const MED::PFamilyInfo& aRightInfo = boost::get<0>(theRight);
 182     return aLeftInfo->GetId() < aRightInfo->GetId();
 183   }
 184
 185
 186   //---------------------------------------------------------------
 187   TEntity2FamilySet
 188   GetEntity2FamilySet(const PWrapper& theWrapper,
 189                       const TEntity2TGeom2ElemInfo& theEntity2TGeom2ElemInfo,
 190                       const TFamilyInfoSet& theFamilyInfoSet)
 191   {
 192     MSG(MYDEBUG,"GetFamiliesByEntity(...)");
 193     TEntity2FamilySet anEntity2FamilySet;
 194
 195     typedef std::map<TInt,PFamilyInfo> TId2Family;
 196     TId2Family anId2Family;
 197     TFamilyInfoSet::const_iterator anIter = theFamilyInfoSet.begin();
 198     for(; anIter != theFamilyInfoSet.end(); anIter++){
 199       const PFamilyInfo& aFamilyInfo = *anIter;
 200       anId2Family.insert(TId2Family::value_type(aFamilyInfo->GetId(),aFamilyInfo));
 201     }
 202
 203     if(!anId2Family.empty()){
 204       typedef std::map<TInt,TInt> TFamilyID2Size;
 205       typedef std::map<EEntiteMaillage,TFamilyID2Size> TEntity2FamilyID;
 206       TEntity2FamilyID anEntity2FamilyID;
 207
 208       if(!theEntity2TGeom2ElemInfo.empty()){
 209         TEntity2TGeom2ElemInfo::const_iterator anIter = theEntity2TGeom2ElemInfo.begin();
 210         for(; anIter != theEntity2TGeom2ElemInfo.end(); anIter++){
 211           const EEntiteMaillage& anEntity = anIter->first;
 212           TFamilyID2Size& aFamilyID2Size = anEntity2FamilyID[anEntity];
 213           const TGeom2ElemInfo& aGeom2ElemInfo = anIter->second;
 214           TGeom2ElemInfo::const_iterator aGeom2ElemInfoIter = aGeom2ElemInfo.begin();
 215           for(; aGeom2ElemInfoIter != aGeom2ElemInfo.end(); aGeom2ElemInfoIter++){
 216             const PElemInfo& aElemInfo = aGeom2ElemInfoIter->second;
 217             if(TInt aNbElem = aElemInfo->GetNbElem()){
 218               for(TInt i = 0; i < aNbElem; i++){
 219                 aFamilyID2Size[aElemInfo->GetFamNum(i)] += 1;
 220               }
 221             }
 222           }
 223         }
 224       }
 225
 226       if(!anEntity2FamilyID.empty()){
 227         TEntity2FamilyID::const_iterator anIter = anEntity2FamilyID.begin();
 228         for(; anIter != anEntity2FamilyID.end(); anIter++){
 229           const EEntiteMaillage& anEntity = anIter->first;
 230           INITMSG(MYDEBUG,"anEntity = "<<anEntity<<":\n");
 231           const TFamilyID2Size& aFamilyID2Size = anIter->second;
 232           TFamilyID2Size::const_iterator anIter2 = aFamilyID2Size.begin();
 233           for(; anIter2 != aFamilyID2Size.end(); anIter2++){
 234             TInt anId = anIter2->first;
 235             TInt aSize = anIter2->second;
 236             TId2Family::const_iterator anIter3 = anId2Family.find(anId);
 237             if(anIter3 != anId2Family.end()){
 238               const PFamilyInfo& aFamilyInfo = anIter3->second;
 239               anEntity2FamilySet[anEntity].insert(TFamilyTSize(aFamilyInfo,aSize));
 240               INITMSG(MYDEBUG,
 241                       "aFamilyName = '"<<aFamilyInfo->GetName()<<
 242                       "' anId = "<<aFamilyInfo->GetId()<<"\n");
 243             }
 244           }
 245         }
 246       }
 247     }
 248     ADDMSG(MYDEBUG,"\n");
 249     return anEntity2FamilySet;
 250   }
 251
 252
 253   //---------------------------------------------------------------
 254   TKey2Gauss
 255   GetKey2Gauss(const PWrapper& theWrapper,
 256                TErr* theErr,
 257                EModeSwitch theMode)
 258   {
 259     INITMSG(MYDEBUG,"GetKey2Gauss - theMode = "<<theMode<<std::endl);
 260     TKey2Gauss aKey2Gauss;
 261     TInt aNbGauss = theWrapper->GetNbGauss(theErr);
 262     for(TInt anId = 1; anId <= aNbGauss; anId++){
 263       TGaussInfo::TInfo aPreInfo = theWrapper->GetGaussPreInfo(anId);
 264       PGaussInfo anInfo = theWrapper->CrGaussInfo(aPreInfo,theMode);
 265       theWrapper->GetGaussInfo(anId,anInfo,theErr);
 266       TGaussInfo::TKey aKey = boost::get<0>(aPreInfo);
 267       aKey2Gauss[aKey] = anInfo;
 268
 269 #ifdef _DEBUG_
 270       const EGeometrieElement& aGeom = boost::get<0>(aKey);
 271       const std::string& aName = boost::get<1>(aKey);
 272       INITMSG(MYDEBUG,
 273               "- aGeom = "<<aGeom<<
 274               "; aName = '"<<aName<<"'"<<
 275               std::endl);
 276 #endif
 277
 278     }
 279     return aKey2Gauss;
 280   }
 281
 282
 283   //---------------------------------------------------------------
 284   PProfileInfo
 285   GetProfileInfo(const PWrapper& theWrapper,
 286                  const std::string& theProfileName,
 287                  TErr* theErr,
 288                  EModeProfil theMode)
 289   {
 290     PProfileInfo anInfo;
 291     TInt aNbProfiles = theWrapper->GetNbProfiles(theErr);
 292     for(TInt anId = 1; anId <= aNbProfiles; anId++){
 293       TProfileInfo::TInfo aPreInfo = theWrapper->GetProfilePreInfo(anId);
 294       const std::string& aName = boost::get<0>(aPreInfo);
 295       if(aName == theProfileName)
 296         return theWrapper->GetPProfileInfo(anId,theMode,theErr);
 297     }
 298     return anInfo;
 299   }
 300
 301
 302   //---------------------------------------------------------------
 303   TMKey2Profile
 304   GetMKey2Profile(const PWrapper& theWrapper,
 305                   TErr* theErr,
 306                   EModeProfil theMode)
 307   {
 308     INITMSG(MYDEBUG,"GetMKey2Profile - theMode = "<<theMode<<std::endl);
 309     TKey2Profile aKey2Profile;
 310     TInt aNbProfiles = theWrapper->GetNbProfiles(theErr);
 311     for(TInt anId = 1; anId <= aNbProfiles; anId++){
 312       TProfileInfo::TInfo aPreInfo = theWrapper->GetProfilePreInfo(anId);
 313       PProfileInfo anInfo = theWrapper->GetPProfileInfo(anId,theMode,theErr);
 314       const std::string& aName = boost::get<0>(aPreInfo);
 315       aKey2Profile[aName] = anInfo;
 316
 317 #ifdef _DEBUG_
 318       INITMSG(MYDEBUG,
 319               "- aName = '"<<aName<<"'"<<
 320               " : "<<
 321               std::endl);
 322       TInt aNbElem = anInfo->GetSize();
 323       for(TInt iElem = 0; iElem < aNbElem; iElem++){
 324         ADDMSG(MYVALUEDEBUG,anInfo->GetElemNum(iElem)<<", ");
 325       }
 326       ADDMSG(MYVALUEDEBUG, std::endl);
 327 #endif
 328
 329     }
 330     return TMKey2Profile(theMode,aKey2Profile);
 331   }
 332
 333   //---------------------------------------------------------------
 334   EEntiteMaillage
 335   GetEntityByFamilyId(PGrilleInfo& theInfo,TInt theId){
 336     TElemNum::iterator aNodeFamIter = (theInfo->myFamNumNode).begin();
 337     for(;aNodeFamIter != (theInfo->myFamNumNode).end(); aNodeFamIter++){
 338       if(theId == *aNodeFamIter)
 339         return eNOEUD;
 340     }
 341     TElemNum::iterator aCellFamIter = (theInfo->myFamNum).begin();
 342     for(;aCellFamIter != (theInfo->myFamNum).end(); aCellFamIter++){
 343       if(theId == *aCellFamIter)
 344         return eMAILLE;
 345     }
 346     EXCEPTION(std::runtime_error, "GetEntityByFamilyId - fails");
 347     return EEntiteMaillage(-1);
 348   }
 349
 350   TFamilyID2NbCells
 351   GetFamilyID2NbCells(PGrilleInfo& theInfo){
 352     TFamilyID2NbCells aFamily2NbCells;
 353     TInt aNbNodes = theInfo->myFamNumNode.size();
 354     TInt aNbCells = theInfo->myFamNum.size();
 355     for(TInt i=0; i<aNbNodes; i++) aFamily2NbCells[theInfo->GetFamNumNode(i)] = 0;
 356     for(TInt i=0; i<aNbCells; i++) aFamily2NbCells[theInfo->GetFamNum(i)] = 0;
 357     for(TInt i=0; i<aNbNodes; i++) aFamily2NbCells[theInfo->GetFamNumNode(i)] += 1;
 358     for(TInt i=0; i<aNbCells; i++) aFamily2NbCells[theInfo->GetFamNum(i)] += 1;
 359     return aFamily2NbCells;
 360   }
 361
 362   EEntiteMaillage ConvertEntity(const EEntiteMaillage& aEntity){
 363     switch( aEntity ){
 364
 365     case eNOEUD_ELEMENT:
 366     case eMAILLE: return eMAILLE; //eNOEUD_ELEMENT it is eMAILLE
 367
 368     case eFACE:
 369     case eARETE:
 370     case eNOEUD: return aEntity; break;
 371     default: return EEntiteMaillage(-1);
 372
 373     }
 374   }
 375 }