709 <<
"WRITE_MED IS EXPERIMENTAL, MAY CONTAIN BUGS, ALWAYS CHECK THE OUTPUT "
713 MEDfileVersionOpen((
char *)file.c_str(), MED_ACC_CREAT, MED_MAJOR_NUM,
714 MED_MINOR_NUM, MED_RELEASE_NUM);
718 "Cannot create MED file");
723 CHKERR MEDfichDesEcr(fid, (
char *)
"MED file generated by MoFEM");
726 char dtUnit[MED_SNAME_SIZE + 1] =
"";
727 char axisName[3 * MED_SNAME_SIZE + 1] =
"";
728 char axisUnit[3 * MED_SNAME_SIZE + 1] =
"";
730 PetscBool is_cubit_meshset = PETSC_TRUE;
731 int med_mesh_name_id = 0;
734 char mesh_name_char[255];
735 std::string mesh_name =
"Mesh";
738 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"MED mesh options",
"");
739 CHKERR PetscOptionsString(
"-med_mesh_name",
"get med mesh name",
"",
740 mesh_name.c_str(), mesh_name_char, 255, PETSC_NULL);
741 ierr = PetscOptionsEnd();
744 mesh_name = mesh_name_char;
747 for (
auto &
m : *meshsets_ptr) {
748 if (
m->getName() == mesh_name) {
749 med_mesh_name_id =
m->getMeshsetId();
750 is_cubit_meshset = PETSC_FALSE;
757 for (
auto &
m : *meshsets_ptr) {
758 if (
m->getMeshsetId() == med_mesh_name_id)
759 mesh_name =
m->getName();
760 max_id = (max_id <
m->getMeshsetId()) ?
m->getMeshsetId() : max_id;
764 CHKERR MEDmeshCr(fid, mesh_name.c_str(), 3, 3, MED_UNSTRUCTURED_MESH,
765 "Mesh created", dtUnit, MED_SORT_DTIT, MED_CARTESIAN,
769 med_int family_id = 0;
770 std::map<std::vector<std::string>, std::tuple<med_int, std::vector<int>>>
772 std::map<EntityHandle, med_int> entityHandle_family_map;
775 shared_meshsets_map[std::vector<string>()] =
776 std::make_tuple(family_id, std::vector<int>());
779 auto get_set_name = [&](
const CubitMeshSets *iit) {
780 std::string bc_type_name;
781 if (iit->getBcTypeULong() &
BLOCKSET) {
783 bc_type_name = iit->getName();
784 if (bc_type_name ==
"NoNameSet") {
785 bc_type_name =
"BLOCKSET_NoNameSet_";
786 bc_type_name += std::to_string(iit->getMeshsetId());
789 }
else if (iit->getBcTypeULong() &
SIDESET ||
790 iit->getBcTypeULong() &
NODESET) {
797 if ((iit->getBcType() & jj_bc_type).any()) {
804 bc_type_name += std::to_string(iit->getMeshsetId());
806 bc_type_name =
"UnknownSet";
816 for (
auto &entity : *write_range_ptr) {
818 std::vector<int> shared_meshsets;
819 std::vector<string> shared_names;
822 auto add_shared_meshset = [&](
auto &other_meshset) {
823 std::string other_name = get_set_name(other_meshset);
824 if (std::find(shared_names.begin(), shared_names.end(), other_name) ==
825 shared_names.end()) {
826 shared_meshsets.push_back(other_meshset->getMeshsetId());
827 shared_names.push_back(other_name);
832 for (
auto &other_meshset : *meshsets_ptr) {
834 if (med_mesh_name_id == other_meshset->getMeshsetId())
837 Range other_entities;
839 CHKERR moab.get_entities_by_handle(other_set, other_entities);
842 EntityType ent_type = moab.type_from_handle(entity);
844 bool is_in_meshset = moab.contains_entities(other_set, &entity, 1);
846 if (ent_type == MBVERTEX) {
849 bool is_in_higher_dim =
false;
850 Range entities_in_higher_dim;
851 for (
int i = 1;
i < 4;
i++) {
852 CHKERR moab.get_entities_by_dimension(other_set,
i,
853 entities_in_higher_dim);
854 if (!entities_in_higher_dim.empty()) {
855 is_in_higher_dim =
true;
859 if (is_in_higher_dim) {
863 add_shared_meshset(other_meshset);
866 add_shared_meshset(other_meshset);
872 auto it = shared_meshsets_map.find(shared_names);
873 if (it == shared_meshsets_map.end()) {
877 it = shared_meshsets_map
879 {shared_names, std::make_tuple(family_id, shared_meshsets)})
883 entityHandle_family_map[entity] = std::get<0>(it->second);
887 for (
const auto &it : shared_meshsets_map) {
889 std::string family_name =
"F_";
890 const auto &[family_id, shared_meshsets] = it.second;
891 family_name += std::to_string(family_id);
893 const auto &shared_meshset_names = it.first;
894 std::string group_name;
895 for (
const auto &name : shared_meshset_names) {
897 std::string meshset_name = name;
898 meshset_name.resize(MED_LNAME_SIZE,
' ');
899 group_name += meshset_name;
903 CHKERR MEDfamilyCr(fid, mesh_name.c_str(), family_name.c_str(), family_id,
904 shared_meshset_names.size(), group_name.c_str());
907 <<
"Creating family " << family_name <<
" with id " << family_id
908 <<
" and " << shared_meshset_names.size() <<
" groups " << std::endl;
913 moab.get_entities_by_type(0, MBVERTEX, verts);
915 std::vector<med_float> coord_med(3 * verts.size());
916 std::vector<med_int> fam;
917 std::vector<med_int> tags;
920 for (Range::iterator it = verts.begin(); it != verts.end(); ++it) {
922 moab.get_coords(&(*it), 1, coords);
923 coord_med[3 * (it - verts.begin())] = coords[0];
924 coord_med[3 * (it - verts.begin()) + 1] = coords[1];
925 coord_med[3 * (it - verts.begin()) + 2] = coords[2];
926 fam.push_back(entityHandle_family_map[*it]);
930 CHKERR MEDmeshNodeWr(fid, mesh_name.c_str(), MED_NO_DT, MED_NO_IT,
931 MED_UNDEF_DT, MED_FULL_INTERLACE, verts.size(),
932 &coord_med[0], MED_FALSE,
"", MED_TRUE, &tags[0],
936 for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
939 moab.get_entities_by_type(0, ent_type, entities,
true);
941 ents_to_write = intersect(entities, *write_range_ptr);
943 if (ents_to_write.empty())
947 std::vector<med_int> tag_number;
948 std::vector<med_int> family_number;
949 std::vector<med_int> connectivity;
952 for (
auto &entity : ents_to_write) {
953 if (ent_type != MBVERTEX) {
955 family_number.push_back(entityHandle_family_map[entity]);
957 tag_number.push_back(entity);
959 std::vector<EntityHandle> conn;
960 moab.get_connectivity(&entity, 1, conn);
961 for (
auto &
c : conn) {
962 connectivity.push_back(
c);
969 auto get_med_element_type = [](EntityType ent_type) {
970 med_geometrie_element
type;
997 med_geometrie_element med_type = get_med_element_type(ent_type);
998 if (ent_type == MBENTITYSET) {
1001 CHKERR MEDmeshElementWr(fid, mesh_name.c_str(), MED_NO_DT, MED_NO_IT, 0.,
1002 MED_CELL, med_type, MED_NODAL, MED_FULL_INTERLACE,
1003 family_number.size(), &connectivity[0], MED_FALSE,
1004 nullptr, MED_TRUE, &tag_number[0], MED_TRUE,
1007 MOFEM_LOG_C(
"MEDWORLD", Sev::inform,
"Writing %i elements of type %i (%s) ",
1008 family_number.size(), med_type,
1009 moab::CN::EntityTypeName(ent_type));