
# Return various pieces of information about the mesh.
#
# Copyright (c) 2001 ICEM CFD Engineering
# Author:  Wayne A. Christopher, wayne@pmac.com
#          Pacific Marketing and Consulting, Berkeley, California
#

# Returns the volume families near the given surface family.

proc mesh_volume_next_to_surface {family} {
    global domain_pid_family full_unstruct
    
    if ![info exists full_unstruct] { return "" }
    
    if ![info exists domain_pid_family] auto_bocos_setpids
    
    set groups [mesh_groups_in_family $family 2]
    set vols [mesh_families_near_groups $groups 3]
    return $vols
}

# Classify all the boundary surface families by which side they are on and
# whether they are touching each other or not.

proc mesh_boundary_surface_families {} {
    global domain_family_pid full_unstruct
    
    if ![info exists full_unstruct] { return "" }
    
    if ![info exists domain_family_pid] auto_bocos_setpids
    
    set conn [grid_connectivity]
    
    set groups ""
    foreach fam [array names domain_family_pid] {
	set gs [mesh_groups_in_family $fam 2]
	if {$gs != ""} {
	    foreach g $gs {
		if [$conn surface_emap_one_sided $g] {
		    lappend groups $g
		    set has_fam($g) $fam
		}
	    }
	}
    }
    set sets [eval $conn classify_groups_adjacent $groups]
    set fsets ""
    foreach set $sets {
	catch {unset got}
	foreach g $set {
	    set got($has_fam($g)) 1
	}
	lappend fsets [array names got]
    }
    return $fsets
}

# What is the minimum and maximum thickness of a volume element that 
# is adjacent to some element in the given surface family?

proc mesh_adjacent_volume_thickness {family} {
    set groups [mesh_groups_in_family $family 2]
    if { [catch {set conn [grid_connectivity]}] } {
    	mess "Mesh not loaded.\n" red
    	return 0
    }
    set min ""; set max ""
    foreach g $groups {
	set mm [$conn adjacent_volume_thickness $g]
	set min1 [lindex $mm 0]; set max1 [lindex $mm 1]
	if {($min == "") || ($min1 < $min)} { set min $min1 }
	if {($max == "") || ($max1 > $max)} { set max $max1 }
    }
    return [list $min $max]
}

# Returns a local coord system -- center, xvec, yvec, zvec -- for the given
# family.  zvec is the normal direction and xvec and yvec are in the plane.

proc mesh_local_coord_system_for_surface {family} {
    set stats [mesh_surface_family_stats $family]
    if {$stats == ""} { return "" }
    set minv [lindex $stats 0]
    set maxv [lindex $stats 1]
    set cent [lindex $stats 2]
    set zvec [lindex $stats 3]
    set diam [lindex $stats 4]
    set flat [lindex $stats 5]
    
    if {$minv == ""} { return "" }
    set diff ""
    set num_non_zero 0
    foreach i {0 1 2} {
	set v1 [lindex $minv $i]; set v2 [lindex $maxv $i]
	set d [expr $v2 - $v1]
	if $d { incr num_non_zero }
	lappend diff $d
    }
    if {$num_non_zero > 2} {
	return "[list $cent] [make_local_coord_system zvec $zvec minv $minv maxv $maxv]"
    } else {
	set d1 [lindex $diff [expr ($flat + 1) % 3]]
	set id1 [expr ($flat + 1) % 3]
	set d2 [lindex $diff [expr ($flat + 2) % 3]]
	set id2 [expr ($flat + 2) % 3]
	if {$d2>$d1} {
	    swap d1 d2
	    swap id1 id2
	}
	set xvec [lreplace "0 0 0" $id1 $id1 1]
	set yvec [lreplace "0 0 0" $id2 $id2 1]
	
	if {[vec_dot [vec_cross $xvec $yvec] $zvec] > 0} {
	    return [list $cent $xvec $yvec $zvec]
	} else {
	    return [list $cent $xvec [vec_minus $yvec] $zvec]
	}
    }
}

# Returns min xyz, max xyz, center, normal direction, diameter, flatdir

proc mesh_surface_family_stats {family {dim 2}} {
    set groups [mesh_groups_in_family $family $dim]
    if 1 { # used to be {$groups == ""}
	set ismesh 0
	set allpts ""
	set allnorms ""
	foreach {sh nums} [get_object_sides_in_family $family $dim] {
	    foreach ii [$sh list_sides 1] {
		set fn [lindex $ii 0]
		set box($fn) [$sh face_bbox [expr $fn - 1]]
		set norms($fn) [lindex [$sh face_dir [expr $fn - 1]] 1]
	    }
	    foreach n $nums {
		if !$n { set n 1 }
		set allpts [concat $allpts $box($n)]
		lappend allnorms $norms($n)
	    }
	}
	if {$allnorms == ""} { return "" }
	set bb [tdv_bbox_of_pts $allpts]
	set norm [lindex $allnorms 0]
	if {$norm == ""} { return "" }
	set norm [vec_minus [vec_normalize $norm]]
	
    } else {
	set ismesh 1
	# these are not images but they do have a get_bbox method ...
	set bb [tdv_bbox_of_images $groups]
    }
    set minv [lindex $bb 0]
    set maxv [lindex $bb 1]
    set cent ""
    set diff ""
    set mindiff ""
    set flatdir ""
    foreach i {0 1 2} {
	set v1 [lindex $minv $i]; set v2 [lindex $maxv $i]
	lappend cent [expr ($v1 + $v2) / 2.0]
	set d [expr $v2 - $v1]
	lappend diff $d
	if {($mindiff == "") || ($d < $mindiff)} {
	    set flatdir $i
	    set mindiff $d
	}
    }
    if $ismesh {
	set conn [grid_connectivity]
	set normdir [lindex [$conn adjacent_volume_thickness [lindex $groups 0]] 2]
	if { $normdir == "" } { 
	    set normdir 1
	    mess "Volume/Internal wall family selected: $family\n"
	}
	set norm [lreplace "0 0 0" $flatdir $flatdir [expr -$normdir]]
	set norm [vec_normalize $norm]
    }

    set d1 [lindex $diff [expr ($flatdir + 1) % 3]]
    set d2 [lindex $diff [expr ($flatdir + 2) % 3]]
    set diam [expr {($d1 > $d2) ? $d1 : $d2}]
    
    return [list $minv $maxv $cent $norm $diam $flatdir]
}

proc mesh_wire_family_diam {fam} {
    set shapes ""
    foreach {sh n} [get_object_sides_in_family $fam 1] { lappend shapes $sh }
    if {$shapes == ""} { return 0 }
    
    set bb [tdv_bbox_of_images $shapes]
    set dist 0
    foreach p1 [lindex $bb 0] p2 [lindex $bb 1] {
	set dist [expr $dist + ($p2 - $p1) * ($p2 - $p1)]
    }
    return [expr sqrt($dist)]
}

proc mesh_display_local_coord_system {family on cent xvec yvec zvec {dim 2}} {
    global tdv_default
    global mesh_local_coord_system
    
    if [info exists mesh_local_coord_system] {
	foreach im $mesh_local_coord_system {
	    $im delete
	}
	unset mesh_local_coord_system
    }
    if !$on return
    
    set mesh_local_coord_system ""
    set dot [tdv_point_im -point $cent -color yellow -size 3]
    $tdv_default add_image $dot
    lappend mesh_local_coord_system $dot
    
    switch $dim {
	2 {
	    set fam_data [mesh_surface_family_stats $family $dim]
	    if { $fam_data == "" } {
		mess "Warning: Local coordinate system could not be displayed for $family.\n" red 
		mess "Please load mesh...\n" red
		return
	    }
	    set diam [expr [lindex $fam_data 4] * .25]
	}
	1 { set diam [expr [mesh_wire_family_diam $family] * .25] }
	default { set diam 1 }
    }
    
    foreach {x color} {x red y green z blue} {
	set vec [set ${x}vec]
	set pt2 ""
	foreach i {0 1 2} {
	    lappend pt2 [expr [lindex $cent $i] + [lindex $vec $i] * $diam]
	}
	set line [tdv_line_im -point1 $cent -point2 $pt2 \
		      -color $color -width 2]
	$tdv_default add_image $line
	lappend mesh_local_coord_system $line
    }
}

proc mesh_show_pml_box {on xmin xmax ymin ymax zmin zmax} {
    global tdv_default
    global domain mesh_pml_box
    
    if [info exists mesh_pml_box] {
	$mesh_pml_box delete
	unset mesh_pml_box
    } 
    if !$on return
    
    set bb [[$domain getval body_shape] get_bbox]
    set xmin [expr [lindex [lindex $bb 0] 0] - $xmin]
    set xmax [expr [lindex [lindex $bb 1] 0] + $xmax]
    set ymin [expr [lindex [lindex $bb 0] 1] - $ymin]
    set ymax [expr [lindex [lindex $bb 1] 1] + $ymax]
    set zmin [expr [lindex [lindex $bb 0] 2] - $zmin]
    set zmax [expr [lindex [lindex $bb 1] 2] + $zmax]
    
    set mesh_pml_box [tdv_prism_geo -point1 [list $xmin $ymin $zmin] \
			  -point2 [list $xmax $ymax $zmax] -linestyle 0xf0]
    $tdv_default add_image $mesh_pml_box
}

proc mesh_show_pml_single {on xmin xmax ymin ymax zmin zmax} {
    global tdv_default
    global domain mesh_pml_single
    
    if [info exists mesh_pml_single] {
	$mesh_pml_single delete
	unset mesh_pml_single
    } 
    if !$on return
    
    set mesh_pml_single [tdv_prism_geo -point1 [list $xmin $ymin $zmin] \
			  -point2 [list $xmax $ymax $zmax] -linestyle 0xf0]
    $tdv_default add_image $mesh_pml_single
}

# ---------------------------------

proc shapes_in_family {fam} {
    set ret ""
    foreach obj [db_list_objects] {
	foreach vv [$obj list_values] {
	    if {[string match fam_* [lindex $vv 0]] &&
		([lindex $vv 1] == $fam)} {
		set sh [$obj getval body_shape]
		if {$sh != ""} {
		    lappend ret $sh
		}
	    }
	}
    }
    return $ret
}

proc shapes_and_sides_in_family {fam} {
    set ret ""
    foreach obj [db_list_objects] {
	foreach vv [$obj list_values] {
	    if {[string match fam_* [lindex $vv 0]] &&
		([lindex $vv 1] == $fam)} {
		set sh [$obj getval body_shape]
		if {$sh != ""} {
		    lappend ret $sh [string range [lindex $vv 0] 4 end]
		}
	    }
	}
    }
    return $ret
}

proc family_is_on_only_2d_objects {fam} {
    foreach sh [shapes_in_family $fam] {
	if [$sh get_volume] { return 0 }
    }
    return 1
}

proc mesh_groups_in_family {family dim} {
    global grid_mapping reduced_group
    
    set ret ""
    foreach dom [array names reduced_group] {
	set obj [lindex $grid_mapping([lindex $dom 0]) 0]

	if {$obj == ""} continue

	set sh [$obj getval body_shape]
	if {$sh == ""} continue
	
	set map $reduced_group($dom)
	if ![$map count] continue

	set side [lindex $dom 1]
	
	if {!$side && ([$obj si_value type] == "fluid")} continue
	if {($dim == 3) && $side} continue
	if {($dim == 2) && !$side} continue
	
	set fams [family_for_side $obj $sh $side]
	if {[lsearch $fams $family] < 0} continue
	
	lappend ret $reduced_group($dom)
    }
    return $ret
}

proc mesh_families_near_groups {groups dim} {
    global full_unstruct domain_pid_family
    
    if ![info exists full_unstruct] { return "" }
    
    set conn [grid_connectivity]
    
    set map [element_map -mesh $full_unstruct]
    foreach eset [$full_unstruct list_sets] {
	if {$dim == 2} {
	    if {[lsearch {QUAD_4 TRI_3} [$eset get -type]] < 0} continue
	} else {
	    if {[lsearch {TETRA_4 PYRA_5 PENTA_6 HEXA_8} [$eset get -type]] < 0} continue
	}
	
	foreach group $groups {
	    $conn near_emap $group $map $eset
	}
    }
    
    set fams ""
    foreach p [$map list_pids] {
	lappend fams $domain_pid_family($p)
    }
    
    $map delete
    return $fams
}

# This one works without any mesh being loaded.  It uses the first segment of
# the first non-point wire in the family.

proc wire_local_coord_system {fam} {
    foreach obj [db_list_objects pipe] {
	if {[$obj getval fam_edge] != $fam} continue
	set sh [$obj getval body_shape]
	if {[$sh getval nverts] < 2} continue
	set p1 [$sh getval vert1]
	set p2 [$sh getval vert2]
	set dir ""
	foreach i {0 1 2} {
	    lappend dir [expr [lindex $p2 $i] - [lindex $p1 $i]]
	}
	return "[list $p1] [make_local_coord_system xvec $dir]"
    }
    return ""
}

# This makes sure that the bars in the domain file are ordered according
# to the wire and cad-curve elements that exist.

proc order_bars_according_to_wires {} {
    global domain_family_pid grid_type order_bars_according_to
    
    if {$grid_type == "hexa"} return ;# pipe meshes don't need this
    
    if ![info exists order_bars_according_to] return
    
    foreach type $order_bars_according_to {
	foreach bar [db_list_objects $type] {
	    set sh [$bar getval body_shape]
	    if {([$sh get -shtype] == "pipe") && \
		    ([$sh getval nverts] < 2)} continue
	    if {[$bar getval fam_edge ""] == ""} continue
	    
	    lappend in_fam([$bar getval fam_edge]) $sh
	}
    }
    
    foreach fam [array names in_fam] {
	set conn [grid_connectivity]
	set pid $domain_family_pid([string toupper $fam])
	set wires ""
	foreach w $in_fam($fam) {
	    if {[$w get -shtype] == "pipe"} {
		set nv [$w getval nverts]
		set wire ""
		for {set i 1} {$i <= $nv} {incr i} {
		    lappend wire [$w si_value vert$i]
		}
	    } else {
		set wire ""
		foreach bsp [$w list_bsplines] {
		    if {[$bsp get -type] == "curve"} {
			lappend wire [$bsp get -number]
		    }
		}
	    }
	    lappend wires $wire
	}
	puts "$conn order_bars_by_wires $wires $pid"
	set ret [$conn order_bars_by_wires $wires $pid]
	if {[lindex $ret 0] != [lindex $ret 1]} {
	    user_error "expected to order [lindex $ret 1] bars but got [lindex $ret 0]"
	}
	grid_invalidate_connectivity
    }
}

# Is a given surface family one-sided or two-sided?  If the surface is mixed
# then it will return -1 and if there is no surface in this family it returns
# ""

proc num_sides_for_surface_family {family} {
    set conn [grid_connectivity]
    set groups [mesh_groups_in_family $family 2]
    set sided ""
    foreach g $groups {
	if [$conn surface_emap_one_sided $g] {
	    if {$sided == ""} { set sided 1 }
	    if {$sided == 2} { set sided -1 }
	} else {
	    if {$sided == ""} { set sided 2 }
	    if {$sided == 1} { set sided -1 }
	}
    }
    return $sided
}

# ------------------------------------------

proc families_with_bocos {bctypes} {
    global famboco_data
    set ret ""
    foreach f [array names famboco_data] {
	foreach bc $famboco_data($f) {
	    if [lmember $bctypes [lindex $bc 1]] {
		lappend ret $f
		break
	    }
	}
    }
    return $ret
}

proc count_nodes_with_bocos {bctypes {not_bctypes ""} {also_families ""}} {
    global full_unstruct
    
    set fams [concat [families_with_bocos $bctypes] $also_families]
    foreach f $fams {
	foreach g [mesh_groups_in_family $f -1] {
	    set found($g) 1
	}
    }
    set groups [array names found]
    $full_unstruct reset_flood_flags
    foreach g $groups {
	$g set_all_flags 1
    }
    
    catch {unset found}
    set not_fams [families_with_bocos $not_bctypes]
    foreach f $not_fams {
	foreach g [mesh_groups_in_family $f -1] {
	    set found($g) 1
	}
    }
    set groups [array names found]
    foreach g $groups {
	$g set_all_flags 0
    }
    return [$full_unstruct num_nodes_on_flagged_element 1]
}

proc mesh_area_of_shape {sh {only_side ""}} {
    global full_unstruct grid_mapping reduced_group

    if ![info exists full_unstruct] {
	error "No mesh exists"
    }
    if ![info exists grid_mapping($sh)] {
	error "No mesh exists for $sh"
    }
    set mnum $grid_mapping($sh)

    set area 0
    foreach side [db_shape_sidenums $sh] {
	if {($only_side != "") && ($side != $only_side)} continue
	
	set varname "reduced_group($mnum $side)"
	if [info exists $varname] {
	    set map [set $varname]
	    set area [expr $area + [$map get_area]]
	}
    }
    return $area
}
