#
# Library functions for the scripting library (scrib)
#
# File: scrib/lib/offset.tcl
#
# $Revision$ on $Date$
# 
# Copyright (c) 1993-2003 ICEM CFD Engineering
#
#--------------------------------------------------------------------
#
# --> if the name of the first proc already exists in interp, return
#
if {"" != [info commands scrib_lib_offset_mesh_byfile]} { return }
#
# --> procs in this file depend on the following scrib/lib/ files--source them
#     I guess we could do this by proc, too, not sure waht's best
#
if {"" != [info commands scrib_lib_setup_source_files]} {
    set required_files "catia geo exec utils"
    scrib_lib_setup_source_files $required_files
}

## needs doc

proc scrib_lib_offset_assemble_mesh_byfile {file exec_args} {
    global proc_scrib_data
    #
    # --> set offset file name
    #     if that file already exists, then continue. it must have been created.
    #
    set proc_scrib_data(offset_mesh) [file root $file]_offset.uns
    if [file exists $proc_scrib_data(offset_mesh)] {
	return 1
    }
    #
    # --> run icemcfd in batch with the scripts/mesh_offset.tcl script
    #
    set script "mesh_offset.tcl"
    set proc_scrib_data(batch_script) $script
    #
    # --> this list of values is given to the batch mode run of med
    #
    #     for a standard scrib/script/*.tcl file, use this proc
    #     to get the args for the call to scrib_lib_exec_medbatch
    #
    set proc_scrib_data(model) $file
    set job_token [scrib_lib_exec_get_unique_job_token [file tail [file root $file]]]
    set med_args ""
    set stat [scrib_lib_exec_medbatch $job_token "$script" $med_args $exec_args]
    if {0 > $stat} {
	;
    }
    
    return $stat
}

## get matno from file with matname

proc scrib_lib_offset_material_mapping_read {matname {filename ""}} {
    
    # 1 is the default material id number
    set matno 1
    
    if {"" == $filename} {
	set wdir [scrib_lib_glue_get_working_directory]
	set filename [file join $wdir material_name_id_mapping.txt]
    }
    if ![file exists $filename] {
	mess "cannot read material name: $filename does not exist...using default of $matno\n" red
	return $matno
    }
    set fp [open $filename r]
    while {![eof $fp]}  {
	set line [string trim [gets $fp]]
	if {$matname == [lindex $line 1]} {
	    set matno [lindex $line 0]
	    break
	}
    }
    close $fp
    
    return $matno
}

## needs doc.

proc scrib_lib_offset_material_mapping_filename {working_dir} {
    
    if {"" == $working_dir} {
	set working_dir [scrib_lib_glue_get_working_directory]
    }
    return [file join $working_dir "material_name_id_mapping.txt"]
}

## needs doc.

proc scrib_lib_offset_material_mapping_write {inputfilename {outputfilename ""}} {
    
    if {"" == $outputfilename} {
	set outputfilename [scrib_lib_offset_material_mapping_filename ""]
    }
    
    if [file exists $outputfilename] {
	return
    }
    set matno 1
    set fp [open $inputfilename r]
    while {![eof $fp]} {
	set line [string trim [gets $fp]]
	if {"#" == [string range $line 0 0]} continue
	set part_num [lindex $line 0]
	if {5 <= [llength $line]} {
	    set matname [lindex $line 2]
	    if {("" != $matname) && (![info exists mat_map_data($matname)])} {
		set mat_map_data($matname) $matno
		incr matno
	    }
	}
    }
    close $fp
    
    foreach el [array names mat_map_data] {
	lappend map_list "$mat_map_data($el) $el"
    }
    set sorted_map_list [lsort -integer -index 0 $map_list]
    
    set fp [open $outputfilename w]
    puts $fp "# Material mapping file\n# ID\tName"
    foreach el $sorted_map_list {
	puts $fp "[lindex $el 0]\t[lindex $el 1]"
    }
    close $fp
}


## This procedure finds curves in layers coming from CAD system representing X-section, and uses curve to offset geometry

proc scrib_lib_offset_find_offset_curve {{fam MESH} {curves ""}} {
    
    #global rad
    set rad [scrib_lib_math_get_val rad]
    
    set min_length 5.0
    set tol .02
    
    if {$curves == ""} {set curves [ic_geo_get_objects curve ""]}
    if {[llength $curves] == 0} {
	mess "no curves given to scrib_lib_offset_do_find_curve_dir\n"
	return 
    }
    foreach curve $curves {
	set pnts [ic_geo_crv_data $curve cps]
	set ends($curve) [list [lindex $pnts 0] [lindex $pnts end]]
	if {[ic_geo_crv_length $curve] < $min_length} {
	    lappend short_curves $curve
	}  
    }
    
    if ![info exists short_curves] return
  
#mess "short_curves are $short_curves\n"  
    set pairs ""
    foreach s_curve $short_curves {
#mess "s_curve is $s_curve\n"
	set p(1) [lindex $ends($s_curve) 0]
	set p(2) [lindex $ends($s_curve) 1]
	set OK 0
	set OK_1 0
	foreach curve $curves {
	    if {[lsearch $short_curves $curve] != -1} continue
	    set pp1 [lindex $ends($curve) 0]
	    set pp2 [lindex $ends($curve) 1]
	    if {([scrib_lib_math_distance $p(1) $pp1] <= $tol) || \
		([scrib_lib_math_distance $p(1) $pp2] <= $tol)} {
		incr OK
		set near_c($OK) $curve
		set near_p($OK) 1
	    } elseif {([scrib_lib_math_distance $p(2) $pp1] <= $tol) || \
		([scrib_lib_math_distance $p(2) $pp2] <= $tol)} {
		incr OK
		set near_c($OK) $curve
		set near_p($OK) 2
	    }
	}
	if {$OK == 1} { 
	    set corner $p($near_p(1))
	    set end $p([expr 3-$near_p(1)])
	    lappend pairs [list $s_curve $near_c(1) $corner $end] 
	}
	if {$OK != 2} continue
	
	set v [scrib_lib_math_vector_from_points $p(1) $p(2)]
	set OK_1 0
	foreach i {1 2} {
	    set curve $near_c($i)
	    set param_a [expr 2.0/[ic_geo_crv_length $curve]]
	    set param_b [expr 1-$param_a]
	    set pnt [ic_geo_cre_pnt_on_crv_at_par $fam $fam/0 0 $curve]
	    set p_loc [ic_geo_get_point_location $pnt]
	    if {[scrib_lib_math_distance $p_loc $p($near_p($i))] <= $tol} {
		set p_near [ic_geo_cre_pnt_on_crv_at_par $fam $fam/0 $param_a $curve]
	    } else {
		set p_near [ic_geo_cre_pnt_on_crv_at_par $fam $fam/0 $param_b $curve]
	    }
	    set p_near_loc [ic_geo_get_point_location $p_near]
	    set v1($i) [scrib_lib_math_vector_from_points $p($near_p($i)) $p_near_loc]
	    set angle [expr [scrib_lib_math_angle_between_vectors $v $v1($i)]/$rad]
	    if {($angle > 105) || ($angle < 75)} break
	    incr OK_1
	}
	if {$OK_1 != 2} continue
	if { [expr [scrib_lib_math_angle_between_vectors $v1(1) $v1(2)]/$rad] > 30 } continue
	
	mess "Side curve recognized as $s_curve \n" blue
	set dir [scrib_lib_offset_find_orient_surface_normal $p(1) $p(2) 0 ] 

	if {[llength $dir] == 1} {     
	    set OK 0
	} else {
	    set thickness [ic_geo_crv_length $s_curve]	    
	    break
	}
    }
    if {[ic_geo_get_objects point $fam] != ""} {ic_delete_geometry point families $fam}

    if ![info exists thickness] {set OK 0}
   
    if !$OK {
	mess "Side curve not found \n" red
	return ""
    } else {
	return [list $thickness $dir]
    }
}

## This procedure was designed to find curves in layers coming from CAD system representing X-section, and finding a curve to use to offset mesh
## This is now replaced with scrib_lib_offset_find_curve which does the same at geometry level

proc scrib_lib_offset_do_find_curve_dir {{fam MESH}} {
    
    #global rad
    set rad [scrib_lib_math_get_val rad]
    
    set min_length 5.0
    set tol .02
    
    set curves [ic_geo_get_objects curve ""]
    foreach curve $curves {
	set pnts [ic_geo_crv_data $curve cps]
	set ends($curve) [list [lindex $pnts 0] [lindex $pnts end]]
	if {[ic_geo_crv_length $curve] < $min_length} {
	    lappend short_curves $curve
	}  
    }
    
    if ![info exists short_curves] return
  
#mess "short_curves are $short_curves\n"  
    set pairs ""
    foreach s_curve $short_curves {
#mess "s_curve is $s_curve\n"
	set p(1) [lindex $ends($s_curve) 0]
	set p(2) [lindex $ends($s_curve) 1]
	set OK 0
	set OK_1 0
	foreach curve $curves {
	    if {[lsearch $short_curves $curve] != -1} continue
	    set pp1 [lindex $ends($curve) 0]
	    set pp2 [lindex $ends($curve) 1]
	    if {([scrib_lib_math_distance $p(1) $pp1] <= $tol) || \
		([scrib_lib_math_distance $p(1) $pp2] <= $tol)} {
		incr OK
		set near_c($OK) $curve
		set near_p($OK) 1
	    } elseif {([scrib_lib_math_distance $p(2) $pp1] <= $tol) || \
		([scrib_lib_math_distance $p(2) $pp2] <= $tol)} {
		incr OK
		set near_c($OK) $curve
		set near_p($OK) 2
	    }
	}
	if {$OK == 1} { 
	    set corner $p($near_p(1))
	    set end $p([expr 3-$near_p(1)])
	    lappend pairs [list $s_curve $near_c(1) $corner $end] 
	}
	if {$OK != 2} continue
	
	set v [scrib_lib_math_vector_from_points $p(1) $p(2)]
	set OK_1 0
	foreach i {1 2} {
	    set curve $near_c($i)
	    set param_a [expr 2.0/[ic_geo_crv_length $curve]]
	    set param_b [expr 1-$param_a]
	    set pnt [ic_geo_cre_pnt_on_crv_at_par $fam $fam/0 0 $curve]
	    set p_loc [ic_geo_get_point_location $pnt]
	    if {[scrib_lib_math_distance $p_loc $p($near_p($i))] <= $tol} {
		set p_near [ic_geo_cre_pnt_on_crv_at_par $fam $fam/0 $param_a $curve]
	    } else {
		set p_near [ic_geo_cre_pnt_on_crv_at_par $fam $fam/0 $param_b $curve]
	    }
	    set p_near_loc [ic_geo_get_point_location $p_near]
	    set v1($i) [scrib_lib_math_vector_from_points $p($near_p($i)) $p_near_loc]
	    set angle [expr [scrib_lib_math_angle_between_vectors $v $v1($i)]/$rad]
	    if {($angle > 105) || ($angle < 75)} break
	    incr OK_1
	}
	if {$OK_1 != 2} continue
	if { [expr [scrib_lib_math_angle_between_vectors $v1(1) $v1(2)]/$rad] > 30 } continue
	
	mess "Side curve recognized as $s_curve \n" blue
	set dir [scrib_lib_offset_get_dir $p(1) $p(2) 0 $fam] 

	if {[llength $dir] == 1} {     
	    set OK 0
	} else {
	    set thickness [ic_geo_crv_length $s_curve]	    
	    break
	}
    }
    if {[ic_geo_get_objects point $fam] != ""} {ic_delete_geometry point families $fam}

    if ![info exists thickness] {set OK 0}
   
    if !$OK {
	mess "Side curve not found \n" red
	return ""
    } else {
	return [list $thickness $dir]
    }
}

proc scib_lib_offset_get_dir {p1 p2} {
mess $p1 
mess $p2
return 
}

## takes two points and finds the nearest element.  It then figures out and orients the mesh, based off nearest point to farther away point
## For example if p1 is closer it will use p1->p2 as the vector and orient the mesh towards p2

proc scrib_lib_offset_find_element_normal {p1 p2 {l 0} {fam MESH}} {
    
    #global rad
    set rad [scrib_lib_math_get_val rad]
    
    set tol 0.5
    #
    # put elements near position "p1" in "subset"
    #
    set subset [ic_uns_subset_create]
    set e_num [ic_uns_subset_add_near_pos $subset $p1]
    if !$e_num {
	mess "No elements near selected curve \r " red
	ic_uns_subset_delete $subset
	return 0
    }
    #
    #
    #
    set elements [ic_uns_subset_list_elements $subset]
    set vertexes [ic_uns_get_element_vertex_positions $subset 1]
    foreach el $elements {
	set type [lindex $el 0]
	set verts [lindex $el 2]
	if {($type != "QUAD_4") && ($type != "TRI_3")} {
	    set n_verts [llength $verts]
	    set vertexes [lreplace $vertexes 0 [expr $n_verts-1]]
	    continue
	}
	if {$type == "QUAD_4"} {
	    for {set i 0} {$i < 4} {incr i} {
		set ep($i) [lindex $vertexes $i]
	    }
	    set OK 0
	    foreach {ep1 ep2} {0 1  0 2  0 3  1 2  1 3  2 3} { 
		if [scrib_lib_math_is_point_on_segment $ep($ep1) $ep($ep2) $p1 $tol] {
		    set OK 1
		    break
		}
	    }
	    if !$OK {
		foreach {ep1 ep2} {0 1  0 2  0 3  1 2  1 3  2 3} { 
		    if [scrib_lib_math_is_point_on_segment $ep($ep1) $ep($ep2) $p2 $tol] {
			set OK 2
			break
		    }
		}
	    }
	    if $OK {
		set corners [list "$ep(0)" "$ep(1)" "$ep(2)" "$ep(3)"]
		break
	    } else {
		set vertexes [lreplace $vertexes 0 3]
	    }
	} else {
	    for {set i 0} {$i < 3} {incr i} {
		set ep($i) [lindex $vertexes $i]
	    }
	    set OK 0
	    foreach {ep1 ep2} {0 1  0 2  1 2} { 
		if [scrib_lib_math_is_point_on_segment $ep($ep1) $ep($ep2) $p1 $tol] {
		    set OK 1
		    break
		}
	    }
	    if !$OK {
		foreach {ep1 ep2} {0 1  0 2 1 2} { 
		    if [scrib_lib_math_is_point_on_segment $ep($ep1) $ep($ep2) $p2 $tol] {
			set OK 2
			break
		    }
		}
	    }
	    if $OK {
		set corners [list "$ep(0)" "$ep(1)" "$ep(2)"]
		break
	    } else {
		set vertexes [lreplace $vertexes 0 2]
	    }
	}
    }
    if ![info exists OK] { set OK 0 }
    if !$OK {
	mess "Can't find appropriate mesh element \n" red
	ic_uns_subset_delete $subset
	return 0
    }
    
    if {$OK == 1} {
	set v [scrib_lib_math_vector_from_points $p1 $p2]
    } else {
	set v [scrib_lib_math_vector_from_points $p2 $p1]
    }
    set x [scrib_lib_math_vector_nrm [scrib_lib_math_vector_from_points $ep(0) $ep(1)]]
    set z [scrib_lib_math_vector_vmult $x [scrib_lib_math_vector_from_points $ep(0) $ep(2)]]
    set a [expr [scrib_lib_math_angle_between_vectors $v $z]/$rad]
    ic_uns_subset_delete $subset
  
    set subset [ic_uns_subset_create]
    ic_uns_subset_add_pick_single element $subset [eval scrib_lib_math_center_point $corners] $z all
    if {$a > 90} {ic_uns_reorient $subset flip}
    ic_uns_reorient all consistent $subset
    ic_uns_subset_delete $subset

#   mess "v is $v;  x is $x;  a is $a\n"   

    return $v
}


## takes two points and finds the nearest surface.  It then figures out and orients the surfaces, based off nearest point to farther away point
## For example if p1 is closer it will use p1->p2 as the vector and orient the mesh towards p2

proc scrib_lib_offset_find_orient_surface_normal {p1 p2 {l 0}} {
    
    #global rad
    set rad [scrib_lib_math_get_val rad]
    
    set tol 0.5
    set l1 [ic_geo_get_point_location $p1]
    set l2 [ic_geo_get_point_location $p2]
    #
    # finds surface nearest position "p1" 
    #
    set surface [scrib_lib_geo_nearest_surface $l1]
    if {$surface == ""} {
	mess "No surface near selected curve \r " red
	return 0
    }
    #
    #Find closest point p1 or p2, then use normal to define whether to use current normal or other.
    #
    set d1 [scrib_lib_geo_distance_from_surfaces "" $l1 $surface]
    set d2  [scrib_lib_geo_distance_from_surfaces "" $l2 $surface]
mess "d1 is $d1; d2 is $d2\n"
    set par [ic_geo_find_nearest_srf_pnt $surface $p1]
    set u_par [lindex $par 0]
    set v_par [lindex $par 1]
    set surf_norm [lindex [ic_geo_get_srf_nrm $u_par $v_par $surface] 0]
mess "surface is $surface; surf_norm is $surf_norm\n"
    if {$d1 < $d2} {
	set v [scrib_lib_math_vector_from_points $p1 $p2]
    } else {
	set v [scrib_lib_math_vector_from_points $p2 $p1]
    }
mess "v is $v\n\n"
#    set x [scrib_lib_math_vector_nrm [scrib_lib_math_vector_from_points $ep(0) $ep(1)]]
#    set z [scrib_lib_math_vector_vmult $x [scrib_lib_math_vector_from_points $ep(0) $ep(2)]]
    set a [expr [scrib_lib_math_angle_between_vectors $v $surf_norm]/$rad]
    #For some reason surf_norm is reporting opposite of what is desired, so this is swapped.
    if {$a > 90} {set reverse 0} else {set reverse 1}
    ic_geo_surface_normals_orient $surface $reverse
    mess "a is $a\n"   

    return $v
}


## Offset the . This returns the thickness if success, -1 otherwise.

proc scrib_lib_offset_from_file {naspid offset_file} {
    
    if ![ic_uns_is_loaded] {
	mess "Error: cannot offset mesh from file without mesh loaded\n" red
	return -1
    }
    
    if {("" == $offset_file) || ![file exists $offset_file]} {
	mess "Error: cannot open offset file \"$offset_file\"\n" red
	return -1
    }
    
    set thickness  1    
    set offset_dir ""
    set mat_name   ""
    
    set fp [open $offset_file r]
    while {![eof $fp]} {
	set line [string trim [gets $fp]]	
	mess "fixme_debug: line $line\n" red
	
	set part_num [lindex $line 0]
	#
	# --> I have no f'n idea what the hell this part_num[1|2] is... --jka
	#	
	set part_num1 0
	set part_num2 0
	if {[string length $part_num] == 7} {
	    set part_num1 [string range $part_num 0 4]
	    set part_num2  [string replace $part_num 4 5]
	    mess "part_num is $part_num, part_num1 is $part_num1, part_num2 is $part_num2\n"
	}
	
	set std_part_num [scrib_lib_assembly_naspid_standardize $part_num]
	
	if { ($part_num == $naspid) || ($std_part_num == $naspid) \
		|| ($part_num1 == $naspid) || ($part_num2 == $naspid) } {
	    
	    mess "line is $line\n" green
	    
	    set mat_name   [lindex $line end-2]
	    set thickness  [lindex $line end-1]
	    set offset_dir [lindex $line end]
	    
	    if {[string first y $proc_scrib_data(offset_dir)] != -1 }  {
		#For examples 12346/5 axis needs to be switched. ?? --jka ?? fixme ??
		
		if {$part_num1 == $proc_scrib_data(naspid) } {
		    if {[string first - $proc_scrib_data(offset_dir)]  != -1}  {
			set proc_scrib_data(offset_dir) y+
		    } else {
			set proc_scrib_data(offset_dir) y-
		    }
		    mess "Warning: switching offset dir to $proc_scrib_data(offset_dir)\n" red
		}
	    }	    
	}
    }
    #
    # --> Offset the mesh
    #
    scrib_lib_offset_mesh $offset_dir $thickness
    #
    # --> Return the thickness
    #    
    return $thickness
}

## Offset the current mesh in direction with thickness (/2)

proc scrib_lib_offset_mesh {direction thickness} {
    
    if ![ic_uns_is_loaded] {
	mess "Error: cannot offset mesh from file without mesh loaded\n" red
	return 0
    }
    
    if {("" == $direction) || ("" == $thickness)} {
	mess "Error: cannot offset mesh in direction \"$direction\" with thickness \"$thickness\"\n" red
	return 0
    }
    
    set offset_distance [expr $thickness / 2.0]
    
    if [catch {ic_uns_offset_mesh map all dir $dir dist $offset_distance} res] {
	mess "Warning: there may have been an error in offsetting the mesh...\n" red
    } else {
	mess "Offset mesh with thickness of \"$thickness\" in direction \"$dir\"\n"
    }
    
    return $stat
}

##

proc scrib_lib_offset_param_material_property_add {param_file mat_id} {
    
    if ![file exists $param_file] {
	mess "Error: cannot find params file \"$param_file\"...\n" red
	return 0
    }
    
    set fp [open $param_file a]
    if {"" == $fp} {
	mess "Error: cannot open params file \"$param_file\" for appending...\n" red	
	return 0
    }
    
    set stat [scrib_lib_offset_param_material_properties_add $fp $mat_id]
    
    close $fp    
    
    return $stat
}

##
# this would use ic_param_*, but that code is f'd.

proc scrib_lib_offset_param_material_property_channel_add {fp mat_id} {
    
    if {"" == [file channels $fp]} {
	mess "Error: bad file channel for adding material property\n" red
	return 0
    }
    
    set mat_name "MAT${mat_id}"
    
    puts $fp "\n"
    puts $fp "begin Isotropic_Material $mat_name"
    puts $fp "    Mass_Density 7.8e-9"
    puts $fp "    Material_ID $mat_id"
    puts $fp "    Poissons_Ratio 0.28"
    puts $fp "    Reference_Temperature \{\}"
    puts $fp "    Shear_Modulus \{\}"
    puts $fp "    Structural_Damping_Coefficient \{\}"
    puts $fp "    Thermal_Expansion_Coefficient \{\}"
    puts $fp "    Youngs_Modulus 2e+5"
    puts $fp "end Isotropic_Material MAT1"
    
    return 1
}

## 

proc scrib_lib_offset_param_shell_properties {param_file material_mapping_file} {

    if ![file exists $param_file] {
	mess "Error: cannot find params file \"$param_file\"...\n" red
	return 0
    }    
    if ![file exists $material_mapping_file] {
	mess "Error: cannot find material mapping file \"$material_mapping_file\"...\n" red
	return 0
    }
    
    # get the array of material types to material id
    
    set fp [open $material_mapping_file r]
    while {![eof $fp]} {
	set line [string trim [gets $fp]]
	if {"#" == [string range $line 0 0]} { continue }
	
	set matno   [lindex $line 0]
	set matname [lindex $line 1]
	
	if {("" != $matno) && ("" != $matname)} {
	    set material_names($matno) $matname
	}    
    }
    close $fp
    unset fp

    # foreach material_id create a params file entry for it

    set fp [open $param_file a]
    if {"" == $fp} {
	mess "Error: cannot open params file \"$param_file\" for appending...\n" red	
	return 0
    }
    
    foreach matno [array names material_names] {
	scrib_lib_offset_param_material_property_channel_add $fp $matno
    }
    
    close $fp
    
    return 1
}

## needs doc

proc scrib_lib_offset_mesh_bydir {} {
    global proc_scrib_data
    #
    #    
    if {[info exists proc_scrib_data(offset_file)] && \
	    [string equal "" $proc_scrib_data(offset_file)]} {
	stdout "cannot offset by file if not specified!\n" red
    } else {
	if ![info exists proc_scrib_data(material_mapping_file)] {
	    set proc_scrib_data(material_mapping_file) ""
	}
	scrib_lib_offset_material_mapping_write $proc_scrib_data(offset_file) \
	    $proc_scrib_data(material_mapping_file)
    }
    #
    if ![info exists proc_scrib_data(wdir)] {
	set proc_scrib_data(wdir) .
    }
    #
    # -->  User might want to skip to a certain file
    #
     if {[info exists proc_scrib_data(skipto)] && $proc_scrib_data(cad_selmethod) == "skipto"} \
	{set skipto $proc_scrib_data(skipto)} else {set skipto ""}
    # 
    # foreach cad file in the directory, call this command
    #
    set command "scrib_lib_offset_assemble_mesh_byfile"
    
    set files [scrib_lib_exec_foreach_file $command domain $proc_scrib_data(wdir) $proc_scrib_data(wdir) $skipto]
    
    stdout "\nFinished offsetting files\n"
    foreach f $files { stdout "\t$f\n" }
}




