
# Code to read in the data from one domain file.
#
# Copyright (c) 1993 ICEM Systems.
# Author:  Wayne A. Christopher, wayne@pmac.com
#          Pacific Marketing and Consulting, Berkeley, California

proc dmess {s} {
    global dfpp
    if ![info exists dfpp] {
	set dfpp [open debug_file w]
    }
    puts $dfpp $s
    flush $dfpp
}

proc load_domain {filename {pref full} {suff ""} {map_group_numbers_to ""}} {
    global ${pref}_unstruct ${pref}_nodeset ${pref}_element_set
    global ${pref}_element_map ${pref}_group grid_type grid_mapping
    global want_bar_elements domain_family_pid domain_node_dimensions
    global grid_tetra_surface grid_tetra_to_quad want_node_elements
    global ${pref}_couplings grid_coupling_walls grid_coupling_cutouts domain
    global dont_report_domain_file_counts subtract_maps_from_couplings
    global grid_coupling_otherparts grid_coupling_otherpart_group
    global solve_optimizing
    
    set verbose 0
    
    #set noisy [expr {![info exists dont_report_domain_file_counts] &&
	#	     !$solve_optimizing}]
    set noisy [expr {![info exists dont_report_domain_file_counts]}]
    if [file exists $filename.matrix] {
	set fp [open $filename.matrix r]
	set mat [gets $fp]
	close $fp
    } else {
	set mat "1 0 0 0   0 1 0 0   0 0 1 0   0 0 0 1"
    }
    
    if {![file exists $filename]} {
	mess "$filename doesn't exist\n" red
	return 0
    }
    if [file isdirectory $filename] {
	mess "$filename is a directory\n" red
	return 0
    }
    if [catch {domainfile -name $filename -mode r -trans_matrix $mat} df] {
	mess "Error: the mesher did not create a valid mesh:\n" red
	mess "$df\n" red
	return 0
    }
    if 0 {
	# This can happen if assemblies completely cover the volume of a cabinet...
	if ![$df nnodes] {
	    mess "Warning: the mesher did not create a valid mesh: 0 nodes\n" red
	    $df close; $df delete
	    return 0
	}
    }
    if $noisy {
	mess "\tNODES: [$df nnodes]\n"
    }
    if {[$df type] != "unstructured"} {error "Wrong type of grid file [$df type]"}
    
    catch {unset domain_family_pid}
    catch {unset domain_node_dimensions}
    
    if $noisy {
	if {[[$domain getval body_shape] get -shtype] == "none"} {
	    mess "Please note: these counts include \"phantom\" cells and nodes which will be deleted later\n"
	}
    }
    
    # Load the nodeset or if there is one, reuse it.
    set nnodes [$df nnodes]
    if {$suff != ""} {
	set ns [set ${pref}_nodeset]
	set node_off [$ns get -numnodes]
	if $nnodes {
	    set nodeptr [$df read_nodes]
	    $ns set_nodes_ptr $node_off $nnodes $nodeptr
	}
    } else {
	set ns [nodeset]
	set ${pref}_nodeset $ns
	set node_off 0
	if $nnodes {
	    set nodeptr [$df read_nodes]
	    $ns set_nodes_ptr 0 $nnodes $nodeptr
	}
    }
    # Make an unstruct mesh.
    if {$suff != ""} {
	set uns [set ${pref}_unstruct]
    } else {
	set uns [unstruct_mesh -nodes $ns]
	set ${pref}_unstruct $uns
    }
 
    # For each section, make an element map and an element set.
    set dname [file tail $filename]
    set nsects [$df nsections]
    
    if {$suff == ""} {
	catch {unset ${pref}_element_set}
	catch {unset ${pref}_element_map}
	set all_em [element_map -mesh $uns]
	set ${pref}_element_map(all_volume) $all_em
    } else {
	set all_em [set ${pref}_element_map(all_volume)]
    }
    
    for {set s 0} {$s < $nsects} {incr s} {
	set si [$df section_info $s]
	set etype [lindex $si 0]
	set estart [lindex $si 1]
	set eend [lindex $si 2]
	set num [expr {$eend - $estart + 1}]
	switch $etype {
	    BAR_2  {
		if ![info exists want_bar_elements] continue
		set tname bar; set nnn BARS
	    }
	    NODE  {
		if ![info exists want_node_elements] continue
		set tname node; set nnn "FIXED NODES"
	    }
	    QUAD_4 {
		set tname quad; set nnn QUADS
	    }
	    HEXA_8 {
		if $grid_tetra_surface continue
		set tname hexa; set nnn HEXAS
	    }
	    TRI_3 {
		set tname tri; set nnn TRIS
	    }
	    PENTA_6 {
		if $grid_tetra_surface continue
		set tname penta; set nnn PENTAS
	    }
	    PYRA_5 {
		if $grid_tetra_surface continue
		set tname pyra; set nnn PYRAS
	    }
	    TETRA_4 {
		if $grid_tetra_surface continue
		set tname tetra; set nnn TETRAS
	    }
	    default continue
	}
	
	if $noisy {
	    mess "\t$nnn: [expr {$eend - $estart + 1}]\n"
	}
	set stuff [$df read_elements $s]
	if [info exists ${pref}_element_set($tname)] {
	    set es [set ${pref}_element_set($tname)]
	    set es_off [$es get -num]
	    $es set_elements_ptr $es_off $num $stuff $node_off
	} else {
	    set es [element_set -mesh $uns -type $etype]
	    set ${pref}_element_set($tname) $es
	    $es set_elements_ptr 0 $num $stuff $node_off
	    set es_off 0
	}
	# If this is a submesh then this will get set to the right thing.
	set range_for_eset($es) [list $estart $eend $es_off]
	
	if [info exists ${pref}_element_map($tname)] {
	    set em [set ${pref}_element_map($tname)]
	} else {
	    set em [element_map -mesh $uns]
	    set ${pref}_element_map($tname) $em
	}
	$em add_range $es $es_off [expr {$num + $es_off - 1}]
	
	$all_em add_range $es 0 [expr {$num + $es_off - 1}]
    }
    
    # For some reason this fails with iceboard meshes.  Since we don't use
    # local refinement yet let's leave it out
    #set ${pref}_couplings [$df read_couplings]
    set ${pref}_couplings {}
    if 0 {
	foreach coup [set ${pref}_couplings] {
	    foreach coupn $coup {
		foreach face $coupn {
		    foreach f $face {
			if {$f > 1000000} {
			    error "bogus node in coupling"
			}
		    }
		    if 0 {
			if {[llength $face] != 4} {
			    set fp [open "error_to_wayne" w]
			    puts $fp [set ${pref}_couplings]
			    close $fp
			    error "non-conformal face must have 4 nodes: it has 3"
			}
		    }
		}
	    }
	}
    }
    set ndomains [$df ndomains]
    if $noisy {
	mess "\tFaces + solids: $ndomains\n"
	set nr [llength [set ${pref}_couplings]]
	set nc 0
	foreach c [set ${pref}_couplings] {
	    incr nc [llength [lindex $c 0]]
	    incr nc [llength [lindex $c 1]]
	}
	if $nc { mess "\tNon-conformal mesh: $nr regions, $nc cells\n" }
    }
    set has_sec_domains 0
    
    if $noisy {
	puts -nonewline "$filename, $ndomains domains: "
    }
    
    # Now figure out the unstructured domains in the file.  These are kept in
    # the variable ${pref}_groups.
    
    array set map_group $map_group_numbers_to
    
    if $verbose { parray map_group }
    if [info exists grid_coupling_walls($suff)] {
	foreach {s m o} $grid_coupling_walls($suff) {
	    set "map_coupling($m $o)" [list $map_group(0) $s]
	    if $verbose { puts "group $m $o goes to $map_group(0) $s" }
	}
    }
    
    if [info exists grid_coupling_cutouts($suff)] {
	foreach n $grid_coupling_cutouts($suff) {
	    set is_cc($n) 1
	}
    }
    if [info exists grid_coupling_otherparts($suff)] {
	foreach {n opname opside} $grid_coupling_otherparts($suff) {
	    set is_cc($n) 1
	    set opdata($n) [list $suff $opname $opside]
	}
    }
    catch {parray opdata}
    #puts "MAP $map_group_numbers_to SUFF $suff"
    #set ffpp [open DEBUG w]

    for {set d 0} {$d < $ndomains} {incr d} {
	set dname [$df dname $d]; set fname $dname
	set utype [$df uns_type $d]
	
	regsub -all _ $dname " " dname
	set dname [string trim $dname]
	
	#dmess "\tregion \"$dname\" $utype\n"
	
	if {($utype == "by_pid") && ![info exists map_coupling($dname)]} {
	    set pid [$df pid_for_name $fname]
	    set domain_family_pid($fname) $pid
	    
	    # Now we have to do the mapping of subgrid groups to overall groups,
	    # for tetra meshes.  If we do it after all the meshes are merged we
	    # won't be able to separate them.
	    
	    #puts "deal with $fname"
	    
	    if {$fname == "DEF_JUNK"} {
		# Just delete them
		foreach es [array names range_for_eset] {
		    $es delete_with_pid $pid -1
		}
		continue
		
	    } elseif [string match CREATED_MATERIAL_* $fname] {
		# These are volume elements that are (maybe) inside some
		# solid object that were cut off from its material point.
		#puts "find_containing_shape_for_created_material $pid"
		set sh [find_containing_shape_for_created_material $pid]
		#puts "sh is $sh"
		if {$sh == ""} {
		    # Just delete them
		    foreach es [array names range_for_eset] {
			$es delete_with_pid $pid -1
		    }
		    continue
		}
		if {($sh == "FLUID") || ([$sh get -object] == $domain)} {
		    foreach es [array names range_for_eset] {
			$es change_pid_in_range $pid $domain_family_pid(_0_0) \
			    [lindex $range_for_eset($es) 2] [$es get -num]
		    }
		    continue
		} else {
		    set obnum $grid_mapping($sh)
		    set sidenum 0
		}
		
	    } else {
		if {[string index $fname 0] != "_"} continue
		if ![string compare $fname __JUNK__] continue
		
		set xx [split $fname _]
		if ![info exists map_group([lindex $xx 1])] continue
		if {([lindex $xx 1] == 0) && ([lindex $xx 2] == 0)} continue ;# don't need explicit fluid groups
		set obnum $map_group([lindex $xx 1])
		set sidenum [lindex $xx 2]
	    }
	    set mvar "${pref}_group($obnum $sidenum)"
	    if [info exists $mvar] {
		set map [set $mvar]
	    } else {
		set map [element_map -mesh [set ${pref}_unstruct]]
		set $mvar $map
	    }
	    $map add_pid $pid NULL
	    $map set_pids 0 -1
	    $map convert_ranges_to_lists
	    
	    # This is where we see this in the tetra case.
	    # What is stored in the opdata array is, for each other assembly, a list of
	    #    this_suff other_suff this_side
	    # that is, the adjacent asm is on this_side side of us.  However the number in
	    # the domain name we are reading is going to be the side of the other assembly.
	    # So deal with that appropriately ...  this of course assumes that assemblies
	    # are hexas.
	    if [info exists opdata([lindex $dname 0])] {
		set opd $opdata([lindex $dname 0])
		set mynum [expr [lindex $opd 2] - 1]
		set othernum [expr 1 + $mynum / 2 + ($mynum % 2 ? 0 : 1)]
		if {[lindex $dname 1] == $othernum} {
		    set grid_coupling_otherpart_group($opd) $map
		    #puts "here2 : grid_coupling_otherpart_group($opd) $map [$map count]"
		}
	    }
	    
	    continue
	}
	switch $dname {
	    "volume nodes" - "surface nodes" - "curve nodes" - "prescribed nodes" {
		# note that this doesn't happen with the auto_mesher and thus
		# won't be a concern with multiple mesh regions
		if {$utype != "range"} continue
		set contents [$df uns_contents $d]
		set domain_node_dimensions($dname) $contents
		continue
	    }
	}
	
	if ![scan $dname %d dummy] continue
	
	set odname $dname
	set this_is_cc 0
	if [info exists map_coupling($dname)] {
	    if $verbose { #puts "map_coupling($dname) exists: $map_coupling($dname)" }
	    set dname $map_coupling($dname)
	} elseif [info exists is_cc([lindex $dname 0])] {
	    if $verbose { #puts "is_cc([lindex $dname 0]) exists" }
	    set this_is_cc 1
	} else {
	    set dn1 [lindex $dname 0]; set dn2 [lindex $dname 1]
	    if [info exists map_group($dn1)] { set dn1 $map_group($dn1) }
	    if $verbose { #puts "switch $dname to \"$dn1 $dn2\"" }
	    set dname "$dn1 $dn2"
	}
	#dmess p1\n
	if $noisy {
#	    puts -nonewline .
        mess "Loading domain [expr $d + 1] of $ndomains domains\n"
	}
	if $this_is_cc {
	    set em [element_map -mesh $uns]
	    lappend subtract_maps_from_couplings $em
	} elseif [info exists ${pref}_group($dname)] {
	    # This means that we loaded up part of this object before and
	    # now we need to load up the rest of it.
	    set em [set ${pref}_group($dname)]
	    #puts "already saw $dname - it has [$em count] now"
	} else {
	    set em [element_map -mesh $uns]
	    set ${pref}_group($dname) $em
	}
	if [info exists opdata([lindex $dname 0])] {
	    set grid_coupling_otherpart_group($opdata([lindex $dname 0])) $em
	    puts "here1 : grid_coupling_otherpart_group($opdata([lindex $dname 0])) = $em [$em count]"
	}
	    
	#dmess p2\n
	switch $utype {
	    range {
		set contents [$df uns_contents $d]
		set low [lfirst $contents]
		set hi [llast $contents]
		if [info exists range_for_eset] {
		    foreach es [array names range_for_eset] {
			set s_low [lindex $range_for_eset($es) 0]
			set s_hi [lindex $range_for_eset($es) 1]
			set off [lindex $range_for_eset($es) 2]
			if {$low < $s_low} {set r_low $s_low} else {set r_low $low}
			if {$hi > $s_hi} {set r_hi $s_hi} else {set r_hi $hi}
			if {$r_hi >= $r_low} {
			    $em add_list_range $es \
				[expr {$r_low - $s_low + $off}] \
				[expr {$r_hi - $s_low + $off}] 
			}
		    }
		}
		set has_sec_domains 1
	    }
	    explicit {
		if [info exists range_for_eset] {
		    set contents [$df uns_contents $d]
		    #dmess p3\n
		    foreach es [array names range_for_eset] {
			set elems {}
			set s_low [lindex $range_for_eset($es) 0]
			set s_hi [lindex $range_for_eset($es) 1]
			set off [lindex $range_for_eset($es) 2]
			foreach l $contents {
			    if {($l >= $s_low) && ($l <= $s_hi)} {
				lappend elems [expr {$l - $s_low + $off}] 
			    }
			}
			#dmess p4\n
			if {$elems != {}} {
			    eval $em add_list $es $elems
			}
			#dmess p5\n
		    }
		    # If this mesh is a auto_mesher mesh then change all 0 pids
		    # to something else
		    $em set_pids_if 0 1 -1
	
		}
		set has_sec_domains 1
	    }
	    by_pid {
		# At this point we are dealing with coupling walls for tetra meshes
		# only
		set pid [$df pid_for_name $fname]
		$em add_pid $pid NULL
		$em set_pids 0 -1
		$em convert_ranges_to_lists
	    }
	}
	#puts "has [$em count] elements; $utype"
    }
    #close $ffpp
    if 0 {
	this is not desirable with mixed meshes
	if $has_sec_domains {
	    set grid_type auto
	} elseif {([info exists ${pref}_element_set(tetra)] ||
		   [info exists ${pref}_element_set(tri)]) &&
		  ($grid_type != "iceboard")} {
	    set grid_type tetra
	} elseif {[info exists ${pref}_element_set(hexa)] ||
		  [info exists ${pref}_element_set(quad)]} {
	    if {($grid_type == "tetra") && $grid_tetra_to_quad &&
		$grid_tetra_surface} {
		# leave it alone
	    } elseif {($grid_type != "hexa") && ($grid_type != "iceboard")} {
		set grid_type global
	    }
	} else {
	    # If this is an empty domain don't worry about it
	    if $nnodes {
		mess "Warning: don't know what kind of grid this is\n" red
	    }
	}
    }
    if $noisy {
	puts " ok"
    }

    # Set all the pids to 0 for the heck of it (also the find_containing_shape stuff will
    # get confused otherwise)
    foreach es [array names range_for_eset] {
	$es set_pids 0 -1 0
    }
    
    $df close
    $df delete
    busy reset
    return 1
}

proc domain_count_elements {file} {
    if {![file exists $file]} {
	mess "$file doesn't exist\n" red
	return 0
    }
    if [file isdirectory $file] {
	mess "$file is a directory\n" red
	return 0
    }
    
    if [catch {domainfile -name $file -mode r} df] {
	mess "Error: the mesher did not create a valid mesh\n" red
	mess "$df\n" red
	return 0
    }
    if {[$df type] != "unstructured"} {
	error "Wrong type of grid file!"
    }
    
    set nsects [$df nsections]
    set num 0
    set nums 0
    for {set s 0} {$s < $nsects} {incr s} {
	set si [$df section_info $s]
	switch [lindex $si 0] {
	    HEXA_8 - PENTA_6 - PYRA_5 - TETRA_4 {
		incr num [expr {[lindex $si 2] - [lindex $si 1] + 1}]
	    }
	    QUAD_4 - TRI_3 {
		incr nums [expr {[lindex $si 2] - [lindex $si 1] + 1}]
	    }
	}
    }
    if !$num {
	if $nums {
	    mess "Warning: the mesher created only surface and no volume elements\n" red
	} else {
	    mess "Warning: the mesher did not create a valid mesh: 0 elements\n" red
	}
    }
    
    $df close
    $df delete
    return $num
}

proc domain_count_nodes {file} {
    if {![file exists $file]} {
	mess "$file doesn't exist\n" red
	return 0
    }
    if [file isdirectory $file] {
	mess "$file is a directory\n" red
	return 0
    }
    
    if [catch {domainfile -name $file -mode r} df] {
	mess "Error: the mesher did not create a valid mesh\n" red
	mess "$df\n" red
	return 0
    }
    if {[$df type] != "unstructured"} {
	error "Wrong type of grid file!"
    }
    
    set num [$df nnodes]
    
    $df close
    $df delete
    return $num
}

proc save_domain {file pref gpref} {
    global ${pref}_unstruct ${pref}_nodeset ${pref}_element_set
    global ${pref}_element_map ${gpref}_group ${pref}_family_name
    global domain_family_pid domain_node_dimensions ${pref}_couplings
    global grid_mesh_parts_separate
    
    puts -nonewline "Writing domain file ..."
    
    set dir [file dirname $file]
    if ![file exists $dir] { cmd_mkdir $dir }
    
    set df [domainfile -name $file -mode w -type unstructured]
    
    set ns [set ${pref}_nodeset]
    set uns [set ${pref}_unstruct]
    set data [$ns get_nodes_ptr 0]
    set nnodes [lindex $data 1]
    $df write_nodes $nnodes [lindex $data 0]
    
    set new_startno($ns) 0
    set new_st [lindex $data 1]
    foreach eset [$uns list_sets] {
	set data [$eset get_elements_ptr 0]
	set nel [lindex $data 0]
	if $nel {
	    set new_startno($eset) $new_st
	    incr new_st $nel
	    $df write_elements [$eset get -type] $nel [lindex $data 1] $nnodes
	}
    }
    
    foreach dn [array names ${gpref}_group] {
	set dom [set ${gpref}_group($dn)]
	if ![$dom count] continue
	
	# We know this is an explicit domain.
	set new ""
	foreach elems [$dom list_elements 0] {
	    set set [lindex $elems 0]
	    set enums [lrange $elems 1 end]
	    if {$enums == ""} continue
	    
	    set st $new_startno($set)
	    foreach n $enums {
		lappend new [expr {$n + $st}]
	    }
	}
	
	if {$new != ""} {
	    $df write_list $dn 2 $new
	}
    }
    
    if !$grid_mesh_parts_separate {
	$df write_couplings_string [set ${pref}_couplings]
    }
    
    set ptr [$uns get_coupling_ptr]
    if {$ptr != 0} {
	$df write_couplings $ptr
    }
    
    foreach n [array names domain_family_pid] {
	if [string match DEF_* $n] continue
	$df set_family_name $domain_family_pid($n) $n
    }
    
    foreach n [array names domain_node_dimensions] {
	set lims $domain_node_dimensions($n)
	$df write_range $n 3 [lindex $lims 0] [lindex $lims 1]
    }
    
    puts " ok"
    $df close
    $df delete
}

proc domain_name_id_mapping {file list} {
    if {![file exists $file]} {
	mess "$file doesn't exist\n" red
	return 0
    }
    if [file isdirectory $file] {
	mess "$file is a directory\n" red
	return 0
    }
    
    if [catch {domainfile -name $file -mode r} df] {
	mess "Error: the mesher did not create a valid mesh\n" red
	mess "$df\n" red
	return 0
    }
    if {[$df type] != "unstructured"} {
	error "Wrong type of grid file!"
    }
    
    set fp [open $list w]
    
    set ndomains [$df ndomains]
    for {set d 0} {$d < $ndomains} {incr d} {
	set dname [$df dname $d]
	set utype [$df uns_type $d]
	
	switch $utype {
	    by_pid {
		puts $fp [list $dname [$df pid_for_name $dname]]
	    }
	}
    }
    close $fp
    
    $df close
    $df delete
}

proc domain_check_structured_couplings {file} {
    if {![file exists $file]} {
	mess "$file doesn't exist\n" red
	return 0
    }
    if [file isdirectory $file] {
	mess "$file is a directory\n" red
	return 0
    }
    
    if [catch {domainfile -name $file -mode r} df] {
	mess "Error: the mesher did not create a valid mesh\n" red
	mess "$df\n" red
	return 0
    }
    if {[$df type] != "unstructured"} {
	error "Wrong type of grid file!"
    }
    
    set nums [$df num_struct_couplings]
    
    $df close
    $df delete
    return $nums
}

