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 (unsigned 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 (unsigned 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();
470 meshDis->setSequentialMEDFilename(mMEDfilename);
472 // get prefix from the original MED filename
473 string strPrefix = removeExtension(mMEDfilename, ".med");
478 for (unsigned itGroup = 0 ; itGroup < mGroups.size() ; itGroup++)
480 Group* currentGroup = mGroups[itGroup];
482 // skip this group if it is a group of nodes
483 if (currentGroup->isGroupOfNodes())
488 char strPartName[256];
489 sprintf(strPartName, "%s_%d", mMeshName, numGroup);
491 char strMEDfilename[256];
492 sprintf(strMEDfilename, "%s_grain%d.med", strPrefix.c_str(), numGroup);
494 Mesh* mesh = createFromGroup(currentGroup, mMeshName);
496 // skip the group which contain all the others groups, even it contains only 1 group
497 if ((mesh->mElements->getNumberOfElements() == mElements->getNumberOfElements()) && (mGroups.size() > 1))
504 MeshDisPart::MULTIPR_WRITE_MESH,
519 Mesh* Mesh::decimate(
520 const char* pFilterName,
522 const char* pNameNewMesh)
524 //---------------------------------------------------------------------
526 //---------------------------------------------------------------------
527 if (pFilterName == NULL) throw NullArgumentException("pFilterName should not be NULL", __FILE__, __LINE__);
528 if (pArgv == NULL) throw NullArgumentException("pArgv should not be NULL", __FILE__, __LINE__);
529 if (pNameNewMesh == NULL) throw NullArgumentException("pNameNewMesh should not be NULL", __FILE__, __LINE__);
531 //---------------------------------------------------------------------
532 // Instanciate filter used for decimation
533 //---------------------------------------------------------------------
534 DecimationFilter* filter = DecimationFilter::create(pFilterName);
536 //---------------------------------------------------------------------
537 // Create new mesh by decimating current one
538 //---------------------------------------------------------------------
539 Mesh* decimatedMesh = filter->apply(this, pArgv, pNameNewMesh);
541 //---------------------------------------------------------------------
543 //---------------------------------------------------------------------
546 return decimatedMesh;
551 void Mesh::getAllPointsOfField(Field* pField, int pTimeStepIt, std::vector<PointOfField>& pPoints)
553 //---------------------------------------------------------------------
555 //---------------------------------------------------------------------
557 if (pField == NULL) throw NullArgumentException("field should not be NULL", __FILE__, __LINE__);
558 if (pTimeStepIt < 1) throw IllegalArgumentException("invalid field iteration; should be >= 1", __FILE__, __LINE__);
560 if (mMeshDim != 3) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
561 if (pField->getType() != MED_FLOAT64) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
562 if (pField->getNumberOfComponents() != 1) throw UnsupportedOperationException("field have more than 1 component (vectorial field, expected scalar field)", __FILE__, __LINE__);
564 //---------------------------------------------------------------------
566 //---------------------------------------------------------------------
568 if (pField->isFieldOnNodes())
570 //-------------------------------------------------------------
571 // Case 1: field of nodes
572 //-------------------------------------------------------------
573 if (mNodes == NULL) throw IllegalStateException("no nodes in the current mesh", __FILE__, __LINE__);
576 for (int itNode = 0, size = mNodes->getNumberOfNodes() ; itNode < size ; itNode++)
578 // collect coordinates and value of the point
579 const med_float* coo = mNodes->getCoordinates(itNode);
581 const med_float* val =
582 reinterpret_cast<const med_float*>(pField->getValue(pTimeStepIt, itNode + 1));
585 pPoints.push_back(PointOfField(coo[0], coo[1], coo[2], val[0]));
590 //-------------------------------------------------------------
591 // Case 2: field of elements
592 //-------------------------------------------------------------
594 if (mElements == NULL) throw IllegalStateException("no elements in the current mesh", __FILE__, __LINE__);
595 if (mElements->getTypeOfPrimitives() != MED_TETRA10) throw UnsupportedOperationException("only support TETRA10 mesh", __FILE__, __LINE__);
597 const string& nameGaussLoc = pField->getNameGaussLoc(pTimeStepIt);
598 GaussLoc* gaussLoc = getGaussLocByName(nameGaussLoc.c_str());
599 if (gaussLoc == NULL) throw IllegalStateException("no Gauss localization for these elements", __FILE__, __LINE__);
601 int numGauss = pField->getNumberOfGaussPointsByElement(pTimeStepIt);
603 int size = gaussLoc->getDim() * gaussLoc->getNumGaussPoints();
604 med_float* cooGaussPts = new med_float[size];
606 int dim = mElements->getTypeOfPrimitives() / 100;
607 int numNodes = mElements->getTypeOfPrimitives() % 100;
608 size = dim * numNodes;
609 med_float* cooNodes = new med_float[size];
612 for (int itElt = 0, size = mElements->getNumberOfElements() ; itElt < size ; itElt++)
614 // get coordinates of nodes of the current elements
615 // OPTIMIZATION: ASSUME TETRA10: ONLY GETS THE 4 FIRST NODES OF EACH ELEMENT
616 mElements->getCoordinates(itElt, mNodes, cooNodes, 4);
618 // compute coordinates of gauss points
619 gaussLoc->getCoordGaussPoints(cooNodes, cooGaussPts);
621 //printArray2D(cooGaussPts, 5, 3, "Gauss pt"); // debug
623 const med_float* val =
624 reinterpret_cast<const med_float*>(pField->getValue(pTimeStepIt, itElt + 1));
626 // for each point of Gauss of the element
627 med_float* srcCoo = cooGaussPts;
628 for (int itPtGauss = 0 ; itPtGauss < numGauss ; itPtGauss++)
630 pPoints.push_back(PointOfField(srcCoo[0], srcCoo[1], srcCoo[2], val[itPtGauss]));
636 delete[] cooGaussPts;
641 float Mesh::evalDefaultRadius(int pN) const
643 if (mFields.size() == 0) return 1.0f;
645 //---------------------------------------------------------------------
646 // Compute default radius
647 //---------------------------------------------------------------------
649 med_float volumeBBox =
650 (mMeshBBoxMax[0] - mMeshBBoxMin[0]) *
651 (mMeshBBoxMax[1] - mMeshBBoxMin[1]) *
652 (mMeshBBoxMax[2] - mMeshBBoxMin[2]);
654 if (isnan(volumeBBox))
659 const med_float k = 0.8; // considered 80% of the volume
661 // get nunmber of gauss points in the field
664 Field* anyField = mFields[mFields.size()-1];
665 int numTimeSteps = anyField->getNumberOfTimeSteps();
667 int numGaussPoints = getNumberOfElements() * anyField->getNumberOfGaussPointsByElement(numTimeSteps-1);
669 med_float radius = med_float(pow( (3.0/4.0) * pN * k * volumeBBox / (3.1415 * numGaussPoints), 1.0/3.0));
671 return float(radius);
680 med_geometrie_element CELL_TYPES[MED_NBR_GEOMETRIE_MAILLE] =
698 char CELL_NAMES[MED_NBR_GEOMETRIE_MAILLE][MED_TAILLE_NOM + 1] =
718 void Mesh::readSequentialMED(const char* pMEDfilename, const char* pMeshName)
722 //---------------------------------------------------------------------
724 //---------------------------------------------------------------------
725 MULTIPR_LOG("Check arguments: ");
726 if (pMEDfilename == NULL) throw NullArgumentException("pMEDfilename should not be NULL", __FILE__, __LINE__);
727 if (pMeshName == NULL) throw NullArgumentException("pMeshName should not be NULL", __FILE__, __LINE__);
730 strncpy(mMEDfilename, pMEDfilename, 256);
731 strncpy(mMeshName, pMeshName, MED_TAILLE_NOM);
733 //---------------------------------------------------------------------
734 // Open MED file (READ_ONLY)
735 //---------------------------------------------------------------------
736 MULTIPR_LOG("Open MED file: ");
737 mMEDfile = MEDouvrir(mMEDfilename, MED_LECTURE); // open MED file for reading
738 if (mMEDfile <= 0) throw FileNotFoundException("MED file not found", __FILE__, __LINE__);
741 //---------------------------------------------------------------------
742 // Check valid HDF format
743 //---------------------------------------------------------------------
744 MULTIPR_LOG("Format HDF: ");
745 if (MEDformatConforme(mMEDfilename) != 0) throw IOException("invalid file", __FILE__, __LINE__);
748 //---------------------------------------------------------------------
750 //---------------------------------------------------------------------
751 MULTIPR_LOG("MED version: ");
752 med_int verMajor, verMinor, verRelease;
753 med_err ret = MEDversionLire(mMEDfile, &verMajor, &verMinor, &verRelease);
754 if (ret != 0) throw IOException("error while reading MED version", __FILE__, __LINE__);
755 MULTIPR_LOG(verMajor << "." << verMinor << "." << verRelease << ": OK\n");
757 //---------------------------------------------------------------------
758 // Check that there is no profil
759 //---------------------------------------------------------------------
760 MULTIPR_LOG("#profils must be 0: ");
761 med_int nbProfils = MEDnProfil(mMEDfile);
762 if (nbProfils != 0) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
765 //---------------------------------------------------------------------
766 // Read all Gauss localizations
767 //---------------------------------------------------------------------
770 //---------------------------------------------------------------------
771 // Read scalars (should be 0)
772 //---------------------------------------------------------------------
773 MULTIPR_LOG("Scalars: ");
774 med_int nbScalars = MEDnScalaire(mMEDfile);
775 if (nbScalars != 0) throw UnsupportedOperationException("not yet implemented", __FILE__, __LINE__);
776 MULTIPR_LOG(nbScalars << ": OK\n");
778 //---------------------------------------------------------------------
780 //---------------------------------------------------------------------
781 // read number of meshes
782 MULTIPR_LOG("Num meshes: ");
783 med_int nbMeshes = MEDnMaa(mMEDfile);
784 if (nbMeshes <= 0) throw IOException("i/o error while reading number of meshes in MED file", __FILE__, __LINE__);
785 MULTIPR_LOG(nbMeshes << ": OK\n");
787 med_int meshIndex = -1;
788 // iteration over mesh to find the mesh we want
789 // for each mesh in the file (warning: first mesh is number 1)
790 for (int itMesh = 1 ; itMesh <= nbMeshes ; itMesh++)
792 char meshName[MED_TAILLE_NOM + 1];
802 if (ret != 0) throw IOException("i/o error while reading mesh information in MED file", __FILE__, __LINE__);
803 MULTIPR_LOG("Mesh: |" << meshName << "|");
805 // test if the current mesh is the mesh we want
806 if (strcmp(pMeshName, meshName) == 0)
808 // *** mesh found ***
809 MULTIPR_LOG(" OK (found)\n");
815 // not the mesh we want: skip this mesh
816 MULTIPR_LOG(" skipped\n");
822 throw IllegalStateException("mesh not found in the given MED file", __FILE__, __LINE__);
825 //---------------------------------------------------------------------
826 // Check mesh validity
827 //---------------------------------------------------------------------
828 // dimension of the mesh must be 3 (= 3D mesh)
829 MULTIPR_LOG("Mesh is 3D: ");
830 if (mMeshDim != 3) throw UnsupportedOperationException("dimension of the mesh should be 3; other dimension not yet implemented", __FILE__, __LINE__);
833 // mesh must not be a grid
834 MULTIPR_LOG("Mesh is not a grid: ");
835 if (mMeshType != MED_NON_STRUCTURE)
836 throw UnsupportedOperationException("grid not supported", __FILE__, __LINE__);
839 // mesh must only contain TETRA10 elements
840 MULTIPR_LOG("Only TETRA10: ");
841 med_connectivite connectivite = MED_NOD; // NODAL CONNECTIVITY ONLY
842 bool onlyTETRA10 = true;
844 for (int itCell = 0 ; itCell < MED_NBR_GEOMETRIE_MAILLE ; itCell++)
846 med_int meshNumCells = MEDnEntMaa(
854 if ((meshNumCells > 0) && (strcmp(CELL_NAMES[itCell], "MED_TETRA10") != 0))
859 if (strcmp(CELL_NAMES[itCell], "MED_TETRA10") == 0)
861 numTetra10 = meshNumCells;
865 if (!onlyTETRA10) throw UnsupportedOperationException("mesh should only contain TETRA10 elements", __FILE__, __LINE__);
866 MULTIPR_LOG(numTetra10 << ": OK\n");
868 // everything is OK...
870 //---------------------------------------------------------------------
871 // Check num joint = 0
872 //---------------------------------------------------------------------
873 MULTIPR_LOG("Num joints: ");
874 med_int numJoints = MEDnJoint(mMEDfile, mMeshName);
875 MULTIPR_LOG(numJoints << ": OK\n");
877 //---------------------------------------------------------------------
878 // Check num equivalence = 0
879 //---------------------------------------------------------------------
880 MULTIPR_LOG("Num equivalences: ");
881 med_int numEquiv = MEDnEquiv(mMEDfile, mMeshName);
882 MULTIPR_LOG(numEquiv << ": OK\n");
884 //---------------------------------------------------------------------
886 //---------------------------------------------------------------------
887 mNodes = new Nodes();
888 mNodes->readMED(mMEDfile, mMeshName, mMeshDim);
889 mNodes->getBBox(mMeshBBoxMin, mMeshBBoxMax);
891 //---------------------------------------------------------------------
893 //---------------------------------------------------------------------
894 mElements = new Elements();
895 mElements->readMED(mMEDfile, mMeshName, mMeshDim, MED_MAILLE, MED_TETRA10);
897 //---------------------------------------------------------------------
899 //---------------------------------------------------------------------
901 finalizeFamiliesAndGroups();
903 //---------------------------------------------------------------------
905 //---------------------------------------------------------------------
908 //---------------------------------------------------------------------
909 // Close the MED file
910 //---------------------------------------------------------------------
911 MULTIPR_LOG("Close MED file: ");
912 ret = MEDfermer(mMEDfile);
913 if (ret != 0) throw IOException("i/o error while closing MED file", __FILE__, __LINE__);
918 void Mesh::writeMED(const char* pMEDfilename)
920 MULTIPR_LOG("Write MED: " << pMEDfilename << endl);
922 if (pMEDfilename == NULL) throw NullArgumentException("pMEDfilename should not be NULL", __FILE__, __LINE__);
923 if (strlen(pMEDfilename) == 0) throw IllegalArgumentException("pMEDfilename size is 0", __FILE__, __LINE__);
925 remove(pMEDfilename);
927 //---------------------------------------------------------------------
928 // Create the new MED file (WRITE_ONLY)
929 //---------------------------------------------------------------------
930 med_idt newMEDfile = MEDouvrir(const_cast<char*>(pMEDfilename), MED_CREATION);
931 if (newMEDfile == -1) throw IOException("", __FILE__, __LINE__);
933 //---------------------------------------------------------------------
935 //---------------------------------------------------------------------
936 // no scalars to write
938 //---------------------------------------------------------------------
939 // Create mesh: must be created first
940 //---------------------------------------------------------------------
941 med_err ret = MEDmaaCr(
948 if (ret != 0) throw IOException("", __FILE__, __LINE__);
949 MULTIPR_LOG(" Create mesh: |" << mMeshName << "|: OK" << endl);
951 //---------------------------------------------------------------------
952 // Write nodes and elements (mesh must exist)
953 //---------------------------------------------------------------------
954 mNodes->writeMED(newMEDfile, mMeshName);
955 MULTIPR_LOG(" Write nodes: ok" << endl);
956 mElements->writeMED(newMEDfile, mMeshName, mMeshDim);
957 MULTIPR_LOG(" write elt: ok" << endl);
959 //---------------------------------------------------------------------
960 // Write families (mesh must exist)
961 //---------------------------------------------------------------------
962 for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
964 Family* fam = mFamilies[itFam];
965 fam->writeMED(newMEDfile, mMeshName);
967 MULTIPR_LOG(" Write families: ok" << endl);
969 //---------------------------------------------------------------------
971 //---------------------------------------------------------------------
974 //---------------------------------------------------------------------
975 // Write Gauss localization (must be written before fields)
976 //---------------------------------------------------------------------
977 for (unsigned itGaussLoc = 0 ; itGaussLoc < mGaussLoc.size() ; itGaussLoc++)
980 GaussLoc* gaussLoc = mGaussLoc[itGaussLoc];
981 gaussLoc->writeMED(newMEDfile);
983 MULTIPR_LOG(" Write Gauss: ok" << endl);
985 //---------------------------------------------------------------------
987 //---------------------------------------------------------------------
988 for (unsigned itField = 0 ; itField < mFields.size() ; itField++)
990 Field* field = mFields[itField];
991 field->writeMED(newMEDfile, mMeshName);
993 MULTIPR_LOG(" Write fields: ok" << endl);
995 //---------------------------------------------------------------------
996 // Close the new MED file
997 //---------------------------------------------------------------------
998 ret = MEDfermer(newMEDfile);
999 if (ret != 0) throw IOException("", __FILE__, __LINE__);
1003 void Mesh::readGaussLoc()
1005 MULTIPR_LOG("Gauss ref: ");
1006 med_int numGauss = MEDnGauss(mMEDfile);
1007 if (numGauss < 0) throw IOException("", __FILE__, __LINE__);
1008 MULTIPR_LOG(numGauss << ": OK\n");
1010 for (int itGauss = 1 ; itGauss <= numGauss ; itGauss++)
1012 GaussLoc* gaussLoc = new GaussLoc();
1013 gaussLoc->readMED(mMEDfile, itGauss);
1015 MULTIPR_LOG((*gaussLoc) << endl);
1017 mGaussLoc.push_back(gaussLoc);
1018 mGaussLocNameToGaussLoc.insert(make_pair(gaussLoc->getName(), gaussLoc));
1023 void Mesh::readFamilies()
1025 med_int numFamilies = MEDnFam(mMEDfile, mMeshName);
1026 if (numFamilies <= 0) throw IOException("", __FILE__, __LINE__);
1028 for (int itFam = 1 ; itFam <= numFamilies ; itFam++)
1030 Family* fam = new Family();
1031 fam->readMED(mMEDfile, mMeshName, itFam);
1032 mFamilies.push_back(fam);
1037 void Mesh::finalizeFamiliesAndGroups()
1039 //---------------------------------------------------------------------
1040 // Build mapping between family id and pointers towards families
1041 //---------------------------------------------------------------------
1042 for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
1044 Family* fam = mFamilies[itFam];
1045 mFamIdToFam.insert(make_pair(fam->getId(), fam));
1048 //---------------------------------------------------------------------
1049 // Fill families of nodes
1050 //---------------------------------------------------------------------
1051 for (int itNode = 1 ; itNode <= mNodes->getNumberOfNodes() ; itNode++)
1053 // get family of the ith nodes
1054 int famIdent = mNodes->getFamIdent(itNode - 1); // MED nodes start at 1
1055 map<med_int, Family*>::iterator itFam = mFamIdToFam.find(famIdent);
1057 if (itFam == mFamIdToFam.end()) throw IllegalStateException("", __FILE__, __LINE__);
1059 Family* fam = (*itFam).second;
1061 // insert the current node to its family
1062 fam->insertElt(itNode);
1063 fam->setIsFamilyOfNodes(true);
1066 //---------------------------------------------------------------------
1067 // Fill families of elements
1068 //---------------------------------------------------------------------
1069 for (int itElt = 1 ; itElt <= mElements->getNumberOfElements() ; itElt++)
1071 // get family of the ith element (MED index start at 1)
1072 int famIdent = mElements->getFamilyIdentifier(itElt - 1);
1073 map<med_int, Family*>::iterator itFam = mFamIdToFam.find(famIdent);
1075 if (itFam == mFamIdToFam.end()) throw IllegalStateException("", __FILE__, __LINE__);
1077 Family* fam = (*itFam).second;
1079 // insert the current node its family
1080 fam->insertElt(itElt);
1081 fam->setIsFamilyOfNodes(false);
1084 //---------------------------------------------------------------------
1086 //---------------------------------------------------------------------
1088 for (unsigned itFam = 0 ; itFam < mFamilies.size() ; itFam++)
1090 mFamilies[itFam]->buildGroups(mGroups, mGroupNameToGroup);
1095 void Mesh::readFields()
1097 //---------------------------------------------------------------------
1098 // Read number of fields
1099 //---------------------------------------------------------------------
1100 MULTIPR_LOG("Read fields: ");
1101 med_int numFields = MEDnChamp(mMEDfile, 0);
1102 if (numFields <= 0) throw IOException("", __FILE__, __LINE__);
1103 MULTIPR_LOG(numFields << ": OK\n");
1105 //---------------------------------------------------------------------
1106 // Iterate over fields
1107 //---------------------------------------------------------------------
1108 // for each field, read number of components and others infos
1109 for (int itField = 1 ; itField <= numFields ; itField++)
1111 Field* field = new Field();
1112 field->readMED(mMEDfile, itField, mMeshName);
1114 // if the nth field does not apply on our mesh => slip it
1115 if (field->isEmpty())
1121 mFields.push_back(field);
1127 ostream& operator<<(ostream& pOs, Mesh& pM)
1129 pOs << "Mesh: " << endl;
1130 pOs << " MED file =|" << pM.mMEDfilename << "|" << endl;
1131 pOs << " Name =|" << pM.mMeshName << "|" << endl;
1132 pOs << " Unv name =|" << pM.mMeshUName << "|" << endl;
1133 pOs << " Desc =|" << pM.mMeshDesc << "|" << endl;
1134 pOs << " Dim =" << pM.mMeshDim << endl;
1135 pOs << " Type =" << ((pM.mMeshType == MED_STRUCTURE)?"STRUCTURE":"NON_STRUCTURE") << endl;
1136 pOs << " BBox =[" << pM.mMeshBBoxMin[0] << " ; " << pM.mMeshBBoxMax[0] << "] x [" << pM.mMeshBBoxMin[1] << " ; " << pM.mMeshBBoxMax[1] << "] x [" << pM.mMeshBBoxMin[2] << " ; " << pM.mMeshBBoxMax[2] << "]" << endl;
1138 if (pM.mFlagPrintAll)
1140 cout << (*(pM.mNodes)) << endl;
1141 cout << (*(pM.mElements)) << endl;
1143 pOs << " Families : #=" << pM.mFamilies.size() << endl;
1144 for (unsigned i = 0 ; i < pM.mFamilies.size() ; i++)
1146 cout << (*(pM.mFamilies[i])) << endl;
1149 pOs << " Groups : #=" << pM.mGroups.size() << endl;
1150 for (unsigned i = 0 ; i < pM.mGroups.size() ; i++)
1152 cout << (*(pM.mGroups[i])) << endl;
1155 pOs << " Gauss loc: #=" << pM.mGaussLoc.size() << endl;
1156 for (unsigned i = 0 ; i < pM.mGaussLoc.size() ; i++)
1158 cout << (*(pM.mGaussLoc[i])) << endl;
1161 pOs << " Fields : #=" << pM.mFields.size() << endl;
1162 for (unsigned i = 0 ; i < pM.mFields.size() ; i++)
1164 cout << (*(pM.mFields[i])) << endl;
1169 pOs << " Nodes : #=" << pM.mNodes->getNumberOfNodes() << endl;
1171 const set<med_int>& setOfNodes = pM.mElements->getSetOfNodes();
1172 if (setOfNodes.size() == 0)
1174 pOs << " Elt : #=" << pM.mElements->getNumberOfElements() << endl;
1178 set<med_int>::iterator itNode = setOfNodes.end();
1180 pOs << " Elt : #=" << pM.mElements->getNumberOfElements() << " node_id_min=" << (*(setOfNodes.begin())) << " node_id_max=" << (*itNode) << endl;
1183 pOs << " Families : #=" << pM.mFamilies.size() << endl;
1184 pOs << " Groups : #=" << pM.mGroups.size() << endl;
1185 pOs << " Gauss loc: #=" << pM.mGaussLoc.size() << endl;
1186 pOs << " Fields : #=" << pM.mFields.size() << endl;
1193 } // namespace multipr