1 // Copyright (C) 2007-2015 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, or (at your option) any later version.
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_MeshCollectionDriver.hxx"
22 #include "MEDPARTITIONER_ConnectZone.hxx"
23 #include "MEDPARTITIONER_MeshCollection.hxx"
24 #include "MEDPARTITIONER_ParaDomainSelector.hxx"
25 #include "MEDPARTITIONER_ParallelTopology.hxx"
26 #include "MEDPARTITIONER_Utils.hxx"
28 #include "MEDCouplingFieldDouble.hxx"
29 #include "MEDCouplingRefCountObject.hxx"
30 #include "MEDCouplingUMesh.hxx"
31 #include "MEDFileField.hxx"
32 #include "MEDFileJoint.hxx"
33 #include "MEDFileMesh.hxx"
34 #include "MEDLoader.hxx"
43 #include <libxml/tree.h>
44 #include <libxml/parser.h>
45 #include <libxml/xpath.h>
46 #include <libxml/xpathInternals.h>
50 using namespace MEDPARTITIONER;
52 MeshCollectionDriver::MeshCollectionDriver(MeshCollection* collection):_collection(collection)
56 /*!reads a unique MED File v>=2.1
57 * and mounts the corresponding mesh in memory
58 *\param filename binary file
59 *\param meshname mesh name in the MED file
61 int MeshCollectionDriver::readSeq(const char* filename, const char* meshname)
63 std::cout << "readSeq" << std::endl;
64 MyGlobals::_File_Names.resize(1);
65 MyGlobals::_File_Names[0]=std::string(filename);
67 ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(filename,meshname);
68 //puts the only mesh in the mesh vector
69 (_collection->getMesh()).push_back(mfm->getLevel0Mesh(false));
70 (_collection->getFaceMesh()).push_back(mfm->getLevelM1Mesh(false));
73 ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy());
74 ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy());
75 (_collection->getCellFamilyIds()).push_back(cellIds);
76 (_collection->getFaceFamilyIds()).push_back(faceIds);
79 (_collection->getFamilyInfo())=mfm->getFamilyInfo();
80 (_collection->getGroupInfo())=mfm->getGroupInfo();
82 (_collection->getCZ()).clear();
84 ParallelTopology* aPT = new ParallelTopology((_collection->getMesh()));
85 _collection->setTopology(aPT, true);
86 _collection->setName(meshname);
87 _collection->setDomainNames(meshname);
92 void MeshCollectionDriver::readMEDFileData(const ParaMEDMEM::MEDFileData* filedata)
94 const int nbDomains = filedata->getMeshes()->getNumberOfMeshes();
95 _collection->getMesh() .resize( nbDomains, 0 );
96 _collection->getFaceMesh() .resize( nbDomains, 0 );
97 _collection->getCellFamilyIds().resize( nbDomains, 0 );
98 _collection->getFaceFamilyIds().resize( nbDomains, 0 );
100 for (int i=0; i<nbDomains; i++)
102 ParaMEDMEM::MEDFileUMesh *mfm = dynamic_cast<ParaMEDMEM::MEDFileUMesh *>(filedata->getMeshes()->getMeshAtPos(i));
104 if ( mfm && mfm->getMeshDimension() > 0 )
105 _collection->setNonEmptyMesh( i );
108 ParallelTopology* aPT = new ParallelTopology(_collection->getMesh());
109 _collection->setTopology(aPT, true);
112 _collection->setName( filedata->getMeshes()->getMeshAtPos(0)->getName() );
113 _collection->setDomainNames( _collection->getName() );
115 if ( ParaDomainSelector* domainSelector = _collection->getParaDomainSelector() )
116 if ( _collection->isParallelMode() )
118 //to know nb of cells on each proc to compute global cell ids from locally global
119 domainSelector->gatherNbOf(_collection->getMesh());
123 void MeshCollectionDriver::readFileData(std::string file,std::string meshname,int idomain) const
125 ParaMEDMEM::MEDFileUMesh* mfm=ParaMEDMEM::MEDFileUMesh::New(file,meshname);
126 readData(mfm,idomain);
130 void MeshCollectionDriver::readData(ParaMEDMEM::MEDFileUMesh* mfm, int idomain) const
132 std::vector<int> nonEmpty=mfm->getNonEmptyLevels();
135 (_collection->getMesh())[idomain]=mfm->getLevel0Mesh(false);
136 //reading families groups
137 ParaMEDMEM::DataArrayInt* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCpy());
138 (_collection->getCellFamilyIds())[idomain]=cellIds;
142 (_collection->getMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want tests;
143 ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
145 (_collection->getCellFamilyIds())[idomain]=empty;
146 std::cout<<"\nNO Level0Mesh (Cells)\n";
150 if (nonEmpty.size()>1 && nonEmpty[1]==-1)
152 (_collection->getFaceMesh())[idomain]=mfm->getLevelM1Mesh(false);
153 //reading families groups
154 ParaMEDMEM::DataArrayInt* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCpy());
155 (_collection->getFaceFamilyIds())[idomain]=faceIds;
156 if (MyGlobals::_Verbose>10)
157 std::cout << "proc " << MyGlobals::_Rank << " : WITH Faces\n";
161 throw INTERP_KERNEL::Exception("no faces");
166 (_collection->getFaceMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want test;
167 ParaMEDMEM::DataArrayInt* empty=ParaMEDMEM::DataArrayInt::New();
168 (_collection->getFaceFamilyIds())[idomain]=empty;
169 if (MyGlobals::_Verbose>10)
170 std::cout << "proc " << MyGlobals::_Rank << " : WITHOUT Faces\n";
173 _collection->getFamilyInfo()=mfm->getFamilyInfo();
174 _collection->getGroupInfo()=mfm->getGroupInfo();
177 void MeshCollectionDriver::readSubdomain(int idomain)
179 std::string meshname=MyGlobals::_Mesh_Names[idomain];
180 std::string file=MyGlobals::_File_Names[idomain];
181 readFileData(file,meshname,idomain);
183 std::vector<std::string> localInformation;
185 localInformation.push_back(str+"ioldDomain="+IntToStr(idomain));
186 localInformation.push_back(str+"meshName="+meshname);
187 MyGlobals::_General_Informations.push_back(SerializeFromVectorOfString(localInformation));
188 std::vector<std::string> localFields=BrowseAllFieldsOnMesh(file, meshname, idomain);
189 if (localFields.size()>0)
190 MyGlobals::_Field_Descriptions.push_back(SerializeFromVectorOfString(localFields));
193 ParaMEDMEM::MEDFileMesh* MeshCollectionDriver::getMesh(int idomain) const
195 ParaMEDMEM::MEDFileUMesh* mfm = ParaMEDMEM::MEDFileUMesh::New();
197 ParaMEDMEM::MEDCouplingUMesh* cellMesh=_collection->getMesh(idomain);
198 ParaMEDMEM::MEDCouplingUMesh* faceMesh=_collection->getFaceMesh(idomain);
199 // std::string cleFilter=Cle1ToStr("filterFaceOnCell",idomain);
200 // ParaMEDMEM::DataArrayInt* filter=0;
201 // if (_collection->getMapDataArrayInt().find(cleFilter)!=_collection->getMapDataArrayInt().end())
203 // filter=_collection->getMapDataArrayInt().find(cleFilter)->second;
204 // int* index=filter->getPointer();
205 // faceMeshFilter=(ParaMEDMEM::MEDCouplingUMesh *) faceMesh->buildPartOfMySelf(index,index+filter->getNbOfElems(),true);
206 // faceMesh=faceMeshFilter;
208 // if (faceMeshFilter!=0)
209 // faceMeshFilter->decrRef();
210 std::string finalMeshName="";
211 if (MyGlobals::_General_Informations.size()!=0)
213 std::size_t found=MyGlobals::_General_Informations[0].find("finalMeshName=");
214 if ((found!=std::string::npos) && (found>0))
216 finalMeshName=ExtractFromDescription(MyGlobals::_General_Informations[0], "finalMeshName=");
219 if (finalMeshName.empty())
221 finalMeshName=_collection->getName();
223 cellMesh->setName(finalMeshName);
224 mfm->setMeshAtLevel( 0, cellMesh );
226 faceMesh->checkCoherency();
227 if (faceMesh->getNumberOfCells()>0)
229 faceMesh->tryToShareSameCoordsPermute(*cellMesh, 1e-10);
230 faceMesh->setName(finalMeshName);
231 mfm->setMeshAtLevel( -1, faceMesh );
234 // ParaMEDMEM::MEDCouplingUMesh* boundaryMesh=0;
235 // if (MyGlobals::_Creates_Boundary_Faces>0)
237 // //try to write Boundary meshes
238 // bool keepCoords=false; //TODO or true
239 // boundaryMesh=(ParaMEDMEM::MEDCouplingUMesh *) cellMesh->buildBoundaryMesh(keepCoords);
240 // boundaryMesh->setName("boundaryMesh");
241 // if (boundaryMesh!=0)
243 // //doing that testMesh becomes second mesh sorted by alphabetical order of name
244 // MEDLoader::WriteUMesh(distfilename, boundaryMesh, false);
245 // boundaryMesh->decrRef();
248 mfm->setFamilyInfo(_collection->getFamilyInfo());
249 mfm->setGroupInfo(_collection->getGroupInfo());
250 std::string key=Cle1ToStr("faceFamily_toArray",idomain);
251 if ( faceMesh->getNumberOfCells()>0 && _collection->getMapDataArrayInt().find(key)!=_collection->getMapDataArrayInt().end())
252 mfm->setFamilyFieldArr(-1,_collection->getMapDataArrayInt().find(key)->second);
253 key=Cle1ToStr("cellFamily_toArray",idomain);
254 if (_collection->getMapDataArrayInt().find(key)!=_collection->getMapDataArrayInt().end())
255 mfm->setFamilyFieldArr(0,_collection->getMapDataArrayInt().find(key)->second);
259 using ParaMEDMEM::MEDCouplingAutoRefCountObjectPtr;
260 using ParaMEDMEM::MEDCouplingSkyLineArray;
261 using ParaMEDMEM::MEDFileJoint;
262 using ParaMEDMEM::MEDFileJointCorrespondence;
263 using ParaMEDMEM::MEDFileJointOneStep;
264 using ParaMEDMEM::MEDFileJoints;
265 using ParaMEDMEM::MEDFileJoints;
267 if ( _collection->getCZ().size() > 0 )
269 MEDCouplingAutoRefCountObjectPtr< MEDFileJoints > joints = MEDFileJoints::New();
271 for ( size_t i = 0; i < _collection->getCZ().size(); ++i )
273 ConnectZone* cz = _collection->getCZ()[i];
275 cz->getLocalDomainNumber() != idomain )
278 std::ostringstream oss;
279 oss << "joint_" << cz->getDistantDomainNumber();
280 cz->setName( oss.str() );
283 std::ostringstream oss;
284 oss << "connect_zone_" << i;
285 cz->setDescription( oss.str() );
288 MEDCouplingAutoRefCountObjectPtr< MEDFileJoint>
289 joint = MEDFileJoint::New( cz->getName(), finalMeshName,
290 finalMeshName, cz->getDistantDomainNumber() );
291 joint->setDescription( cz->getDescription() );
292 joints->pushJoint( joint );
294 MEDCouplingAutoRefCountObjectPtr< MEDFileJointOneStep> j1st = MEDFileJointOneStep::New();
295 joint->pushStep( j1st );
297 const MEDCouplingSkyLineArray * nodeCorr = cz->getNodeCorresp();
300 MEDCouplingAutoRefCountObjectPtr< MEDFileJointCorrespondence >
301 corr = MEDFileJointCorrespondence::New( nodeCorr->getValueArray() );
302 j1st->pushCorrespondence( corr );
305 std::vector< std::pair< int,int > > types = cz->getEntities();
306 INTERP_KERNEL::NormalizedCellType t1, t2;
307 for ( size_t it = 0; it < types.size(); ++it )
309 const MEDCouplingSkyLineArray * cellCorr =
310 cz->getEntityCorresp( types[it].first, types[it].second );
311 if ( cellCorr && cellCorr->getNumberOf() > 0 )
313 t1 = INTERP_KERNEL::NormalizedCellType( types[it].first );
314 t2 = INTERP_KERNEL::NormalizedCellType( types[it].second );
315 MEDCouplingAutoRefCountObjectPtr< MEDFileJointCorrespondence>
316 corr = MEDFileJointCorrespondence::New( cellCorr->getValueArray(), t1, t2 );
317 j1st->pushCorrespondence( corr );
321 mfm->setJoints( joints );
327 ParaMEDMEM::MEDCouplingFieldDouble* MeshCollectionDriver::getField(std::string key, std::string description, ParaMEDMEM::DataArrayDouble* data, ParaMEDMEM::MEDFileMesh* mfm, int idomain) const
329 std::string desc=description;
330 if (MyGlobals::_Verbose>20)
331 std::cout << "proc " << MyGlobals::_Rank << " : write field " << desc << std::endl;
332 std::string meshName, fieldName;
333 int typeField, DT, IT, entity;
334 FieldShortDescriptionToData(desc, fieldName, typeField, entity, DT, IT);
335 double time=StrToDouble(ExtractFromDescription(desc, "time="));
336 int typeData=StrToInt(ExtractFromDescription(desc, "typeData="));
337 std::string entityName=ExtractFromDescription(desc, "entityName=");
338 ParaMEDMEM::MEDCouplingFieldDouble* field=0;
341 std::cout << "WARNING : writeMedFile : typeData " << typeData << " not implemented for fields\n";
343 if (entityName=="MED_CELL")
345 //there is a field of idomain to write
346 field=ParaMEDMEM::MEDCouplingFieldDouble::New(ParaMEDMEM::ON_CELLS,ParaMEDMEM::ONE_TIME);
348 if (entityName=="MED_NODE_ELEMENT")
350 //there is a field of idomain to write
351 field=ParaMEDMEM::MEDCouplingFieldDouble::New(ParaMEDMEM::ON_GAUSS_NE,ParaMEDMEM::ONE_TIME);
355 std::cout << "WARNING : writeMedFile : entityName " << entityName << " not implemented for fields\n";
357 if (field && typeData==6)
359 field->setName(fieldName);
360 field->setMesh(mfm->getGenMeshAtLevel(0));
361 ParaMEDMEM::DataArrayDouble *da=data;
362 //get information for components etc..
363 std::vector<std::string> r1;
364 r1=SelectTagsInVectorOfString(MyGlobals::_General_Informations,"fieldName="+fieldName);
365 r1=SelectTagsInVectorOfString(r1,"typeField="+IntToStr(typeField));
366 r1=SelectTagsInVectorOfString(r1,"DT="+IntToStr(DT));
367 r1=SelectTagsInVectorOfString(r1,"IT="+IntToStr(IT));
368 //not saved in file? field->setDescription(ExtractFromDescription(r1[0], "fieldDescription="));
369 int nbc=StrToInt(ExtractFromDescription(r1[0], "nbComponents="));
370 if (nbc==da->getNumberOfComponents())
372 for (int i=0; i<nbc; i++)
373 da->setInfoOnComponent(i,ExtractFromDescription(r1[0], "componentInfo"+IntToStr(i)+"="));
377 std::cerr << "Problem On field " << fieldName << " : number of components unexpected " << da->getNumberOfComponents() << std::endl;
380 field->setTime(time,DT,IT);
381 field->checkCoherency();
386 void MeshCollectionDriver::writeMedFile(int idomain, const std::string& distfilename) const
388 ParaMEDMEM::MEDFileMesh* mfm = getMesh( idomain );
389 mfm->write(distfilename,2);
391 std::string key="/inewFieldDouble="+IntToStr(idomain)+"/";
392 std::map<std::string,ParaMEDMEM::DataArrayDouble*>::iterator it;
394 for (it=_collection->getMapDataArrayDouble().begin() ; it!=_collection->getMapDataArrayDouble().end(); it++)
396 size_t found=(*it).first.find(key);
397 if (found==std::string::npos)
399 ParaMEDMEM::MEDCouplingFieldDouble* field=0;
400 field = getField(key, (*it).first, (*it).second, mfm, idomain);
404 MEDLoader::WriteField(distfilename,field,false);
406 catch(INTERP_KERNEL::Exception& e)
408 //cout trying rewrite all data, only one field defined
409 std::string tmp,newName=distfilename;
410 std::string fieldName;
411 fieldName=field->getName();
412 tmp+="_"+fieldName+"_"+IntToStr(nbfFieldFound)+".med";
413 newName.replace(newName.find(".med"),4,tmp);
414 std::cout << "WARNING : writeMedFile : create a new file name with only one field because MEDLoader::WriteField throw:" << newName << std::endl;
415 MEDLoader::WriteField(newName,field,true);
421 ParaMEDMEM::MEDFileData* MeshCollectionDriver::getMEDFileData()
423 ParaMEDMEM::MEDFileData* newdata = ParaMEDMEM::MEDFileData::New();
425 ParaMEDMEM::MEDCouplingAutoRefCountObjectPtr<ParaMEDMEM::MEDFileMeshes> meshes;
426 ParaMEDMEM::MEDCouplingAutoRefCountObjectPtr<ParaMEDMEM::MEDFileFields> fields;
427 meshes = ParaMEDMEM::MEDFileMeshes::New();
428 fields = ParaMEDMEM::MEDFileFields::New();
430 for (size_t i=0; i<_collection->getMesh().size(); i++)
432 ParaMEDMEM::MEDFileMesh* mfm = getMesh( i );
433 meshes->pushMesh(mfm);
435 std::string key="/inewFieldDouble="+IntToStr(i)+"/";
436 std::map<std::string,ParaMEDMEM::DataArrayDouble*>::iterator it;
437 ParaMEDMEM::MEDFileFieldMultiTS* fieldsMTS = ParaMEDMEM::MEDFileFieldMultiTS::New();
438 for (it=_collection->getMapDataArrayDouble().begin() ; it!=_collection->getMapDataArrayDouble().end(); it++)
440 size_t found=(*it).first.find(key);
441 if (found==std::string::npos)
443 ParaMEDMEM::MEDCouplingFieldDouble* field=0;
444 field=getField(key, (*it).first, (*it).second, mfm, i);
445 ParaMEDMEM::MEDFileField1TS* f1ts = ParaMEDMEM::MEDFileField1TS::New();
446 f1ts->setFieldNoProfileSBT(field);
447 fieldsMTS->pushBackTimeStep(f1ts);
452 fields->pushField(fieldsMTS);
454 fieldsMTS->decrRef();
457 newdata->setMeshes(meshes);
458 newdata->setFields(fields);