1 // Copyright (C) 2007-2023 CEA, EDF
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 "MEDCouplingSkyLineArray.hxx"
31 #include "MEDCouplingUMesh.hxx"
32 #include "MEDFileData.hxx"
33 #include "MEDFileField.hxx"
34 #include "MEDFileJoint.hxx"
35 #include "MEDFileMesh.hxx"
36 #include "MEDLoader.hxx"
45 #include <libxml/tree.h>
46 #include <libxml/parser.h>
47 #include <libxml/xpath.h>
48 #include <libxml/xpathInternals.h>
52 using namespace MEDPARTITIONER;
54 MeshCollectionDriver::MeshCollectionDriver(MeshCollection* collection):_collection(collection)
58 /*!reads a unique MED File v>=2.1
59 * and mounts the corresponding mesh in memory
60 *\param filename binary file
61 *\param meshname mesh name in the MED file
63 int MeshCollectionDriver::readSeq(const char* filename, const char* meshname)
65 std::cout << "readSeq" << std::endl;
66 MyGlobals::_File_Names.resize(1);
67 MyGlobals::_File_Names[0]=std::string(filename);
69 MEDCoupling::MEDFileUMesh* mfm=MEDCoupling::MEDFileUMesh::New(filename,meshname);
70 //puts the only mesh in the mesh vector
71 (_collection->getMesh()).push_back(mfm->getLevel0Mesh(false));
72 (_collection->getFaceMesh()).push_back(mfm->getLevelM1Mesh(false));
75 MEDCoupling::DataArrayIdType* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCopy());
76 MEDCoupling::DataArrayIdType* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCopy());
77 (_collection->getCellFamilyIds()).push_back(cellIds);
78 (_collection->getFaceFamilyIds()).push_back(faceIds);
81 (_collection->getFamilyInfo())=mfm->getFamilyInfo();
82 (_collection->getGroupInfo())=mfm->getGroupInfo();
84 (_collection->getCZ()).clear();
86 ParallelTopology* aPT = new ParallelTopology((_collection->getMesh()));
87 _collection->setTopology(aPT, true);
88 _collection->setName(meshname);
89 _collection->setDomainNames(meshname);
94 void MeshCollectionDriver::readMEDFileData(const MEDCoupling::MEDFileData* filedata)
96 const int nbDomains = filedata->getMeshes()->getNumberOfMeshes();
97 _collection->getMesh() .resize( nbDomains, 0 );
98 _collection->getFaceMesh() .resize( nbDomains, 0 );
99 _collection->getCellFamilyIds().resize( nbDomains, 0 );
100 _collection->getFaceFamilyIds().resize( nbDomains, 0 );
102 for (int i=0; i<nbDomains; i++)
104 MEDCoupling::MEDFileUMesh *mfm = dynamic_cast<MEDCoupling::MEDFileUMesh *>(filedata->getMeshes()->getMeshAtPos(i));
106 if ( mfm && mfm->getMeshDimension() > 0 )
107 _collection->setNonEmptyMesh( i );
110 ParallelTopology* aPT = new ParallelTopology(_collection->getMesh());
111 _collection->setTopology(aPT, true);
114 _collection->setName( filedata->getMeshes()->getMeshAtPos(0)->getName() );
115 _collection->setDomainNames( _collection->getName() );
117 if ( ParaDomainSelector* domainSelector = _collection->getParaDomainSelector() )
118 if ( _collection->isParallelMode() )
120 //to know nb of cells on each proc to compute global cell ids from locally global
121 domainSelector->gatherNbOf(_collection->getMesh());
125 void MeshCollectionDriver::readFileData(std::string file,std::string meshname,int idomain) const
127 MEDCoupling::MEDFileUMesh* mfm=MEDCoupling::MEDFileUMesh::New(file,meshname);
128 readData(mfm,idomain);
132 void MeshCollectionDriver::readData(MEDCoupling::MEDFileUMesh* mfm, int idomain) const
134 std::vector<int> nonEmpty=mfm->getNonEmptyLevels();
137 (_collection->getMesh())[idomain]=mfm->getLevel0Mesh(false);
138 //reading families groups
139 MEDCoupling::DataArrayIdType* cellIds(mfm->getFamilyFieldAtLevel(0)->deepCopy());
140 (_collection->getCellFamilyIds())[idomain]=cellIds;
144 (_collection->getMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want tests;
145 MEDCoupling::DataArrayIdType* empty=MEDCoupling::DataArrayIdType::New();
147 (_collection->getCellFamilyIds())[idomain]=empty;
148 std::cout<<"\nNO Level0Mesh (Cells)\n";
152 if (nonEmpty.size()>1 && nonEmpty[1]==-1)
154 (_collection->getFaceMesh())[idomain]=mfm->getLevelM1Mesh(false);
155 //reading families groups
156 MEDCoupling::DataArrayIdType* faceIds(mfm->getFamilyFieldAtLevel(-1)->deepCopy());
157 (_collection->getFaceFamilyIds())[idomain]=faceIds;
158 if (MyGlobals::_Verbose>10)
159 std::cout << "proc " << MyGlobals::_Rank << " : WITH Faces\n";
163 throw INTERP_KERNEL::Exception("no faces");
168 (_collection->getFaceMesh())[idomain]=CreateEmptyMEDCouplingUMesh(); // or 0 if you want test;
169 MEDCoupling::DataArrayIdType* empty=MEDCoupling::DataArrayIdType::New();
170 (_collection->getFaceFamilyIds())[idomain]=empty;
171 if (MyGlobals::_Verbose>10)
172 std::cout << "proc " << MyGlobals::_Rank << " : WITHOUT Faces\n";
175 _collection->getFamilyInfo()=mfm->getFamilyInfo();
176 _collection->getGroupInfo()=mfm->getGroupInfo();
179 void MeshCollectionDriver::readSubdomain(int idomain)
181 std::string meshname=MyGlobals::_Mesh_Names[idomain];
182 std::string file=MyGlobals::_File_Names[idomain];
183 readFileData(file,meshname,idomain);
185 std::vector<std::string> localInformation;
187 localInformation.push_back(str+"ioldDomain="+IntToStr(idomain));
188 localInformation.push_back(str+"meshName="+meshname);
189 MyGlobals::_General_Informations.push_back(SerializeFromVectorOfString(localInformation));
190 std::vector<std::string> localFields=BrowseAllFieldsOnMesh(file, meshname, idomain);
191 if (localFields.size()>0)
192 MyGlobals::_Field_Descriptions.push_back(SerializeFromVectorOfString(localFields));
195 MEDCoupling::MEDFileMesh* MeshCollectionDriver::getMesh(int idomain) const
197 MEDCoupling::MEDFileUMesh* mfm = MEDCoupling::MEDFileUMesh::New();
199 MEDCoupling::MEDCouplingUMesh* cellMesh=_collection->getMesh(idomain);
200 MEDCoupling::MEDCouplingUMesh* faceMesh=_collection->getFaceMesh(idomain);
201 // std::string cleFilter=Cle1ToStr("filterFaceOnCell",idomain);
202 // MEDCoupling::DataArrayInt* filter=0;
203 // if (_collection->getMapDataArrayInt().find(cleFilter)!=_collection->getMapDataArrayInt().end())
205 // filter=_collection->getMapDataArrayInt().find(cleFilter)->second;
206 // int* index=filter->getPointer();
207 // faceMeshFilter=(MEDCoupling::MEDCouplingUMesh *) faceMesh->buildPartOfMySelf(index,index+filter->getNbOfElems(),true);
208 // faceMesh=faceMeshFilter;
210 // if (faceMeshFilter!=0)
211 // faceMeshFilter->decrRef();
212 std::string finalMeshName="";
213 if (MyGlobals::_General_Informations.size()!=0)
215 std::size_t found=MyGlobals::_General_Informations[0].find("finalMeshName=");
216 if ((found!=std::string::npos) && (found>0))
218 finalMeshName=ExtractFromDescription(MyGlobals::_General_Informations[0], "finalMeshName=");
221 if (finalMeshName.empty())
223 finalMeshName=_collection->getName();
225 cellMesh->setName(finalMeshName);
226 mfm->setMeshAtLevel( 0, cellMesh );
228 faceMesh->checkConsistencyLight();
229 if (faceMesh->getNumberOfCells()>0)
231 faceMesh->tryToShareSameCoordsPermute(*cellMesh, 1e-10);
232 faceMesh->setName(finalMeshName);
233 mfm->setMeshAtLevel( -1, faceMesh );
236 // MEDCoupling::MEDCouplingUMesh* boundaryMesh=0;
237 // if (MyGlobals::_Create_Boundary_Faces>0)
239 // //try to write Boundary meshes
240 // bool keepCoords=false; //TODO or true
241 // boundaryMesh=(MEDCoupling::MEDCouplingUMesh *) cellMesh->buildBoundaryMesh(keepCoords);
242 // boundaryMesh->setName("boundaryMesh");
243 // if (boundaryMesh!=0)
245 // //doing that testMesh becomes second mesh sorted by alphabetical order of name
246 // WriteUMesh(distfilename, boundaryMesh, false);
247 // boundaryMesh->decrRef();
250 mfm->setFamilyInfo(_collection->getFamilyInfo());
251 mfm->setGroupInfo(_collection->getGroupInfo());
252 std::string key=Cle1ToStr("faceFamily_toArray",idomain);
253 if ( faceMesh->getNumberOfCells()>0 && _collection->getMapDataArrayInt().find(key)!=_collection->getMapDataArrayInt().end())
254 mfm->setFamilyFieldArr(-1,_collection->getMapDataArrayInt().find(key)->second);
255 key=Cle1ToStr("cellFamily_toArray",idomain);
256 if (_collection->getMapDataArrayInt().find(key)!=_collection->getMapDataArrayInt().end())
257 mfm->setFamilyFieldArr(0,_collection->getMapDataArrayInt().find(key)->second);
261 using MEDCoupling::MCAuto;
262 using MEDCoupling::MEDCouplingSkyLineArray;
263 using MEDCoupling::MEDFileJoint;
264 using MEDCoupling::MEDFileJointCorrespondence;
265 using MEDCoupling::MEDFileJointOneStep;
266 using MEDCoupling::MEDFileJoints;
267 using MEDCoupling::MEDFileJoints;
269 if ( _collection->getCZ().size() > 0 )
271 MCAuto< MEDFileJoints > joints = MEDFileJoints::New();
273 for ( size_t i = 0; i < _collection->getCZ().size(); ++i )
275 ConnectZone* cz = _collection->getCZ()[i];
277 cz->getLocalDomainNumber() != idomain )
280 std::ostringstream oss;
281 oss << "joint_" << cz->getDistantDomainNumber();
282 cz->setName( oss.str() );
285 std::ostringstream oss;
286 oss << "connect_zone_" << i;
287 cz->setDescription( oss.str() );
290 MCAuto< MEDFileJoint>
291 joint = MEDFileJoint::New( cz->getName(), finalMeshName,
292 finalMeshName, cz->getDistantDomainNumber() );
293 joint->setDescription( cz->getDescription() );
294 joints->pushJoint( joint );
296 MCAuto< MEDFileJointOneStep> j1st = MEDFileJointOneStep::New();
297 joint->pushStep( j1st );
299 const MEDCouplingSkyLineArray * nodeCorr = cz->getNodeCorresp();
302 MCAuto< MEDFileJointCorrespondence >
303 corr = MEDFileJointCorrespondence::New( nodeCorr->getValuesArray() );
304 j1st->pushCorrespondence( corr );
307 std::vector< std::pair< mcIdType,mcIdType > > types = cz->getEntities();
308 INTERP_KERNEL::NormalizedCellType t1, t2;
309 for ( size_t it = 0; it < types.size(); ++it )
311 const MEDCouplingSkyLineArray * cellCorr =
312 cz->getEntityCorresp( types[it].first, types[it].second );
313 if ( cellCorr && cellCorr->getNumberOf() > 0 )
315 t1 = INTERP_KERNEL::NormalizedCellType( types[it].first );
316 t2 = INTERP_KERNEL::NormalizedCellType( types[it].second );
317 MCAuto< MEDFileJointCorrespondence>
318 corr = MEDFileJointCorrespondence::New( cellCorr->getValuesArray(), t1, t2 );
319 j1st->pushCorrespondence( corr );
323 mfm->setJoints( joints );
329 MEDCoupling::MEDCouplingFieldDouble* MeshCollectionDriver::getField(std::string key, std::string description, MEDCoupling::DataArrayDouble* data, MEDCoupling::MEDFileMesh* mfm, int idomain) const
331 std::string desc=description;
332 if (MyGlobals::_Verbose>20)
333 std::cout << "proc " << MyGlobals::_Rank << " : write field " << desc << std::endl;
334 std::string meshName, fieldName;
335 int typeField, DT, IT, entity;
336 FieldShortDescriptionToData(desc, fieldName, typeField, entity, DT, IT);
337 double time=StrToDouble(ExtractFromDescription(desc, "time="));
338 int typeData=StrToInt(ExtractFromDescription(desc, "typeData="));
339 std::string entityName=ExtractFromDescription(desc, "entityName=");
340 MEDCoupling::MEDCouplingFieldDouble* field=0;
343 std::cout << "WARNING : writeMedFile : typeData " << typeData << " not implemented for fields\n";
345 if (entityName=="MED_CELL")
347 //there is a field of idomain to write
348 field=MEDCoupling::MEDCouplingFieldDouble::New(MEDCoupling::ON_CELLS,MEDCoupling::ONE_TIME);
350 if (entityName=="MED_NODE_ELEMENT")
352 //there is a field of idomain to write
353 field=MEDCoupling::MEDCouplingFieldDouble::New(MEDCoupling::ON_GAUSS_NE,MEDCoupling::ONE_TIME);
357 std::cout << "WARNING : writeMedFile : entityName " << entityName << " not implemented for fields\n";
359 if (field && typeData==6)
361 field->setName(fieldName);
362 field->setMesh(mfm->getMeshAtLevel(0));
363 MEDCoupling::DataArrayDouble *da=data;
364 //get information for components etc..
365 std::vector<std::string> r1;
366 r1=SelectTagsInVectorOfString(MyGlobals::_General_Informations,"fieldName="+fieldName);
367 r1=SelectTagsInVectorOfString(r1,"typeField="+IntToStr(typeField));
368 r1=SelectTagsInVectorOfString(r1,"DT="+IntToStr(DT));
369 r1=SelectTagsInVectorOfString(r1,"IT="+IntToStr(IT));
370 //not saved in file? field->setDescription(ExtractFromDescription(r1[0], "fieldDescription="));
371 std::size_t nbc=StrToInt(ExtractFromDescription(r1[0], "nbComponents="));
372 if (nbc==da->getNumberOfComponents())
374 for (unsigned int i=0; i<nbc; i++)
375 da->setInfoOnComponent(i,ExtractFromDescription(r1[0], "componentInfo"+IntToStr(i)+"="));
379 std::cerr << "Problem On field " << fieldName << " : number of components unexpected " << da->getNumberOfComponents() << std::endl;
382 field->setTime(time,DT,IT);
383 field->checkConsistencyLight();
388 void MeshCollectionDriver::writeMedFile(int idomain, const std::string& distfilename) const
390 MEDCoupling::MEDFileMesh* mfm = getMesh( idomain );
391 mfm->write(distfilename,2);
393 std::string key="/inewFieldDouble="+IntToStr(idomain)+"/";
394 std::map<std::string,MEDCoupling::DataArrayDouble*>::iterator it;
396 for (it=_collection->getMapDataArrayDouble().begin() ; it!=_collection->getMapDataArrayDouble().end(); it++)
398 size_t found=(*it).first.find(key);
399 if (found==std::string::npos)
401 MEDCoupling::MEDCouplingFieldDouble* field=0;
402 field = getField(key, (*it).first, (*it).second, mfm, idomain);
406 WriteField(distfilename,field,false);
408 catch(INTERP_KERNEL::Exception& e)
410 //cout trying rewrite all data, only one field defined
411 std::string tmp,newName=distfilename;
412 std::string fieldName;
413 fieldName=field->getName();
414 tmp+="_"+fieldName+"_"+IntToStr(nbfFieldFound)+".med";
415 newName.replace(newName.find(".med"),4,tmp);
416 std::cout << "WARNING : writeMedFile : create a new file name with only one field because WriteField throw:" << newName << std::endl;
417 WriteField(newName,field,true);
423 MEDCoupling::MEDFileData* MeshCollectionDriver::getMEDFileData()
425 MEDCoupling::MEDFileData* newdata = MEDCoupling::MEDFileData::New();
427 MEDCoupling::MCAuto<MEDCoupling::MEDFileMeshes> meshes;
428 MEDCoupling::MCAuto<MEDCoupling::MEDFileFields> fields;
429 meshes = MEDCoupling::MEDFileMeshes::New();
430 fields = MEDCoupling::MEDFileFields::New();
432 for (unsigned int i=0; i<_collection->getMesh().size(); i++)
434 MEDCoupling::MEDFileMesh* mfm = getMesh( i );
435 meshes->pushMesh(mfm);
437 std::string key="/inewFieldDouble="+IntToStr(i)+"/";
438 std::map<std::string,MEDCoupling::DataArrayDouble*>::iterator it;
439 MEDCoupling::MEDFileFieldMultiTS* fieldsMTS = MEDCoupling::MEDFileFieldMultiTS::New();
440 for (it=_collection->getMapDataArrayDouble().begin() ; it!=_collection->getMapDataArrayDouble().end(); it++)
442 size_t found=(*it).first.find(key);
443 if (found==std::string::npos)
445 MEDCoupling::MEDCouplingFieldDouble* field=0;
446 field=getField(key, (*it).first, (*it).second, mfm, i);
447 MEDCoupling::MEDFileField1TS* f1ts = MEDCoupling::MEDFileField1TS::New();
448 f1ts->setFieldNoProfileSBT(field);
449 fieldsMTS->pushBackTimeStep(f1ts);
454 fields->pushField(fieldsMTS);
456 fieldsMTS->decrRef();
459 newdata->setMeshes(meshes);
460 newdata->setFields(fields);