1 // Project MULTIPR, IOLS WP1.2.1 - EDF/CS
2 // Partitioning/decimation module for the SALOME v3.2 platform
5 * \file MULTIPR_Mesh.cxx
7 * \brief see MULTIPR_Mesh.hxx
9 * \author Olivier LE ROUX - CS, Virtual Reality Dpt
14 //*****************************************************************************
16 //*****************************************************************************
18 #include "MULTIPR_Mesh.hxx"
19 #include "MULTIPR_Nodes.hxx"
20 #include "MULTIPR_Elements.hxx"
21 #include "MULTIPR_Family.hxx"
22 #include "MULTIPR_Profil.hxx"
23 #include "MULTIPR_GaussLoc.hxx"
24 #include "MULTIPR_Field.hxx"
25 #include "MULTIPR_MeshDis.hxx"
26 #include "MULTIPR_PointOfField.hxx"
27 #include "MULTIPR_DecimationFilter.hxx"
28 #include "MULTIPR_Utils.hxx"
29 #include "MULTIPR_Exceptions.hxx"
30 #include "MULTIPR_Globals.hxx"
31 #include "MULTIPR_API.hxx"
40 //*****************************************************************************
41 // Class Mesh implementation
42 //*****************************************************************************
61 mMEDfilename[0] = '\0';
68 mMeshType = MED_NON_STRUCTURE;
70 for (int itDim = 0 ; itDim < 3 ; itDim++)
72 mMeshBBoxMin[itDim] = numeric_limits<med_float>::quiet_NaN();
73 mMeshBBoxMax[itDim] = numeric_limits<med_float>::quiet_NaN();
76 if (mNodes != NULL) { delete mNodes; mNodes = NULL; }
77 if (mElements != NULL) { delete mElements; mElements = NULL; }
79 for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
81 delete mFamilies[itFam];
86 for (unsigned itGroup = 0 ; itGroup < mGroups.size() ; itGroup++)
88 delete mGroups[itGroup];
91 mGroupNameToGroup.clear();
93 for (unsigned itGaussLoc = 0 ; itGaussLoc < mGaussLoc.size() ; itGaussLoc++)
95 delete mGaussLoc[itGaussLoc];
98 mGaussLocNameToGaussLoc.clear();
100 for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
102 delete mFields[itField];
106 for (unsigned itProfil = 0 ; itProfil < mProfils.size() ; itProfil++)
108 delete mProfils[itProfil];
112 mFlagPrintAll = false;
116 vector<string> Mesh::getNameScalarFields() const
120 for (int itField = 0 ; itField < mFields.size() ; itField++)
122 Field* currentField = mFields[itField];
124 // only get scalar fields, not vectorial fields
125 // (because, currently, decimation can only be performed on scalar fields)
126 if (currentField->getNumberOfComponents() == 1)
128 res.push_back(currentField->getName());
136 int Mesh::getTimeStamps(const char* pFieldName) const
138 for (int itField = 0 ; itField < mFields.size() ; itField++)
140 Field* currentField = mFields[itField];
141 if (strcmp(currentField->getName(), pFieldName) == 0)
143 return currentField->getNumberOfTimeSteps();
151 Field* Mesh::getFieldByName(const char* pFieldName) const
153 if (pFieldName == NULL) throw NullArgumentException("pFieldName should not be NULL", __FILE__, __LINE__);
155 Field* retField = NULL;
158 for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
160 Field* currentField = mFields[itField];
161 if (strcmp(pFieldName, currentField->getName()) == 0)
164 retField = currentField;
173 GaussLoc* Mesh::getGaussLocByName(const char* pGaussLocName) const
175 if (pGaussLocName == NULL) throw NullArgumentException("pGaussLocName should not be NULL", __FILE__, __LINE__);
177 map<string, GaussLoc*>::const_iterator itGaussLoc = mGaussLocNameToGaussLoc.find(pGaussLocName);
178 GaussLoc* retGaussLoc = NULL;
180 if (itGaussLoc != mGaussLocNameToGaussLoc.end())
182 retGaussLoc = (*itGaussLoc).second;
189 int Mesh::getNumberOfElements() const
191 if (mElements == NULL) throw IllegalStateException("", __FILE__, __LINE__);
193 return mElements->getNumberOfElements();
197 Mesh* Mesh::createFromSetOfElements(const std::set<med_int>& pSetOfElements, const char* pNewMeshName)
199 if (pNewMeshName == NULL) throw NullArgumentException("pNewMeshName", __FILE__, __LINE__);
201 //---------------------------------------------------------------------
203 //---------------------------------------------------------------------
204 Mesh* mesh = new Mesh();
206 //---------------------------------------------------------------------
207 // Build name of the new mesh
208 //---------------------------------------------------------------------
209 strcpy(mesh->mMeshName, pNewMeshName);
211 MULTIPR_LOG("Mesh name=|" << mesh->mMeshName << "|" << endl);
213 //---------------------------------------------------------------------
214 // Fill general infos
215 //---------------------------------------------------------------------
216 strcpy(mesh->mMeshUName, mMeshUName);
217 strcpy(mesh->mMeshDesc, mMeshDesc);
219 mesh->mMeshDim = mMeshDim;
220 mesh->mMeshType = mMeshType;
222 MULTIPR_LOG("Mesh u. name=|" << mesh->mMeshUName << "|" << endl);
223 MULTIPR_LOG("Mesh desc=|" << mesh->mMeshDesc << "|" << endl);
224 MULTIPR_LOG("Mesh dim=" << mesh->mMeshDim << endl);
225 MULTIPR_LOG("Mesh Type=" << mesh->mMeshType << endl);
227 //---------------------------------------------------------------------
228 // Build nodes and elements
229 //---------------------------------------------------------------------
230 // get all elements involved
231 mesh->mElements = mElements->extractSubSet(pSetOfElements);
232 MULTIPR_LOG((*(mesh->mElements)) << endl);
234 // get all nodes involved
235 const set<med_int> setOfNodes = mesh->mElements->getSetOfNodes();
236 mesh->mNodes = mNodes->extractSubSet(setOfNodes);
237 MULTIPR_LOG((*(mesh->mNodes)) << endl);
239 //---------------------------------------------------------------------
241 //---------------------------------------------------------------------
242 mesh->mElements->remap();
243 MULTIPR_LOG((*(mesh->mElements)) << endl);
245 //---------------------------------------------------------------------
247 //---------------------------------------------------------------------
248 MULTIPR_LOG("Build fam.:" << endl);
249 // get families of nodes
251 set<med_int> famOfNodes = mesh->mNodes->getSetOfFamilies();
252 for (set<med_int>::iterator itFam = famOfNodes.begin() ; itFam != famOfNodes.end() ; itFam++)
254 Family* famSrc = mFamIdToFam[*itFam];
255 cout << (*famSrc) << endl;
256 Family* famDest = famSrc->extractGroup(NULL);
257 mesh->mFamilies.push_back(famDest);
261 // get families of elements
263 set<med_int> famOfElt = mesh->mElements->getSetOfFamilies();
264 for (set<med_int>::iterator itFam = famOfElt.begin() ; itFam != famOfElt.end() ; itFam++)
266 Family* famSrc = mFamIdToFam[*itFam];
267 Family* famDest = famSrc->extractGroup(NULL);
268 mesh->mFamilies.push_back(famDest);
272 MULTIPR_LOG("Finalize:");
274 // fill families with elements and build groups
275 mesh->finalizeFamiliesAndGroups();
279 //---------------------------------------------------------------------
280 // Create new fields and collect Gauss
281 //---------------------------------------------------------------------
283 set<string> newSetOfGauss;
284 for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
286 Field* currentField = mFields[itField];
289 if (currentField->isFieldOnNodes())
291 newField = currentField->extractSubSet(setOfNodes);
295 newField = currentField->extractSubSet(pSetOfElements);
298 if (!newField->isEmpty())
300 mesh->mFields.push_back(newField);
301 newField->getSetOfGaussLoc(newSetOfGauss);
304 MULTIPR_LOG("Collect fields: ok: #gauss=" << newSetOfGauss.size() << endl);
306 //---------------------------------------------------------------------
308 //---------------------------------------------------------------------
309 for (set<string>::iterator itSet = newSetOfGauss.begin() ; itSet != newSetOfGauss.end(); itSet++)
311 const string& keyName = (*itSet);
313 GaussLoc* gaussLoc = getGaussLocByName(keyName.c_str());
314 if (gaussLoc != NULL)
316 GaussLoc* copyGaussLoc = new GaussLoc(*gaussLoc);
317 mesh->mGaussLoc.push_back(copyGaussLoc);
318 mesh->mGaussLocNameToGaussLoc.insert(make_pair(copyGaussLoc->getName(), copyGaussLoc));
322 //---------------------------------------------------------------------
324 //---------------------------------------------------------------------
325 mesh->mNodes->getBBox(mesh->mMeshBBoxMin, mesh->mMeshBBoxMax);
331 Mesh* Mesh::createFromGroup(const Group* pGroup, const char* pNewMeshName)
333 if (pGroup == NULL) throw NullArgumentException("pGroup should not be NULL", __FILE__, __LINE__);
334 if (pNewMeshName == NULL) throw NullArgumentException("pNewMeshName should not be NULL", __FILE__, __LINE__);
335 if (strlen(pNewMeshName) > MED_TAILLE_NOM) throw IllegalArgumentException("pNewMeshName length too long", __FILE__, __LINE__);
337 //---------------------------------------------------------------------
339 //---------------------------------------------------------------------
340 Mesh* mesh = new Mesh();
342 //---------------------------------------------------------------------
343 // Build name of the new mesh
344 //---------------------------------------------------------------------
345 strcpy(mesh->mMeshName, pNewMeshName);
347 MULTIPR_LOG("Mesh name=|" << mesh->mMeshName << "|" << endl);
349 //---------------------------------------------------------------------
350 // Fill general infos
351 //---------------------------------------------------------------------
352 strcpy(mesh->mMeshUName, mMeshUName);
353 strcpy(mesh->mMeshDesc, mMeshDesc);
355 mesh->mMeshDim = mMeshDim;
356 mesh->mMeshType = mMeshType;
358 MULTIPR_LOG("Mesh u. name=|" << mesh->mMeshUName << "|" << endl);
359 MULTIPR_LOG("Mesh desc=|" << mesh->mMeshDesc << "|" << endl);
360 MULTIPR_LOG("Mesh dim=" << mesh->mMeshDim << endl);
361 MULTIPR_LOG("Mesh Type=" << mesh->mMeshType << endl);
363 //---------------------------------------------------------------------
364 // Build nodes and elements
365 //---------------------------------------------------------------------
366 // get all elements involved
367 const set<med_int> setOfElt = pGroup->getSetOfElt();
368 mesh->mElements = mElements->extractSubSet(setOfElt);
369 MULTIPR_LOG((*(mesh->mElements)) << endl);
371 // get all nodes involved
372 const set<med_int> setOfNodes = mesh->mElements->getSetOfNodes();
373 mesh->mNodes = mNodes->extractSubSet(setOfNodes);
374 MULTIPR_LOG((*(mesh->mNodes)) << endl);
376 //---------------------------------------------------------------------
378 //---------------------------------------------------------------------
379 mesh->mElements->remap();
380 MULTIPR_LOG((*(mesh->mElements)) << endl);
382 //---------------------------------------------------------------------
384 //---------------------------------------------------------------------
385 MULTIPR_LOG("Build fam.:" << endl);
386 // get families of nodes
388 set<med_int> famOfNodes = mesh->mNodes->getSetOfFamilies();
389 for (set<med_int>::iterator itFam = famOfNodes.begin() ; itFam != famOfNodes.end() ; itFam++)
391 Family* famSrc = mFamIdToFam[*itFam];
392 Family* famDest = famSrc->extractGroup(pGroup->getName().c_str());
393 mesh->mFamilies.push_back(famDest);
397 // get families of elements
399 set<med_int> famOfElt = mesh->mElements->getSetOfFamilies();
400 for (set<med_int>::iterator itFam = famOfElt.begin() ; itFam != famOfElt.end() ; itFam++)
402 Family* famSrc = mFamIdToFam[*itFam];
403 Family* famDest = famSrc->extractGroup(pGroup->getName().c_str());
404 mesh->mFamilies.push_back(famDest);
408 MULTIPR_LOG("Finalize:");
410 // fill families with elements and build groups
411 mesh->finalizeFamiliesAndGroups();
415 //---------------------------------------------------------------------
416 // Create new fields and collect Gauss
417 //---------------------------------------------------------------------
419 set<string> newSetOfGauss;
420 for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
422 Field* currentField = mFields[itField];
425 if (currentField->isFieldOnNodes())
427 newField = currentField->extractSubSet(setOfNodes);
431 newField = currentField->extractSubSet(setOfElt);
434 if (!newField->isEmpty())
436 mesh->mFields.push_back(newField);
437 newField->getSetOfGaussLoc(newSetOfGauss);
440 MULTIPR_LOG("Collect fields: ok: #gauss=" << newSetOfGauss.size() << endl);
442 //---------------------------------------------------------------------
444 //---------------------------------------------------------------------
445 for (set<string>::iterator itSet = newSetOfGauss.begin() ; itSet != newSetOfGauss.end(); itSet++)
447 const string& keyName = (*itSet);
449 GaussLoc* gaussLoc = getGaussLocByName(keyName.c_str());
450 if (gaussLoc != NULL)
452 GaussLoc* copyGaussLoc = new GaussLoc(*gaussLoc);
453 mesh->mGaussLoc.push_back(copyGaussLoc);
454 mesh->mGaussLocNameToGaussLoc.insert(make_pair(copyGaussLoc->getName(), copyGaussLoc));
458 //---------------------------------------------------------------------
460 //---------------------------------------------------------------------
461 mesh->mNodes->getBBox(mesh->mMeshBBoxMin, mesh->mMeshBBoxMax);
467 MeshDis* Mesh::splitGroupsOfElements()
469 MeshDis* meshDis = new MeshDis();
471 // get prefix from the original MED filename
472 string strPrefix = removeExtension(mMEDfilename, ".med");
477 for (unsigned itGroup = 0 ; itGroup < mGroups.size() ; itGroup++)
479 Group* currentGroup = mGroups[itGroup];
481 // skip this group if it is a group of nodes
482 if (currentGroup->isGroupOfNodes())
487 char strPartName[256];
488 sprintf(strPartName, "%s_%d", mMeshName, numGroup);
490 char strMEDfilename[256];
491 sprintf(strMEDfilename, "%s_grain%d.med", strPrefix.c_str(), numGroup);
493 Mesh* mesh = createFromGroup(currentGroup, mMeshName);
495 // skip the group which contain all the others groups, even it contains only 1 group
496 if ((mesh->mElements->getNumberOfElements() == mElements->getNumberOfElements()) && (mGroups.size() > 1))
503 MeshDisPart::MULTIPR_WRITE_MESH,
518 Mesh* Mesh::decimate(
519 const char* pFilterName,
521 const char* pNameNewMesh)
523 //---------------------------------------------------------------------
525 //---------------------------------------------------------------------
526 if (pFilterName == NULL) throw NullArgumentException("pFilterName should not be NULL", __FILE__, __LINE__);
527 if (pArgv == NULL) throw NullArgumentException("pArgv should not be NULL", __FILE__, __LINE__);
528 if (pNameNewMesh == NULL) throw NullArgumentException("pNameNewMesh should not be NULL", __FILE__, __LINE__);
530 //---------------------------------------------------------------------
531 // Instanciate filter used for decimation
532 //---------------------------------------------------------------------
533 DecimationFilter* filter = DecimationFilter::create(pFilterName);
535 //---------------------------------------------------------------------
536 // Create new mesh by decimating current one
537 //---------------------------------------------------------------------
538 Mesh* decimatedMesh = filter->apply(this, pArgv, pNameNewMesh);
540 //---------------------------------------------------------------------
542 //---------------------------------------------------------------------
545 return decimatedMesh;
550 void Mesh::getAllPointsOfField(Field* pField, int pTimeStepIt, std::vector<PointOfField>& pPoints)
552 //---------------------------------------------------------------------
554 //---------------------------------------------------------------------
556 if (pField == NULL) throw NullArgumentException("field should not be NULL", __FILE__, __LINE__);
557 if (pTimeStepIt < 1) throw IllegalArgumentException("invalid field iteration; should be >= 1", __FILE__, __LINE__);
559 if (mMeshDim != 3) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
560 if (pField->getType() != MED_FLOAT64) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
561 if (pField->getNumberOfComponents() != 1) throw UnsupportedOperationException("field have more than 1 component (vectorial field, expected scalar field)", __FILE__, __LINE__);
563 //---------------------------------------------------------------------
565 //---------------------------------------------------------------------
567 if (pField->isFieldOnNodes())
569 //-------------------------------------------------------------
570 // Case 1: field of nodes
571 //-------------------------------------------------------------
572 if (mNodes == NULL) throw IllegalStateException("no nodes in the current mesh", __FILE__, __LINE__);
575 for (int itNode = 0, size = mNodes->getNumberOfNodes() ; itNode < size ; itNode++)
577 // collect coordinates and value of the point
578 const med_float* coo = mNodes->getCoordinates(itNode);
580 const med_float* val =
581 reinterpret_cast<const med_float*>(pField->getValue(pTimeStepIt, itNode + 1));
584 pPoints.push_back(PointOfField(coo[0], coo[1], coo[2], val[0]));
589 //-------------------------------------------------------------
590 // Case 2: field of elements
591 //-------------------------------------------------------------
593 if (mElements == NULL) throw IllegalStateException("no elements in the current mesh", __FILE__, __LINE__);
594 if (mElements->getTypeOfPrimitives() != MED_TETRA10) throw UnsupportedOperationException("only support TETRA10 mesh", __FILE__, __LINE__);
596 const string& nameGaussLoc = pField->getNameGaussLoc(pTimeStepIt);
597 GaussLoc* gaussLoc = getGaussLocByName(nameGaussLoc.c_str());
598 if (gaussLoc == NULL) throw IllegalStateException("no Gauss localization for these elements", __FILE__, __LINE__);
600 int numGauss = pField->getNumberOfGaussPointsByElement(pTimeStepIt);
602 int size = gaussLoc->getDim() * gaussLoc->getNumGaussPoints();
603 med_float* cooGaussPts = new med_float[size];
605 int dim = mElements->getTypeOfPrimitives() / 100;
606 int numNodes = mElements->getTypeOfPrimitives() % 100;
607 size = dim * numNodes;
608 med_float* cooNodes = new med_float[size];
611 for (int itElt = 0, size = mElements->getNumberOfElements() ; itElt < size ; itElt++)
613 // get coordinates of nodes of the current elements
614 // OPTIMIZATION: ASSUME TETRA10: ONLY GETS THE 4 FIRST NODES OF EACH ELEMENT
615 mElements->getCoordinates(itElt, mNodes, cooNodes, 4);
617 // compute coordinates of gauss points
618 gaussLoc->getCoordGaussPoints(cooNodes, cooGaussPts);
620 //printArray2D(cooGaussPts, 5, 3, "Gauss pt");
622 const med_float* val =
623 reinterpret_cast<const med_float*>(pField->getValue(pTimeStepIt, itElt + 1));
625 // for each point of Gauss of the element
626 med_float* srcCoo = cooGaussPts;
627 for (int itPtGauss = 0 ; itPtGauss < numGauss ; itPtGauss++)
629 pPoints.push_back(PointOfField(srcCoo[0], srcCoo[1], srcCoo[2], val[itPtGauss]));
635 delete[] cooGaussPts;
640 float Mesh::evalDefaultRadius(int pN) const
642 if (mFields.size() == 0) return 1.0f;
644 //---------------------------------------------------------------------
645 // Compute default radius
646 //---------------------------------------------------------------------
648 med_float volumeBBox =
649 (mMeshBBoxMax[0] - mMeshBBoxMin[0]) *
650 (mMeshBBoxMax[1] - mMeshBBoxMin[1]) *
651 (mMeshBBoxMax[2] - mMeshBBoxMin[2]);
653 if (isnan(volumeBBox))
658 const med_float k = 0.8; // considered 80% of the volume
660 // get nunmber of gauss points in the field
663 Field* anyField = mFields[mFields.size()-1];
664 int numTimeSteps = anyField->getNumberOfTimeSteps();
666 int numGaussPoints = getNumberOfElements() * anyField->getNumberOfGaussPointsByElement(numTimeSteps-1);
668 med_float radius = med_float(pow( (3.0/4.0) * pN * k * volumeBBox / (3.1415 * numGaussPoints), 1.0/3.0));
670 return float(radius);
679 med_geometrie_element CELL_TYPES[MED_NBR_GEOMETRIE_MAILLE] =
697 char CELL_NAMES[MED_NBR_GEOMETRIE_MAILLE][MED_TAILLE_NOM + 1] =
717 void Mesh::readSequentialMED(const char* pMEDfilename, const char* pMeshName)
721 //---------------------------------------------------------------------
723 //---------------------------------------------------------------------
724 MULTIPR_LOG("Check arguments: ");
725 if (pMEDfilename == NULL) throw NullArgumentException("pMEDfilename should not be NULL", __FILE__, __LINE__);
726 if (pMeshName == NULL) throw NullArgumentException("pMeshName should not be NULL", __FILE__, __LINE__);
729 strncpy(mMEDfilename, pMEDfilename, 256);
730 strncpy(mMeshName, pMeshName, MED_TAILLE_NOM);
732 //---------------------------------------------------------------------
733 // Open MED file (READ_ONLY)
734 //---------------------------------------------------------------------
735 MULTIPR_LOG("Open MED file: ");
736 mMEDfile = MEDouvrir(mMEDfilename, MED_LECTURE); // open MED file for reading
737 if (mMEDfile <= 0) throw FileNotFoundException("MED file not found", __FILE__, __LINE__);
740 //---------------------------------------------------------------------
741 // Check valid HDF format
742 //---------------------------------------------------------------------
743 MULTIPR_LOG("Format HDF: ");
744 if (MEDformatConforme(mMEDfilename) != 0) throw IOException("invalid file", __FILE__, __LINE__);
747 //---------------------------------------------------------------------
749 //---------------------------------------------------------------------
750 MULTIPR_LOG("MED version: ");
751 med_int verMajor, verMinor, verRelease;
752 med_err ret = MEDversionLire(mMEDfile, &verMajor, &verMinor, &verRelease);
753 if (ret != 0) throw IOException("error while reading MED version", __FILE__, __LINE__);
754 MULTIPR_LOG(verMajor << "." << verMinor << "." << verRelease << ": OK\n");
756 //---------------------------------------------------------------------
757 // Check that there is no profil
758 //---------------------------------------------------------------------
759 MULTIPR_LOG("#profils must be 0: ");
760 med_int nbProfils = MEDnProfil(mMEDfile);
761 if (nbProfils != 0) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
764 //---------------------------------------------------------------------
765 // Read all Gauss localizations
766 //---------------------------------------------------------------------
769 //---------------------------------------------------------------------
770 // Read scalars (should be 0)
771 //---------------------------------------------------------------------
772 MULTIPR_LOG("Scalars: ");
773 med_int nbScalars = MEDnScalaire(mMEDfile);
774 if (nbScalars != 0) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
775 MULTIPR_LOG(nbScalars << ": OK\n");
777 //---------------------------------------------------------------------
779 //---------------------------------------------------------------------
780 // read number of meshes
781 MULTIPR_LOG("Num meshes: ");
782 med_int nbMeshes = MEDnMaa(mMEDfile);
783 if (nbMeshes <= 0) throw IOException("i/o error while reading number of meshes in MED file", __FILE__, __LINE__);
784 MULTIPR_LOG(nbMeshes << ": OK\n");
786 med_int meshIndex = -1;
787 // iteration over mesh to find the mesh we want
788 // for each mesh in the file (warning: first mesh is number 1)
789 for (int itMesh = 1 ; itMesh <= nbMeshes ; itMesh++)
791 char meshName[MED_TAILLE_NOM + 1];
801 if (ret != 0) throw IOException("i/o error while reading mesh information in MED file", __FILE__, __LINE__);
802 MULTIPR_LOG("Mesh: |" << meshName << "|");
804 // test if the current mesh is the mesh we want
805 if (strcmp(pMeshName, meshName) == 0)
807 // *** mesh found ***
808 MULTIPR_LOG(" OK (found)\n");
814 // not the mesh we want: skip this mesh
815 MULTIPR_LOG(" skipped\n");
821 throw IllegalStateException("mesh not found in the given MED file", __FILE__, __LINE__);
824 //---------------------------------------------------------------------
825 // Check mesh validity
826 //---------------------------------------------------------------------
827 // dimension of the mesh must be 3 (= 3D mesh)
828 MULTIPR_LOG("Mesh is 3D: ");
829 if (mMeshDim != 3) throw UnsupportedOperationException("dimension of the mesh should be 3; other dimension not yet implemented", __FILE__, __LINE__);
832 // mesh must not be a grid
833 MULTIPR_LOG("Mesh is not a grid: ");
834 if (mMeshType != MED_NON_STRUCTURE)
835 throw UnsupportedOperationException("grid not supported", __FILE__, __LINE__);
838 // mesh must only contain TETRA10 elements
839 MULTIPR_LOG("Only TETRA10: ");
840 med_connectivite connectivite = MED_NOD; // NODAL CONNECTIVITY ONLY
841 bool onlyTETRA10 = true;
843 for (int itCell = 0 ; itCell < MED_NBR_GEOMETRIE_MAILLE ; itCell++)
845 med_int meshNumCells = MEDnEntMaa(
853 if ((meshNumCells > 0) && (strcmp(CELL_NAMES[itCell], "MED_TETRA10") != 0))
858 if (strcmp(CELL_NAMES[itCell], "MED_TETRA10") == 0)
860 numTetra10 = meshNumCells;
864 if (!onlyTETRA10) throw UnsupportedOperationException("mesh should only contain TETRA10 elements", __FILE__, __LINE__);
865 MULTIPR_LOG(numTetra10 << ": OK\n");
867 // everything is OK...
869 //---------------------------------------------------------------------
870 // Check num joint = 0
871 //---------------------------------------------------------------------
872 MULTIPR_LOG("Num joints: ");
873 med_int numJoints = MEDnJoint(mMEDfile, mMeshName);
874 MULTIPR_LOG(numJoints << ": OK\n");
876 //---------------------------------------------------------------------
877 // Check num equivalence = 0
878 //---------------------------------------------------------------------
879 MULTIPR_LOG("Num equivalences: ");
880 med_int numEquiv = MEDnEquiv(mMEDfile, mMeshName);
881 MULTIPR_LOG(numEquiv << ": OK\n");
883 //---------------------------------------------------------------------
885 //---------------------------------------------------------------------
886 mNodes = new Nodes();
887 mNodes->readMED(mMEDfile, mMeshName, mMeshDim);
888 mNodes->getBBox(mMeshBBoxMin, mMeshBBoxMax);
890 //---------------------------------------------------------------------
892 //---------------------------------------------------------------------
893 mElements = new Elements();
894 mElements->readMED(mMEDfile, mMeshName, mMeshDim, MED_MAILLE, MED_TETRA10);
896 //---------------------------------------------------------------------
898 //---------------------------------------------------------------------
900 finalizeFamiliesAndGroups();
902 //---------------------------------------------------------------------
904 //---------------------------------------------------------------------
907 //---------------------------------------------------------------------
908 // Close the MED file
909 //---------------------------------------------------------------------
910 MULTIPR_LOG("Close MED file: ");
911 ret = MEDfermer(mMEDfile);
912 if (ret != 0) throw IOException("i/o error while closing MED file", __FILE__, __LINE__);
917 void Mesh::writeMED(const char* pMEDfilename)
919 MULTIPR_LOG("Write MED: " << pMEDfilename << endl);
921 if (pMEDfilename == NULL) throw NullArgumentException("pMEDfilename should not be NULL", __FILE__, __LINE__);
922 if (strlen(pMEDfilename) == 0) throw IllegalArgumentException("pMEDfilename size is 0", __FILE__, __LINE__);
924 remove(pMEDfilename);
926 //---------------------------------------------------------------------
927 // Create the new MED file (WRITE_ONLY)
928 //---------------------------------------------------------------------
929 med_idt newMEDfile = MEDouvrir(const_cast<char*>(pMEDfilename), MED_CREATION);
930 if (newMEDfile == -1) throw IOException("", __FILE__, __LINE__);
932 //---------------------------------------------------------------------
934 //---------------------------------------------------------------------
935 // no scalars to write
937 //---------------------------------------------------------------------
938 // Create mesh: must be created first
939 //---------------------------------------------------------------------
940 med_err ret = MEDmaaCr(
947 if (ret != 0) throw IOException("", __FILE__, __LINE__);
948 MULTIPR_LOG(" Create mesh: |" << mMeshName << "|: OK" << endl);
950 //---------------------------------------------------------------------
951 // Write nodes and elements (mesh must exist)
952 //---------------------------------------------------------------------
953 mNodes->writeMED(newMEDfile, mMeshName);
954 MULTIPR_LOG(" Write nodes: ok" << endl);
955 mElements->writeMED(newMEDfile, mMeshName, mMeshDim);
956 MULTIPR_LOG(" write elt: ok" << endl);
958 //---------------------------------------------------------------------
959 // Write families (mesh must exist)
960 //---------------------------------------------------------------------
961 for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
963 Family* fam = mFamilies[itFam];
964 fam->writeMED(newMEDfile, mMeshName);
966 MULTIPR_LOG(" Write families: ok" << endl);
968 //---------------------------------------------------------------------
970 //---------------------------------------------------------------------
973 //---------------------------------------------------------------------
974 // Write Gauss localization (must be written before fields)
975 //---------------------------------------------------------------------
976 for (unsigned itGaussLoc = 0 ; itGaussLoc < mGaussLoc.size() ; itGaussLoc++)
979 GaussLoc* gaussLoc = mGaussLoc[itGaussLoc];
980 gaussLoc->writeMED(newMEDfile);
982 MULTIPR_LOG(" Write Gauss: ok" << endl);
984 //---------------------------------------------------------------------
986 //---------------------------------------------------------------------
987 for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
989 Field* field = mFields[itField];
990 field->writeMED(newMEDfile, mMeshName);
992 MULTIPR_LOG(" Write fields: ok" << endl);
994 //---------------------------------------------------------------------
995 // Close the new MED file
996 //---------------------------------------------------------------------
997 ret = MEDfermer(newMEDfile);
998 if (ret != 0) throw IOException("", __FILE__, __LINE__);
1002 void Mesh::readGaussLoc()
1004 MULTIPR_LOG("Gauss ref: ");
1005 med_int numGauss = MEDnGauss(mMEDfile);
1006 if (numGauss < 0) throw IOException("", __FILE__, __LINE__);
1007 MULTIPR_LOG(numGauss << ": OK\n");
1009 for (int itGauss = 1 ; itGauss <= numGauss ; itGauss++)
1011 GaussLoc* gaussLoc = new GaussLoc();
1012 gaussLoc->readMED(mMEDfile, itGauss);
1014 MULTIPR_LOG((*gaussLoc) << endl);
1016 mGaussLoc.push_back(gaussLoc);
1017 mGaussLocNameToGaussLoc.insert(make_pair(gaussLoc->getName(), gaussLoc));
1022 void Mesh::readFamilies()
1024 med_int numFamilies = MEDnFam(mMEDfile, mMeshName);
1025 if (numFamilies <= 0) throw IOException("", __FILE__, __LINE__);
1027 for (int itFam = 1 ; itFam <= numFamilies ; itFam++)
1029 Family* fam = new Family();
1030 fam->readMED(mMEDfile, mMeshName, itFam);
1031 mFamilies.push_back(fam);
1036 void Mesh::finalizeFamiliesAndGroups()
1038 //---------------------------------------------------------------------
1039 // Build mapping between family id and pointers towards families
1040 //---------------------------------------------------------------------
1041 for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
1043 Family* fam = mFamilies[itFam];
1044 mFamIdToFam.insert(make_pair(fam->getId(), fam));
1047 //---------------------------------------------------------------------
1048 // Fill families of nodes
1049 //---------------------------------------------------------------------
1050 for (int itNode = 1 ; itNode <= mNodes->getNumberOfNodes() ; itNode++)
1052 // get family of the ith nodes
1053 int famIdent = mNodes->getFamIdent(itNode - 1); // MED nodes start at 1
1054 map<med_int, Family*>::iterator itFam = mFamIdToFam.find(famIdent);
1056 if (itFam == mFamIdToFam.end()) throw IllegalStateException("", __FILE__, __LINE__);
1058 Family* fam = (*itFam).second;
1060 // insert the current node to its family
1061 fam->insertElt(itNode);
1062 fam->setIsFamilyOfNodes(true);
1065 //---------------------------------------------------------------------
1066 // Fill families of elements
1067 //---------------------------------------------------------------------
1068 for (int itElt = 1 ; itElt <= mElements->getNumberOfElements() ; itElt++)
1070 // get family of the ith element (MED index start at 1)
1071 int famIdent = mElements->getFamilyIdentifier(itElt - 1);
1072 map<med_int, Family*>::iterator itFam = mFamIdToFam.find(famIdent);
1074 if (itFam == mFamIdToFam.end()) throw IllegalStateException("", __FILE__, __LINE__);
1076 Family* fam = (*itFam).second;
1078 // insert the current node its family
1079 fam->insertElt(itElt);
1080 fam->setIsFamilyOfNodes(false);
1083 //---------------------------------------------------------------------
1085 //---------------------------------------------------------------------
1087 for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
1089 mFamilies[itFam]->buildGroups(mGroups, mGroupNameToGroup);
1094 void Mesh::readFields()
1096 //---------------------------------------------------------------------
1097 // Read number of fields
1098 //---------------------------------------------------------------------
1099 MULTIPR_LOG("Read fields: ");
1100 med_int numFields = MEDnChamp(mMEDfile, 0);
1101 if (numFields <= 0) throw IOException("", __FILE__, __LINE__);
1102 MULTIPR_LOG(numFields << ": OK\n");
1104 //---------------------------------------------------------------------
1105 // Iterate over fields
1106 //---------------------------------------------------------------------
1107 // for each field, read number of components and others infos
1108 for (int itField = 1 ; itField <= numFields ; itField++)
1110 Field* field = new Field();
1111 field->readMED(mMEDfile, itField, mMeshName);
1113 // if the nth field does not apply on our mesh => slip it
1114 if (field->isEmpty())
1120 mFields.push_back(field);
1126 ostream& operator<<(ostream& pOs, Mesh& pM)
1128 pOs << "Mesh: " << endl;
1129 pOs << " MED file =|" << pM.mMEDfilename << "|" << endl;
1130 pOs << " Name =|" << pM.mMeshName << "|" << endl;
1131 pOs << " Unv name =|" << pM.mMeshUName << "|" << endl;
1132 pOs << " Desc =|" << pM.mMeshDesc << "|" << endl;
1133 pOs << " Dim =" << pM.mMeshDim << endl;
1134 pOs << " Type =" << ((pM.mMeshType == MED_STRUCTURE)?"STRUCTURE":"NON_STRUCTURE") << endl;
1135 pOs << " BBox =[" << pM.mMeshBBoxMin[0] << " ; " << pM.mMeshBBoxMax[0] << "] x [" << pM.mMeshBBoxMin[1] << " ; " << pM.mMeshBBoxMax[1] << "] x [" << pM.mMeshBBoxMin[2] << " ; " << pM.mMeshBBoxMax[2] << "]" << endl;
1137 if (pM.mFlagPrintAll)
1139 cout << (*(pM.mNodes)) << endl;
1140 cout << (*(pM.mElements)) << endl;
1142 pOs << " Families : #=" << pM.mFamilies.size() << endl;
1143 for (unsigned i = 0 ; i < pM.mFamilies.size() ; i++)
1145 cout << (*(pM.mFamilies[i])) << endl;
1148 pOs << " Groups : #=" << pM.mGroups.size() << endl;
1149 for (unsigned i = 0 ; i < pM.mGroups.size() ; i++)
1151 cout << (*(pM.mGroups[i])) << endl;
1154 pOs << " Gauss loc: #=" << pM.mGaussLoc.size() << endl;
1155 for (unsigned i = 0 ; i < pM.mGaussLoc.size() ; i++)
1157 cout << (*(pM.mGaussLoc[i])) << endl;
1160 pOs << " Fields : #=" << pM.mFields.size() << endl;
1161 for (unsigned i = 0 ; i < pM.mFields.size() ; i++)
1163 cout << (*(pM.mFields[i])) << endl;
1168 pOs << " Nodes : #=" << pM.mNodes->getNumberOfNodes() << endl;
1170 const set<med_int>& setOfNodes = pM.mElements->getSetOfNodes();
1171 if (setOfNodes.size() == 0)
1173 pOs << " Elt : #=" << pM.mElements->getNumberOfElements() << endl;
1177 set<med_int>::iterator itNode = setOfNodes.end();
1179 pOs << " Elt : #=" << pM.mElements->getNumberOfElements() << " node_id_min=" << (*(setOfNodes.begin())) << " node_id_max=" << (*itNode) << endl;
1182 pOs << " Families : #=" << pM.mFamilies.size() << endl;
1183 pOs << " Groups : #=" << pM.mGroups.size() << endl;
1184 pOs << " Gauss loc: #=" << pM.mGaussLoc.size() << endl;
1185 pOs << " Fields : #=" << pM.mFields.size() << endl;
1192 } // namespace multipr