Back to index

salome-med  6.5.0
MEDPARTITIONER_MeshCollectionDriver.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "MEDPARTITIONER_ParallelTopology.hxx"
00021 #include "MEDPARTITIONER_MeshCollectionDriver.hxx"
00022 #include "MEDPARTITIONER_MeshCollection.hxx"
00023 #include "MEDPARTITIONER_ParaDomainSelector.hxx"
00024 #include "MEDPARTITIONER_Utils.hxx"
00025 
00026 #include "MEDCouplingUMesh.hxx"
00027 #include "MEDCouplingFieldDouble.hxx"
00028 #include "MEDLoader.hxx"
00029 #include "MEDFileMesh.hxx"
00030 
00031 #include <map>
00032 #include <set>
00033 #include <vector>
00034 #include <string>
00035 #include <fstream>
00036 #include <iostream>
00037 
00038 #include <libxml/tree.h>
00039 #include <libxml/parser.h>
00040 #include <libxml/xpath.h>
00041 #include <libxml/xpathInternals.h>
00042 
00043 #include "med.h"
00044 
00045 using namespace MEDPARTITIONER;
00046 
00047 MeshCollectionDriver::MeshCollectionDriver(MeshCollection* collection):_collection(collection)
00048 {
00049 }
00050 
00056 int MeshCollectionDriver::readSeq(const char* filename, const char* meshname)
00057 {
00058   std::cout << "readSeq" << std::endl;
00059   MyGlobals::_File_Names.resize(1);
00060   MyGlobals::_File_Names[0]=std::string(filename);
00061 
00062   ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(filename,meshname);
00063   //puts the only mesh in the mesh vector
00064   (_collection->getMesh()).push_back(mfm->getLevel0Mesh(false));
00065   (_collection->getFaceMesh()).push_back(mfm->getLevelM1Mesh(false));
00066 
00067   //reading family ids
00068   ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy());
00069   ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy());
00070   (_collection->getCellFamilyIds()).push_back(cellIds);
00071   (_collection->getFaceFamilyIds()).push_back(faceIds); 
00072 
00073   //reading groups
00074   (_collection->getFamilyInfo())=mfm->getFamilyInfo();
00075   (_collection->getGroupInfo())=mfm->getGroupInfo();
00076 
00077   (_collection->getCZ()).clear();
00078   
00079   ParallelTopology* aPT = new ParallelTopology((_collection->getMesh()));
00080   _collection->setTopology(aPT);
00081   _collection->setName(meshname);
00082   _collection->setDomainNames(meshname);
00083   return 0;
00084 }
00085 
00086 
00087 //================================================================================
00091 //================================================================================
00092 
00093 void MeshCollectionDriver::readSubdomain(std::vector<int*>& cellglobal,
00094                                          std::vector<int*>& faceglobal,
00095                                          std::vector<int*>& nodeglobal, int idomain)
00096 {
00097   std::string meshname=MyGlobals::_Mesh_Names[idomain];
00098   std::string file=MyGlobals::_File_Names[idomain];
00099 
00100   ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(file.c_str(),meshname.c_str());
00101   std::vector<int> nonEmpty=mfm->getNonEmptyLevels();
00102   
00103   try 
00104     { 
00105       (_collection->getMesh())[idomain]=mfm->getLevel0Mesh(false); 
00106       //reading families groups
00107       ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy());
00108       (_collection->getCellFamilyIds())[idomain]=cellIds;
00109     }
00110   catch(...)
00111     { 
00112       (_collection->getMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want tests;
00113       ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
00114       empty->alloc(0,1);
00115       (_collection->getCellFamilyIds())[idomain]=empty;
00116       std::cout << "\nNO Level0Mesh (Cells)\n";
00117     }
00118   try 
00119     { 
00120       if (nonEmpty.size()>1 && nonEmpty[1]==-1)
00121         {
00122           (_collection->getFaceMesh())[idomain]=mfm->getLevelM1Mesh(false);
00123           //reading families groups
00124           ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy());
00125           (_collection->getFaceFamilyIds())[idomain]=faceIds;
00126           if (MyGlobals::_Verbose>10)
00127             std::cout << "proc " << MyGlobals::_Rank << " : WITH Faces\n";
00128 
00129         }
00130       else
00131         {
00132           throw INTERP_KERNEL::Exception("no faces");
00133         }
00134     }
00135   catch(...)
00136     {
00137       (_collection->getFaceMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want test;
00138       ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
00139       (_collection->getFaceFamilyIds())[idomain]=empty;
00140       if (MyGlobals::_Verbose>10)
00141         std::cout << "proc " << MyGlobals::_Rank << " : WITHOUT Faces\n";
00142     }
00143   
00144   //reading groups
00145   _collection->getFamilyInfo()=mfm->getFamilyInfo();
00146   _collection->getGroupInfo()=mfm->getGroupInfo();
00147 
00148   mfm->decrRef();
00149   
00150   std::vector<std::string> localInformation;
00151   std::string str;
00152   localInformation.push_back(str+"ioldDomain="+IntToStr(idomain));
00153   localInformation.push_back(str+"meshName="+meshname);
00154   MyGlobals::_General_Informations.push_back(SerializeFromVectorOfString(localInformation));
00155   std::vector<std::string> localFields=BrowseAllFieldsOnMesh(file, meshname, idomain);
00156   if (localFields.size()>0) 
00157     MyGlobals::_Field_Descriptions.push_back(SerializeFromVectorOfString(localFields));
00158 }
00159 
00160 
00161 void MeshCollectionDriver::readSubdomain(int idomain)
00162 {
00163   std::string meshname=MyGlobals::_Mesh_Names[idomain];
00164   std::string file=MyGlobals::_File_Names[idomain];
00165 
00166   ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(file.c_str(),meshname.c_str());
00167   std::vector<int> nonEmpty=mfm->getNonEmptyLevels();
00168   
00169   try 
00170     { 
00171       (_collection->getMesh())[idomain]=mfm->getLevel0Mesh(false); 
00172       //reading families groups
00173       ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy());
00174       (_collection->getCellFamilyIds())[idomain]=cellIds;
00175     }
00176   catch(...)
00177     { 
00178       (_collection->getMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want tests;
00179       ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
00180       empty->alloc(0,1);
00181       (_collection->getCellFamilyIds())[idomain]=empty;
00182       std::cout<<"\nNO Level0Mesh (Cells)\n";
00183     }
00184   try 
00185     { 
00186       if (nonEmpty.size()>1 && nonEmpty[1]==-1)
00187         {
00188           (_collection->getFaceMesh())[idomain]=mfm->getLevelM1Mesh(false);
00189           //reading families groups
00190           ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy());
00191           (_collection->getFaceFamilyIds())[idomain]=faceIds;
00192           if (MyGlobals::_Verbose>10)
00193             std::cout << "proc " << MyGlobals::_Rank << " : WITH Faces\n";
00194         }
00195       else
00196         {
00197           throw INTERP_KERNEL::Exception("no faces");
00198         }
00199     }
00200   catch(...)
00201     {
00202       (_collection->getFaceMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want test;
00203       ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
00204       (_collection->getFaceFamilyIds())[idomain]=empty;
00205       if (MyGlobals::_Verbose>10)
00206         std::cout << "proc " << MyGlobals::_Rank << " : WITHOUT Faces\n";
00207     }
00208   
00209   //reading groups
00210   _collection->getFamilyInfo()=mfm->getFamilyInfo();
00211   _collection->getGroupInfo()=mfm->getGroupInfo();
00212 
00213   mfm->decrRef();
00214   
00215   std::vector<std::string> localInformation;
00216   std::string str;
00217   localInformation.push_back(str+"ioldDomain="+IntToStr(idomain));
00218   localInformation.push_back(str+"meshName="+meshname);
00219   MyGlobals::_General_Informations.push_back(SerializeFromVectorOfString(localInformation));
00220   std::vector<std::string> localFields=BrowseAllFieldsOnMesh(file, meshname, idomain);
00221   if (localFields.size()>0) 
00222     MyGlobals::_Field_Descriptions.push_back(SerializeFromVectorOfString(localFields));
00223 }
00224 
00225 
00226 void MeshCollectionDriver::writeMedFile(int idomain, const std::string& distfilename) const
00227 {
00228   std::vector<const ParaMEDMEM::MEDCouplingUMesh*> meshes;
00229   ParaMEDMEM::MEDCouplingUMesh* cellMesh=_collection->getMesh(idomain);
00230   ParaMEDMEM::MEDCouplingUMesh* faceMesh=_collection->getFaceMesh(idomain);
00231   ParaMEDMEM::MEDCouplingUMesh* faceMeshFilter=0;
00232   
00233   std::string finalMeshName=ExtractFromDescription(MyGlobals::_General_Informations[0], "finalMeshName=");
00234   std::string cleFilter=Cle1ToStr("filterFaceOnCell",idomain);
00235   ParaMEDMEM::DataArrayInt* filter=0;
00236   if (_collection->getMapDataArrayInt().find(cleFilter)!=_collection->getMapDataArrayInt().end())
00237     {
00238       filter=_collection->getMapDataArrayInt().find(cleFilter)->second;
00239       int* index=filter->getPointer();
00240       faceMeshFilter=(ParaMEDMEM::MEDCouplingUMesh *) faceMesh->buildPartOfMySelf(index,index+filter->getNbOfElems(),true);
00241       faceMesh=faceMeshFilter;
00242     }
00243   cellMesh->setName(finalMeshName.c_str());
00244   meshes.push_back(cellMesh);
00245   
00246   faceMesh->checkCoherency();
00247   if (faceMesh->getNumberOfCells()>0)
00248     {
00249       faceMesh->tryToShareSameCoordsPermute(*cellMesh, 1e-10);
00250       meshes.push_back(faceMesh);
00251     }
00252   
00253   ParaMEDMEM::MEDCouplingUMesh* boundaryMesh=0;
00254   if (MyGlobals::_Creates_Boundary_Faces>0)
00255     {
00256       //try to write Boundary meshes
00257       bool keepCoords=false; //TODO or true
00258       boundaryMesh=(ParaMEDMEM::MEDCouplingUMesh *) cellMesh->buildBoundaryMesh(keepCoords);
00259       boundaryMesh->setName("boundaryMesh");
00260     }
00261   
00262   MEDLoader::WriteUMeshes(distfilename.c_str(), meshes, true);
00263   if (faceMeshFilter!=0)
00264     faceMeshFilter->decrRef();
00265   
00266   if (boundaryMesh!=0)
00267     {
00268       //doing that testMesh becomes second mesh sorted by alphabetical order of name
00269       MEDLoader::WriteUMesh(distfilename.c_str(), boundaryMesh, false);
00270       boundaryMesh->decrRef();
00271     }
00272   ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(distfilename.c_str(), _collection->getMesh(idomain)->getName());
00273   
00274   mfm->setFamilyInfo(_collection->getFamilyInfo());
00275   mfm->setGroupInfo(_collection->getGroupInfo());
00276   
00277   std::string key=Cle1ToStr("faceFamily_toArray",idomain);
00278   if (_collection->getMapDataArrayInt().find(key)!=_collection->getMapDataArrayInt().end())
00279     {
00280       ParaMEDMEM::DataArrayInt *fam=_collection->getMapDataArrayInt().find(key)->second;
00281       ParaMEDMEM::DataArrayInt *famFilter=0;
00282       if (filter!=0)
00283         {
00284           int* index=filter->getPointer();
00285           int nbTuples=filter->getNbOfElems();
00286           //not the good one...buildPartOfMySelf do not exist for DataArray 
00287           //Filter=fam->renumberAndReduce(index, filter->getNbOfElems());
00288           famFilter=ParaMEDMEM::DataArrayInt::New();
00289           famFilter->alloc(nbTuples,1);
00290           int* pfamFilter=famFilter->getPointer();
00291           int* pfam=fam->getPointer();
00292           for (int i=0; i<nbTuples; i++)
00293             pfamFilter[i]=pfam[index[i]];
00294           fam=famFilter;
00295           mfm->setFamilyFieldArr(-1,fam);
00296           famFilter->decrRef();
00297         }
00298     }
00299   
00300   key=Cle1ToStr("cellFamily_toArray",idomain);
00301   if (_collection->getMapDataArrayInt().find(key)!=_collection->getMapDataArrayInt().end())
00302     mfm->setFamilyFieldArr(0,_collection->getMapDataArrayInt().find(key)->second);
00303   
00304   mfm->write(distfilename.c_str(),0);
00305   key="/inewFieldDouble="+IntToStr(idomain)+"/";
00306     
00307   std::map<std::string,ParaMEDMEM::DataArrayDouble*>::iterator it;
00308   int nbfFieldFound=0;
00309   for (it=_collection->getMapDataArrayDouble().begin() ; it!=_collection->getMapDataArrayDouble().end(); it++)
00310     {
00311       std::string desc=(*it).first;
00312       size_t found=desc.find(key);
00313       if (found==std::string::npos)
00314         continue;
00315       if (MyGlobals::_Verbose>20)
00316         std::cout << "proc " << MyGlobals::_Rank << " : write field " << desc << std::endl;
00317       std::string meshName, fieldName;
00318       int typeField, DT, IT, entity;
00319       FieldShortDescriptionToData(desc, fieldName, typeField, entity, DT, IT);
00320       double time=StrToDouble(ExtractFromDescription(desc, "time="));
00321       int typeData=StrToInt(ExtractFromDescription(desc, "typeData="));
00322       std::string entityName=ExtractFromDescription(desc, "entityName=");
00323       ParaMEDMEM::MEDCouplingFieldDouble* field=0;
00324       if (typeData!=6)
00325         {
00326           std::cout << "WARNING : writeMedFile : typeData " << typeData << " not implemented for fields\n";
00327           continue;
00328         }
00329       if (entityName=="MED_CELL")
00330         {
00331           //there is a field of idomain to write
00332           field=ParaMEDMEM::MEDCouplingFieldDouble::New(ParaMEDMEM::ON_CELLS,ParaMEDMEM::ONE_TIME);
00333         }
00334       if (entityName=="MED_NODE_ELEMENT")
00335         {
00336           //there is a field of idomain to write
00337           field=ParaMEDMEM::MEDCouplingFieldDouble::New(ParaMEDMEM::ON_GAUSS_NE,ParaMEDMEM::ONE_TIME);
00338         }
00339       if (!field)
00340         {
00341           std::cout << "WARNING : writeMedFile : entityName " << entityName << " not implemented for fields\n";
00342           continue;
00343         }
00344       nbfFieldFound++;
00345       field->setName(fieldName.c_str());
00346       field->setMesh(mfm->getLevel0Mesh(false));
00347       ParaMEDMEM::DataArrayDouble *da=(*it).second;
00348     
00349       //get information for components etc..
00350       std::vector<std::string> r1;
00351       r1=SelectTagsInVectorOfString(MyGlobals::_General_Informations,"fieldName="+fieldName);
00352       r1=SelectTagsInVectorOfString(r1,"typeField="+IntToStr(typeField));
00353       r1=SelectTagsInVectorOfString(r1,"DT="+IntToStr(DT));
00354       r1=SelectTagsInVectorOfString(r1,"IT="+IntToStr(IT));
00355       //not saved in file? field->setDescription(ExtractFromDescription(r1[0], "fieldDescription=").c_str());
00356       int nbc=StrToInt(ExtractFromDescription(r1[0], "nbComponents="));
00357       if (nbc==da->getNumberOfComponents())
00358         {
00359           for (int i=0; i<nbc; i++) 
00360             da->setInfoOnComponent(i,ExtractFromDescription(r1[0], "componentInfo"+IntToStr(i)+"=").c_str());
00361         }
00362       else
00363         {
00364           std::cerr << "Problem On field " << fieldName << " : number of components unexpected " << da->getNumberOfComponents() << std::endl;
00365         }
00366     
00367       field->setArray(da);
00368       field->setTime(time,DT,IT);
00369       field->checkCoherency();
00370       try
00371         {
00372           MEDLoader::WriteField(distfilename.c_str(),field,false);
00373         }
00374       catch(INTERP_KERNEL::Exception& e)
00375         {
00376           //cout trying rewrite all data, only one field defined
00377           std::string tmp,newName=distfilename;
00378           tmp+="_"+fieldName+"_"+IntToStr(nbfFieldFound)+".med";
00379           newName.replace(newName.find(".med"),4,tmp);
00380           std::cout << "WARNING : writeMedFile : create a new file name with only one field because MEDLoader::WriteField throw:" << newName << std::endl;
00381           MEDLoader::WriteField(newName.c_str(),field,true);
00382         }
00383     }
00384   mfm->decrRef();
00385 }