513void msh2geo (istream& in, ostream& out,
string sys_coord_name,
bool do_upgrade)
520 "msh2geo: input stream does not contains a gmsh mesh file ($MeshFormat not found).");
521 double gmsh_fmt_version;
522 size_t file_type, float_data_size;
523 in >> gmsh_fmt_version >> file_type >> float_data_size;
525 "gmsh format version " << gmsh_fmt_version <<
" founded ; expect version 2.2 (HINT: use gmsh -format msh2)");
526 check_macro (file_type == 0,
"msh2geo: unsupported gmsh non-ascii format");
527 check_macro (
scatch(in,
"$EndMeshFormat",
true),
"msh2geo: gmsh input error: $EndMeshFormat not found.");
531 bool full_match =
true;
532 bool partial_match = !full_match;
533 check_macro (
scatch(in,
"$",partial_match),
"msh2geo: gmsh input error: no more label found.");
535 map<size_t, string> phys;
539 if (label ==
"PhysicalNames") {
541 for (
size_t i = 0; i < nphys; i++) {
543 size_t dom_dim, dom_idx;
544 in >> dom_dim >> dom_idx;
548 if (c !=
'"') name.push_back(c);
552 }
while (c !=
'"' && c !=
'\n');
554 size_t start = 0, end = name.length();
555 if (name[start] ==
'"') start++;
556 if (name[end-1] ==
'"') end--;
557 name = name.substr(start,end-start);
559 for (
size_t i = 0; i < name.size(); i++) {
560 if (name[i] ==
' ' || name[i] ==
'\t') name[i] =
'_';
562 phys[dom_idx] = name.substr(start,end-start);
564 check_macro (
scatch(in,
"$EndPhysicalNames",
true),
"msh2geo: gmsh input error ($EndPhysicalNames not found).");
565 check_macro (
scatch(in,
"$",partial_match),
"msh2geo: gmsh input error: no more label found.");
571 if (label !=
"Nodes") {
572 check_macro (
scatch(in,
"$Nodes",
true),
"msh2geo: $Nodes not found in gmsh file");
576 vector<point> node (nnode);
577 double infty = numeric_limits<double>::max();
578 point xmin ( infty, infty, infty);
579 point xmax (-infty, -infty, -infty);
580 for (
size_t k = 0; k < node.size(); k++) {
583 for (
size_t i = 0; i < 3; ++i) {
586 for (
size_t j = 0 ; j < 3; j++) {
587 xmin[j] = min(node[k][j], xmin[j]);
588 xmax[j] = max(node[k][j], xmax[j]);
595 if (xmax[2] == xmin[2]) {
597 if (xmax[1] == xmin[1]) {
599 if (xmax[0] == xmin[0]) dim = 0;
605 check_macro (
scatch(in,
"$Elements",
true),
"msh2geo: $Elements not found in gmsh file");
608 vector<geo_element> element (n_element);
609 vector<size_t> node_subgeo_variant (node.size(), std::numeric_limits<size_t>::max());
610 array<map<size_t, list<size_t> >,4> domain_map;
613 for (
size_t i = 0; i < element.size(); i++) {
614 size_t id, dummy, gmshtype;
615 in >>
id >> gmshtype;
618 size_t domain_idx = 0;
619 for (
size_t j = 0 ; j < n_tag_gmsh; j++) {
626 domain_idx = tag_dummy;
630 "msh2geo: element #" <<
id <<
": unexpected gmsh type '" << gmshtype <<
"'");
632 "msh2geo: element #" <<
id <<
": unsupported gmsh type '" << gmshtype <<
"'");
634 element[i].setname (gmsh_table[gmshtype].name);
635 element[i].setorder(gmsh_table[gmshtype].order);
636 element[i].setgmshtype(gmshtype);
637 size_t nv = gmsh_table[gmshtype].nv;
638 size_t nn = gmsh_table[gmshtype].nn_tot;
639 size_t dim_i = element[i].dimension();
640 size_t variant = element[i].variant();
641 map_dim = max(map_dim,dim_i);
643 order = element[i].order();
645 check_macro (order == element[i].order(),
"msh2geo: unexpected order="<<element[i].order());
647 element[i].resize(nn);
648 for (
size_t j = 0; j < nn; j++) {
652 for (
size_t subgeo_variant = 0; subgeo_variant <= variant; subgeo_variant++) {
653 for (
size_t loc_inod = reference_element_first_inod_by_variant(variant,element[i].order(),subgeo_variant),
654 loc_nnod = reference_element_last_inod_by_variant(variant,element[i].order(),subgeo_variant);
655 loc_inod < loc_nnod; loc_inod++) {
656 node_subgeo_variant [element[i][loc_inod]] = subgeo_variant;
661 if (domain_idx != 0) {
662 domain_map[dim_i][domain_idx].push_back(i);
665 array<size_t,reference_element__max_variant> loc_ndof_by_variant;
666 reference_element_init_local_nnode_by_variant (order, loc_ndof_by_variant);
668 size_t last_index = 0;
669 for (
size_t variant = reference_element_first_variant_by_dimension(map_dim);
670 variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
671 for (
size_t ie = 0; ie < element.size(); ++ie) {
672 if (element[ie].variant() != variant)
continue;
673 element[ie].setindex (last_index);
675 case 0: element[ie].resize(1); element[ie][0] = last_index;
break;
676 case 1: element[ie].edge_indirect(0).setindex(last_index);
break;
677 case 2: element[ie].face_indirect(0).setindex(last_index);
break;
687 vector<size_t> old2new_inode (node.size(), std::numeric_limits<size_t>::max());
688 if (
true || !do_upgrade) {
690 size_t new_inode = 0;
691 for (
size_t subgeo_variant = 0; subgeo_variant < reference_element_last_variant_by_dimension(map_dim); subgeo_variant++) {
692 for (
size_t old_inode = 0; old_inode < node.size(); old_inode++) {
693 if (node_subgeo_variant [old_inode] != subgeo_variant)
continue;
694 old2new_inode[old_inode] = new_inode++;
699 for (
size_t inode = 0; inode < node.size(); inode++) {
700 old2new_inode[inode] = inode;
703 for (
size_t inode = 0; inode < node.size(); inode++) {
704 check_macro (old2new_inode[inode] != std::numeric_limits<size_t>::max(),
"msh2geo: invalid permutation");
707 vector<size_t> new2old_inode (node.size(), std::numeric_limits<size_t>::max());
708 for (
size_t inode = 0; inode < node.size(); inode++) {
709 new2old_inode[old2new_inode[inode]] = inode;
711 for (
size_t inode = 0; inode < node.size(); inode++) {
712 check_macro (new2old_inode[inode] != std::numeric_limits<size_t>::max(),
"msh2geo: invalid inverse permutation for inode="<<inode);
716 vector<point> new_node (node.size());
717 for (
size_t old_inode = 0; old_inode < node.size(); old_inode++) {
718 new_node [old2new_inode[old_inode]] = node [old_inode];
720 for (
size_t inode = 0; inode < node.size(); inode++) {
721 node [inode] = new_node [inode];
724 for (
size_t ie = 0; ie < element.size(); ++ie) {
725 for (
size_t i = 0; i < element[ie].size(); ++i) {
726 element[ie][i] = old2new_inode[element[ie][i]];
732 vector<index_set> edge_ball (node.size());
733 vector<geo_element_side>
edge;
734 if (do_upgrade && map_dim >= 2) {
735 list <geo_element_side> edge_list;
737 for (
size_t ie = 0; ie < element.size(); ++ie) {
738 if (element[ie].
dimension() != map_dim)
continue;
739 for (
size_t iedg = 0; iedg < element[ie].n_edge(); ++iedg) {
740 size_t iv0 = element[ie].edge_local_vertex(iedg, 0);
741 size_t iv1 = element[ie].edge_local_vertex(iedg, 1);
742 size_t inod0 = element[ie][iv0];
743 size_t inod1 = element[ie][iv1];
747 "msh2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<
")");
748 if (iedge_set.size() == 1)
continue;
749 edge_ball[inod0].insert (nedg);
750 edge_ball[inod1].insert (nedg);
751 geo_element_side new_edge;
752 new_edge.setname(
'e');
755 new_edge.setindex(nedg);
756 edge_list.push_back (new_edge);
760 edge.resize (edge_list.size());
761 std::copy (edge_list.begin(), edge_list.end(),
edge.begin());
764 if (do_upgrade && map_dim >= 2) {
765 for (
size_t ie = 0; ie < element.size(); ++ie) {
766 if (element[ie].
dimension() != map_dim)
continue;
767 for (
size_t iedg = 0; iedg < element[ie].n_edge(); ++iedg) {
768 size_t iv0 = element[ie].edge_local_vertex(iedg, 0);
769 size_t iv1 = element[ie].edge_local_vertex(iedg, 1);
770 size_t inod0 = element[ie][iv0];
771 size_t inod1 = element[ie][iv1];
775 "msh2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<
")");
776 size_t iedge = *(iedge_set.begin());
777 element[ie].edge_indirect(iedg).setindex (iedge);
778 element[ie].edge_indirect(iedg).setorientation((
edge[iedge][0] == inod0) ? 1 : -1);
785 vector<index_set> face_ball (node.size());
786 vector<geo_element_side>
face;
787 if (do_upgrade && map_dim >= 3) {
788 list <geo_element_side> face_list;
790 for (
size_t ie = 0; ie < element.size(); ++ie) {
791 if (element[ie].
dimension() != map_dim)
continue;
792 for (
size_t loc_ifac = 0; loc_ifac < element[ie].n_face(); ++loc_ifac) {
793 size_t iv0 = element[ie].face_local_vertex(loc_ifac, 0);
794 size_t inod0 = element[ie][iv0];
796 for (
size_t j = 1; j < element[ie].face_n_vertex(loc_ifac); ++j) {
797 size_t ivj = element[ie].face_local_vertex(loc_ifac, j);
798 size_t inodj = element[ie][ivj];
802 "msh2geo: error: connectivity problem (iface.size="<<iface_set.size()<<
")");
803 if (iface_set.size() == 1)
continue;
804 geo_element_side new_face;
805 char face_name = (element[ie].face_n_vertex(loc_ifac) == 3) ?
't' :
'q';
806 new_face.setname(face_name);
807 new_face.setindex(nfac);
808 for (
size_t j = 0; j < element[ie].face_n_vertex(loc_ifac); ++j) {
809 size_t ivj = element[ie].face_local_vertex(loc_ifac, j);
810 size_t inodj = element[ie][ivj];
811 face_ball[inodj].insert (nfac);
814 face_list.push_back (new_face);
818 face.resize (face_list.size());
819 std::copy (face_list.begin(), face_list.end(),
face.begin());
822 if (do_upgrade && map_dim >= 3) {
823 for (
size_t ie = 0; ie < element.size(); ++ie) {
824 if (element[ie].
dimension() != map_dim)
continue;
825 for (
size_t loc_ifac = 0; loc_ifac < element[ie].n_face(); ++loc_ifac) {
826 size_t iv0 = element[ie].face_local_vertex(loc_ifac, 0);
827 size_t inod0 = element[ie][iv0];
829 for (
size_t j = 1; j < element[ie].face_n_vertex(loc_ifac); ++j) {
830 size_t ivj = element[ie].face_local_vertex(loc_ifac, j);
831 size_t inodj = element[ie][ivj];
835 "msh2geo: error: connectivity problem (iface.size="<<iface_set.size()<<
")");
836 size_t iface = *(iface_set.begin());
837 orientation_type orient;
839 if (element[ie].face_n_vertex(loc_ifac) == 3) {
840 geo_element_get_orientation_and_shift (
842 element[ie][element[ie].face_local_vertex(loc_ifac, 0)],
843 element[ie][element[ie].face_local_vertex(loc_ifac, 1)],
844 element[ie][element[ie].face_local_vertex(loc_ifac, 2)],
847 geo_element_get_orientation_and_shift (
849 element[ie][element[ie].face_local_vertex(loc_ifac, 0)],
850 element[ie][element[ie].face_local_vertex(loc_ifac, 1)],
851 element[ie][element[ie].face_local_vertex(loc_ifac, 2)],
852 element[ie][element[ie].face_local_vertex(loc_ifac, 3)],
855 element[ie].face_indirect(loc_ifac).setindex (iface);
856 element[ie].face_indirect(loc_ifac).setorientation (orient);
857 element[ie].face_indirect(loc_ifac).setshift (shift);
864 array<size_t,reference_element__max_variant> size_by_variant;
865 array<size_t,reference_element__max_variant+1> first_by_variant;
866 size_by_variant.fill (0);
867 first_by_variant.fill (0);
870 for (
size_t inode = 0; inode < node.size(); inode++) {
871 if (node_subgeo_variant [inode] == reference_element__p)
n_vertex++;
873 size_by_variant [reference_element__p] =
n_vertex;
874 if (map_dim >= 2 &&
edge.size() != 0) {
875 size_by_variant [reference_element__e] =
edge.size();
877 if (map_dim >= 3 &&
face.size() != 0) {
878 for (
size_t loc_ifac = 0; loc_ifac <
face.size(); ++loc_ifac) {
879 size_by_variant [
face[loc_ifac].variant()]++;
882 for (
size_t ie = 0; ie < element.size(); ++ie) {
883 if (element[ie].
dimension() != map_dim)
continue;
884 size_t variant = element[ie].variant();
885 size_by_variant [variant]++;
887 for (
size_t dim = 0; dim <= 3; ++dim) {
888 for (
size_t variant = reference_element_first_variant_by_dimension(dim);
889 variant < reference_element_last_variant_by_dimension(dim); variant++) {
890 first_by_variant [variant+1] = size_by_variant [variant];
899 std::fill (old2new_inode.begin(), old2new_inode.end(), std::numeric_limits<size_t>::max());
900 std::fill (new2old_inode.begin(), new2old_inode.end(), std::numeric_limits<size_t>::max());
901 for (
size_t variant = reference_element_first_variant_by_dimension(map_dim);
902 variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
903 for (
size_t ie = 0; ie < element.size(); ++ie) {
904 if (element[ie].variant() != variant)
continue;
906 std::vector<size_t> new_K;
908 order, size_by_variant, first_by_variant, loc_ndof_by_variant, old_K, new_K);
909 for (
size_t loc_inod = 0; loc_inod < old_K.
size(); loc_inod++) {
910 size_t old_inod = old_K [loc_inod];
911 size_t new_inod = new_K [loc_inod];
912 if (old2new_inode [old_inod] != std::numeric_limits<size_t>::max())
continue;
913 old2new_inode [old_inod] = new_inod;
919 for (
size_t inode = 0; inode < node.size(); inode++) {
920 old2new_inode[inode] = inode;
923 for (
size_t inode = 0; inode < node.size(); inode++) {
924 check_macro (old2new_inode[inode] != std::numeric_limits<size_t>::max(),
"msh2geo: invalid permutation");
927 for (
size_t inode = 0; inode < node.size(); inode++) {
928 new2old_inode[old2new_inode[inode]] = inode;
930 for (
size_t inode = 0; inode < node.size(); inode++) {
931 check_macro (new2old_inode[inode] != std::numeric_limits<size_t>::max(),
"msh2geo: invalid inverse permutation for inode="<<inode);
940 static const char* geo_variant_name [reference_element__max_variant] = {
949 out << setprecision(numeric_limits<double>::digits10)
955 <<
" dimension\t" << dim << endl;
956 if (sys_coord_name !=
"cartesian") {
957 out <<
" coordinate_system " << sys_coord_name << endl;
960 out <<
" order\t\t" << order << endl;
962 out <<
" nodes\t\t" << node.size() << endl;
965 for (
size_t variant = reference_element_first_variant_by_dimension(map_dim);
966 variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
967 if (size_by_variant[variant] > 0) {
968 out <<
" " << geo_variant_name [variant] <<
"\t" << size_by_variant[variant] << endl;
972 if (map_dim >= 3 && size_by_variant[reference_element__t] != 0) {
973 out <<
" triangles\t" << size_by_variant[reference_element__t] << endl;
975 if (map_dim >= 3 && size_by_variant[reference_element__q] != 0) {
976 out <<
" quadrangles\t" << size_by_variant[reference_element__q] << endl;
978 if (map_dim >= 2 && size_by_variant[reference_element__e] != 0) {
979 out <<
" edges\t" << size_by_variant[reference_element__e] << endl;
981 out <<
"end header" << endl
984 for (
size_t inode = 0; inode < node.size(); inode++) {
985 for (
size_t i = 0; i < dim; ++i) {
986 out << node[new2old_inode[inode]][i];
987 if (i+1 != dim) out <<
" ";
993 for (
size_t variant = reference_element_first_variant_by_dimension(map_dim);
994 variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
995 for (
size_t ie = 0; ie < element.size(); ie++) {
996 if (element[ie].variant() != variant)
continue;
998 if (element[ie].name() !=
'e' || element[ie].order() > 1) {
999 out << element[ie].name() <<
"\t";
1001 if (element[ie].order() > 1) {
1002 out <<
"p" << element[ie].order() <<
" ";
1004 for (
size_t iloc = 0, nloc = element[ie].size(); iloc < nloc; iloc++) {
1005 out << element[ie][iloc];
1006 if (iloc+1 != nloc) out <<
" ";
1010 out << element[ie].name() <<
"\t";
1011 for (
size_t iloc = 0, nloc = element[ie].
n_vertex(); iloc < nloc; iloc++) {
1012 out << element[ie][iloc];
1013 if (iloc+1 != nloc) out <<
" ";
1021 for (
size_t variant = reference_element_first_variant_by_dimension(2);
1022 variant < reference_element_last_variant_by_dimension(2); variant++) {
1023 for (
size_t ie = 0; ie <
face.size(); ie++) {
1024 if (
face[ie].variant() != variant)
continue;
1025 out <<
face[ie].name() <<
"\t";
1026 for (
size_t iloc = 0, nloc =
face[ie].
n_vertex(); iloc < nloc; iloc++) {
1027 out <<
face[ie][iloc];
1028 if (iloc+1 != nloc) out <<
" ";
1034 for (
size_t ie = 0, ne =
edge.size(); ie < ne; ++ie) {
1035 out <<
"e\t" << old2new_inode[
edge[ie][0]]
1036 <<
" " << old2new_inode[
edge[ie][1]] << endl;
1042 for (
size_t d = 0;
d <= 3; ++
d) {