#=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~=~
#                    numgeo - GiD interface
#             Copyright (C) 2023-today Jan Machacek, Rahul Dogra
#=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~=~
#
# History
# - 26.11.2023, J. Machacek - Initial version
# - 03.03.2025, R. Dogra - Added truss elements
# - 07.03.2025, R. Dogra - Added beam elements
# - 21.06.2025, R. Dogra - Removed repeation for contacts
# - 24.06.2025, R. Dogra - Corrected 20 and 27 hexahedral element
# - 16.09.2025, R. Dogra - Added surfaces for flow condition (in/out-flow)
# - 01.01.2026, R. Dogra - Added errors for the element types
#
# Memo:
# - get_domnode_attribute returns the current value of "n" (name) of $group
#
#
#=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~=~

proc numgeo::WriteMesh { filename } {

  append filename "-mesh.inp"

  customlib::InitWriteFile $filename
  
  # xml document to get some tree data
  set root [$::gid_groups_conds::doc documentElement] ;
  
  customlib::WriteString "**=================================================================="
  customlib::WriteString "**            numgeo input file"  
  customlib::WriteString "**=================================================================="
  customlib::WriteString "**"  
  customlib::WriteString "** Input file for the finite element code numgeo (Machacek & Staubach)"
  customlib::WriteString "**"   
  customlib::WriteString "** The input file was created by numgeo-GiD (Machacek & Dogra)"   
  customlib::WriteString "**"   
  customlib::WriteString "** numgeo is available from www.numgeo.de"   
  customlib::WriteString "** numgeo-GiD can be downloaded from github XXX" 
  customlib::WriteString "**"  
  customlib::WriteString "** Contact:"
  customlib::WriteString "** - Jan Machacek, jan.machacek@tu-darmstadt.de"
  customlib::WriteString "** - Rahul Dogra, rahul.dogra@tu-darmstadt.de"
  customlib::WriteString "**"
  customlib::WriteString "**=================================================================="
  
  
  #
  # get the dimensions of the model: 2D (plane-strain), 2D: Axissymmetric, 3D
  # 
  set dim_path {string(//container[@n="_dimensions_units_"]/value[@n="_dimensions_"]/@v)}
  set dim [$root selectNodes $dim_path]
  
  #
  # Nodes
  # 
  
  set units_mesh [gid_groups_conds::give_mesh_unit]
  customlib::WriteString ""
  customlib::WriteString "**"
  customlib::WriteString "** Nodes"
  customlib::WriteString "**"
  customlib::WriteString ""
  customlib::WriteString "*Node"
  if {$dim == "2D (plane-strain)"} {
    customlib::WriteCoordinates "%5d, %14.5e, %14.5e%.0s\n"
  } elseif {$dim == "2D: Axissymmetric"} {
    customlib::WriteCoordinates "%5d, %14.5e, %14.5e%.0s\n"
  } elseif {$dim == "3D"} {
    customlib::WriteCoordinates "%5d, %14.5e, %14.5e, %14.5e\n"
  }

  #
  # Elements
  #
  customlib::WriteString ""
  customlib::WriteString "**"
  customlib::WriteString "** Elements"
  customlib::WriteString "**"
  customlib::WriteString ""
  
  set document [$::gid_groups_conds::doc documentElement]
  set xp_solid [format_xpath {container[@n="_element_type_"]/condition[@n="_element_types_"]/group}]
  set xp_truss [format_xpath {container[@n="_element_type_"]/condition[@n="_truss_element_types_"]/group}]
  set solid_groups [$document selectNodes $xp_solid]
  set truss_groups [$document selectNodes $xp_truss]
  if {[llength $solid_groups] == 0 && [llength $truss_groups] == 0} {
    error [= "No element type assigned. Go to: Assign element type (or Assign truss/beam)"]
  }

  # detect presence of 2D element types in the mesh
  set has_tri  [expr {[llength [GiD_Mesh list -element_type Triangle element]] > 0}]
  set has_quad [expr {[llength [GiD_Mesh list -element_type Quadrilateral element]] > 0}]
  # detect presence of 3D element types in the mesh
  set has_hex  [expr {[llength [GiD_Mesh list -element_type Hexahedra element]] > 0}]
  set has_tet  [expr {[llength [GiD_Mesh list -element_type Tetrahedra element]] > 0}]
  # 2D/Axisymmetric selected, but 3D elements exist
  if {$dim ne "3D" && ($has_hex || $has_tet)} {
    set shapes {}
    if {$has_hex} { lappend shapes "Hexahedra" }
    if {$has_tet} { lappend shapes "Tetrahedra" }
    error [= "Mesh contains 3D elements ([join $shapes {, }]) but the selected dimension is '$dim'. Please set Dimensions/Units → Dimensions = 3D, or remesh using 2D elements (Triangle/Quadrilateral)."]
  }
  # 3D selected, but 2D elements exist
  if {$dim eq "3D" && ($has_tri || $has_quad)} {
    set shapes {}
    if {$has_tri}  { lappend shapes "Triangle" }
    if {$has_quad} { lappend shapes "Quadrilateral" }
    error [= "Mesh contains 2D elements ([join $shapes {, }]) but the selected dimension is '3D'. Please set Dimensions/Units → Dimensions = 2D (plane-strain) or 2D: Axissymmetric, or remesh using 3D elements (Hexahedra/Tetrahedra)."]
  }

  # assemble all triangular elements (Linear and quadratic)
  set elements_tri [GiD_Mesh list -element_type Triangle element]
  set elements_tri3 {}
  set elements_tri6 {}
  foreach i $elements_tri {
    set nodes [GiD_Mesh get element $i connectivities]
    if {[llength $nodes] == 3} {
      lappend elements_tri3 $i
    } elseif {[llength $nodes] == 6} {
      lappend elements_tri6 $i
    }
  }

  # assemble all rectangular elements (Linear and quadratic)
  set elements_quad [GiD_Mesh list -element_type Quadrilateral element]
  set elements_quad4 {}
  set elements_quad8 {}
  set elements_quad20 {}
  foreach i $elements_quad {
    set nodes [GiD_Mesh get element $i connectivities]
    if {[llength $nodes] == 4} {
      lappend elements_quad4 $i
    } elseif {[llength $nodes] == 8} {
      lappend elements_quad8 $i
    } 
  }
  
  # assemble all hexahedral elements (Linear and quadratic)
  set elements_hex [GiD_Mesh list -element_type Hexahedra element]
  set elements_hex8 {}
  set elements_hex20 {}
  set elements_hex27 {}
  foreach i $elements_hex {
    set nodes [GiD_Mesh get element $i connectivities]
    if {[llength $nodes] == 8} {
      lappend elements_hex8 $i
    } elseif {[llength $nodes] == 20} {
      lappend elements_hex20 $i
    } elseif {[llength $nodes] == 27} {
      lappend elements_hex27 $i
    }
  }
  
   # assemble all tetrahedral elements (Linear and quadratic)
  set elements_tet [GiD_Mesh list -element_type Tetrahedra element]
  set elements_tet4 {}
  set elements_tet10 {}
  foreach i $elements_tet {
    set nodes [GiD_Mesh get element $i connectivities]
    if {[llength $nodes] == 4} {
      lappend elements_tet4 $i
    } elseif {[llength $nodes] == 10} {
      lappend elements_tet10 $i
    }
  }

  set document [$::gid_groups_conds::doc documentElement]
  set xp [format_xpath {container[@n="_element_type_"]/condition[@n="_element_types_"]/group}]
  # loop over all groups to which element types were assigned
  foreach group [$document selectNodes $xp] {
    set group_name [get_domnode_attribute $group n]
    set node [$group selectNodes "./value\[@n = '_element_formulation_'\]"]
    set eltype [get_domnode_attribute $node v]
    set node [$group selectNodes "./value\[@n = '_integration_'\]"]
    set integration [get_domnode_attribute $node v]
    set group_element_ids [GiD_EntitiesGroups get $group_name elements]

    #
    # 2D elements
    #
    if {$dim == "2D (plane-strain)"} {
	
      #
      # 3 node triangle
      #
      if {[llength $elements_tri3] > 0} {
        set keyword [return_eltype $dim $eltype $integration {Triangular} {3}]
        set print_eltype 0
        # loop over all group members
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          # get the element id and check if it is the correct element type
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_tri3 $element_id] >= 0} {
            if {$print_eltype == 0} {
              # print header, only for first member of the group
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            # print connectivity of i-th element of group
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3
            customlib::WriteString "$element_id, $c1, $c2, $c3"
          }
        }
        customlib::WriteString ""
      }
      #
      # 6 node triangle
      #
      if {[llength $elements_tri6] > 0} {
        set keyword [return_eltype $dim $eltype $integration {Triangular} {6}]
        set print_eltype 0
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_tri6 $element_id] >= 0} {
            if {$print_eltype == 0} {
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3 c4 c5 c6
            customlib::WriteString "$element_id, $c1, $c2, $c3, $c4, $c5, $c6"
          }
        }
        customlib::WriteString ""
      }
      #
      # 4 node quadrilateral
      #
      if {[llength $elements_quad4] > 0} {
        set keyword [return_eltype $dim $eltype $integration {Quadrilateral} {4}]
        set print_eltype 0
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_quad4 $element_id] >= 0} {
            if {$print_eltype == 0} {
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3 c4
            customlib::WriteString "$element_id, $c1, $c2, $c3, $c4"
          }
        }
        customlib::WriteString ""
      }
      #
      # 4 node quadrilateral
      #
      if {[llength $elements_quad8] > 0} {
        set keyword [return_eltype $dim $eltype $integration {Quadrilateral} {8}]
		    set print_eltype 0
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_quad8 $element_id] >= 0} {
            if {$print_eltype == 0} {
              customlib::WriteString "*Element, type=$keyword"
			        set print_eltype 1
            }
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3 c4 c5 c6 c7 c8
            customlib::WriteString "$element_id, $c1, $c2, $c3, $c4, $c5, $c6, $c7, $c8"
          }
        }
        customlib::WriteString ""
      }
      #
      # 2D axisymmetric elements
      #
    } elseif {$dim == "2D: Axissymmetric"} {
      #
      # 3 node triangle - axisymmetric
      #
      if {[llength $elements_tri3] > 0} {
        set keyword [return_eltype $dim $eltype $integration {Triangular} {3}]
        set print_eltype 0
        # loop over all group members
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          # get the element id and check if it is the correct element type
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_tri3 $element_id] >= 0} {
            if {$print_eltype == 0} {
              # print header, only for first member of the group
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            # print connectivity of i-th element of group
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3
            customlib::WriteString "$element_id, $c1, $c2, $c3"
          }
        }
        customlib::WriteString ""
      }
      #
      # 6 node triangle - axisymmetric
      #
      if {[llength $elements_tri6] > 0} {
        set keyword [return_eltype $dim $eltype $integration {Triangular} {6}]
        set print_eltype 0
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_tri6 $element_id] >= 0} {
            if {$print_eltype == 0} {
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3 c4 c5 c6
            customlib::WriteString "$element_id, $c1, $c2, $c3, $c4, $c5, $c6"
          }
        }
        customlib::WriteString ""
      }
      #
      # 4 node quadrilateral - axisymmetric
      #
      if {[llength $elements_quad4] > 0} {
        set keyword [return_eltype $dim $eltype $integration {Quadrilateral} {4}]
        set print_eltype 0
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_quad4 $element_id] >= 0} {
            if {$print_eltype == 0} {
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3 c4
            customlib::WriteString "$element_id, $c1, $c2, $c3, $c4"
          }
        }
        customlib::WriteString ""
      }
      #
      # 8 node quadrilateral - axisymmetric
      #
      if {[llength $elements_quad8] > 0} {
        set keyword [return_eltype $dim $eltype $integration {Quadrilateral} {8}]
        set print_eltype 0
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_quad8 $element_id] >= 0} {
            if {$print_eltype == 0} {
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3 c4 c5 c6 c7 c8
            customlib::WriteString "$element_id, $c1, $c2, $c3, $c4, $c5, $c6, $c7, $c8"
          }
        }
        customlib::WriteString ""
      }
      #
      # 3D elements
      #
    } elseif {$dim == "3D"} {
      #
      # 8 node hexahedra
      #
      if {[llength $elements_hex8] > 0} {
        set keyword [return_eltype $dim $eltype $integration {Hexahedra} {8}]
        set print_eltype 0
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_hex8 $element_id] >= 0} {
            if {$print_eltype == 0} {
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3 c4 c5 c6 c7 c8
            customlib::WriteString "$element_id, $c1, $c2, $c3, $c4, $c5, $c6, $c7, $c8"
          }
        }
        customlib::WriteString ""
      }
      #
      # 20 node hexahedra
      #
      if {[llength $elements_hex20] > 0} {
        set keyword [return_eltype $dim $eltype $integration {Hexahedra} {20}]
        set print_eltype 0
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_hex20 $element_id] >= 0} {
            if {$print_eltype == 0} {
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13 c14 c15 c16 c17 c18 c19 c20
            customlib::WriteString "$element_id, $c1, $c2, $c3, $c4, $c5, $c6, $c7, $c8, $c9, $c10, $c11, $c12, $c17, $c18, $c19, $c20, $c13, $c14, $c15, $c16"
          }
        }
        customlib::WriteString ""
      }
      #
      # 27 node hexahedra
      #
      if {[llength $elements_hex27] > 0} {
        set keyword [return_eltype $dim $eltype $integration {Hexahedra} {27}]
        set print_eltype 0
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_hex27 $element_id] >= 0} {
            if {$print_eltype == 0} {
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12 c13 c14 c15 c16 c17 c18 c19 c20 c21 c22 c23 c24 c25 c26 c27
            customlib::WriteString "$element_id, $c1, $c2, $c3, $c4, $c5, $c6, $c7, $c8, $c9, $c10, $c11, $c12, $c17, $c18, $c19, $c20, $c13, $c14, $c15, $c16, $c25, $c23, $c22, $c24, $c21, $c26, $c27"
          }
        }
        customlib::WriteString ""
      }
      #
      # 4 node tetrahedra
      #
      if {[llength $elements_tet4] > 0} {
        set keyword [return_eltype $dim $eltype $integration {Tetrahedra} {4}]
        set print_eltype 0
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_tet4 $element_id] >= 0} {
            if {$print_eltype == 0} {
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3 c4
            customlib::WriteString "$element_id, $c1, $c2, $c3, $c4"
          }
        }
        customlib::WriteString ""
      }
      #
      # 10 node tetrahedra
      #
      if {[llength $elements_tet10] > 0} {
        set keyword [return_eltype $dim $eltype $integration {Tetrahedra} {10}]
		    set print_eltype 0
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_tet10 $element_id] >= 0} {
            if {$print_eltype == 0} {
              customlib::WriteString "*Element, type=$keyword"
			        set print_eltype 1
            }
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3 c4 c5 c6 c7 c8 c9 c10
            customlib::WriteString "$element_id, $c1, $c2, $c3, $c4, $c5, $c6, $c7, $c8, $c9, $c10"
          }
        }
        customlib::WriteString ""
      }
    }
  }
  
  #
  # Truss elements (Linear and quadratic)
  #
  set elements_lin [GiD_Mesh list -element_type Linear element]
  set elements_lin2 {}
  set elements_lin3 {}
  foreach i $elements_lin {
    set nodes [GiD_Mesh get element $i connectivities]
    if {[llength $nodes] == 2} {
      lappend elements_lin2 $i
    } elseif {[llength $nodes] == 3} {
      lappend elements_lin3 $i
    }
  }
  
  set document [$::gid_groups_conds::doc documentElement]
  set xp [format_xpath {container[@n="_element_type_"]/condition[@n="_truss_element_types_"]/group}]
  # loop over all groups to which element types were assigned
  foreach group [$document selectNodes $xp] {
    set group_name [get_domnode_attribute $group n]
    set node [$group selectNodes "./value\[@n = '_truss_beam_'\]"]
    set t_b_type [get_domnode_attribute $node v]
    set group_element_ids [GiD_EntitiesGroups get $group_name elements]
    if {$dim == "2D (plane-strain)"} {
  	  #
      # 2 node Linear
      #
      if {[llength $elements_lin2] > 0} {
        set keyword [return_eltype2 $dim $t_b_type {Linear} {2}]
        set print_eltype 0
        # loop over all group members
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          # get the element id and check if it is the correct element type
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_lin2 $element_id] >= 0} {
            if {$print_eltype == 0} {
              # print header, only for first member of the group
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            # print connectivity of i-th element of group
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2
            customlib::WriteString "$element_id, $c1, $c2"
          }
        }
        customlib::WriteString ""
      }
      # 3 node Linear
      if {[llength $elements_lin3] > 0} {
        set keyword [return_eltype2 $dim $t_b_type  {Linear} {3}]
        set print_eltype 0
        # loop over all group members
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          # get the element id and check if it is the correct element type
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_lin3 $element_id] >= 0} {
            if {$print_eltype == 0} {
              # print header, only for first member of the group
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            # print connectivity of i-th element of group
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3
            customlib::WriteString "$element_id, $c1, $c2, $c3"
          }
        }
        customlib::WriteString ""
      }
    } elseif {$dim == "2D: Axissymmetric"} {
      # 2 node Linear
      if {[llength $elements_lin2] > 0} {
        set keyword [return_eltype2 $dim $t_b_type {Linear} {2}]
        set print_eltype 0
        # loop over all group members
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          # get the element id and check if it is the correct element type
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_lin2 $element_id] >= 0} {
            if {$print_eltype == 0} {
              # print header, only for first member of the group
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            # print connectivity of i-th element of group
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2
            customlib::WriteString "$element_id, $c1, $c2"
          }
        }
        customlib::WriteString ""
      }
      # 3 node Linear
      if {[llength $elements_lin3] > 0} {
        set keyword [return_eltype2 $dim $t_b_type {Linear} {3}]
        set print_eltype 0
        # loop over all group members
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          # get the element id and check if it is the correct element type
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_lin3 $element_id] >= 0} {
            if {$print_eltype == 0} {
              # print header, only for first member of the group
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            # print connectivity of i-th element of group
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3
            customlib::WriteString "$element_id, $c1, $c2, $c3"
          }
        }
        customlib::WriteString ""
      }
    } elseif {$dim == "3D"} {
      # 2 node Linear
      if {[llength $elements_lin2] > 0} {
        set keyword [return_eltype2 $dim $t_b_type {Linear} {2}]
        set print_eltype 0
        # loop over all group members
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          # get the element id and check if it is the correct element type
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_lin2 $element_id] >= 0} {
            if {$print_eltype == 0} {
              # print header, only for first member of the group
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            # print connectivity of i-th element of group
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2
            customlib::WriteString "$element_id, $c1, $c2"
          }
        }
        customlib::WriteString ""
      }
      # 3 node Linear
      if {[llength $elements_lin3] > 0} {
        set keyword [return_eltype2 $dim $t_b_type {Linear} {3}]
        set print_eltype 0
        # loop over all group members
        for {set i 0} {$i < [llength $group_element_ids]} {incr i} {
          # get the element id and check if it is the correct element type
          set element_id [lindex $group_element_ids $i]
          if {[lsearch -exact $elements_lin3 $element_id] >= 0} {
            if {$print_eltype == 0} {
              # print header, only for first member of the group
              customlib::WriteString "*Element, type=$keyword"
              set print_eltype 1
            }
            # print connectivity of i-th element of group
            lassign [GiD_Mesh get element $element_id connectivities] c1 c2 c3
            customlib::WriteString "$element_id, $c1, $c2, $c3"
          }
        }
        customlib::WriteString ""
      }
    }
  }
  #
  # Node sets
  #
  
  customlib::WriteString ""
  customlib::WriteString "**"
  customlib::WriteString "** Node sets"
  customlib::WriteString "**"
  customlib::WriteString ""
  
  set node_set_list {}
  
  set groups [GiD_Groups list]
  foreach group_name $groups {
    set group_entity_ids [GiD_EntitiesGroups get $group_name nodes]
    set entity_count [llength $group_entity_ids]
    if {$entity_count > 0} {
      lappend node_set_list $group_name
      set variables_per_line 16
      customlib::WriteString ""
      customlib::WriteString "*NSet, nset=$group_name"
      for {set i 0} {$i < $entity_count} {incr i $variables_per_line} {
        set end_index [expr {$i + $variables_per_line - 1}]
        set end_index [expr {$end_index >= $entity_count ? $entity_count - 1 : $end_index}]
        customlib::WriteString [join [lrange $group_entity_ids $i $end_index] ", "]
      }
    }
  }
  
  #
  # Element sets
  #

  customlib::WriteString ""
  customlib::WriteString "**"
  customlib::WriteString "** Element sets"
  customlib::WriteString "**"
  customlib::WriteString ""
  
  set groups [GiD_Groups list]
  foreach group_name $groups {
    set group_entity_ids [GiD_EntitiesGroups get $group_name elements]
    set entity_count [llength $group_entity_ids]
    if {$entity_count > 0} {
      set variables_per_line 16
      customlib::WriteString ""
      customlib::WriteString "*Elset, elset=$group_name"
      for {set i 0} {$i < $entity_count} {incr i $variables_per_line} {
        set end_index [expr {$i + $variables_per_line - 1}]
        set end_index [expr {$end_index >= $entity_count ? $entity_count - 1 : $end_index}]
        customlib::WriteString [join [lrange $group_entity_ids $i $end_index] ", "]
      }
    }
  }
  
  #
  # Surfaces
  #

  customlib::WriteString ""
  customlib::WriteString "**"
  customlib::WriteString "** Surfaces (generated from node sets)"
  customlib::WriteString "**"
  
  # we create surfaces from node sets
  # to keep the cost low we only create surfaces for node sets used to assign surface loads or contact conditions.
  
  # Warning: this seems to be a rather naive approach - cant this be done better?
  
  set group_names {}
  set xp [format_xpath {container[@n="_steps_"]/blockdata}]
  set list [$root selectNodes $xp]
  # loop over all steps
  set step_num 0
  foreach istep $list {
    set step_num [expr $step_num + 1]
    # all groups to which DSload "Drain" is applied
    set xp "./container\[@n='_loads_'\]/condition\[@n='_drain_'\]/group"
    set groups [$istep selectNodes $xp]
    foreach group $groups {
      set group_name [get_domnode_attribute $group n]
      if {[lsearch -exact $group_names $group_name] == -1} {
        lappend group_names $group_name
      }
    }
    
    # all groups to which DSload "Distributed load" is applied
    set xp "./container\[@n='_loads_'\]/condition\[@n='_dis_load_'\]/group"
    set groups [$istep selectNodes $xp]
    foreach group $groups {
      set group_name [get_domnode_attribute $group n]
      if {[lsearch -exact $group_names $group_name] == -1} {
        lappend group_names $group_name
      }
    }
	
  # all groups to which DSload "Pressure" is applied
    set xp "./container\[@n='_loads_'\]/condition\[@n='_surface_pressure_'\]/group"
    set groups [$istep selectNodes $xp]
    foreach group $groups {
      set group_name [get_domnode_attribute $group n]
      if {[lsearch -exact $group_names $group_name] == -1} {
        lappend group_names $group_name
      }
    }
	
  # all groups to which DSload "Pressure" is applied
    set xp "./container\[@n='_loads_'\]/condition\[@n='_flow_condition_'\]/group"
    set groups [$istep selectNodes $xp]
    foreach group $groups {
      set group_name [get_domnode_attribute $group n]
      if {[lsearch -exact $group_names $group_name] == -1} {
        lappend group_names $group_name
      }
    }
	
  # all groups to which DSload "Water Pressure" is applied
    set xp "./container\[@n='_loads_'\]/condition\[@n='_surface_water_pressure_'\]/group"
    set groups [$istep selectNodes $xp]
    foreach group $groups {
      set group_name [get_domnode_attribute $group n]
      if {[lsearch -exact $group_names $group_name] == -1} {
        lappend group_names $group_name
      }
    }
    
    # all groups to which DSload "Compliant Base" is applied
    set xp "./container\[@n='_loads_'\]/condition\[@n='_compliant_base_'\]/group"
    set groups [$istep selectNodes $xp]
    foreach group $groups {
      set group_name [get_domnode_attribute $group n]
      if {[lsearch -exact $group_names $group_name] == -1} {
        lappend group_names $group_name
      }
    }
  }
    
  # all groups to which contacts have been assigned
  set xp [format_xpath {container[@n="_contact1_"]/blockdata}]
  set list [$root selectNodes $xp]
  foreach icon $list {
    # NEW XPATH to reach groupList -> group
    set xp [format_xpath {condition[@n='_contact_pair_']//group}]
    set groups [$icon selectNodes $xp]
    foreach group $groups {
      set group_name [get_domnode_attribute $group n]
      if {[lsearch -exact $group_names $group_name] == -1} {
        lappend group_names $group_name
      }
    }
  }

  foreach name $group_names {
    customlib::WriteString ""
    customlib::WriteString "*Surface, name=$name, type = node"
    customlib::WriteString "$name"
  }
  
  #
  # Close file
  #
  
  customlib::EndWriteFile
}