# Library functions for the scripting library (scrib)
#
# File: scrib/lib/quad.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_quad_mesh]} { 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 0 {
    # can't do this because tetra includes quad includes tetra...
    if {"" != [info commands scrib_lib_setup_source_files]} {
	set required_files "tetra utils quality"
	scrib_lib_setup_source_files $required_files
    }
}

## needs doc. xxx

proc scrib_lib_quad_doit {{extra ""} {try_twice 1} } {
    global proc_scrib_data
    #ic_set_meshing_params global 0 gnat $proc_scrib_data(ignore_tol)
    if ![info exists proc_scrib_data(ignore_tol) ] {set proc_scrib_data(ignore_tol) 1}
    if ![info exists proc_scrib_data(mapped_mesh) ] {set proc_scrib_data(mapped_mesh) 0}
    if ![info exists proc_scrib_data(number_rings) ] {set proc_scrib_data(number_rings) 0}

    set cmdargs "element 2 iters 25 proj 1 improvement 2 \
	    ele_tol $proc_scrib_data(ignore_tol) geo_tol [expr $proc_scrib_data(ignore_tol) * 0.25] \
	    time_max 720 error_family ERROR try_twice $try_twice \
	    adjust_nodes 1 "
    
    if {($proc_scrib_data(mapped_mesh) != 0) && ([lsearch $extra block] == -1)} {
	append cmdargs "block $proc_scrib_data(mapped_mesh) "
    }
    
    if $proc_scrib_data(number_rings)  {
	#if {[lsearch $extra bunch ] == -1} {
	append cmdargs " curve_family {} offset [expr $proc_scrib_data(number_rings) * -1] \
		offset_h $proc_scrib_data(ratio_rings)"
	#}
    }
    if {"" != $extra} {
	append cmdargs " $extra"
    }    
    
    # this is necessary to get the pt-loops (shadow-weld-pts) to work --jka
    scrib_lib_geo_unload_reload
    
#ic_save_tetin "pre_quad.tin"
#puts "ic_quad $cmdargs\n"
    ic_uns_create_empty
    ic_uns_set_modified 0
    if [catch {eval ic_quad $cmdargs} res] {
	mess "$res\n"
    }
    
    set failed_loops [ic_geo_get_objects loop ERROR]
    if {$failed_loops != ""} {
    	mess "failed loops are $failed_loops\n" blue
    }
    return $failed_loops
}

## loop per each surface

proc scrib_lib_quad_create_loops_for_each_surface {{fam MESH} {surfs ""} {points ""}} {
    
    if {$surfs == ""} {
	set surfs [ic_geo_get_objects surface ""]
    }
    if {"" == $surfs} {
	mess "got no surfs from geomobj in scrib_lib_quad_mesh!\n" red
	return 1
    }
    
    set return_loops ""
    
    foreach surf $surfs {
	if [catch {ic_geo_boundary surface $surf 0 1} bounds] {
	    mess "Error: $bounds\n"
	    continue
	}
	if {"" == $bounds} {
	    return ""
	}
	#mess "for surface $surf got bounds: $bounds\n" green
    	set new_loop [ic_geo_create_loop $fam..$surf $fam names $bounds 0 $surf "" $points]
	lappend return_loops $new_loop
    }

    return $return_loops
}

##

proc scrib_lib_quad_create_loops_for_surface {surf {fam MESH} {points ""}} {
    
    if [catch {ic_geo_boundary surface $surf 0 1} bounds] {
	mess "Error: $bounds\n"
	return ""
    }
    if {"" == $bounds} {
	return ""
    }
    #mess "In creating loop for surface $surf got bounds: $bounds\n" green
    
    set new_loop [ic_geo_create_loop $fam..$surf $fam names $bounds 0 $surf "" $points]
    
    return $new_loop
}

## Loop per set of curves

proc scrib_lib_quad_create_loops_from_curves {{fam MESH} {curves ""} {surfs ""}} {
    if {"" == $curves} {
	mess "got no curves from geomobj in scrib_lib_quad_mesh!\n" red
	return 1
    }

    return [ic_geo_create_loop $fam..FLANGE $fam names $curves 0 $surfs]
}

## loop for collected regions of surfaces.  For now this is being used for beads.. BEN

proc scrib_lib_quad_create_loop_for_group {fam surfs {points ""} {corners ""} } {
 
    if {"" == $surfs} {
	mess "got no surfs from geomobj in scrib_lib_quad_create_loops_for_grouped_surfaces!\n" red
	return 1
    }
    
    set curves ""
    foreach surf $surfs { 
	#be careful of multiple edge surface case, especially surfaces that have only one boundary
	set new_curves [scrib_lib_utils_list_remove_duplicated [ic_geo_boundary surface $surf 0 1]]
	if {$new_curves != ""} { set curves [concat $curves $new_curves] }
	if ![info exists surflist] {
	    set surflist $surf..
	} else {
	    append surflist $surf..
	} 
    }
    #mess "$surflist\n"
    set bounds [scrib_lib_utils_list_remove_duplicated $curves]
    #set loop_name  [ic_geo_create_loop $fam..$surflist $fam names $bounds 0 $surfs]
    set loop_name  [ic_geo_create_loop "" $fam names $bounds 0 $surfs "" $points $corners]
    return $loop_name
}

## get surfaces from a list of loops where loops have names fam..surf1..surf2.. etc or fam..surf

proc scrib_lib_quad_find_surfs_from_loops { loop_list } {
    set surf_list ""
    foreach loop $loop_list {
	set surfs ""
	set loc [string first .. $loop]
	set fam [string range $loop 0 $loc ]
	set s [string range $loop [expr $loc + 2] end]

	while {[string first .. $s] != -1} {
	    set loc [string first .. $s]
	    lappend surfs [string range $s 0 [expr $loc -1]]
	    set s [string range $s [expr $loc + 2] end]
	}	
	if {$s != ""} {lappend surfs $s	 }   
	lappend surf_list $surfs
    }
    return $surf_list
}


## quad mesher. needs doc

proc scrib_lib_quad_mesh { {fam MESH} {surfs ""} } {
    global proc_scrib_data    
    
    set q_error 0
    if {$surfs == ""} {
	set surfs [ic_geo_get_objects surface ""]
    }
    if {"" == $surfs} {
	mess "got no surfs from geomobj in scrib_lib_quad_mesh!\n" red
	return 1
    }
      
    foreach surf $surfs {
	set bounds [ic_geo_boundary surface $surf 0 1]
	#mess "for surface $surf got bounds: $bounds\n" green
    	ic_geo_create_loop $fam..$surf $fam names $bounds 0 $surf
    }
    ic_uns_create_empty
    ic_uns_set_modified 0
    
    #mess "created loops of surfs $surfs: [ic_geo_get_objects loop ""]\n" green
    
    set q_error [scrib_lib_quad_doit]
    
    set vis_loops [ic_geo_get_objects loop ERROR]
    #the return from the above does not seem to work properly.
    #instead we have to get the failed loops as done below:    
    if {$vis_loops != ""} {
	mess "failed loops are $vis_loops\n"
	set q_error 1
    }
    if {$q_error == 1 } {
	if {"" == $vis_loops} {
	    mess "error in quad but vis_loops is null!\n" red
	} else {
	    if {$proc_scrib_data(failure_mode) == 1}  {
   		switch [ tk_dialog .fail_quad "Quad failed" \
	    		"Do you want to continue?" \
	    		"" -1 Yes No ] {
			0 { 
	    			# yes
			   }
			1 { return 1 }
    		}
	    }
	    mess "will check for folded surfaces for $vis_loops"
	    set q_error [scrib_lib_quad_mesh_try2 $vis_loops $fam $surfs]
	}
    }
    if !$q_error {
	catch {ic_change_tri_quad skew 160 warp 20 iters 0 \
		conver 0.25 improvement 0 quadrization 0}
	ic_uns_set_node_dimension all same
    }
    
    return $q_error
}

## try to segment folded surfaces to fix these failed surfaces. needs doc.
proc scrib_lib_quad_mesh_try2 {failed_loops {fam MESH} {surfs ""}} {
    global proc_scrib_data
    
    if {"" == $failed_loops} {
	mess "Warning: called scrib_lib_quad_mesh_try2 with no failed_loops!\n" red
	# return 1
    }
    #set surfs ""
    foreach loop $failed_loops {
	
	if {[ic_geo_check_objects_exist loop $loop] == "" } { continue }
	ic_unload_mesh
	set loc [string first .. $loop]
	set surf [string range $loop [expr $loc + 2] end]
	lappend surfs $surf
    }
    ic_delete_geometry loop all ""
    ic_geo_trm_srfs_by_curvature $surfs 120
    ic_build_topo $proc_scrib_data(build_topo_tol)
    ic_geo_delete_unattached

    set q_error 0
    if {$surfs == ""} {
	set surfs [ic_geo_get_objects surface ""]
    }
    if {"" == $surfs} {
	mess "got no surfs from geomobj in scrib_lib_quad_mesh!\n" red
	return 1
    }
    foreach surf $surfs {
	set bounds [ic_geo_boundary surface $surf 0 1]
	#mess "for surface $surf got bounds: $bounds\n" green
    	ic_geo_create_loop $fam..$surf $fam names $bounds 0 $surf
    }
    ic_uns_create_empty
    ic_uns_set_modified 0
    
    set q_error [scrib_lib_quad_doit]
    
    set vis_loops [ic_geo_get_objects loop ERROR]
    #the return from the above does not seem to work properly.
    #instead we have to get the failed loops as done below:    
    if {$vis_loops != ""} {
	mess "failed loops are $vis_loops\n"
	set q_error 1
    }
    if {$q_error == 1 } {
	if {"" == $vis_loops} {
	    mess "error in quad but vis_loops is null!\n" red
	} else {
	    if { $proc_scrib_data(failure_mode) == 1} {
   		switch [ tk_dialog .fail_quad2 "Quad failed second time" \
	    		"Do you want to continue?" \
	    		"" -1 Yes No ] {
			0 { 
	    			# yes
			   }
			1 { return 1 }
    		}
	    }
	    mess "will try to join bad loops with $vis_loops"
	    set q_error [scrib_lib_quad_mesh_try3 $vis_loops $fam $surfs]
	}
    }
    return $q_error    
}

## try to group loops to resolve. needs doc.
proc scrib_lib_quad_mesh_try3 {failed_loops {fam MESH} {surfs ""}} {
    global proc_scrib_data
    
    set failed 1
    set used_surfs ""
    
    if {"" == $failed_loops} {
	mess "Warning: called scrib_lib_quad_mesh_try3 with no failed_loops!\n" red
	# return 1
    }
    
    foreach loop $failed_loops {
	
	if {[ic_geo_check_objects_exist loop $loop] == "" } { continue }
	ic_unload_mesh
	set loc [string first .. $loop]
	set surf [string range $loop [expr $loc + 2] end]
	mess "loop $loop attached to $surf will be deleted\n"
	set curves [lsort -unique [ic_geo_boundary surface $surf 0 1]]
	
	#======================
	
	set outer_curves [lsort -unique [ic_geo_boundary surface $surf 1 1]]
	set inner_curves [scrib_lib_utils_list_remove_duplicate \
		[concat $curves $outer_curves]]
#======== INNER METHOD of loop adding
	if {$inner_curves != ""} {
	    set jsurfs $surf
	    #
	    # just add missing regions
	    #
	    mess "Using inner method. Will add all internal surfaces\n"
	    lappend used_surfs $surf
	    set used_curves $curves
	    set done 1
	    while $done {
		set new_surfs ""
		foreach c $inner_curves {
		    set isrf [ic_geo_incident curve $c 1]
		    foreach srf $isrf {
			if {[lsearch -exact $used_surfs $srf] != -1} {continue}	  	
		   	lappend new_surfs $srf
			lappend used_surfs $srf
			lappend jsurfs $srf
		    }
		}
		mess "new_surfs $new_surfs will be added\n"
		if {$new_surfs == ""} {
		    set done 0
		    break
		}
		set new_curves ""
		foreach ss $new_surfs {
		    set bcrvs [ic_geo_boundary surface $ss 0 1]
		    foreach bcrv $bcrvs {
			if {[lsearch -exact $used_curves $bcrv] != -1} {continue}	  	
			lappend new_curves $bcrv
			lappend used_curves $bcrv			
		    }
		}
		mess "new_curves are $new_curves\n"
		set inner_curves $new_curves
	    }
#============ BREAK CYLINDRICAL SURFACES ===========
	# elseif {[ic_surface_curvature $surf] > 90} {
	#    set new_surfs [ic_geo_trm_srfs_by_curvature $surf 90]
	#
	#}
#============ ADD OUTER LAYER OF SURFACES ==========
	} else {
	    #
	    # add one layer of surfaces around bad loop, and remesh these.
	    #
	    mess "Using outer method.  Will add one layer of surfaces.\n"
	    set ok_curves ""
	    #
	    # check attached curves to see if they are valid to use neighbor surface $isrf
	    #
	    foreach curve $curves {	
		set isrf [lsort -unique [ic_geo_incident curve $curve 1]]
		#(![check_normals $isrf 400 70])
		if {([llength $isrf] != 2)} {continue}
		set ind [lsearch -exact $isrf $surf]
		set neighbor_surf [lreplace $isrf $ind $ind]
		mess "failure surf $surf has \
			neighbor $neighbor_surf of curve $curve\n"
		if {[lsearch -exact $used_surfs $neighbor_surf] != -1} {
		    mess "neighbor $neighbor_surf is used.\n"
		    continue
		}
		lappend ok_curves [list $curve [ic_geo_crv_length $curve]] 
	    }
	    #
	    # ok_curves are now valid (ideally) to attach surface, 
	    # listed as curve_name/length pairs
	    #
	    if {$ok_curves == ""} { continue }
	    set ok_curves [lsort -index 1 -decreasing $ok_curves]
	    set jsurfs ""
	    foreach pair $ok_curves {
		set crv [lindex $pair 0]
		set jsurfs [concat [ic_geo_incident curve $crv 1] $jsurfs] 
	    }
	    set jsurfs [lsort -unique $jsurfs]
	}
	set loop_crvs ""
	foreach s $jsurfs {
	    foreach c [ic_geo_boundary surface $s 0 1] {
		lappend loop_crvs $c
	    }
	}
	set loop_crvs [scrib_lib_utils_list_remove_duplicate $loop_crvs]
	mess "loop will be created from curves $loop_crvs\n"
	ic_geo_create_loop $fam..failed$failed $fam names $loop_crvs 0 $jsurfs
	foreach srf $jsurfs {
	    ic_delete_geometry loop names $fam..$srf 0
	}
	set failure($fam..failed$failed) $jsurfs
	set used_surfs [lsort -unique [concat $used_surfs $jsurfs]]
  	incr failed
    }
    if 0 {
	set q_error 0
	if [ic_quad element 2 iters 25 proj 1 improvement 2 \
		ele_tol $proc_scrib_data(ignore_tol) geo_tol 0.05 time_max 720 \
		error_family ERROR try_twice -10 curve_family HOLE offset 2 offset_h 0.6] {
	    set q_error 1
	}
    } else {
	set q_error [scrib_lib_quad_doit]
	#the return from the above does not seem to work properly.
	#instead we have to get the failed loops as done below:
    	set vis_loops [ic_geo_get_objects loop ERROR]
    	if {$vis_loops != ""} {
	    mess "failed loops are $vis_loops\n"
	    set q_error 1
    	}
    }
    if {$q_error == 1 } {    	
	#Now we want to save the good mesh, and mesh the bad surfaces with tetra.
	#before doing this, we need to delete unattached bar elements.
	ic_save_unstruct partial_quad.uns
	set bad_loops [ic_geo_get_objects loop ERROR]
	if {"" == $bad_loops} {
	    mess "There was an error but no bad loops returned!\n" red	    
	    #return $q_error
	}
	
	mess "loops that failed on second try are $bad_loops\n"
	ic_delete_geometry loop names $bad_loops	
	
	set bad_surfs ""
	foreach l $bad_loops {
	    set loc [string first .. $l]
	    if {[string range $l [expr $loc + 2] [expr $loc + 7]] == "failed"} {
		foreach s $failure($l) {
		    lappend bad_surfs $s
		}
	    } else {
		lappend bad_surfs [string range $l [expr $loc + 2] end]
	    }
	}
	if {"" == $bad_surfs} {
	    mess "There are no surfs to mesh in Tetra!\n" red
	} else {
	    if  {$proc_scrib_data(failure_mode) == 1} {
   		switch [ tk_dialog .fail_quad3 "Quad failed third time" \
	    		"Do you want to continue?" \
	    		"" -1 Yes No ] {
			0 { 
	    			# yes
			   }
			1 { return 1 }
    		}
	    }
	    mess "Meshing surfaces in Tetra: $bad_surfs\n"
	    mess "need to mesh with partial tet/quad method"
	    set q_error [scrib_lib_tetra_partial_mesh $bad_surfs $fam]
	}
    }
    
    return $q_error
}

# ( this is just for debuggin the quad_clean proc, saving domains along the way. --jka )

proc scrib_lib_quad_clean_debug_save {} {
    global scrib_lib_quad_clean_debug_save_db_qmp
    
    #xx
    set debug_it 0
    
    if !$debug_it { return }
    
    if ![info exists scrib_lib_quad_clean_debug_save_db_qmp] {
	set scrib_lib_quad_clean_debug_save_db_qmp 1
    } else {
	incr scrib_lib_quad_clean_debug_save_db_qmp
    }
    
    set name quad_clean_debug${scrib_lib_quad_clean_debug_save_db_qmp}.uns
    ic_save_unstruct $name
}

## quad cleanup. needs doc

proc scrib_lib_quad_clean {{bolt_hole_bar_fam HOLE} {outer_ring_bars OUTER_RING} } {
    global proc_scrib_data
    
    if 0 {
	# first try to combine to quads    
	#   ic_change_tri_quad skew 160 warp 20 iters 6 conver 0.25 \
		#	    improvement 0 quadrization 0        
	# then merge nodes to try to get rid of slivers    
	#set merge_nodes [ic_uns_subset_create]
    }
    
    
    # --> D E L E T E   A L L   D U P L I C A T E   E L E M E N T S
    scrib_lib_quad_clean_debug_save
    
    scrib_lib_quality_duplicate_elements_delete
    
    
    # ======== PUT BOLT QUAD_4S INTO HOLE FAMILY ========
    set orig_fam 	[scrib_lib_assembly_famname_get $proc_scrib_data(naspid)]
    set boltquads [scrib_lib_quad_find_offset_quad_layers $proc_scrib_data(number_rings) $bolt_hole_bar_fam]
    #if {[lsearch [ic_uns_subset_list] $boltquads] != -1} { } this below is better --jka
    if {("" != $boltquads) && [ic_uns_subset_exists $boltquads]} {
    	if [ic_uns_subset_count_elements $boltquads] {
	    ic_uns_change_family $boltquads $bolt_hole_bar_fam
	}
    }
    
    scrib_lib_quad_clean_debug_save
    
    # ----> M E R G E   N O D E S   O F   F I  R S T   L A Y E R
    # ----> A R O U N D   R I N G S   O F   Q U A D S (bolt_hole_bar_fam)
   
    set tol $proc_scrib_data(ignore_tol)
    set stat [scrib_lib_quad_merge_ring_nodes $bolt_hole_bar_fam $tol]    
    
    
    # ----> the following sets families not to smooth.
    #       basically all QUADS that should not be edited
    #       are thrown into $bolt_hole_bar_fam            --BEN
    
    
    set fams_for_smoothing [ic_uns_list_families]
    foreach fam [list $bolt_hole_bar_fam ORFN] {
    	set ind [lsearch $fams_for_smoothing $fam]
    	set fams_for_smoothing [lreplace $fams_for_smoothing $ind $ind]
    }
    
    scrib_lib_quad_clean_debug_save
    
    # ----> C H E C K   A N D   F I X   O V E R L A P S
    
    scrib_lib_quad_remesh_by_diagnostic overlap $fams_for_smoothing {TRI_3 QUAD_4} 2
    

    # ----> D E L E T E { disconnected bars, hanging bar elements, disconnected vertices }
    
    scrib_lib_quad_clean_debug_save
    
    scrib_lib_quality_disconnected_bars_delete    
    scrib_lib_quality_hanging_bars_delete    
    ic_delete_disconnected_vertices
    
    scrib_lib_quad_clean_debug_save
    
    
    # ----> C R E A T E   B A R S   O N   O U T E R   R I N G
    
    
    # foreach quad in the return of subset_for_quad, create a bar element at top nodes 
    # top nodes are those not in $nodes
if 0 {   
    if [ic_uns_count_in_families $bolt_hole_bar_fam] {
    	foreach item [ic_uns_subset_list_elements $boltquads] {
	    if {[lindex $item 0] != "QUAD_4"} continue
	    set quad_num [lindex $item 1]
	    set top_nodes ""
	    foreach quad_node [lindex $item 2] {
	    	if {[lsearch -exact $nodes $quad_node] == -1} {
		    lappend top_nodes $quad_node
	    	}
	    }
    	    if {[llength $top_nodes] != 2} {
	    	mess "There are more or less than 2 nodes on outer boundary of quad $quad_num, \n" red
	    	mess "Not an O-grid, so not creating bar element\n" red
    	    } else {
	    	#set tmp [ic_uns_subset_create]
		#ic_uns_subset_add_node_numbers $tmp $top_nodes
		#mess "Bar element will be created between $top_nodes in family $outer_ring_bars\n" blue
		#if ![ic_uns_check_family $outer_ring_bars] {ic_uns_new_family $outer_ring_bars}
	    	#ic_uns_create_element BAR_2 $outer_ring_bars $tmp 1
	    }
	}
    }
}
    if {[info exists boltquads] && ("" != $boltquads)} {
	ic_uns_subset_delete $boltquads
    }
    scrib_lib_quad_clean_debug_save
    
    
    # ---> Reset fams_for_smoothing to include {bolt_hole_bar_fam , outer_ring_bars}
    
    
    set fams_for_smoothing [ic_uns_list_families]
    foreach fam [list $bolt_hole_bar_fam $outer_ring_bars ORFN] {
    	set ind [lsearch $fams_for_smoothing $fam]
    	set fams_for_smoothing [lreplace $fams_for_smoothing $ind $ind]
    }
    mess "Cleaning families: $fams_for_smoothing...\n"


    # ----> R E M O V E   B A R S   A T   G I V E N   F E A T U R E   A N G L E
    
    
    if {[info exists proc_scrib_data(ignore_angle)] && (0 < $proc_scrib_data(ignore_angle))} {
	
	scrib_lib_quad_remove_bars_by_feature_angle $fams_for_smoothing {NODE BAR_2 LINE_2 TRI_3 QUAD_4} \
		$proc_scrib_data(ignore_angle)
    }
    scrib_lib_quad_clean_debug_save
    
    scrib_lib_uns_unload_reload ;#xxx this helps it along... --jka
    
    if {0 >= [ic_count_nodes]} {
	mess "Error mesh $thismesh has no nodes...returning from quad cleaning process\n" red
	return 1
    }
    
    set fams_for_smoothing [ic_uns_list_families]
    foreach fam [list $bolt_hole_bar_fam $outer_ring_bars ORFN] {
    	set ind [lsearch $fams_for_smoothing $fam]
    	set fams_for_smoothing [lreplace $fams_for_smoothing $ind $ind]
    }
    mess "Cleaning families: $fams_for_smoothing...\n"
    
    ic_uns_set_node_dimension all same    
    
    #xx set all_types {NODE BAR_2 TRI_3 QUAD_4}
    
    
    # ----> M E R G E   N O D E S   A T   T O L E R A N C E (not yet we don't) xxx
        
    
    # ----> S M O O T H   A L L   T R I S   A N D   Q U A D S
    
    
    set stat [scrib_lib_quad_smooth_by_quality $fams_for_smoothing \
	    {TRI_3 QUAD_4} {Custom quality}]
    scrib_lib_quad_clean_debug_save
    
    
    # ----> R E M E S H I N G   O N   T Y P E S  { TRI_3 QUAD_4 }
    
    if 1 {
	#scrib_lib_uns_unload_reload
	set stat [scrib_lib_quad_remesh_by_quality $fams_for_smoothing \
		{TRI_3 QUAD_4} {Custom quality}]
	scrib_lib_quad_clean_debug_save
    }
    
    # ----> R E M E S H   B A D   E L E M E N T S
    
    
    set stat [scrib_lib_quad_remesh_bad_elements_by_quality $fams_for_smoothing \
	    {TRI_3 QUAD_4} {Custom quality}]
    scrib_lib_quad_clean_debug_save
    
    # ----> S M O O T H   A L L   T R I S   A N D   Q U A D S  ( A G A I N )
    
    
    set stat [scrib_lib_quad_smooth_by_quality $fams_for_smoothing \
	    {TRI_3 QUAD_4} {Custom quality}]
    scrib_lib_quad_clean_debug_save
    
    
    # ----> M O V E   Q U A D _ 4 s   I N   B O L T   F A M I L Y   T O   O R I G _ F A M
    
    
    if {[ic_uns_count_in_families $bolt_hole_bar_fam] && [info exists orig_fam]} {
	set quadsubset [ic_uns_subset_create [ic_uns_subset_unused_name bsm_moveboltsback]]
	ic_uns_subset_add_families_and_types $quadsubset $bolt_hole_bar_fam {TRI_3 QUAD_4}
    	ic_uns_change_family $quadsubset $orig_fam
	ic_uns_subset_delete $quadsubset
	scrib_lib_quad_clean_debug_save
    }
    
    
    # ----> A G A I N   D E L E T E { disconnected bars, hanging bar elements, disconnected vertices }
    
    
    scrib_lib_quality_disconnected_bars_delete    
    scrib_lib_quality_hanging_bars_delete    
    ic_delete_disconnected_vertices
    scrib_lib_quad_clean_debug_save
    
    
    return 0
}

## This procedure will look at bars in fam $fam, and if they are attached to valid O-ring 
## return a subset of quads

proc scrib_lib_quad_check_ogrid {nodes fam} {
    
   set good_quads ""
   set att_nodes ""
   set remaining_nodes ""
   foreach node $nodes {
       set tmpsubset [ic_uns_subset_create]
       ic_uns_subset_add_near_nodes $tmpsubset $node
       set elem_list ""
       set node_list ""
       foreach el [ic_uns_subset_list_elements $tmpsubset] {
	   if {[lindex $el 0] == "QUAD_4" } {
	       lappend elem_list [lindex $el 1] 
	       set node_list [concat $node_list [lindex $el 2]]
	   }
       } 
       ic_uns_subset_delete $tmpsubset
       #
       #
       #       
       if {[llength $elem_list] == 2} {
	   append good_quads $elem_list
	   set att_nodes [concat  $att_nodes $node_list]
	   continue
       } else {
	   lappend remaining_nodes $node
	   mess "node $node was found to not be connected to Oring\n" red
       }
       
   }
   set att_nodes [lsort -unique $att_nodes]
   set outer_nodes "" 
   foreach n $att_nodes {
	if {[lsearch $nodes $n] == -1} {lappend outer_nodes $n}
   }
   foreach node $outer_nodes {
	set tmpsubset [ic_uns_subset_create]
    	ic_uns_subset_add_near_nodes $tmpsubset $node
	set elem_list ""
	set node_list ""
    	foreach el [ic_uns_subset_list_elements $tmpsubset] {
	    if {[lindex $el 0] == "QUAD_4" } {
	    	lappend elem_list [lindex $el 1]
		#append node_list [lindex $el 2]
	    }
    	} 
	ic_uns_subset_delete $tmpsubset
	if {[llength $elem_list] == 4} {
	    set good_quads [concat $good_quads $elem_list]
	    #append att_nodes $node_list
	    continue
	} else {
	    mess "There are not 4 elements attached to $node\n" red
	    mess "Can not add $elem_list\n" red
	}

   }  
   set good_quads [lsort -unique $good_quads] 
   mess "good_quads are $good_quads\n" blue
   return [list $good_quads $remaining_nodes]
}

## needs doc. xxx

proc scrib_lib_quad_connected_to_nodes {nodes {make_ogrid 0} {num_layers 2}} {
    
    set nsubset [ic_uns_subset_create]
    
    foreach node $nodes {
	if !$make_ogrid {
  	    ic_uns_subset_add_near_nodes $nsubset $node
	    continue
	}	
	set tmpsubset [ic_uns_subset_create]
    	ic_uns_subset_add_near_nodes $tmpsubset $node
	set OK 1
    	set list [ic_uns_subset_list_elements $tmpsubset]
    	#mess "DEBUG: list is $list\n"    
	#
    	# find number attached (2 quads would be ideal, however we will assume 
	# 2 tris will also work)
	#
    	set elem_list ""
    	foreach el $list {
	    if {([lindex $el 0] == "QUAD_4") || ([lindex $el 0] == "TRI_3")} {
	    	set num [lindex $el 1]
	   	lappend elem_list $num 
	    	set etype($num) [lindex $el 0]
	    	set enodes($num) [lindex $el 2]
	    	#mess "elem_list is $elem_list \n"
	    }
    	}
    	set num_connected [llength $elem_list]
	if {$num_connected == 2} {
	    ic_uns_subset_add_near_nodes $nsubset $node
	    continue
	}
	set OK 0
	#
	# --> if there are more than 2 elements attached, it is not an O-grid ------
	#
	mess "node is $node \n" green
    	set diagonals ""
	set diag_nodes ""
    	set base_elements ""
	set base_nodes ""
    	foreach el $elem_list {
	    set nodes_on_hole($el) ""
	    foreach n $enodes($el) { 
	    	if  {[lsearch $nodes $n] != -1} { 
		    lappend nodes_on_hole($el) $n
		    if {$n != $node} {set away_base_node($el) $n}
	    	}
	    }
	    if {[llength $nodes_on_hole($el)] < 2} {
		lappend diagonals $el
		append diag_nodes $enodes($el) 
	    } else {
	  	lappend base_elements $el
		#append base_nodes $enodes($el)
	    }
	}
	mess "base elements are $base_elements\n" blue
	mess "diagonals are $diagonals\n" blue
	if {[llength $base_elements] != 2} {
	    mess "The bars are not correctly attached to the mesh.\n" red
	    continue
	}
	#For simple cases
	#------Find  key nodes, 
	#
	#                                      quad_fourth_node
	#                                             /    \                        
	#                                           /        \                   
	#                      X____tip_node       tip_node____X
	#                      |                   \          /                  |
	#                      |       Base       \      /       Base      |
	#         away_base_node ____ node _____ away_base_node
	
	#Compare the nodes of diagonals with base------------------------------
	set order_list_dim ""
	foreach base $base_elements { 
	    foreach n $enodes($base) {
		if {$n == $node} continue
		if {[lsearch $diag_nodes $n] == -1} continue
		set tip_node($base) $n

		#lappend order_list_dim "$base $tip_node_dim($base)"	
	    } 
	}
	#reorder base_elements according to dimension (lowest dim first)
	#set ordered [lsort -index 1 $order_list_dim]
	#set base_elements ""
	#foreach o $ordered {lappend base_elements [lindex $o 0]}
	#mess "tip node for $base is $tip_node($base), dim is $tip_node_dim($base); \
	#	position is $tip_node_pos($base) \n" blue
	#------------------------------------------
	set base_tri ""
	set base_quad ""
	foreach base $base_elements {
	    if {$etype($base) == "TRI_3"} { 
		lappend base_tri $base} else {lappend base_quad $base}
	}
	set problem 0
	foreach d $diagonals {
	    if {$etype($d) == "TRI_3"} {
		set quad_fourth_node($d) ""
		continue
	    } else {
		set att_nodes ""
		set quad_fourth_node($d) ""
		foreach e $base_elements {append att_nodes $enodes($e)}
		foreach e $diagonals {
		    if {$e == $d} {continue} else {append att_nodes $enodes($e)}
		}
		foreach n $enodes($d) {
		    if {[lsearch $att_nodes $n] == -1} {lappend quad_fourth_node($d) $n}
		}
		if {[llength $quad_fourth_node($d)] > 1} {
		    mess "There should only be one node not attached to other elements. \
			quad_fourth_node($d) is $quad_fourth_node($d)\n" red
		     set problem 1
		}
	    }
	}
	if $problem {continue}
	#===============================================
	#Find base\diag, and operate off these pairs
	#===============================================
 	if {[llength $diagonals] == 1 } {
	    if {$base_tri != ""} {
		set base [lindex $base_tri 0] 
		if {[lindex $base_tri 1] != ""} {set base2 [lindex $base_tri 1]} else {set base2 $base_quad}
	    } else {
		set base [lindex $base_quad 0]
		set base2 [lindex $base_quad 1]
	    }
	    set OK [scrib_lib_quad_reorganize_connections $base $etype($base) $diagonals \
		$etype($diagonals) $node $away_base_node($base) $tip_node($base) \
		$tip_node($base2) $quad_fourth_node($diagonals)]

	} elseif {[llength $diagonals] == 2} {
	    #================ relate diag with its base
	    set diag(1) [lindex $diagonals 0]
	    set diag(2) [lindex $diagonals 1]
	    set base_el(1) ""
	    set base_el(2) ""
	    foreach base $base_elements {
		set common 1
		foreach n $enodes($base) {
		    if {($n == $node) || ([lsearch $enodes($diag(1)) $n] == -1)} {continue
		    } else {set common 0}
		}
		if $common {lappend base_el(2) $base} else {lappend base_el(1) $base}
	    }
	    if {([llength $base_el(1)] > 1) || ([llength $base_el(2) ] > 1) } {
		mess "base_el(1) is $base_el(1); base_el(2) is $base_el(2) ..... ERROR\n" red
		set OK 0
		continue
	    }
	    foreach n $enodes($diag(1)) {
		if {([lsearch $enodes($diag(2)) $n] != -1) && ($n != $node)} {set common_node $n}
	    }
	    foreach ind {1 2} {
		set OK_$ind [scrib_lib_quad_reorganize_connections \
			$base_el($ind) $etype($base_el($ind)) \
			$diag($ind) $etype($diag($ind)) $node $away_base_node($base_el($ind)) \
			$tip_node($base_el($ind)) $common_node $quad_fourth_node($diag($ind))]
	    }
	    if {$OK_1 && $OK_2 } {set OK 1} else {set OK 0}
	} else {
	    mess "Does not currently support cases with [llength $diagonals] diagonals\n" red
	}
	#============================================
 	if $OK {
	    ic_uns_subset_add_near_nodes $nsubset $node
	}
    }
    return $nsubset
}

##  needs doc. xxx

proc scrib_lib_quad_reorganize_connections { base base_type diag diag_type node \
	away_base_node tip_node {away_diag_node ""} {quad_fourth_node ""}} {

    mess "base/diag is $base / $diag \n" green
    #find dimensions/positions of nodes
    foreach n [list $node $away_base_node $tip_node $away_diag_node \
		$quad_fourth_node] {
	if {$n == ""} {continue}
	set tmp [ic_uns_subset_create]
	ic_uns_subset_add_node_numbers $tmp $n
	set item [lindex [ic_uns_subset_list_elements $tmp] 0]
    	set tip_node_dim($n) [lindex $item 3]
   	set tip_node_pos($n) [lindex $item 2]
	ic_uns_subset_delete $tmp	
    }

    if {($base_type == "TRI_3") && ($diag_type == "QUAD_4")} {
	if {($tip_node_dim($tip_node) > 1) && 0} {
	    set tmp [ic_uns_subset_create]
	    ic_uns_subset_add_node_numbers $tmp "$away_base_node $tip_node"
	    ic_uns_merge_nodes $tmp "" {0 0 0} {0 0 0} 0
	    ic_uns_subset_delete $tmp
	    return 1
	} else {
	    #
	    # delete diagonal and base tri
	    #
	    set elset [ic_uns_subset_create]
	    ic_uns_subset_add_number $elset "$base" 1 TRI_3	
	    ic_uns_subset_add_number $elset "$diag" 1 QUAD_4
	    set fam [lindex [ic_uns_subset_list_families $elset] 0]
	    ic_delete_elements subset $elset 
	    ic_uns_subset_delete $elset		
	    #
	    # recreate elements to fill hole
	    #
	    set nodesubset [ic_uns_subset_create]
	    ic_uns_subset_add_node_numbers $nodesubset $node 
	    ic_uns_subset_add_node_numbers $nodesubset $away_diag_node 
	    ic_uns_subset_add_node_numbers $nodesubset $away_base_node
	    ic_uns_subset_add_node_numbers $nodesubset $tip_node
	    ic_uns_create_element QUAD_4 $fam $nodesubset 0
	    ic_uns_subset_delete $nodesubset
	    #
	    #
	    #
	    set node2subset [ic_uns_subset_create]
	    ic_uns_subset_add_node_numbers $node2subset "$tip_node	$away_diag_node \
		$quad_fourth_node"
	    ic_uns_create_element TRI_3 $fam $node2subset 0
	    ic_uns_subset_delete $node2subset
	    #
	    return 1
	}
    }
    if {($base_type == "TRI_3") && ($diag_type == "TRI_3")} {
	set elset [ic_uns_subset_create]
 	ic_uns_subset_add_number $elset "$diag" 1 TRI_3
	ic_uns_subset_add_number $elset "$base" 1 TRI_3
	set error_trap [ic_uns_subset_list_elements $elset]
	mess $error_trap
	if {[llength $error_trap] != 2} {
	    mess "There are not two elements in diag/quad pair?  \
		Maybe two elements have the same number, or one element doesn't exist.\n" red
	    return 0
	}
	set nodes1 [lindex [lindex $error_trap 0] 2]
	set nodes2 [lindex [lindex $error_trap 1] 2]
	set num_con_nodes 0 
	foreach n $nodes1 {
	    if {[lsearch $nodes2 $n] != -1} {incr num_con_nodes}
	}
	if {$num_con_nodes != 2} {
	    mess "There is a messed up attachment at this node.\n" red
	    return 0
	}
	ic_uns_quad_remesh $elset element 4 proj 1 
	ic_uns_subset_delete $elset
	return 1
    }
    if {$base_type == "QUAD_4"} {
	#if both the common node and the base_tip are the same dim, move to middle.
	if { $tip_node_dim($away_diag_node) == $tip_node_dim($tip_node)} {
	    foreach i { 0 1 2 } {
		set ii 0
		foreach b [list $away_diag_node $tip_node] {
		    set vv($ii) [lindex $tip_node_pos($b) $i]
		    incr ii
		}
		 set v($i) [expr (($vv(0) + $vv(1) ) / 2 )]
	    }
	    set nodeset [ic_uns_subset_create]
	    ic_uns_subset_add_node_numbers $nodeset "$away_diag_node $tip_node "
	    ic_uns_move_nodes_exact $nodeset 1 {1 1 1} "$v(0) $v(1) $v(2)"
	    ic_uns_set_node_dimension $nodeset same
	    ic_uns_subset_delete $nodeset
	}
	#order nodes by dimension
	set unsorted ""	
  	foreach n [list $away_diag_node $tip_node $quad_fourth_node] {
	   if {$n != ""} {lappend unsorted [list $n $tip_node_dim($n)]}
	}
	set sorted [lsort -index 1 $unsorted]
	set nodeset [ic_uns_subset_create]
	foreach nn $sorted {
	    set n [lindex $nn 0]
	    if {$n != ""} {ic_uns_subset_add_node_numbers $nodeset $n}
	}
	ic_uns_merge_nodes $nodeset "" {0 0 0} {0 0 0} 0
	ic_uns_subset_delete $nodeset
	return 1
    }
}

## needs doc.  xxx

proc scrib_lib_quad_remesh_by_quality {fams types metric {quality_limit 0.6} {eltype 2} {ignore_dim 0} } {
    
    if ![ic_uns_is_loaded] {
	mess "cannot remesh quads by quality without a mesh loaded!\n" red
	return -1
    }
    
    # --> I *think* this makes is safer...just in case tetin is not loaded --jka
    if [ic_geo_is_loaded] {
	set projval 1
	set geo_activeval 1
    } else {
	set projval 0
	set geo_activeval 0
    }
    mess "Remeshing parts of the quad mesh with[expr {$projval? "" : "out"}] geometry loaded...\n"
    
    # --> Create a map of the given families and types
    set quality_limit [expr $quality_limit * 1.0]
    set tmpmap [ic_uns_subset_create [ic_uns_subset_unused_name bsm_tmpmap]]
    ic_uns_subset_add_families_and_types $tmpmap $fams $types	
    scrib_lib_quad_clean_debug_save
    
    # --> Set the given metric on that map

    set metric_weights [scrib_lib_quality_standard_weights $metric]
    ic_uns_set_metric_weights $metric_weights 1
    set stat [ic_uns_metric $tmpmap $metric]
    if ![lindex $stat 0] {	
	ic_uns_subset_delete $tmpmap
	return -1 
    }
    
    # --> Create a map into which we put the elements from above map in rangle -1000->.6

    set badmap [ic_uns_subset_create [ic_uns_subset_unused_name bsm_quad_remesh_by_quality]]
    ic_uns_set_with_diagnostic $badmap $tmpmap [list [list -1000 $quality_limit]]
    
    set nbads [ic_uns_subset_count_elements $badmap]
    if $nbads {
	
	mess "found $nbads elements which need to be remeshed..."
	
	#xx Ben, look at this proc in med/lib/med_batch/unstruct_funcs.tcl and make sure you *really*
	#xx want to be using it. note that it first *clear*s the given map! --jka
	#xx ic_uns_subset_add_families_and_types $badmap $fams {NODE BAR_2 TRI_3 QUAD_4}
	
	ic_uns_subset_add_layer $badmap 2 $fams all ;#$types
	mess "adding [expr {[ic_uns_subset_count_elements $badmap] - $nbads}] elements of types $types in the surrounding layers of fams $fams...\n"
	
	scrib_lib_quad_clean_debug_save
	#xx if you want to see what the badmap looks like, uncomment next line
	#xx return
	
	#xx isn't this (how=7) just "ic_remesh_bad_elements" elsewhere in quad_clean you have:
	#xx    ic_remesh_bad_elements $badmap 2 1 0 1 20 (twice)
	#xx ...
	
	if [catch {ic_uns_quad_remesh $badmap all_dim $ignore_dim element $eltype proj $projval ele_tol 0 \
		how 7 geo_active $geo_activeval iters 20} res] {
	    mess "\n$res\n" red
	}
	
	mess "...done\n"
	
    } else {
	mess "...all element qualities acceptable!\n"
    }
    
    ic_uns_subset_delete $badmap
    ic_uns_subset_delete $tmpmap
    
    return 0
}

## needs doc.  xxx

proc scrib_lib_quad_remesh_bad_elements_by_quality {fams types metric {quality_limit 0.6} {eltype 2} {ignore_dim 0} } {
        
    if ![ic_uns_is_loaded] {
	mess "cannot remesh bad elements by quality without a mesh loaded!\n" red
	return -1
    }
    
    mess "\nRemeshing bad elements by quality metric $metric...\n"
    
    set tmpmap [ic_uns_subset_create [ic_uns_subset_unused_name bsm_remesh_bad_elements]]   
    ic_uns_subset_add_families_and_types $tmpmap $fams $types
    
    # --> G E T   M E T R I C S
    set quality_limit [expr $quality_limit * 1.0]    
    set metric_weights [scrib_lib_quality_standard_weights $metric]
    ic_uns_set_metric_weights $metric_weights 1
    set stat [ic_uns_metric $tmpmap $metric]
    if ![lindex $stat 0] { return -1 }
    
    set badmap [ic_uns_subset_create [ic_uns_subset_unused_name bsm_quad_remesh_by_quality]]
    ic_uns_set_with_diagnostic $badmap $tmpmap [list [list -1000 $quality_limit]]
    ic_uns_subset_delete $tmpmap
    
    set nbads [ic_uns_subset_count_elements $badmap]
    if $nbads {
	
	mess "found $nbads elements which need to be remeshed\n"
	# watch out for all_dim of 1---it has to be 0 elsewhere!! xxx --jka
	if [catch {ic_remesh_bad_elements $badmap $eltype 1 0 $ignore_dim 20} res] {
	    mess "\n$res\n"
	}
	mess "...done\n"
    } else {
	mess "...all element qualities acceptable!\n"
    }
    
    ic_uns_subset_delete $badmap

    return 0
}

## needs doc. types is ignored!!!!xxxxx

proc scrib_lib_quad_smooth_by_quality {fams types metric {quality_limit 0.5} {additional_args ""}} {
    
    if ![ic_uns_is_loaded] {
	mess "cannot smooth by quality without a mesh loaded!\n" red
	return -1
    }
    
    mess "Smoothing elements of type $types in families $fams according to $metric...\n"
    
    set tmpmap [ic_uns_subset_create [ic_uns_subset_unused_name bsm_smooth_by_quality]]    
    ic_uns_subset_add_families_and_types $tmpmap $fams $types    
    
    # --> G E T   M E T R I C S
    
    set metric_weights [scrib_lib_quality_standard_weights $metric]
    ic_uns_set_metric_weights $metric_weights 1
    
    set stat [ic_uns_metric $tmpmap $metric]
    ic_uns_subset_delete $tmpmap
    
    if {0 < [lindex $stat 0]} {
	set args " upto $quality_limit iterations 10 metric [list $metric] family $fams"
	foreach t $types {
	    lappend args smooth
	    lappend args $t
	}
	if {$args != ""} {set args [concat $args $additional_args]}
	mess "args are $args\n"
	set cmd  ic_smooth_elements 
	if [catch {eval $cmd $args} res] {
	    mess "\n$res\n" red
	}
    }
    
    mess "...done\n"
    
    return 0
}

##  needs doc.  xxx

proc scrib_lib_quad_remesh_by_diagnostic {diagnostic fams types nlayers {eltype 2} {ignore_dim 0}} {
    return [ic_quad_remesh_by_diagnostic $diagnostic $fams $types $nlayers $eltype $ignore_dim]
}

## needs doc. xxx

proc scrib_lib_quad_find_offset_quad_layers { num_layers bar_fam } {
    
    if {0 >= [ic_uns_count_in_families $bar_fam]} { return "" }
    
    set barsubset [ic_uns_subset_create [ic_uns_subset_unused_name bsm_boltbars]]
    ic_uns_subset_add_families_and_types $barsubset $bar_fam {BAR_2 LINE_2}
    set ring_of_quads ""
    
    set subsets [scrib_lib_uns_get_connected_bar_regions $barsubset]
    
    foreach sub $subsets {
	#first find if first layer of elements is a ring of quads. 
	#mess "[ic_uns_subset_list_elements $sub]\n" green
	set els_to_add ""
	set nogo 0
	set used_nodes [scrib_lib_uns_subset_get_nodes $sub]
	foreach node $used_nodes {
	    set nsubset [ic_uns_subset_create]
	    ic_uns_subset_add_near_nodes $nsubset $node
	    if {[scrib_lib_uns_number_type $nsubset QUAD_4] != 2} { 
		set nogo 1
	    }
	    ic_uns_subset_delete $nsubset
	    if {1 == $nogo} { break }
	}
	if $nogo {
	    mess "stopping layer adding on nodes $used_nodes\n" green
	    continue
	} else {
	    mess "found good first layer, continuing adding\n" blue
	    mess "num_layers is $num_layers\n"
	}
	for {set i 1} {$i < $num_layers} {incr i} {
	    ic_uns_subset_add_layer $sub 1 all all 
	    set new_nodes [scrib_lib_uns_subset_get_nodes $sub]
	    set node_dup [concat $new_nodes $used_nodes] 
	    set top_nodes [scrib_lib_utils_list_remove_duplicated $node_dup]	
	    #mess "used_nodes are $used_nodes\n new_nodes are $new_nodes\n top_nodes are $top_nodes\n" grey
	    set used_nodes $new_nodes
	    set nogo 0
	    foreach node $top_nodes {
		set nsubset [ic_uns_subset_create]
		ic_uns_subset_add_near_nodes $nsubset $node
		if {[scrib_lib_uns_number_type $nsubset QUAD_4] != 4} { 
		    set nogo 1
		}
		ic_uns_subset_delete $nsubset
		if {1 == $nogo} { break }
	    }
	    if $nogo {
		mess "stopping layer adding on nodes $top_nodes\n" green
		break
	    }
	}
	
    }
    set subset_for_quads [ic_uns_subset_create]
    foreach sub $subsets {
	ic_uns_subset_add_from $subset_for_quads $sub
	ic_uns_subset_delete $sub
    }
    return $subset_for_quads
}

## needs doc.  xxx

proc scrib_lib_quad_merge_ring_nodes {boltfam tol} {
    
    if [ic_uns_count_in_families $boltfam] {
	
	mess "Merging nodes of first layer around rings of quads on family $boltfam with tolerance $tol...\n"
	
	set nodes_to_merge [scrib_lib_uns_get_nodes_between_fams all]
	if ![llength $nodes_to_merge] { return 1 }
	set name [ic_uns_subset_unused_name "bsm_merge_ring_nodes"]
	set quadsubset [ic_uns_subset_create $name]
	
	ic_uns_subset_add_near_nodes $quadsubset $nodes_to_merge
	
	foreach i [ic_uns_subset_list_elements $quadsubset] {
	    
	    if {[lindex $i 3] == $boltfam} {
		continue
	    }	    
	    append nodes_to_merge [lindex $i 2]
	}
	set nodes_to_merge [lsort -unique $nodes_to_merge]
	
	ic_uns_subset_clear $quadsubset
	ic_uns_subset_add_node_numbers $quadsubset $nodes_to_merge
	ic_uns_merge_nodes $quadsubset $tol
	ic_uns_subset_delete $quadsubset
	
	mess "...done\n"
    }
    
    return 0
}

## needs doc

proc scrib_lib_quad_remove_bars_by_feature_angle {fams types angle} {

    if {0 >= $angle} {
	mess "cannot remove bars by feature angle $angle" red
	return -1
    }
     
    if ![ic_uns_is_loaded] {
	mess "cannot remove bars at feature angle without a mesh loaded!\n" red
	return -1
    }
    
    mess "Removing bars below feature angle $angle...\n"        
    
    if {$fams == ""} { set fams [ic_uns_list_families] }

    set bt_set [ic_uns_subset_create [ic_uns_subset_unused_name bsm_remove_bars_angle]]
    
    ic_uns_subset_add_families_and_types $bt_set $fams $types
    
    ic_uns_build_mesh_topo $bt_set -filter_bars -angle $angle -filter_points -angle $angle
    ic_uns_subset_delete $bt_set
    
    mess "...done\n"

if 0 {
    scrib_lib_uns_unload_reload ;#xxx this helps it along... --jka
    
    if {0 >= [ic_count_nodes]} {
	mess "Error mesh $thismesh has no nodes...returning from quad cleaning process\n" red
	return 1
    }

    #is the below still necessary?  I think this should be able to be removed if the ic_uns_build_topo
    #updates nodes properly...  BEN
    
    ic_uns_set_node_dimension all same    
    
    #xx set all_types {NODE BAR_2 TRI_3 QUAD_4}
    
}        
    return 0
}


## needs doc.  xxx

proc scrib_lib_quad_convert_tri_quad {fams {frozenfams ""}} {
    
    if ![ic_uns_is_loaded] {
	mess "cannot convert tri to quad without a mesh loaded!\n" red
	return -1
    }
    
    mess "\nConverting elements from tri to quad in families $fams.\n"
    
    set types {TRI_3 QUAD_4}

    set tmpmap [ic_uns_subset_create [ic_uns_subset_unused_name bsm_tri_quad]] 
    ic_uns_subset_add_families_and_types $tmpmap $fams $types    

   #the below is necessary due to a bug in ic_uns_subset_add_families.  Somehow the frozenfams are also added?
   if {$frozenfams != ""} { 
	set tmp2map [ic_uns_subset_create]
	ic_uns_subset_add_families_and_types $tmp2map $frozenfams $types
	ic_uns_subset_subtract_from $tmpmap $tmp2map
	ic_uns_subset_delete $tmp2map
    }
  
    if [catch {ic_change_tri_quad_2 $tmpmap 1 0 0 1 160 20 0
	} res] {
	    mess "\n$res\n" red
    }
    ic_uns_subset_delete $tmpmap
    #mess "...done\n"
    
    return 0
}


proc scrib_lib_quad_remesh_bad_elements {fams types metric {quality_limit 0.6} \
	{eltype 2} {ignore_dim 0} {frozenfams ""} } {
        
    if ![ic_uns_is_loaded] {
	mess "cannot remesh bad elements by quality without a mesh loaded!\n" red
	return -1
    }
    
    mess "Remeshing bad elements by quality metric $metric...\n"

    set tmpmap [ic_uns_subset_create [ic_uns_subset_unused_name bsm_remesh_bad_elements]]   
    ic_uns_subset_add_families_and_types $tmpmap $fams $types

   #the below is necessary due to a bug in ic_uns_subset_add_families.  Somehow the frozenfams are also added?
   if {$frozenfams != ""} {
	set tmp2map [ic_uns_subset_create]
	ic_uns_subset_add_families_and_types $tmp2map $frozenfams $types 
 	ic_uns_subset_subtract_from $tmpmap $tmp2map
	ic_uns_subset_delete $tmp2map
    }
 
    # --> G E T   M E T R I C S
    set quality_limit [expr $quality_limit * 1.0]    
    set metric_weights [scrib_lib_quality_standard_weights $metric]
    ic_uns_set_metric_weights $metric_weights 1
    set stat [ic_uns_metric $tmpmap $metric]
    if ![lindex $stat 0] { return -1 }
    
    set badmap [ic_uns_subset_create [ic_uns_subset_unused_name bsm_quad_remesh_by_quality]]
    ic_uns_set_with_diagnostic $badmap $tmpmap [list [list -1000 $quality_limit]]
    #ic_uns_subset_delete $tmpmap
    
    set nbads [ic_uns_subset_count_elements $badmap]
    if $nbads {
	
	mess "found $nbads elements which need to be remeshed\n"
	# watch out for all_dim of 1---it has to be 0 elsewhere!! xxx --jka
	if [catch { ic_remesh_bad_elements $tmpmap $eltype 1 0 $ignore_dim 20 $quality_limit} res] {
	    mess "\n$res\n"
	}
	mess "...done\n"
    } else {
	mess "...all element qualities acceptable!\n"
    }
    ic_uns_subset_delete $badmap    
    ic_uns_subset_delete $tmpmap

    return 0
}


proc scrib_lib_quad_laplace_smooth {fams {upto 1.0} } {
    if ![ic_uns_is_loaded] {
	mess "cannot laplace smooth elements without a mesh loaded!\n" red
	return -1
    }
    mess "\nMeshing fams $fams with laplace smoother...\n"
    set types {QUAD_4 TRI_3}    
    set tmpmap [ic_uns_subset_create [ic_uns_subset_unused_name bsm_laplace_smooth]]   
    ic_uns_subset_add_families_and_types $tmpmap $fams $types
    ic_uns_subset_remove_layer $tmpmap 1 all

    if [catch {ic_uns_quad_remesh $tmpmap quality $upto element 3 proj 1 max_steps 100 \
		   how 5 angle 0 iters 20 surf 0 element_new TRI_3 conver 0.025 \
		   geo_active 1 merge_tol 0.10 keep_fixed_nodes 1} res] {
	mess "\n$res\n"
    }

    ic_uns_subset_delete $tmpmap
}

## this will take a set of points within a max distance and attach the points to the loops of the surfaces

proc scrib_lib_quad_add_attached_points_to_loops {pt_list max_dist} {
    foreach ll [ic_geo_get_objects loop ""] {
	foreach s [ic_geo_incident loop $ll] {
	    lappend surf_loop($s) $ll
	}
    }
    foreach ps [scrib_lib_geo_nearest_surface_point_pair $pt_list $max_dist] {
	set pp [lindex $ps 0]
	set ss [lindex $ps 1]  
	if ![info exists surf_loop($ss)] {
	    mess "surface $ss does not seem to have a loop bounding it.\n" red
	    continue
	}
	set loop [lindex $surf_loop($ss) 0] 
	set pps [concat [ic_geo_loop_get_objects $loop point] $pp ]
	set pps [lsort -unique $pps]
	ic_geo_modify_loop $loop [ic_geo_boundary loop $loop] [ic_geo_incident loop $loop] $pps
    }

}


























