v0.8.23
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 (const MOFEMuuid &uuid, 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
template<class IFACE >
MoFEMErrorCode registerInterface (const MOFEMuuid &uuid, bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE , bool VERIFY = false>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, IFACE *&iface) const
 Get interface by uuid and return reference to pointer of 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 ()
 
virtual MoFEMErrorCode getLibVersion (Version &version) const
 Get library version. More...
 
virtual const MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version) const
 Get database major version. More...
 
virtual MoFEMErrorCode getInterfaceVersion (Version &version) const
 Get database major version. More...
 
template<>
MoFEMErrorCode getInterface (const MOFEMuuid &uuid, UnknownInterface *&iface) const
 

Public Attributes

std::vector< std::string > meshNames
 list of meshes in MED file More...
 
std::vector< EntityHandlemeshMeshsets
 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

- Protected Member Functions inherited from MoFEM::UnknownInterface
boost::typeindex::type_index getClassIdx (const MOFEMuuid &uid) const
 Get type name for interface Id. More...
 
MOFEMuuid getUId (const boost::typeindex::type_index &class_idx) const
 Get interface Id for class name. More...
 

Detailed Description

Interface for load MED files.

Definition at line 36 of file MedInterface.hpp.

Constructor & Destructor Documentation

◆ MedInterface()

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

Definition at line 62 of file MedInterface.cpp.

63  : cOre(const_cast<Core &>(core)), flgFile(PETSC_FALSE) {}
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 65 of file MedInterface.cpp.

65  {
66  Interface &m_field = cOre;
67  char mesh_file_name[255];
69  ierr = PetscOptionsBegin(m_field.get_comm(), "", "MED Interface", "none");
70  CHKERRG(ierr);
71  ierr = PetscOptionsString("-med_file", "med file name", "", "mesh.med",
72  mesh_file_name, 255, &flgFile);
73  CHKERRG(ierr);
74  ierr = PetscOptionsEnd();
75  CHKERRG(ierr);
76  if (flgFile == PETSC_TRUE) {
77  medFileName = std::string(mesh_file_name);
78  } else {
79  medFileName = std::string("mesh.med");
80  }
82 }
MoFEM::Core & cOre
core database
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:500
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:543
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507
char mesh_file_name[255]
PetscBool flgFile
true if file name given in command line
std::string medFileName
MED file name.

◆ getFlgFile()

PetscBool MoFEM::MedInterface::getFlgFile ( ) const

Check if file name is given in command line.

Definition at line 52 of file MedInterface.hpp.

52 { return flgFile; }
PetscBool flgFile
true if file name given in command line

◆ 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 590 of file MedInterface.cpp.

591  {
592 
593  Interface &m_field = cOre;
595  MeshsetsManager *meshsets_manager_ptr;
596  CHKERR m_field.getInterface(meshsets_manager_ptr);
597 
598  int max_id = 0;
600  max_id = (max_id < cit->getMeshsetId()) ? cit->getMeshsetId() : max_id;
601  }
602  max_id++;
603 
604  // cerr << group_elem_map.size() << endl;
605  for (std::map<string, Range>::const_iterator git = group_elem_map.begin();
606  git != group_elem_map.end(); git++) {
607  // cerr << "AAA\n";
608  CHKERR meshsets_manager_ptr->addMeshset(BLOCKSET, max_id, git->first);
609  CubitMeshSet_multiIndex::index<
610  Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator cit;
611  cit =
612  meshsets_manager_ptr->getMeshsetsMultindex()
613  .get<Composite_Cubit_msId_And_MeshSetType_mi_tag>()
614  .find(boost::make_tuple(max_id, CubitBCType(BLOCKSET).to_ulong()));
615  EntityHandle meshsets = cit->getMeshset();
616  if (!git->second.empty()) {
617  CHKERR m_field.get_moab().add_entities(meshsets, git->second);
618  }
619  max_id++;
620  // cerr << git->second << endl;
621  }
622 
623  // if(verb>0) {
625  PetscPrintf(
626  m_field.get_comm(), "%s\n",
627  static_cast<std::ostringstream &>(std::ostringstream().seekp(0) << *cit)
628  .str()
629  .c_str());
630  }
631  // }
632 
634 }
MoFEM::Core & cOre
core database
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
#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.
#define CHKERR
Inline error check.
Definition: definitions.h:595
std::bitset< 32 > CubitBCType
Definition: Types.hpp:63
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ 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 84 of file MedInterface.cpp.

84  {
85  Interface &m_field = cOre;
87  med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
88  if (fid < 0) {
89  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
90  "Unable to open file '%s'", file.c_str());
91  }
92  med_int num_fields = MEDnField(fid);
93  for (int index = 0; index < num_fields; index++) {
94 
95  med_int num_comp = MEDfieldnComponent(fid, index + 1);
96  if (num_comp <= 0) {
97  SETERRQ(m_field.get_comm(), MOFEM_IMPOSIBLE_CASE,
98  "Could not get number of components for MED field");
99  }
100 
101  char name[MED_NAME_SIZE + 1], mesh_name[MED_NAME_SIZE + 1];
102  char dt_unit[MED_SNAME_SIZE + 1];
103  std::vector<char> comp_name(num_comp * MED_SNAME_SIZE + 1);
104  std::vector<char> comp_unit(num_comp * MED_SNAME_SIZE + 1);
105  med_int num_steps = 0;
106  med_type_champ type;
107  med_bool local_mesh;
108  if (MEDfieldInfo(fid, index + 1, name, mesh_name, &local_mesh, &type,
109  &comp_name[0], &comp_unit[0], dt_unit, &num_steps) < 0) {
110  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
111  "Could not get MED field info");
112  }
113 
114  std::string field_name = std::string(name);
115  fieldNames[field_name] = FieldData();
116  fieldNames[field_name].fieldName = field_name;
117  fieldNames[field_name].meshName = std::string(mesh_name);
118  fieldNames[field_name].localMesh = (local_mesh == MED_TRUE);
119  for (int ff = 0; ff != num_comp; ff++) {
120  fieldNames[field_name].componentNames.push_back(
121  std::string(&comp_name[ff * MED_SNAME_SIZE], MED_SNAME_SIZE));
122  fieldNames[field_name].componentUnits.push_back(
123  std::string(&comp_unit[ff * MED_SNAME_SIZE], MED_SNAME_SIZE));
124  }
125  fieldNames[field_name].dtUnit = std::string(&dt_unit[0]);
126  fieldNames[field_name].ncSteps = num_steps;
127 
128  if (verb > 0) {
129  std::ostringstream ss;
130  ss << fieldNames[name] << std::endl;
131  CHKERR PetscPrintf(m_field.get_comm(), ss.str().c_str());
132  }
133  }
134  if (MEDfileClose(fid) < 0) {
135  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
136  "Unable to close file '%s'", file.c_str());
137  }
139 }
std::map< std::string, FieldData > fieldNames
MoFEM::Core & cOre
core database
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
#define CHKERR
Inline error check.
Definition: definitions.h:595
double FieldData
Field data type.
Definition: Types.hpp:36
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ 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 141 of file MedInterface.cpp.

141  {
143  if (medFileName.empty()) {
145  }
148 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
MoFEMErrorCode getFileNameFromCommandLine(int verb=1)
Get MED file name from command line.
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406
MoFEMErrorCode medGetFieldNames(const string &file, int verb=1)
Get field names in MED file.
std::string medFileName
MED file name.

◆ query_interface()

MoFEMErrorCode MoFEM::MedInterface::query_interface ( const MOFEMuuid uuid,
UnknownInterface **  iface 
) const
virtual

Implements MoFEM::UnknownInterface.

Definition at line 50 of file MedInterface.cpp.

51  {
53  *iface = NULL;
54  if (uuid == IDD_MOFEMMedInterface) {
55  *iface = const_cast<MedInterface *>(this);
57  }
58  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
60 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:500
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507
static const MOFEMuuid IDD_MOFEMMedInterface

◆ 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 493 of file MedInterface.cpp.

495  {
496  //
497  Interface &m_field = cOre;
499 
500  med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
501  if (fid < 0) {
502  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
503  "Unable to open file '%s'", file.c_str());
504  }
505  med_int v[3], vf[3];
506  MEDlibraryNumVersion(&v[0], &v[1], &v[2]);
507  MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
508  // if(verb>1) {
509  // PetscPrintf(
510  // m_field.get_comm(),
511  // "Reading MED file V%d.%d.%d using MED library V%d.%d.%d\n",
512  // vf[0], vf[1], vf[2], v[0], v[1], v[2]
513  // );
514  // }
515  if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
516  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
517  "Cannot read MED file older than V2.2");
518  }
519 
520  // clear groups
521  group_elem_map.clear();
522 
523  med_int num_families = MEDnFamily(fid, meshNames[index].c_str());
524  if (num_families < 0) {
525  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
526  "Could not read MED families");
527  }
528  for (int i = 0; i < num_families; i++) {
529  med_int num_attrib =
530  (vf[0] == 2)
531  ? MEDnFamily23Attribute(fid, meshNames[index].c_str(), i + 1)
532  : 0;
533  med_int num_groups = MEDnFamilyGroup(fid, meshNames[index].c_str(), i + 1);
534  if (num_attrib < 0 || num_groups < 0) {
535  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
536  "Could not read MED groups or attributes");
537  }
538 
539  std::vector<med_int> attribId(num_attrib + 1);
540  std::vector<med_int> attrib_val(num_attrib + 1);
541  std::vector<char> attrib_des(MED_COMMENT_SIZE * num_attrib + 1);
542  std::vector<char> group_names(MED_LNAME_SIZE * num_groups + 1);
543  char family_name[MED_NAME_SIZE + 1];
544  med_int family_num;
545 
546  if (vf[0] == 2) { // MED2 file
547  if (MEDfamily23Info(fid, meshNames[index].c_str(), i + 1, family_name,
548  &attribId[0], &attrib_val[0], &attrib_des[0],
549  &family_num, &group_names[0]) < 0) {
550  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
551  "Could not read info for MED2 family %d", i + 1);
552  }
553  } else {
554  if (MEDfamilyInfo(fid, meshNames[index].c_str(), i + 1, family_name,
555  &family_num, &group_names[0]) < 0) {
556  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
557  "Could not read info for MED3 family %d", i + 1);
558  }
559  }
560 
561  // cerr << family_name << " " << family_num << " " << num_groups << endl;
562  for (int g = 0; g != num_groups; g++) {
563  std::string name =
564  std::string(&group_names[MED_LNAME_SIZE * g], MED_LNAME_SIZE - 1);
565  name.resize(NAME_TAG_SIZE - 1);
566  if (family_elem_map.find(family_num) == family_elem_map.end()) {
567  PetscPrintf(
568  PETSC_COMM_SELF,
569  "Warring: \n Family %d not read, likely type of element is not "
570  "added "
571  "to moab database. Currently only triangle, quad, tetrahedral and "
572  "hexahedral elements are read to moab database\n",
573  family_num);
574  } else {
575  group_elem_map[name].merge(family_elem_map.at(family_num));
576  // cerr << string(&group_names[MED_LNAME_SIZE*g]) << endl;
577  }
578  }
579  }
580 
581  if (MEDfileClose(fid) < 0) {
582  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
583  "Unable to close file '%s'", file.c_str());
584  }
585 
587 }
MoFEM::Core & cOre
core database
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
std::vector< std::string > meshNames
list of meshes in MED file
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ 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 652 of file MedInterface.cpp.

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

◆ 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 150 of file MedInterface.cpp.

150  {
151  Interface &m_field = cOre;
153 
154  med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
155  if (fid < 0) {
156  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
157  "Unable to open file '%s'", file.c_str());
158  }
159 
160  med_int v[3], vf[3];
161  MEDlibraryNumVersion(&v[0], &v[1], &v[2]);
162  MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
163 
164  if (verb > 0) {
165  PetscPrintf(m_field.get_comm(),
166  "Reading MED file V%d.%d.%d using MED library V%d.%d.%d\n",
167  vf[0], vf[1], vf[2], v[0], v[1], v[2]);
168  }
169  if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
170  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
171  "Cannot read MED file older than V2.2");
172  }
173 
174  for (int i = 0; i < MEDnMesh(fid); i++) {
175  char mesh_name[MED_NAME_SIZE + 1], mesh_desc[MED_COMMENT_SIZE + 1];
176  bzero(mesh_name, MED_NAME_SIZE);
177  bzero(mesh_desc, MED_COMMENT_SIZE);
178  med_int space_dim;
179  med_mesh_type mesh_type;
180  med_int mesh_dim, n_step;
181  char dt_unit[MED_SNAME_SIZE + 1];
182  char axis_name[3 * MED_SNAME_SIZE + 1], axis_unit[3 * MED_SNAME_SIZE + 1];
183  med_sorting_type sorting_type;
184  med_axis_type axis_type;
185  if (MEDmeshInfo(fid, i + 1, mesh_name, &space_dim, &mesh_dim, &mesh_type,
186  mesh_desc, dt_unit, &sorting_type, &n_step, &axis_type,
187  axis_name, axis_unit) < 0) {
188  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
189  "Unable to read mesh information");
190  }
191  meshNames.push_back(std::string(mesh_name));
192  if (verb > 0) {
193  PetscPrintf(m_field.get_comm(), "Check mesh %s nsteps %d\n", mesh_name,
194  n_step);
195  }
196  }
197 
198  std::map<int, Range> family_elem_map;
199  std::map<string, Range> group_elem_map;
200 
201  for (unsigned int ii = 0; ii != meshNames.size(); ii++) {
202  CHKERR readMesh(file, ii, family_elem_map, verb);
203  CHKERR readFamily(file, ii, family_elem_map, group_elem_map, verb);
204  CHKERR makeBlockSets(group_elem_map, verb);
205  }
206 
207  if (MEDfileClose(fid) < 0) {
208  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
209  "Unable to close file '%s'", file.c_str());
210  }
211 
213 }
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
MoFEM::Core & cOre
core database
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
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
std::vector< std::string > meshNames
list of meshes in MED file
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406

◆ 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 636 of file MedInterface.cpp.

636  {
638  if (medFileName.empty()) {
640  }
641  CHKERR readMed(medFileName, verb);
643 }
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
MoFEMErrorCode getFileNameFromCommandLine(int verb=1)
Get MED file name from command line.
MoFEMErrorCode readMed(const string &file, int verb=1)
read MED file
#define CHKERR
Inline error check.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406
std::string medFileName
MED file name.

◆ 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 238 of file MedInterface.cpp.

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

◆ 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 645 of file MedInterface.cpp.

645  {
646  Interface &m_field = cOre;
648  SETERRQ(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
650 }
MoFEM::Core & cOre
core database
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:500
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507

Member Data Documentation

◆ cOre

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

core database

Definition at line 128 of file MedInterface.hpp.

◆ fieldNames

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

Definition at line 124 of file MedInterface.hpp.

◆ flgFile

PetscBool MoFEM::MedInterface::flgFile
private

true if file name given in command line

Definition at line 130 of file MedInterface.hpp.

◆ medFileName

std::string MoFEM::MedInterface::medFileName

MED file name.

Definition at line 125 of file MedInterface.hpp.

◆ meshMeshsets

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

meshset for each mesh

Definition at line 123 of file MedInterface.hpp.

◆ meshNames

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

list of meshes in MED file

Definition at line 122 of file MedInterface.hpp.


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