v0.13.1
Classes | Public Member Functions | Public Attributes | Private Member Functions | Private Attributes | List of all members
MoFEM::MedInterface Struct Reference

Interface for load MED files. More...

#include <src/interfaces/MedInterface.hpp>

Inheritance diagram for MoFEM::MedInterface:
[legend]
Collaboration diagram for MoFEM::MedInterface:
[legend]

Classes

struct  FieldData
 

Public Member Functions

MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const
 
 MedInterface (const MoFEM::Core &core)
 
MoFEMErrorCode getFileNameFromCommandLine (int verb=1)
 Get MED file name from command line. More...
 
PetscBool getFlgFile () const
 Check if file name is given in command line. More...
 
MoFEMErrorCode medGetFieldNames (const string &file, int verb=1)
 Get field names in MED file. More...
 
MoFEMErrorCode medGetFieldNames (int verb=1)
 get field names in MED file More...
 
MoFEMErrorCode readMed (const string &file, int verb=1)
 read MED file More...
 
MoFEMErrorCode readMed (int verb=1)
 read MED file More...
 
MoFEMErrorCode writeMed (const string &file, int verb=1)
 write MED file More...
 
MoFEMErrorCode readFields (const std::string &file_name, const std::string &field_name, const bool load_series=false, const int only_step=-1, int verb=1)
 
- Public Member Functions inherited from MoFEM::UnknownInterface
virtual MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
 
template<class IFACE >
MoFEMErrorCode registerInterface (bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()=default
 

Public Attributes

std::vector< std::string > meshNames
 list of meshes in MED file More...
 
std::vector< EntityHandle > meshMeshsets
 meshset for each mesh More...
 
std::map< std::string, FieldDatafieldNames
 
std::string medFileName
 MED file name. More...
 

Private Member Functions

MoFEMErrorCode readMesh (const string &file, const int index, std::map< int, Range > &family_elem_map, int verb=1)
 read mesh from MED file More...
 
MoFEMErrorCode readFamily (const string &file, const int index, const std::map< int, Range > &family_elem_mapint, std::map< string, Range > &group_elem_map, int verb=1)
 read family and groups More...
 
MoFEMErrorCode makeBlockSets (const std::map< string, Range > &group_elem_map, int verb=1)
 make from groups meshsets More...
 

Private Attributes

MoFEM::CorecOre
 core database More...
 
PetscBool flgFile
 true if file name given in command line More...
 

Additional Inherited Members

- Static Public Member Functions inherited from MoFEM::UnknownInterface
static MoFEMErrorCode getLibVersion (Version &version)
 Get library version. More...
 
static MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version)
 Get database major version. More...
 
static MoFEMErrorCode setFileVersion (moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
 Get database major version. More...
 
static MoFEMErrorCode getInterfaceVersion (Version &version)
 Get database major version. More...
 

Detailed Description

Interface for load MED files.

Definition at line 23 of file MedInterface.hpp.

Constructor & Destructor Documentation

◆ MedInterface()

MoFEM::MedInterface::MedInterface ( const MoFEM::Core core)

Definition at line 47 of file MedInterface.cpp.

48 : cOre(const_cast<Core &>(core)), flgFile(PETSC_FALSE) {
49
50 if (!LogManager::checkIfChannelExist("MEDWORLD")) {
51 auto core_log = logging::core::get();
52 core_log->add_sink(
54 LogManager::setLog("MEDWORLD");
55 MOFEM_LOG_TAG("MEDWORLD", "MED");
56 }
57
58 MOFEM_LOG("MEDWORLD", Sev::noisy) << "Mashset manager created";
59}
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
CoreTmp< 0 > Core
Definition: Core.hpp:1086
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
static bool checkIfChannelExist(const std::string channel)
Check if channel exist.
Definition: LogManager.cpp:374
MoFEM::Core & cOre
core database
PetscBool flgFile
true if file name given in command line

Member Function Documentation

◆ getFileNameFromCommandLine()

MoFEMErrorCode MoFEM::MedInterface::getFileNameFromCommandLine ( int  verb = 1)

Get MED file name from command line.

Parameters
verbverbosity level
Returns
error code

Definition at line 61 of file MedInterface.cpp.

61 {
62 Interface &m_field = cOre;
63 char mesh_file_name[255];
65 ierr = PetscOptionsBegin(m_field.get_comm(), "", "MED Interface", "none");
67 ierr = PetscOptionsString("-med_file", "med file name", "", "mesh.med",
68 mesh_file_name, 255, &flgFile);
70 ierr = PetscOptionsEnd();
72 if (flgFile == PETSC_TRUE) {
73 medFileName = std::string(mesh_file_name);
74 } else {
75 medFileName = std::string("mesh.med");
76 }
78}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
std::string medFileName
MED file name.

◆ getFlgFile()

PetscBool MoFEM::MedInterface::getFlgFile ( ) const

Check if file name is given in command line.

Definition at line 39 of file MedInterface.hpp.

39{ return flgFile; }

◆ makeBlockSets()

MoFEMErrorCode MoFEM::MedInterface::makeBlockSets ( const std::map< string, Range > &  group_elem_map,
int  verb = 1 
)
private

make from groups meshsets

Parameters
group_elem_mapmap of groups and elements
verbverbosity level
Returns
error code

Definition at line 612 of file MedInterface.cpp.

613 {
614
615 Interface &m_field = cOre;
617 MeshsetsManager *meshsets_manager_ptr;
618 CHKERR m_field.getInterface(meshsets_manager_ptr);
619
620 int max_id = 0;
622 max_id = (max_id < cit->getMeshsetId()) ? cit->getMeshsetId() : max_id;
623 }
624 max_id++;
625
626 for (std::map<string, Range>::const_iterator git = group_elem_map.begin();
627 git != group_elem_map.end(); git++) {
628 CHKERR meshsets_manager_ptr->addMeshset(BLOCKSET, max_id, git->first);
630 Composite_Cubit_msId_And_MeshsetType_mi_tag>::type::iterator cit;
631 cit =
632 meshsets_manager_ptr->getMeshsetsMultindex()
633 .get<Composite_Cubit_msId_And_MeshsetType_mi_tag>()
634 .find(boost::make_tuple(max_id, CubitBCType(BLOCKSET).to_ulong()));
635 EntityHandle meshsets = cit->getMeshset();
636 if (!git->second.empty()) {
637 CHKERR m_field.get_moab().add_entities(meshsets, git->second);
638 }
639 max_id++;
640 }
641
643 MOFEM_LOG("MEDWORLD", Sev::verbose) << *cit;
644
646}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
std::bitset< 32 > CubitBCType
Definition: Types.hpp:52

◆ medGetFieldNames() [1/2]

MoFEMErrorCode MoFEM::MedInterface::medGetFieldNames ( const string &  file,
int  verb = 1 
)

Get field names in MED file.

Parameters
filefile name
verbverbosity level
Returns
error code

Definition at line 80 of file MedInterface.cpp.

80 {
81 Interface &m_field = cOre;
83 med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
84 if (fid < 0) {
85 SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
86 "Unable to open file '%s'", file.c_str());
87 }
88 med_int num_fields = MEDnField(fid);
89 for (int index = 0; index < num_fields; index++) {
90
91 med_int num_comp = MEDfieldnComponent(fid, index + 1);
92 if (num_comp <= 0) {
93 SETERRQ(m_field.get_comm(), MOFEM_IMPOSIBLE_CASE,
94 "Could not get number of components for MED field");
95 }
96
97 char name[MED_NAME_SIZE + 1], mesh_name[MED_NAME_SIZE + 1];
98 char dt_unit[MED_SNAME_SIZE + 1];
99 std::vector<char> comp_name(num_comp * MED_SNAME_SIZE + 1);
100 std::vector<char> comp_unit(num_comp * MED_SNAME_SIZE + 1);
101 med_int num_steps = 0;
102 med_type_champ type;
103 med_bool local_mesh;
104 if (MEDfieldInfo(fid, index + 1, name, mesh_name, &local_mesh, &type,
105 &comp_name[0], &comp_unit[0], dt_unit, &num_steps) < 0) {
106 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
107 "Could not get MED field info");
108 }
109
110 std::string field_name = std::string(name);
112 fieldNames[field_name].fieldName = field_name;
113 fieldNames[field_name].meshName = std::string(mesh_name);
114 fieldNames[field_name].localMesh = (local_mesh == MED_TRUE);
115 for (int ff = 0; ff != num_comp; ff++) {
116 fieldNames[field_name].componentNames.push_back(
117 std::string(&comp_name[ff * MED_SNAME_SIZE], MED_SNAME_SIZE));
118 fieldNames[field_name].componentUnits.push_back(
119 std::string(&comp_unit[ff * MED_SNAME_SIZE], MED_SNAME_SIZE));
120 }
121 fieldNames[field_name].dtUnit = std::string(&dt_unit[0]);
122 fieldNames[field_name].ncSteps = num_steps;
123
124 if (verb > 0)
125 MOFEM_LOG("MEDWORLD", Sev::inform) << fieldNames[name];
126 }
127 if (MEDfileClose(fid) < 0) {
128 SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
129 "Unable to close file '%s'", file.c_str());
130 }
132}
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:35
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
double FieldData
Field data type.
Definition: Types.hpp:25
constexpr auto field_name
std::map< std::string, FieldData > fieldNames

◆ medGetFieldNames() [2/2]

MoFEMErrorCode MoFEM::MedInterface::medGetFieldNames ( int  verb = 1)

get field names in MED file

File name is get form command line

Parameters
verbverbosity level
Returns
error code

Definition at line 134 of file MedInterface.cpp.

134 {
136 if (medFileName.empty()) {
138 }
141}
MoFEMErrorCode getFileNameFromCommandLine(int verb=1)
Get MED file name from command line.
MoFEMErrorCode medGetFieldNames(const string &file, int verb=1)
Get field names in MED file.

◆ query_interface()

MoFEMErrorCode MoFEM::MedInterface::query_interface ( boost::typeindex::type_index  type_index,
UnknownInterface **  iface 
) const
virtual

Implements MoFEM::UnknownInterface.

Definition at line 41 of file MedInterface.cpp.

42 {
43 *iface = const_cast<MedInterface *>(this);
44 return 0;
45}
MedInterface(const MoFEM::Core &core)

◆ readFamily()

MoFEMErrorCode MoFEM::MedInterface::readFamily ( const string &  file,
const int  index,
const std::map< int, Range > &  family_elem_mapint,
std::map< string, Range > &  group_elem_map,
int  verb = 1 
)
private

read family and groups

Parameters
filefile name
indexmesh index
family_elem_mapintmap of families and elements
group_elem_mapmap of groups and elements
verbverbosity level
Returns
error code

Definition at line 517 of file MedInterface.cpp.

519 {
520 //
521 Interface &m_field = cOre;
523
524 med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
525 if (fid < 0) {
526 SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
527 "Unable to open file '%s'", file.c_str());
528 }
529 med_int v[3], vf[3];
530 MEDlibraryNumVersion(&v[0], &v[1], &v[2]);
531 MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
532
533 if (verb > 1) {
534 MOFEM_LOG_C("MEDWORLD", Sev::noisy,
535 "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
536 vf[1], vf[2], v[0], v[1], v[2]);
537 }
538 if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
539 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
540 "Cannot read MED file older than V2.2");
541 }
542
543 // clear groups
544 group_elem_map.clear();
545
546 med_int num_families = MEDnFamily(fid, meshNames[index].c_str());
547 if (num_families < 0) {
548 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
549 "Could not read MED families");
550 }
551 for (int i = 0; i < num_families; i++) {
552 med_int num_attrib =
553 (vf[0] == 2)
554 ? MEDnFamily23Attribute(fid, meshNames[index].c_str(), i + 1)
555 : 0;
556 med_int num_groups = MEDnFamilyGroup(fid, meshNames[index].c_str(), i + 1);
557 if (num_attrib < 0 || num_groups < 0) {
558 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
559 "Could not read MED groups or attributes");
560 }
561
562 std::vector<med_int> attribId(num_attrib + 1);
563 std::vector<med_int> attrib_val(num_attrib + 1);
564 std::vector<char> attrib_des(MED_COMMENT_SIZE * num_attrib + 1);
565 std::vector<char> group_names(MED_LNAME_SIZE * num_groups + 1);
566 char family_name[MED_NAME_SIZE + 1];
567 med_int family_num;
568
569 if (vf[0] == 2) { // MED2 file
570 if (MEDfamily23Info(fid, meshNames[index].c_str(), i + 1, family_name,
571 &attribId[0], &attrib_val[0], &attrib_des[0],
572 &family_num, &group_names[0]) < 0) {
573 SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
574 "Could not read info for MED2 family %d", i + 1);
575 }
576 } else {
577 if (MEDfamilyInfo(fid, meshNames[index].c_str(), i + 1, family_name,
578 &family_num, &group_names[0]) < 0) {
579 SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
580 "Could not read info for MED3 family %d", i + 1);
581 }
582 }
583
584 // cerr << family_name << " " << family_num << " " << num_groups << endl;
585 for (int g = 0; g != num_groups; g++) {
586 std::string name =
587 std::string(&group_names[MED_LNAME_SIZE * g], MED_LNAME_SIZE - 1);
588 name.resize(NAME_TAG_SIZE - 1);
589 if (family_elem_map.find(family_num) == family_elem_map.end()) {
591 "MEDWORLD", Sev::warning,
592 "Family %d not read, likely type of element is not added "
593 "to moab database. Currently only triangle, quad, tetrahedral and "
594 "hexahedral elements are read to moab database",
595 family_num);
596 } else {
597 group_elem_map[name].merge(family_elem_map.at(family_num));
598 // cerr << string(&group_names[MED_LNAME_SIZE*g]) << endl;
599 }
600 }
601 }
602
603 if (MEDfileClose(fid) < 0) {
604 SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
605 "Unable to close file '%s'", file.c_str());
606 }
607
609}
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
constexpr double g
std::vector< std::string > meshNames
list of meshes in MED file

◆ readFields()

MoFEMErrorCode MoFEM::MedInterface::readFields ( const std::string &  file_name,
const std::string &  field_name,
const bool  load_series = false,
const int  only_step = -1,
int  verb = 1 
)

Read fields

Parameters
file_namefile name
file_indexmesh index
loadSeriesload field into series
onlyStepread only one step
Returns
error code

Definition at line 664 of file MedInterface.cpp.

667 {
668
669 Interface &m_field = cOre;
671 med_idt fid = MEDfileOpen((char *)file_name.c_str(), MED_LECTURE);
672 if (fid < 0) {
673 SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
674 "Unable to open file '%s'", file_name.c_str());
675 }
676
677 med_int num_comp = MEDfieldnComponentByName(fid, field_name.c_str());
678 if (num_comp <= 0) {
679 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
680 "Could not get number of components for MED field");
681 }
682
683 char meshName[MED_NAME_SIZE + 1];
684 char dtUnit[MED_SNAME_SIZE + 1];
685 std::vector<char> compName(num_comp * MED_SNAME_SIZE + 1);
686 std::vector<char> compUnit(num_comp * MED_SNAME_SIZE + 1);
687 med_int numSteps = 0;
688 med_type_champ type;
689 med_bool localMesh;
690 if (MEDfieldInfoByName(fid, field_name.c_str(), meshName, &localMesh, &type,
691 &compName[0], &compUnit[0], dtUnit, &numSteps) < 0) {
692 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
693 "Could not get MED field info");
694 }
695
696 // Get meshset
697 MeshsetsManager *meshsets_manager_ptr;
698 CHKERR m_field.getInterface(meshsets_manager_ptr);
699 const CubitMeshSets *cubit_meshset_ptr;
700 CHKERR meshsets_manager_ptr->getCubitMeshsetPtr(meshName, &cubit_meshset_ptr);
701 EntityHandle meshset = cubit_meshset_ptr->getMeshset();
702
703 // Paraview can only plot field which have 1, 3, or 9 components. If field has
704 // more that 9 comonents, it is stored on MOAB mesh but not viable in
705 // ParaView.
706 int num_comp_msh = (num_comp <= 1) ? 1
707 : (num_comp <= 3) ? 3
708 : (num_comp <= 9) ? 9
709 : num_comp;
710
711 // Create tag to store nodal or cell values read form med file. Note that tag
712 // has prefix MED to avoid collision with other tags.
713 Tag th;
714 std::string tag_name = "MED_" + field_name;
715 {
716 std::vector<double> def_val(num_comp_msh, 0);
717 CHKERR m_field.get_moab().tag_get_handle(
718 tag_name.c_str(), num_comp_msh, MB_TYPE_DOUBLE, th,
719 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val[0]);
720 }
721
722 // Warning! The ordering of the elements in the last two lists is
723 // important: it should match the ordering of the MSH element types
724 // (when elements are saved without tags, this governs the order
725 // with which we implicitly index them in GModel::readMED)
726 const med_entity_type entType[] = {MED_NODE, MED_CELL, MED_NODE_ELEMENT};
727 const med_geometrie_element eleType[] = {
728 MED_NONE, MED_SEG2, MED_TRIA3, MED_QUAD4, MED_TETRA4, MED_HEXA8,
729 MED_PENTA6, MED_PYRA5, MED_SEG3, MED_TRIA6, MED_QUAD9, MED_TETRA10,
730 MED_HEXA27, MED_POINT1, MED_QUAD8, MED_HEXA20, MED_PENTA15, MED_PYRA13};
731 // const int nodesPerEle[] = {0, 2, 3, 4, 4, 8, 6, 5, 3, 6, 9, 10, 27, 1, 8,
732 // 20, 15, 13};
733
734 std::vector<std::pair<int, int>> pairs;
735 for (unsigned int i = 0; i < sizeof(entType) / sizeof(entType[0]); i++) {
736 for (unsigned int j = 0; j < sizeof(eleType) / sizeof(eleType[0]); j++) {
737 if ((!i && !j) || j) {
738 med_int n = numSteps;
739 if (n > 0) {
740 pairs.push_back(std::pair<int, int>(i, j));
741 numSteps = std::max(numSteps, n);
742 }
743 if (!i && !j)
744 break; // MED_NOEUD does not care about eleType
745 }
746 }
747 }
748
749 if (numSteps < 1 || pairs.empty()) {
750 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
751 "Nothing to import from MED file");
752 }
753
754 if (load_series) {
755 SETERRQ(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
756 }
757
758 for (int step = (only_step == -1) ? 0 : only_step; step < numSteps; step++) {
759
760 if (!load_series && only_step != step)
761 continue;
762
763 // cerr << only_step << " " << step << endl;
764
765 // FIXME: in MED3 we might want to loop over all profiles instead
766 // of relying of the default one
767
768 // FIXME: MED3 allows to store multi-step meshes; we should
769
770 for (unsigned int pair = 0; pair < pairs.size(); pair++) {
771
772 // get step info
773 med_entite_maillage ent = entType[pairs[pair].first];
774 med_geometrie_element ele = eleType[pairs[pair].second];
775 med_int numdt, numit, ngauss;
776 med_float dt;
777 if (MEDfieldComputingStepInfo(fid, field_name.c_str(), step + 1, &numdt,
778 &numit, &dt) < 0) {
779 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
780 "Could not read step info");
781 }
782
783 char locName[MED_NAME_SIZE + 1], profileName[MED_NAME_SIZE + 1];
784
785 // get number of values in the field (numVal takes the number of
786 // Gauss points or the number of nodes per element into account,
787 // but not the number of components)
788 med_int profileSize;
789 med_int numVal = MEDfieldnValueWithProfile(
790 fid, field_name.c_str(), numdt, numit, ent, ele, 1,
791 MED_COMPACT_STMODE, profileName, &profileSize, locName, &ngauss);
792 numVal *= ngauss;
793
794 if (numVal <= 0)
795 continue;
796
797 // int mult = 1;
798 // if(ent == MED_NODE_ELEMENT) {
799 // mult = nodesPerEle[pairs[pair].second];
800 //}
801 // else if(ngauss != 1){
802 // mult = ngauss;
803 //}
804
805 // read field data
806 std::vector<double> val(numVal * num_comp);
807 if (MEDfieldValueWithProfileRd(fid, field_name.c_str(), numdt, numit, ent,
808 ele, MED_COMPACT_STMODE, profileName,
809 MED_FULL_INTERLACE, MED_ALL_CONSTITUENT,
810 (unsigned char *)&val[0]) < 0) {
811 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
812 "Could not read field values");
813 }
814
815 // if (verb > 2) {
816 // // FIXME: This not looks ok for me
817 // cerr << ent << " " << ele << endl;
818 // cerr << string(meshName) << " : " << string(profileName) << " : "
819 // << string(locName) << " : " << profileSize << " : " << ngauss
820 // << endl;
821 // }
822
823 switch (ent) {
824 case MED_CELL: {
825 EntityType ent_type = MBMAXTYPE;
826 switch (ele) {
827 case MED_TETRA4:
828 case MED_TETRA10:
829 ent_type = MBTET;
830 break;
831 case MED_HEXA8:
832 ent_type = MBHEX;
833 break;
834 default:
835 MOFEM_LOG_C("MEDWORLD", Sev::warning,
836 "Not yet implemented for this cell %d", ele);
837 }
838 if (ent_type != MBMAXTYPE) {
839 if (ngauss == 1) {
840 Range ents;
841 CHKERR m_field.get_moab().get_entities_by_type(meshset, ent_type,
842 ents, true);
843 double e_vals[num_comp_msh];
844 bzero(e_vals, sizeof(double) * num_comp_msh);
845 std::vector<double>::iterator vit = val.begin();
846 for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
847 for (int ii = 0; ii != num_comp; ii++, vit++) {
848 e_vals[ii] = *vit;
849 }
850 CHKERR m_field.get_moab().tag_set_data(th, &*eit, 1, e_vals);
851 }
852 } else {
853 Range ents;
854 CHKERR m_field.get_moab().get_entities_by_type(meshset, ent_type,
855 ents, true);
856 if (ents.size() * ngauss * num_comp != val.size()) {
857 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
858 "data inconsistency");
859 }
860 // FIXME simply average gauss values, far from perfect need fix
861 double e_vals[num_comp_msh];
862 std::vector<double>::iterator vit = val.begin();
863 for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
864 bzero(e_vals, sizeof(double) * num_comp_msh);
865 for (int gg = 0; gg != ngauss; gg++) {
866 for (int ii = 0; ii != num_comp; ii++, vit++) {
867 e_vals[ii] += *vit / ngauss;
868 }
869 }
870 CHKERR m_field.get_moab().tag_set_data(th, &*eit, 1, e_vals);
871 }
872 }
873 }
874 // SETERRQ1(
875 // m_field.get_comm(),
876 // MOFEM_NOT_IMPLEMENTED,
877 // "Not implemented for more gauss pts ngauss = %d",
878 // ngauss
879 // );
880 } break;
881 case MED_NODE:
882 case MED_NODE_ELEMENT: {
883 EntityType ent_type = MBVERTEX;
884 Range ents;
885 CHKERR m_field.get_moab().get_entities_by_type(meshset, ent_type, ents,
886 true);
887 double e_vals[num_comp_msh];
888 bzero(e_vals, sizeof(double) * num_comp_msh);
889 std::vector<double>::iterator vit = val.begin();
890 for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
891 for (int ii = 0; ii != num_comp; ii++, vit++) {
892 e_vals[ii] = *vit;
893 }
894 CHKERR m_field.get_moab().tag_set_data(th, &*eit, 1, e_vals);
895 }
896 } break;
897 default:
898 MOFEM_LOG_C("MEDWORLD", Sev::inform, "Entity type %d not implemented",
899 ent);
900 }
901 }
902 }
903
905}
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
FTensor::Index< 'n', SPACE_DIM > n
double dt
Definition: heat_method.cpp:26
FTensor::Index< 'j', 3 > j

◆ readMed() [1/2]

MoFEMErrorCode MoFEM::MedInterface::readMed ( const string &  file,
int  verb = 1 
)

read MED file

Parameters
filefiled name
verbverbosity level
Returns
error code

Definition at line 143 of file MedInterface.cpp.

143 {
144 Interface &m_field = cOre;
146
147 med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
148 if (fid < 0) {
149 SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
150 "Unable to open file '%s'", file.c_str());
151 }
152
153 med_int v[3], vf[3];
154 MEDlibraryNumVersion(&v[0], &v[1], &v[2]);
155 MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
156
157 if (verb > 0)
158 MOFEM_LOG_C("MEDWORLD", Sev::inform,
159 "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
160 vf[1], vf[2], v[0], v[1], v[2]);
161
162 if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
163 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
164 "Cannot read MED file older than V2.2");
165 }
166
167 for (int i = 0; i < MEDnMesh(fid); i++) {
168 char mesh_name[MED_NAME_SIZE + 1], mesh_desc[MED_COMMENT_SIZE + 1];
169 bzero(mesh_name, MED_NAME_SIZE);
170 bzero(mesh_desc, MED_COMMENT_SIZE);
171 med_int space_dim;
172 med_mesh_type mesh_type;
173 med_int mesh_dim, n_step;
174 char dt_unit[MED_SNAME_SIZE + 1];
175 char axis_name[3 * MED_SNAME_SIZE + 1], axis_unit[3 * MED_SNAME_SIZE + 1];
176 med_sorting_type sorting_type;
177 med_axis_type axis_type;
178 if (MEDmeshInfo(fid, i + 1, mesh_name, &space_dim, &mesh_dim, &mesh_type,
179 mesh_desc, dt_unit, &sorting_type, &n_step, &axis_type,
180 axis_name, axis_unit) < 0) {
181 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
182 "Unable to read mesh information");
183 }
184 meshNames.push_back(std::string(mesh_name));
185 if (verb > 0) {
186 MOFEM_LOG_C("MEDWORLD", Sev::inform, "Check mesh %s nsteps %d", mesh_name,
187 n_step);
188 }
189 }
190
191 std::map<int, Range> family_elem_map;
192 std::map<string, Range> group_elem_map;
193
194 for (unsigned int ii = 0; ii != meshNames.size(); ii++) {
195 CHKERR readMesh(file, ii, family_elem_map, verb);
196 CHKERR readFamily(file, ii, family_elem_map, group_elem_map, verb);
197 CHKERR makeBlockSets(group_elem_map, verb);
198 }
199
200 if (MEDfileClose(fid) < 0) {
201 SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
202 "Unable to close file '%s'", file.c_str());
203 }
204
206}
MoFEMErrorCode readFamily(const string &file, const int index, const std::map< int, Range > &family_elem_mapint, std::map< string, Range > &group_elem_map, int verb=1)
read family and groups
MoFEMErrorCode makeBlockSets(const std::map< string, Range > &group_elem_map, int verb=1)
make from groups meshsets
MoFEMErrorCode readMesh(const string &file, const int index, std::map< int, Range > &family_elem_map, int verb=1)
read mesh from MED file

◆ readMed() [2/2]

MoFEMErrorCode MoFEM::MedInterface::readMed ( int  verb = 1)

read MED file

File name is form command line

Parameters
verbverbosity level
Returns
error code

Definition at line 648 of file MedInterface.cpp.

648 {
650 if (medFileName.empty()) {
652 }
655}
MoFEMErrorCode readMed(const string &file, int verb=1)
read MED file

◆ readMesh()

MoFEMErrorCode MoFEM::MedInterface::readMesh ( const string &  file,
const int  index,
std::map< int, Range > &  family_elem_map,
int  verb = 1 
)
private

read mesh from MED file

Parameters
filefile name
indexindex of mesh
family_elem_mapmap of families and elements
verbverbosity level
Returns
error code

Definition at line 247 of file MedInterface.cpp.

249 {
250
251 Interface &m_field = cOre;
253
254 med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
255 if (fid < 0) {
256 SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
257 "Unable to open file '%s'", file.c_str());
258 }
259 med_int v[3], vf[3];
260 MEDlibraryNumVersion(&v[0], &v[1], &v[2]);
261 MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
262 if (verb > 1)
263 MOFEM_LOG_C("MEDWORLD", Sev::noisy,
264 "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
265 vf[1], vf[2], v[0], v[1], v[2]);
266
267 if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
268 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
269 "Cannot read MED file older than V2.2");
270 }
271
272 char mesh_name[MED_NAME_SIZE + 1], mesh_desc[MED_COMMENT_SIZE + 1];
273 bzero(mesh_name, MED_NAME_SIZE + 1);
274 bzero(mesh_desc, MED_COMMENT_SIZE + 1);
275 med_int space_dim;
276 med_mesh_type mesh_type;
277 med_int mesh_dim, n_step;
278 char dt_unit[MED_SNAME_SIZE + 1];
279 char axis_name[3 * MED_SNAME_SIZE + 1], axis_unit[3 * MED_SNAME_SIZE + 1];
280 med_sorting_type sorting_type;
281 med_axis_type axis_type;
282 if (MEDmeshInfo(fid, index + 1, mesh_name, &space_dim, &mesh_dim, &mesh_type,
283 mesh_desc, dt_unit, &sorting_type, &n_step, &axis_type,
284 axis_name, axis_unit) < 0) {
285 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
286 "Unable to read mesh information");
287 }
288 if (verb > 0)
289 MOFEM_LOG_C("MEDWORLD", Sev::inform, "Reading mesh %s nsteps %d", mesh_name,
290 n_step);
291
292 switch (axis_type) {
293 case MED_CARTESIAN:
294 break;
295 case MED_SPHERICAL:
296 case MED_CYLINDRICAL:
297 default:
298 SETERRQ(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED,
299 "Curvilinear coordinate system implemented");
300 }
301 if (space_dim < 2) {
302 SETERRQ1(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED,
303 "Not implemented for space dim %d", space_dim);
304 }
305
306 EntityHandle mesh_meshset;
307 {
308 MeshsetsManager *meshsets_manager_ptr;
309 CHKERR m_field.getInterface(meshsets_manager_ptr);
310 int max_id = 0;
312 max_id = (max_id < cit->getMeshsetId()) ? cit->getMeshsetId() : max_id;
313 }
314 max_id++;
315 CHKERR meshsets_manager_ptr->addMeshset(BLOCKSET, max_id,
316 std::string(mesh_name));
318 Composite_Cubit_msId_And_MeshsetType_mi_tag>::type::iterator cit;
319 cit =
320 meshsets_manager_ptr->getMeshsetsMultindex()
321 .get<Composite_Cubit_msId_And_MeshsetType_mi_tag>()
322 .find(boost::make_tuple(max_id, CubitBCType(BLOCKSET).to_ulong()));
323 max_id++;
324 mesh_meshset = cit->getMeshset();
325 meshMeshsets.push_back(mesh_meshset);
326 }
327
328 med_bool change_of_coord, geo_transform;
329 med_int num_nodes = MEDmeshnEntity(
330 fid, mesh_name, MED_NO_DT, MED_NO_IT, MED_NODE, MED_NO_GEOTYPE,
331 MED_COORDINATE, MED_NO_CMODE, &change_of_coord, &geo_transform);
332 if (num_nodes < 0) {
333 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
334 "Could not read number of MED nodes");
335 }
336 if (num_nodes == 0) {
337 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
338 "No nodes in MED mesh");
339 }
340 if (verb > 0)
341 MOFEM_LOG_C("MEDWORLD", Sev::inform, "Read number of nodes %d", num_nodes);
342
343 std::vector<med_float> coord_med(space_dim * num_nodes);
344 if (MEDmeshNodeCoordinateRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
345 MED_NO_INTERLACE, &coord_med[0]) < 0) {
346 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
347 "Could not read MED node coordinates");
348 }
349
350 // Add vertices to moab
351 ReadUtilIface *iface;
352 vector<double *> arrays_coord;
353 EntityHandle startv;
354 CHKERR m_field.get_moab().query_interface(iface);
355 CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays_coord);
356 Range verts(startv, startv + num_nodes - 1);
357 std::copy(&coord_med[0 * num_nodes], &coord_med[1 * num_nodes],
358 arrays_coord[0]);
359 std::copy(&coord_med[1 * num_nodes], &coord_med[2 * num_nodes],
360 arrays_coord[1]);
361 if (space_dim == 2) {
362 std::fill(arrays_coord[2], &arrays_coord[2][num_nodes], 0.);
363 } else {
364 std::copy(&coord_med[2 * num_nodes], &coord_med[3 * num_nodes],
365 arrays_coord[2]);
366 }
367 CHKERR m_field.get_moab().add_entities(mesh_meshset, verts);
368 family_elem_map.clear();
369
370 // get family for vertices
371 {
372 std::vector<med_int> nodes_tags(num_nodes, 0);
373 if (MEDmeshEntityFamilyNumberRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
374 MED_NODE, MED_NO_GEOTYPE,
375 &nodes_tags[0]) < 0) {
376 nodes_tags.clear();
377 // SETERRQ(
378 // m_field.get_comm(),
379 // MOFEM_OPERATION_UNSUCCESSFUL,
380 // "No family number for elements: using 0 as default family number"
381 // );
382 }
383 for (int i = 0; i < num_nodes; i++) {
384 // cerr << verts[i] << " " /*<< ele_tags[j] << " "*/ << nodes_tags[i] <<
385 // endl;
386 family_elem_map[nodes_tags.empty() ? i : nodes_tags[i]].insert(verts[i]);
387 }
388 }
389
390 // read elements (loop over all possible MSH element types)
391 for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
392
393 auto types = moab2med_element_type(ent_type);
394
395 for (auto type : types) {
396
397 // get number of cells
398 med_bool change_of_coord;
399 med_bool geo_transform;
400 med_int num_ele = MEDmeshnEntity(
401 fid, mesh_name, MED_NO_DT, MED_NO_IT, MED_CELL, type,
402 MED_CONNECTIVITY, MED_NODAL, &change_of_coord, &geo_transform);
403
404 // get connectivity
405 if (num_ele <= 0)
406 continue;
407
408 int num_nod_per_ele = type % 100;
409 if (verb > 0)
410 MOFEM_LOG("MEDWORLD", Sev::inform)
411 << "Reading elements " << num_ele << " of type "
412 << moab::CN::EntityTypeName(ent_type) << " number of nodes "
413 << num_nod_per_ele;
414
415 std::vector<med_int> conn_med(num_ele * num_nod_per_ele);
416 if (MEDmeshElementConnectivityRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
417 MED_CELL, type, MED_NODAL,
418 MED_FULL_INTERLACE, &conn_med[0]) < 0) {
419 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
420 "Could not read MED elements");
421 }
422
423 // cerr << "type " << ent_type << " ";
424 // cerr << "num_ele " << num_ele << " " << num_nod_per_ele << endl;;
425
426 Range ents;
427
428 if (ent_type != MBVERTEX) {
429 EntityHandle *conn_moab;
430 EntityHandle starte;
431 CHKERR iface->get_element_connect(num_ele, num_nod_per_ele, ent_type, 0,
432 starte, conn_moab);
433 switch (ent_type) {
434 // FIXME: Some connectivity could not work, need to run and test
435 case MBTET: {
436 int ii = 0;
437 for (int ee = 0; ee != num_ele; ee++) {
438 EntityHandle n[4];
439 for (int nn = 0; nn != num_nod_per_ele; nn++) {
440 // conn_moab[ii] = verts[conn_med[ii]-1];
441 n[nn] = verts[conn_med[ii + nn] - 1];
442 }
443
444 std::array<int, 10> nodes_map{
445
446 1, 0, 2, 3,
447
448 4, 6, 5, 8,
449 7, 9
450
451 };
452
453 for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
454 conn_moab[ii] = n[nodes_map[nn]];
455 }
456 }
457 } break;
458 default: {
459 int ii = 0;
460 for (int ee = 0; ee != num_ele; ee++) {
461 for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
462 // cerr << conn_med[ii] << " ";
463 conn_moab[ii] = verts[conn_med[ii] - 1];
464 }
465 // cerr << endl;
466 }
467 }
468 }
469 CHKERR iface->update_adjacencies(starte, num_ele, num_nod_per_ele,
470 conn_moab);
471 ents = Range(starte, starte + num_ele - 1);
472 } else {
473 // This is special case when in med vertices are defined as elements
474 int ii = 0;
475 std::vector<EntityHandle> conn_moab(num_ele * num_nod_per_ele);
476 for (int ee = 0; ee != num_ele; ++ee)
477 for (int nn = 0; nn != num_nod_per_ele; ++nn, ++ii)
478 conn_moab[ii] = verts[conn_med[ii] - 1];
479 ents.insert_list(conn_moab.begin(), conn_moab.end());
480 }
481
482 // Add elements to family meshset
483 CHKERR m_field.get_moab().add_entities(mesh_meshset, ents);
484
485 // get family for cells
486 {
487 std::vector<med_int> fam(num_ele, 0);
488 if (MEDmeshEntityFamilyNumberRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
489 MED_CELL, type, &fam[0]) < 0) {
490 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
491 "No family number for elements: using 0 as default family "
492 "number");
493 }
494 // std::vector<med_int> ele_tags(num_ele);
495 // if(MEDmeshEntityNumberRd(
496 // fid, mesh_name,MED_NO_DT,MED_NO_IT,MED_CELL,type,&ele_tags[0]) < 0
497 // ) {
498 // ele_tags.clear();
499 // }
500 for (int j = 0; j < num_ele; j++) {
501 // cerr << ents[j] << " " /*<< ele_tags[j] << " "*/ << fam[j] << endl;
502 family_elem_map[fam[j]].insert(ents[j]);
503 }
504 }
505 }
506 }
507
508 if (MEDfileClose(fid) < 0) {
509 SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
510 "Unable to close file '%s'", file.c_str());
511 }
512
514}
static std::vector< med_geometrie_element > moab2med_element_type(const EntityType type)
std::vector< EntityHandle > meshMeshsets
meshset for each mesh

◆ writeMed()

MoFEMErrorCode MoFEM::MedInterface::writeMed ( const string &  file,
int  verb = 1 
)

write MED file

Parameters
filefile name
verbverbosity level
Returns
error code

Definition at line 657 of file MedInterface.cpp.

657 {
658 Interface &m_field = cOre;
660 SETERRQ(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
662}

Member Data Documentation

◆ cOre

MoFEM::Core& MoFEM::MedInterface::cOre
private

core database

Definition at line 115 of file MedInterface.hpp.

◆ fieldNames

std::map<std::string, FieldData> MoFEM::MedInterface::fieldNames

Definition at line 111 of file MedInterface.hpp.

◆ flgFile

PetscBool MoFEM::MedInterface::flgFile
private

true if file name given in command line

Definition at line 117 of file MedInterface.hpp.

◆ medFileName

std::string MoFEM::MedInterface::medFileName

MED file name.

Definition at line 112 of file MedInterface.hpp.

◆ meshMeshsets

std::vector<EntityHandle> MoFEM::MedInterface::meshMeshsets

meshset for each mesh

Definition at line 110 of file MedInterface.hpp.

◆ meshNames

std::vector<std::string> MoFEM::MedInterface::meshNames

list of meshes in MED file

Definition at line 109 of file MedInterface.hpp.


The documentation for this struct was generated from the following files: