711 <<
"WRITE_MED IS EXPERIMENTAL, MAY CONTAIN BUGS, ALWAYS CHECK THE OUTPUT "
715 MEDfileVersionOpen((
char *)file.c_str(), MED_ACC_CREAT, MED_MAJOR_NUM,
716 MED_MINOR_NUM, MED_RELEASE_NUM);
720 "Cannot create MED file");
725 CHKERR MEDfichDesEcr(fid, (
char *)
"MED file generated by MoFEM");
728 char dtUnit[MED_SNAME_SIZE + 1] =
"";
729 char axisName[3 * MED_SNAME_SIZE + 1] =
"";
730 char axisUnit[3 * MED_SNAME_SIZE + 1] =
"";
732 PetscBool is_cubit_meshset = PETSC_TRUE;
733 int med_mesh_name_id = 0;
736 char mesh_name_char[255];
737 std::string mesh_name =
"Mesh";
740 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"MED mesh options",
"");
741 CHKERR PetscOptionsString(
"-med_mesh_name",
"get med mesh name",
"",
742 mesh_name.c_str(), mesh_name_char, 255, PETSC_NULL);
743 ierr = PetscOptionsEnd();
746 mesh_name = mesh_name_char;
749 for (
auto &
m : *meshsets_ptr) {
750 if (
m->getName() == mesh_name) {
751 med_mesh_name_id =
m->getMeshsetId();
752 is_cubit_meshset = PETSC_FALSE;
759 for (
auto &
m : *meshsets_ptr) {
760 if (
m->getMeshsetId() == med_mesh_name_id)
761 mesh_name =
m->getName();
762 max_id = (max_id <
m->getMeshsetId()) ?
m->getMeshsetId() : max_id;
766 CHKERR MEDmeshCr(fid, mesh_name.c_str(), 3, 3, MED_UNSTRUCTURED_MESH,
767 "Mesh created", dtUnit, MED_SORT_DTIT, MED_CARTESIAN,
771 med_int family_id = 0;
772 std::map<std::vector<std::string>, std::tuple<med_int, std::vector<int>>>
774 std::map<EntityHandle, med_int> entityHandle_family_map;
777 shared_meshsets_map[std::vector<string>()] =
778 std::make_tuple(family_id, std::vector<int>());
781 auto get_set_name = [&](
const CubitMeshSets *iit) {
782 std::string bc_type_name;
783 if (iit->getBcTypeULong() &
BLOCKSET) {
785 bc_type_name = iit->getName();
786 if (bc_type_name ==
"NoNameSet") {
787 bc_type_name =
"BLOCKSET_";
788 bc_type_name += std::to_string(iit->getMeshsetId());
791 }
else if (iit->getBcTypeULong() &
SIDESET ||
792 iit->getBcTypeULong() &
NODESET) {
799 if ((iit->getBcType() & jj_bc_type).any()) {
801 std::string temp_bc_type_name = string(
CubitBCNames[jj + 1]);
802 if (temp_bc_type_name.length() > 10) {
803 temp_bc_type_name = temp_bc_type_name.substr(0, 4);
805 bc_type_name += temp_bc_type_name;
810 bc_type_name += std::to_string(iit->getMeshsetId());
812 bc_type_name =
"UnknownSet";
822 for (
auto &entity : *write_range_ptr) {
824 std::vector<int> shared_meshsets;
825 std::vector<string> shared_names;
828 auto add_shared_meshset = [&](
auto &other_meshset) {
829 std::string other_name = get_set_name(other_meshset);
830 if (std::find(shared_names.begin(), shared_names.end(), other_name) ==
831 shared_names.end()) {
832 shared_meshsets.push_back(other_meshset->getMeshsetId());
833 shared_names.push_back(other_name);
838 for (
auto &other_meshset : *meshsets_ptr) {
840 if (med_mesh_name_id == other_meshset->getMeshsetId())
843 Range other_entities;
845 moab.get_entities_by_handle(other_set, other_entities);
846 if (other_entities.empty())
849 bool is_in_meshset = moab.contains_entities(other_set, &entity, 1);
852 add_shared_meshset(other_meshset);
857 auto it = shared_meshsets_map.find(shared_names);
858 if (it == shared_meshsets_map.end()) {
862 it = shared_meshsets_map
864 {shared_names, std::make_tuple(family_id, shared_meshsets)})
868 entityHandle_family_map[entity] = std::get<0>(it->second);
872 for (
const auto &it : shared_meshsets_map) {
874 std::string family_name =
"F_";
875 const auto &[family_id, shared_meshsets] = it.second;
876 family_name += std::to_string(family_id);
878 const auto &shared_meshset_names = it.first;
879 std::string group_name;
880 for (
const auto &name : shared_meshset_names) {
882 std::string meshset_name = name;
883 meshset_name.resize(MED_LNAME_SIZE,
' ');
884 group_name += meshset_name;
888 CHKERR MEDfamilyCr(fid, mesh_name.c_str(), family_name.c_str(), family_id,
889 shared_meshset_names.size(), group_name.c_str());
892 <<
"Creating family " << family_name <<
" with id " << family_id
893 <<
" and " << shared_meshset_names.size() <<
" groups " << std::endl;
898 verts = write_range_ptr->subset_by_type(MBVERTEX);
901 std::vector<med_float> coord_med(3 * verts.size());
902 std::vector<med_int> fam;
903 std::vector<med_int> tags;
904 std::map<EntityHandle, med_int> entityHandle_tag_map;
905 med_int med_node_num = 1;
908 for (Range::iterator it = verts.begin(); it != verts.end(); ++it) {
910 moab.get_coords(&(*it), 1, coords);
911 coord_med[3 * (it - verts.begin())] = coords[0];
912 coord_med[3 * (it - verts.begin()) + 1] = coords[1];
913 coord_med[3 * (it - verts.begin()) + 2] = coords[2];
914 fam.push_back(entityHandle_family_map[*it]);
917 entityHandle_tag_map[*it] = med_node_num;
921 CHKERR MEDmeshNodeWr(fid, mesh_name.c_str(), MED_NO_DT, MED_NO_IT,
922 MED_UNDEF_DT, MED_FULL_INTERLACE, verts.size(),
923 &coord_med[0], MED_FALSE,
"", MED_TRUE, &tags[0],
927 double last_tag = tags.back();
929 for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
932 moab.get_entities_by_type(0, ent_type, entities,
true);
934 ents_to_write = intersect(entities, *write_range_ptr);
936 if (ents_to_write.empty())
940 std::vector<med_int> tag_number;
941 std::vector<med_int> family_number;
942 std::vector<med_int> connectivity;
945 for (
auto &entity : ents_to_write) {
946 if (ent_type != MBVERTEX) {
949 family_number.push_back(entityHandle_family_map[entity]);
951 tag_number.push_back(last_tag);
953 std::vector<EntityHandle> conn;
954 moab.get_connectivity(&entity, 1, conn);
956 for (
auto &
c : conn) {
957 connectivity.push_back(entityHandle_tag_map[
c]);
964 auto get_med_element_type = [](EntityType ent_type) {
965 med_geometrie_element
type;
992 med_geometrie_element med_type = get_med_element_type(ent_type);
993 if (ent_type == MBENTITYSET) {
996 CHKERR MEDmeshElementWr(fid, mesh_name.c_str(), MED_NO_DT, MED_NO_IT, 0.,
997 MED_CELL, med_type, MED_NODAL, MED_FULL_INTERLACE,
998 family_number.size(), &connectivity[0], MED_FALSE,
999 nullptr, MED_TRUE, &tag_number[0], MED_TRUE,
1002 MOFEM_LOG_C(
"MEDWORLD", Sev::inform,
"Writing %i elements of type %i (%s) ",
1003 family_number.size(), med_type,
1004 moab::CN::EntityTypeName(ent_type));