1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #include "MEDPARTITIONER_ParallelTopology.hxx"
21 #include "MEDPARTITIONER_MeshCollectionDriver.hxx"
22 #include "MEDPARTITIONER_MeshCollection.hxx"
23 #include "MEDPARTITIONER_ParaDomainSelector.hxx"
24 #include "MEDPARTITIONER_Utils.hxx"
26 #include "MEDCouplingUMesh.hxx"
27 #include "MEDCouplingFieldDouble.hxx"
28 #include "MEDLoader.hxx"
29 #include "MEDFileMesh.hxx"
38 #include <libxml/tree.h>
39 #include <libxml/parser.h>
40 #include <libxml/xpath.h>
41 #include <libxml/xpathInternals.h>
45 using namespace MEDPARTITIONER;
47 MeshCollectionDriver::MeshCollectionDriver(MeshCollection* collection):_collection(collection)
51 /*!reads a unique MED File v>=2.1
52 * and mounts the corresponding mesh in memory
53 *\param filename binary file
54 *\param meshname mesh name in the MED file
56 int MeshCollectionDriver::readSeq(const char* filename, const char* meshname)
58 std::cout << "readSeq" << std::endl;
59 MyGlobals::_File_Names.resize(1);
60 MyGlobals::_File_Names[0]=std::string(filename);
62 ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(filename,meshname);
63 //puts the only mesh in the mesh vector
64 (_collection->getMesh()).push_back(mfm->getLevel0Mesh(false));
65 (_collection->getFaceMesh()).push_back(mfm->getLevelM1Mesh(false));
68 ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy());
69 ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy());
70 (_collection->getCellFamilyIds()).push_back(cellIds);
71 (_collection->getFaceFamilyIds()).push_back(faceIds);
74 (_collection->getFamilyInfo())=mfm->getFamilyInfo();
75 (_collection->getGroupInfo())=mfm->getGroupInfo();
77 (_collection->getCZ()).clear();
79 ParallelTopology* aPT = new ParallelTopology((_collection->getMesh()));
80 _collection->setTopology(aPT);
81 _collection->setName(meshname);
82 _collection->setDomainNames(meshname);
87 //================================================================================
89 * \brief Return mesh dimension from distributed med file had being read
91 //================================================================================
93 void MeshCollectionDriver::readSubdomain(std::vector<int*>& cellglobal,
94 std::vector<int*>& faceglobal,
95 std::vector<int*>& nodeglobal, int idomain)
97 std::string meshname=MyGlobals::_Mesh_Names[idomain];
98 std::string file=MyGlobals::_File_Names[idomain];
100 ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(file.c_str(),meshname.c_str());
101 std::vector<int> nonEmpty=mfm->getNonEmptyLevels();
105 (_collection->getMesh())[idomain]=mfm->getLevel0Mesh(false);
106 //reading families groups
107 ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy());
108 (_collection->getCellFamilyIds())[idomain]=cellIds;
112 (_collection->getMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want tests;
113 ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
115 (_collection->getCellFamilyIds())[idomain]=empty;
116 std::cout << "\nNO Level0Mesh (Cells)\n";
120 if (nonEmpty.size()>1 && nonEmpty[1]==-1)
122 (_collection->getFaceMesh())[idomain]=mfm->getLevelM1Mesh(false);
123 //reading families groups
124 ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy());
125 (_collection->getFaceFamilyIds())[idomain]=faceIds;
126 if (MyGlobals::_Verbose>10)
127 std::cout << "proc " << MyGlobals::_Rank << " : WITH Faces\n";
132 throw INTERP_KERNEL::Exception("no faces");
137 (_collection->getFaceMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want test;
138 ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
139 (_collection->getFaceFamilyIds())[idomain]=empty;
140 if (MyGlobals::_Verbose>10)
141 std::cout << "proc " << MyGlobals::_Rank << " : WITHOUT Faces\n";
145 _collection->getFamilyInfo()=mfm->getFamilyInfo();
146 _collection->getGroupInfo()=mfm->getGroupInfo();
150 std::vector<std::string> localInformation;
152 localInformation.push_back(str+"ioldDomain="+IntToStr(idomain));
153 localInformation.push_back(str+"meshName="+meshname);
154 MyGlobals::_General_Informations.push_back(SerializeFromVectorOfString(localInformation));
155 std::vector<std::string> localFields=BrowseAllFieldsOnMesh(file, meshname, idomain);
156 if (localFields.size()>0)
157 MyGlobals::_Field_Descriptions.push_back(SerializeFromVectorOfString(localFields));
161 void MeshCollectionDriver::readSubdomain(int idomain)
163 std::string meshname=MyGlobals::_Mesh_Names[idomain];
164 std::string file=MyGlobals::_File_Names[idomain];
166 ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(file.c_str(),meshname.c_str());
167 std::vector<int> nonEmpty=mfm->getNonEmptyLevels();
171 (_collection->getMesh())[idomain]=mfm->getLevel0Mesh(false);
172 //reading families groups
173 ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy());
174 (_collection->getCellFamilyIds())[idomain]=cellIds;
178 (_collection->getMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want tests;
179 ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
181 (_collection->getCellFamilyIds())[idomain]=empty;
182 std::cout<<"\nNO Level0Mesh (Cells)\n";
186 if (nonEmpty.size()>1 && nonEmpty[1]==-1)
188 (_collection->getFaceMesh())[idomain]=mfm->getLevelM1Mesh(false);
189 //reading families groups
190 ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy());
191 (_collection->getFaceFamilyIds())[idomain]=faceIds;
192 if (MyGlobals::_Verbose>10)
193 std::cout << "proc " << MyGlobals::_Rank << " : WITH Faces\n";
197 throw INTERP_KERNEL::Exception("no faces");
202 (_collection->getFaceMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want test;
203 ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
204 (_collection->getFaceFamilyIds())[idomain]=empty;
205 if (MyGlobals::_Verbose>10)
206 std::cout << "proc " << MyGlobals::_Rank << " : WITHOUT Faces\n";
210 _collection->getFamilyInfo()=mfm->getFamilyInfo();
211 _collection->getGroupInfo()=mfm->getGroupInfo();
215 std::vector<std::string> localInformation;
217 localInformation.push_back(str+"ioldDomain="+IntToStr(idomain));
218 localInformation.push_back(str+"meshName="+meshname);
219 MyGlobals::_General_Informations.push_back(SerializeFromVectorOfString(localInformation));
220 std::vector<std::string> localFields=BrowseAllFieldsOnMesh(file, meshname, idomain);
221 if (localFields.size()>0)
222 MyGlobals::_Field_Descriptions.push_back(SerializeFromVectorOfString(localFields));
226 void MeshCollectionDriver::writeMedFile(int idomain, const std::string& distfilename) const
228 std::vector<const ParaMEDMEM::MEDCouplingUMesh*> meshes;
229 ParaMEDMEM::MEDCouplingUMesh* cellMesh=_collection->getMesh(idomain);
230 ParaMEDMEM::MEDCouplingUMesh* faceMesh=_collection->getFaceMesh(idomain);
231 //ParaMEDMEM::MEDCouplingUMesh* faceMeshFilter=0;
233 std::string finalMeshName=ExtractFromDescription(MyGlobals::_General_Informations[0], "finalMeshName=");
234 // std::string cleFilter=Cle1ToStr("filterFaceOnCell",idomain);
235 // ParaMEDMEM::DataArrayInt* filter=0;
236 // if (_collection->getMapDataArrayInt().find(cleFilter)!=_collection->getMapDataArrayInt().end())
238 // filter=_collection->getMapDataArrayInt().find(cleFilter)->second;
239 // int* index=filter->getPointer();
240 // faceMeshFilter=(ParaMEDMEM::MEDCouplingUMesh *) faceMesh->buildPartOfMySelf(index,index+filter->getNbOfElems(),true);
241 // faceMesh=faceMeshFilter;
243 cellMesh->setName(finalMeshName.c_str());
244 meshes.push_back(cellMesh);
246 faceMesh->checkCoherency();
247 if (faceMesh->getNumberOfCells()>0)
249 faceMesh->tryToShareSameCoordsPermute(*cellMesh, 1e-10);
250 meshes.push_back(faceMesh);
253 //ParaMEDMEM::MEDCouplingUMesh* boundaryMesh=0;
254 // if (MyGlobals::_Creates_Boundary_Faces>0)
256 // //try to write Boundary meshes
257 // bool keepCoords=false; //TODO or true
258 // boundaryMesh=(ParaMEDMEM::MEDCouplingUMesh *) cellMesh->buildBoundaryMesh(keepCoords);
259 // boundaryMesh->setName("boundaryMesh");
262 MEDLoader::WriteUMeshes(distfilename.c_str(), meshes, true);
263 // if (faceMeshFilter!=0)
264 // faceMeshFilter->decrRef();
266 // if (boundaryMesh!=0)
268 // //doing that testMesh becomes second mesh sorted by alphabetical order of name
269 // MEDLoader::WriteUMesh(distfilename.c_str(), boundaryMesh, false);
270 // boundaryMesh->decrRef();
272 ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(distfilename.c_str(), _collection->getMesh(idomain)->getName().c_str());
274 mfm->setFamilyInfo(_collection->getFamilyInfo());
275 mfm->setGroupInfo(_collection->getGroupInfo());
277 std::string key=Cle1ToStr("faceFamily_toArray",idomain);
278 if ( meshes.size() == 2 &&
279 _collection->getMapDataArrayInt().find(key)!=_collection->getMapDataArrayInt().end())
281 ParaMEDMEM::DataArrayInt *fam=_collection->getMapDataArrayInt().find(key)->second;
282 mfm->setFamilyFieldArr(-1,fam);
285 key=Cle1ToStr("cellFamily_toArray",idomain);
286 if (_collection->getMapDataArrayInt().find(key)!=_collection->getMapDataArrayInt().end())
287 mfm->setFamilyFieldArr(0,_collection->getMapDataArrayInt().find(key)->second);
289 mfm->write(distfilename.c_str(),0);
290 key="/inewFieldDouble="+IntToStr(idomain)+"/";
292 std::map<std::string,ParaMEDMEM::DataArrayDouble*>::iterator it;
294 for (it=_collection->getMapDataArrayDouble().begin() ; it!=_collection->getMapDataArrayDouble().end(); it++)
296 std::string desc=(*it).first;
297 size_t found=desc.find(key);
298 if (found==std::string::npos)
300 if (MyGlobals::_Verbose>20)
301 std::cout << "proc " << MyGlobals::_Rank << " : write field " << desc << std::endl;
302 std::string meshName, fieldName;
303 int typeField, DT, IT, entity;
304 FieldShortDescriptionToData(desc, fieldName, typeField, entity, DT, IT);
305 double time=StrToDouble(ExtractFromDescription(desc, "time="));
306 int typeData=StrToInt(ExtractFromDescription(desc, "typeData="));
307 std::string entityName=ExtractFromDescription(desc, "entityName=");
308 ParaMEDMEM::MEDCouplingFieldDouble* field=0;
311 std::cout << "WARNING : writeMedFile : typeData " << typeData << " not implemented for fields\n";
314 if (entityName=="MED_CELL")
316 //there is a field of idomain to write
317 field=ParaMEDMEM::MEDCouplingFieldDouble::New(ParaMEDMEM::ON_CELLS,ParaMEDMEM::ONE_TIME);
319 if (entityName=="MED_NODE_ELEMENT")
321 //there is a field of idomain to write
322 field=ParaMEDMEM::MEDCouplingFieldDouble::New(ParaMEDMEM::ON_GAUSS_NE,ParaMEDMEM::ONE_TIME);
326 std::cout << "WARNING : writeMedFile : entityName " << entityName << " not implemented for fields\n";
330 field->setName(fieldName.c_str());
331 field->setMesh(mfm->getLevel0Mesh(false));
332 ParaMEDMEM::DataArrayDouble *da=(*it).second;
334 //get information for components etc..
335 std::vector<std::string> r1;
336 r1=SelectTagsInVectorOfString(MyGlobals::_General_Informations,"fieldName="+fieldName);
337 r1=SelectTagsInVectorOfString(r1,"typeField="+IntToStr(typeField));
338 r1=SelectTagsInVectorOfString(r1,"DT="+IntToStr(DT));
339 r1=SelectTagsInVectorOfString(r1,"IT="+IntToStr(IT));
340 //not saved in file? field->setDescription(ExtractFromDescription(r1[0], "fieldDescription=").c_str());
341 int nbc=StrToInt(ExtractFromDescription(r1[0], "nbComponents="));
342 if (nbc==da->getNumberOfComponents())
344 for (int i=0; i<nbc; i++)
345 da->setInfoOnComponent(i,ExtractFromDescription(r1[0], "componentInfo"+IntToStr(i)+"=").c_str());
349 std::cerr << "Problem On field " << fieldName << " : number of components unexpected " << da->getNumberOfComponents() << std::endl;
353 field->setTime(time,DT,IT);
354 field->checkCoherency();
357 MEDLoader::WriteField(distfilename.c_str(),field,false);
359 catch(INTERP_KERNEL::Exception& e)
361 //cout trying rewrite all data, only one field defined
362 std::string tmp,newName=distfilename;
363 tmp+="_"+fieldName+"_"+IntToStr(nbfFieldFound)+".med";
364 newName.replace(newName.find(".med"),4,tmp);
365 std::cout << "WARNING : writeMedFile : create a new file name with only one field because MEDLoader::WriteField throw:" << newName << std::endl;
366 MEDLoader::WriteField(newName.c_str(),field,true);