# 
# shockinfo is a function which takes as argument the result of a
# PDEPlot calculation, and a jump condition, to calculate the paths
# of any shocks which occur in the solution.
#
# These paths are returned as 3D plot curves, so that they can be
# superimposed on the plot of the solution.
#
# The new solution is also returned, with characteristics ending at
# the shocks.
#
# The jump condition is a function with four arguments:
# the values of the solution on both sides of the shock, x and t. 
# It must return the derivative of the shock path. Example:
#   jumpcond := (u1,u2,x,t)->(u1^2 - u2^2)/(u1-u2);
#
# Assumptions about the characterics:
#    - initial locations of the characteristics are monotonically increasing
#      (this can be relaxed by changing the way the algorithm finds the next char
#	that can collide with a shock, as well as some other small changes)
#    - time steps are equal for all characteristics and all times (this can 
#       be relaxed by evaluating the "interval" variable separately for each
#	calculation)
#
#
# by Josef Drexler jdrexler@julian.uwo.ca
#

shockinfo := proc(solution,jumpcond)

local mesh,charset,nchars,ntimes,charlist,shocklist,nshocks,breakingtime,
	curtime,mintime,eventt,i,j,k,l,char1,char2,delta,getcharinfo,
	collt,collx,eventc,eventx,leftval,rightval,interval,shocknum,
	slope,spath,sinfo,lininterp,timenow,gettime,shockof,getcharloc,
	findcharcharcollision,findcharshockcollision,findshockshockcollision;


# OK, now for the actual algorithm. 
# Repeat these steps for all time intervals in the solution:
# 1) Find out if any of the following events have occured
#    a) Two (or more) characteristics collide
#    b) One (or more) characteristic collides with a shock
#    c) Two (or more) shocks merge into a single shock
#    All of the events must be treated in sequence, because a) can
#    create a new shock, and c) remove one, so that it's imperative
#    to have the right order.  This is accomplished here by just finding
#    which event occurs first, and treats this one only, and then 
#    see if after this there are any events remaining in the interval.
#
#    The collision detection works by a linear interpolation in the 
#    time interval, and checks whether the two lines intersect
#
# 2) Treat the event as appropriate
#    a)	Mark the first and last characteristics that collide. These might
#	not be adjacent, if several characteristics collide at the same
#	location at the same time. Remove these characteristics from all 
#	following calculations - they have ended at the shock.
#	Then, create a new shock with the data from these chars
#    b)	Remove the chars as well, and modify the shock path
#    c)	Remove all but one shock, join the chars into the remaining
#	shock and calculate the new shock path


# Representation of shocks internally
#
# The information about all shocks is stored in two variables:
#
# spath is a list of 3D paths of the shocks, with a new point being added
#    every time the slope of the shock path changes.  Therefore this is not
#    usually an equally spaced (in time) curve.
#
# sinfo is a list of properties for each shock. For shock number j we have
#    sinfo[j][1] is the set of characteristics that have ended at this shock
#    sinfo[j][2] is the most recent slope of the shock path
#    sinfo[j][3] is the value of the solution on the "left" side of the shock
#    sinfo[j][4] is the value on the "right" side.
#    Here, "left" means the side with lower-numbered characteristics, conversely
#    "right" is the side where characteristics have higher numbers.
#    [2] is calculated from the jump condition, and the values in [3] and [4].
#    [3] and [4] are calculated by interpolating the outer-most characteristics
#    that have collided with the shock to the shock path on either side.


# Function output
#
# Depending on the setting of infolevel[shockinfo], shockinfo will return the 
# following information:
#    0  No output
#    1  The breaking time (first occurence of a shock) will be printed
#    2  Show each time step (might be useful for long calculations) as well
#       as the time and location of all new shocks
#    3  Show information about all collisions
#    4..7  Increasingly verbose information about the collision detection
#

#  ================================
#
#     Supporting functions
#
#  ================================

# a simple function for linear interpolation
# interpolates the char (two elements of a MESH for consecutive times) to
# the given time t.  For ind=1 it returns the location, for ind=3 the solution.

lininterp := (char,ind,t)->t/interval*(char[2][ind]-char[1][ind])+char[1][ind];


# Return two time steps of the mesh

getcharloc := proc(i,t)
[mesh[i,t],mesh[i,t+1]];
end:


# find the value of the time at a certain time step. Some chars may be undefined
# which makes this nontrivial

gettime := proc(timestep)
local i,t;
for i from 1 to nchars do
	t:=mesh[i,timestep][2];
	if (t<>undefined) then break; fi:
od:
t;
end:


# find out which shock a char is part of. This is only valid for
# chars that have ended at a shock.  In addition, no char should
# be part of more than one shock.

shockof := proc(charnum,shocklist)
local shocknum,k;
shocknum := -1;
for k in shocklist do
	userinfo(7,shockinfo,`Shock`,k,`has chars`,sinfo[k][1]);
	if member(charnum,sinfo[k][1]) then
		if (shocknum<>-1) then
			ERROR(`Char`,charnum,`has no unique shock.`);
		fi:
		shocknum := k;
	fi:
od:
if (shocknum=-1) then
	for k in shocklist do
		userinfo(6,shockinfo,`Shock`,k,`has chars`,sinfo[k][1]);
	od:
	ERROR(`Char`,charnum,`has no associated shock.`);
fi:
shocknum;
end:


# find out the location of a char. If it is now part of a shock,
# take the data of the shock instead.  For shockval=3 it returns
# the solution on the "left" side of the shock (see above), and
# for shockval=4 the solution on the right side.

getcharinfo := proc(charnum,shocklist,shockval)
local shocknum;
if member(charnum,charlist) then
	getcharloc(charnum,curtime);		# it's still a valid char
else
	userinfo(6,shockinfo,`Char`,charnum,`has ended, using shock info`);

	shocknum := shockof(charnum,shocklist);
	# return the shock data
	[ [eventx, gettime(curtime), sinfo[shocknum][shockval] ],
	  [eventx, gettime(curtime+1), sinfo[shocknum][shockval] ] ];
fi:
end;


# The following aren't really functions, but rather procedures still
# using the variables of shockinfo.  They are separate functions to
# make the code less confusing.


#  ============================
#
#  Char/char collision detection
#
#  ============================

findcharcharcollision := proc()
local lastvalid;

lastvalid:=-1:	# the most recent valid characteristic found so far

# iterate through all remaining characteristics 
for i from 1 to nchars-1 do
	userinfo(6,shockinfo,`Is char`,i,`alive:`,member(i,charlist));
	if not member(i,charlist) then next; fi;  # char has ended

	char1:=getcharloc(i,curtime);
	userinfo(6,shockinfo,`Char1 is at`,char1[1][1],char1[2][1]);

	if ((char1[1][1]=undefined) or (char1[2][1]=undefined) ) then
		userinfo(6,shockinfo,`Char is undefined.`);
		next;	   # the char is undefined at this point (e.g. out of range)
	fi;

	j:=lastvalid;
	lastvalid:=i;

	userinfo(5,shockinfo,`The two chars are`,i,j);
	if (j=-1) then next; fi;	# it's the first valid char, but we
					# need two for a collision

	char2:=getcharloc(j,curtime);
	userinfo(6,shockinfo,`Char2 is at`,char2[1][1],char2[2][1]);

	# we have two valid chars, do they collide? 
	# Find out by linear interpolation

	delta:=(char1[2][1]-char1[1][1]-char2[2][1]+char2[1][1]);
	userinfo(6,shockinfo,`Delta is`,delta);	# difference in slopes
	if (delta=0) then next; fi;		# they are parallel

	collt:=(char2[1][1]-char1[1][1])*interval/delta;	# collision time
	collx:=lininterp(char1,1,collt);			# and coordinate

	userinfo(4,shockinfo,`Chars`,i,j,`collide at`,collt,collx,mintime);

	# chars that intersect at t=0 are presumed to be in a fan and not a shock
	if ((collt=0) and (char1[1][2]=0)) then next; fi;

	# if collt is in [0,interval) they collide in the current interval at collx
	if (collt<0) then next; fi;		# have already collided
	if (collt>mintime) then next; fi;	# another collision happens earlier

	# if both time and location are the same, add the current
	# char to that event. This happens if more than two chars
	# cross at the same time.
	if ( (collt=eventt) and (collx=eventx) ) then
		userinfo(3,shockinfo,`Chars`,i,j,`join this event.`);
		eventc:=eventc union {i,j};
		next;
	fi;

	userinfo(3,shockinfo,`Found char/char collision`,collt,collx,i,j);
	eventt:=collt;	# time of event
	eventx:=collx;	# location
	eventc:={i,j};	# set of chars involved
	mintime:=collt;

od:	# end of the char/char loop
end:


#  ============================
#
#  Shock/char collision detection
#
#  ============================

findcharshockcollision := proc()
local i,j,k,twochars;

# now look if there are new chars that cross existing shocks
for i from 1 to nshocks do
	userinfo(6,shockinfo,`Is shock`,i,`alive:`,member(i,shocklist));
	if not member(i,shocklist) then next; fi;  # shock has ended

	# find the set of two chars that might intersect with the shock
	# these are the ones adjacent to the first and last char in the shock
	twochars:=[min(op(sinfo[i][1]))-1,max(op(sinfo[i][1]))+1];

	# check on both sides of the shock, whether the next char has hit
	for j from 1 to 2 do:
		if not member(twochars[j],charlist) then next; fi;	# char not valid
		char1:=getcharloc(twochars[j],curtime);
		userinfo(6,shockinfo,`Shock char`,j,`is`,twochars[j],
					`at`,char1[1][1],char1[2][1]);

		if ((char1[1][1]=undefined) or (char1[2][1]=undefined) ) then
			userinfo(6,shockinfo,`Char is undefined.`);
			next;   # the char is undefined at this point (e.g. out of range)
		fi;

		# Find out by linear interpolation

		delta:=(char1[2][1]-char1[1][1])/interval;	# slope of the char
		userinfo(6,shockinfo,`Char slope is`,delta);
		if (delta=sinfo[i][2]) then next; fi;		# they are parallel

		k:=nops(spath[i]);	# k=last shock path segment
		userinfo(6,shockinfo,`Shock segment`,k,` has x,t,m = `,
			spath[i][k][1],spath[i][k][2],sinfo[i][2]);

		collt:=(spath[i][k][1]-char1[1][1]
			-sinfo[i][2]*(spath[i][k][2]-char1[1][2]))
			/(delta-sinfo[i][2]);
		collx:=lininterp(char1,1,collt);

		userinfo(4,shockinfo,`Char`,twochars[j],`collides with shock`,i,
				`at`,collt,collx,mintime);

		if (collt<0) then next; fi;		# have already collided
		if (collt>mintime) then next; fi;	# another event happens earlier

		# if both time and location are the same, add the current
		# char to that event. This happens if more than two chars
		# cross at the same time.
		if ( (collt=eventt) and (collx=eventx) ) then
			userinfo(3,shockinfo,`Shock`,i,`and char`,twochars[j],
						`join this event.`);
			eventc:=eventc union {-i,twochars[j]};
			next;
		fi;

		userinfo(3,shockinfo,`Found char/shock collision`,
					collt,collx,-i,twochars[j]);
		eventt:=collt;	# time of event
		eventx:=collx;	# location
		eventc:={-i,twochars[j]};	# set of chars involved
		mintime:=collt;

	od:	# end of loop over two sides
od:	# end of shock/char loop
end:

#  ============================
#
#  Shock/shock collision detection
#
#  ============================

findshockshockcollision := proc()

for i in shocklist do
    for j in {$i+1..nshocks} intersect shocklist do
	userinfo(7,shockinfo,`Do shocks`,i,j,`collide?`);

	# Find out by linear interpolation

	if (sinfo[i][2]=sinfo[j][2]) then next; fi;  # they are parallel

	k:=nops(spath[i]);	# last shock path segment
	userinfo(6,shockinfo,`Shock`,i,`segment`,k,` has x,t,m = `,
		spath[i][k][1],spath[i][k][2],sinfo[i][2]);

	l:=nops(spath[j]);
	userinfo(6,shockinfo,`Shock`,j,`segment`,l,` has x,t,m = `,
		spath[j][l][1],spath[j][l][2],sinfo[j][2]);

	collt:=(spath[j][l][1]-spath[i][k][1]
		+sinfo[i][2]*spath[i][k][2]-sinfo[j][2]*spath[j][l][2])
			/(sinfo[i][2]-sinfo[j][2]);
	collx:=spath[i][k][1]+sinfo[i][2]*(collt-spath[i][k][2]);

		# convert to interval time...
	collt:=collt-timenow;

	userinfo(4,shockinfo,`Shocks`,i,j,`collide at`,
		collt,collx,mintime);

	if (collt<0) then next; fi;		# have already collided
	if (collt>mintime) then next; fi;	# another event happens earlier

	# if both time and location are the same, add the current
	# shocks to that event
	if ( (collt=eventt) and (collx=eventx) ) then
		userinfo(3,shockinfo,`Shocks`,i,j,`join this event.`);
		eventc:=eventc union {-i,-j};
		next;
	fi;

	userinfo(3,shockinfo,`Found shock/shock collision`,collt,collx,-i,-j);
	eventt:=collt;		# time of event
	eventx:=collx;		# location
	eventc:={-i,-j};	# set of shocks involved
	mintime:=collt;
    od:
od:	# end of shock/shock loop
end:

#  ==================================
#
#    End of supporting functions
#
#  ==================================




#  ==================================
#
#    Shock finding algorithm
#
#  ==================================


# Extract the actual characteristic curves from the solution
mesh:=op(1,op(1,solution)):
nchars:=nops(mesh);	# number of characteristics
ntimes:=nops(mesh[1]);	# number of time steps

# Convert it to an array (a list doesn't work well for many characteristics)
# (it's now an array of lists for each characteristic)
mesh:=array(1..nchars,1..ntimes,mesh);

if (nops(solution)=6+nchars) then	# PDEPlot with basechar=true option 
	charset:=array(1..nchars,1..ntimes-1);
	for i from 1 to nchars do
		k := [op(op(i+2,solution))];
		for j from 1 to ntimes-1 do
			charset[i,j]:=k[j];
		od:
	od:
		# solution element 1 is mesh, 2 is IC, 3.. is base chars
else
	charset:=undefined;
fi:

interval := gettime(2) - gettime(1);

charlist:={$1..nchars}:		# set of "surviving" characteristics, i.e.
				# those that have not ended at a shock.
shocklist:={}:			# the same thing for shocks
nshocks:=0:			# total number of shocks found so far
breakingtime:=undefined:	# will contain the first occurence of a shock

spath := [];

userinfo(3,shockinfo,`I have`,nchars,`characteristics and`,ntimes,`time steps`);

# Now do the time loop. Iterate only to ntimes-1 because we consider
# intervals. The inner `while' loop is over all events in the interval, and
# is terminated explicitly if no more events are found.

for curtime from 1 to ntimes-1 do while true do

	# find current time
	timenow := gettime(curtime);
	userinfo(2,shockinfo,`The time is now`,curtime,timenow);
	userinfo(5,shockinfo,`The surviving characteristics are`,charlist);
	userinfo(5,shockinfo,`The active shocks are`,shocklist);

	# continue the dead chars along their shocks with undefined z
	i := {$1..nchars} minus charlist;
	for j in shocklist do
		mintime := gettime(curtime+1);
		char1 := spath[j][nops(spath[j])];
		k := {min(op(sinfo[j][1])),max(op(sinfo[j][1]))};
		i := i minus k;
		for l from 1 to 2 do
			mesh[k[l],curtime+1][1]:=
				char1[1]+sinfo[j][2]*(mintime-char1[2]);
			mesh[k[l],curtime+1][3]:=sinfo[j][2+l];
		od:
	od:
	for l in i do
		mesh[l,curtime+1][1]:=undefined;
		mesh[l,curtime+1][3]:=undefined;
	od:

	mintime:=interval:
			# first occurence of an intersection in this interval
			# note that all times are given as fraction of 
			# the time interval

	eventt:=-1:	# the time of the event which occurs next

	# Go through all possible events and find the earliest one

	findcharcharcollision();
	findcharshockcollision();
	findshockshockcollision();


	#  ============================
	#
	#  Event processing
	#
	#  ============================

	if (eventt=-1) then break; fi;	# no more events found, terminate loop
	shocknum:=-1;
	mintime:=eventt+char1[1][2];

	userinfo(4,shockinfo,`Chars (+ve) and shocks (-ve) in the event:`,eventc);
	i := select(a->a<0,eventc);	# shocks are all negative elements
	if (i<>{}) then
		# shocks are involved
		userinfo(4,shockinfo,`The event has shocks`);

		# what to do:
		# - find leftmost and rightmost chars from all shocks combined
		# - remove all but one shock
		# - replace this shock by the one generated from the two chars

		# i has shock numbers with negative sign. Correct this.
		j := {};
		for k in i do
			j := j union {-k};
		od:
		k := i;
		i := j;
		j := eventc minus k;	# j: the chars from the list
		eventc := j;

		for k in i do
			userinfo(6,shockinfo,`Adding chars `,sinfo[k][1]);
			eventc:=eventc union sinfo[k][1];
		od:
		userinfo(5,shockinfo,`All contributing chars are`,eventc);

		char1 := min(op(eventc));
		char2 := max(op(eventc));

		userinfo(6,shockinfo,`The two important chars are`,char1,char2);

		char1 := getcharinfo(char1,i,3);
		char2 := getcharinfo(char2,i,4);
		userinfo(6,shockinfo,`Char1 is `,char1);
		userinfo(6,shockinfo,`Char2 is `,char2);

		# remove all but first shock
		shocknum := i[1];
		i := i minus {shocknum};

		for k in i do
			spath[k] := [op(spath[k]),[eventx,mintime,0]]
		od:
		shocklist := shocklist minus i;
	else
		# only chars -> generate new shock
		userinfo(2,shockinfo,`Starting new shock at`,eventx,eventt);

		nshocks:=nshocks+1;
		shocknum:=nshocks;

		shocklist := shocklist union {nshocks};

		char1 := getcharloc(min(op(eventc)),curtime);
		char2 := getcharloc(max(op(eventc)),curtime);

		spath := [op(spath),[]];
	fi:

	# generate the new shock data
	if (breakingtime=undefined) then
		breakingtime:=mintime;
		userinfo(1,shockinfo,`The breaking time is`,breakingtime);
	fi;

	# find the values by linear interpolation
	leftval := lininterp(char1,3,eventt);
	rightval := lininterp(char2,3,eventt);

	userinfo(5,shockinfo,`Jump at shock is`,leftval,rightval,
			`at x`,eventx,`t`,mintime);

	slope := jumpcond(leftval,rightval,eventx,mintime);

	userinfo(4,shockinfo,`Jump at shock is`,leftval,rightval,
		`and the slope is`,slope);
	userinfo(5,shockinfo,`The char slopes are`,
		(char1[2][1]-char1[1][1])/interval,
		(char2[2][1]-char2[1][1])/interval);

	# spath is an irregular 3D curve of the shock paths
	spath[shocknum] := [op(spath[shocknum]),[eventx,mintime,0]];

	# sinfo is [chars, slope, solution left, solution right];
	sinfo[shocknum]:=[eventc, slope, leftval, rightval];

	# if basechars were in the PDEPlot result, end the chars at the shock
	if (charset<>undefined) then
		for i in eventc intersect charlist do
			charset[i,curtime][2]:=[eventx,mintime,0];

			for j from curtime+1 to ntimes-1 do
				charset[i,j]:=
				    [[undefined,undefined,0],[undefined,undefined,0]];
			od:
		od:
	fi:

	# remove chars from further consideration.
	charlist := charlist minus eventc;

	# end the chars at the shock
	mesh[min(op(eventc)),curtime+1]:=[eventx,mintime,leftval];
	mesh[max(op(eventc)),curtime+1]:=[eventx,mintime,rightval];
	for i in eventc minus {min(op(eventc)),max(op(eventc))} do
		j := lininterp(getcharloc(i,curtime),3,mintime);
		mesh[i,curtime+1]:=[eventx,mintime,j];
	od:

	# repeat this until nothing else happens
	od:	# end of event loop

od:	# end of the time loop


# continue all shocks to the end of the last interval
timenow := gettime(ntimes);

for i in shocklist do
	char1 := spath[i][nops(spath[i])];
	j := char1[1] + (timenow - char1[2]) * sinfo[i][2];
	userinfo(6,shockinfo,`Continuing shock`,i,`to`,[j,timenow,0]);
	spath[i] := [op(spath[i]),[j,timenow,0]];
od:

# return the following list:
# - the list of shockpaths
# - new PLOT3D structure, with the modified MESH as well as new characteristic
#   curves if basechar=true was given to PDEPlot

mesh:=convert(mesh,listlist);

if (charset<>undefined) then
	# convert charset to a list of CURVES
	charset:=array(1..nchars,convert(charset,listlist));
	for i from 1 to nchars do
		charset[i]:=CURVES(op(charset[i]));
	od:
	charset:=convert(charset,listlist);
	[spath,PLOT3D(MESH(mesh),op(2,solution),op(charset),op(2+nchars..nops(solution),solution))];
else
	[spath,PLOT3D(MESH(mesh),op(2..nops(solution),solution))];
fi;

end: