/*
 $Id: nonlin.cc,v 1.2 1996/10/04 15:07:10 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "nonlin.h"
#include "numerics.h"

#include "int.h"
#include "nodeco.h"

#include "cmdpars.h"
extern CmdPars Cmd;

//-------------------------------------------------------------------------

static const Real MachMax = machMax(Real(0.0));

//-------------------------------------------------------------------------

NonLinearity:: NonLinearity() : u(1), H(1), areaWeight(1), 
				uppObstacle(1), lowObstacle(1)
{ } 

NonLinearity:: ~NonLinearity() { }

//-------------------------------------------------------------------------

Real NonLinearity:: lowerObstacle(Vector<Real>& /*x*/) const { return -MachMax; }
Real NonLinearity:: upperObstacle(Vector<Real>& /*x*/) const { return  MachMax; }

//-------------------------------------------------------------------------

void NonLinearity:: addAreaWeight(const Vector<Real>& LNorm,
				  const Vector<int>& nodes)
{
    int i;
    FORALL(nodes,i) { if (nodes[i]) areaWeight[nodes[i]] += LNorm[i]; }
}
//-------------------------------------------------------------------------


void NonLinearity:: update(const Interface& interface)
{
    int i;

    areaWeight.resize(interface.Dim());
    FORALL(areaWeight,i)  areaWeight[i] = 0.0;

    updateObstacles(uppObstacle, lowObstacle, interface);
}
//-------------------------------------------------------------------------


void NonLinearity:: defaultObstacleUpdate(Vector<Real>& uppO,Vector<Real>& lowO, 
					  const Interface& interface) const
{
    int i;

    uppO.resize(interface.Dim());
    lowO.resize(interface.Dim());

    FORALL(uppO,i)
    {
	lowO[i] = -MachMax;
	uppO[i] =  MachMax;
    }
}
//-------------------------------------------------------------------------


void NonLinearity:: fctObstacleUpdate(Vector<Real>& uppO, Vector<Real>& lowO, 
				      const Interface& interface) const
{
    int i;
    const int dim 	   = interface.Dim();
    const int spaceDim = interface.SpaceDim();

    uppO.resize(dim);
    lowO.resize(dim);

    Vector<Real> x(spaceDim);
    NodeCoordinates nc(spaceDim, 1, dim);

    interface.getNodeCoordinates(nc);

    FORALL(uppO,i)
    {
	nc.getCoordinate(i,x);
	lowO[i] = lowerObstacle(x);
	uppO[i] = upperObstacle(x);
   }
}
//-------------------------------------------------------------------------


Real NonLinearity:: HValue(int node, Real uVal)  const
{
    if (node <= lowObstacle.high())
    {
	if (uVal < lowObstacle[node])
	{
	    cout << "\n** NonLinearity::HValue(...): value " << uVal
	         << " of node " << node << " below lower obstacle: " 
	    	 << lowObstacle[node] << "\n";
	    //cout.flush(); abort();
	}
	if (uVal > uppObstacle[node])
	{
	    cout << "\n** NonLinearity::HValue(...): value " << uVal
	    	 << " of node " << node << " above upper obstacle: " 
	    	 << uppObstacle[node] << "\n";
	    //cout.flush(); abort();
	}
    }

    return HValue(uVal);
}
//-------------------------------------------------------------------------


Real NonLinearity:: HValue(Real uVal)  const
{
    int i;

    if      (uVal <= u[1]) return linearInterp(uVal, u[1], u[2], H[1], H[2]);
    else if (uVal >= u[M]) return linearInterp(uVal, u[M-1], u[M], H[M-1], H[M]);
    else
    {
	for (i=2; i<=M; ++i)
	{
	    if (uVal <= u[i])
	    {
		if (equal(uVal,u[i]) && equal(u[i-1],u[i]))
		   return 0.5 * (H[i-1] + H[i]);
		
		else return linearInterp(uVal, u[i-1], u[i], H[i-1], H[i]);
	    }
	}
    }

    cout << "\n*** NonLinearity::HValue(...): uVal " << uVal << " not located\n";
    cout.flush(); abort();
    return 0.0;
}
//-------------------------------------------------------------------------

// --		 interpolate H = H(u),  H1 = H(u1), H2 = H(u2)

Real NonLinearity:: linearInterp(Real ui, Real u1, Real u2, Real H1, Real H2) 
									   const
{
    Real x = (ui-u1)/(u2-u1);	
    return H1 + x*(H2-H1);
}
//-------------------------------------------------------------------------


Real PorousMedia:: HValue(Real uVal)  const
{
    int i;

    if (uVal <= u[1])
    {
	if (uVal > 0.0)  return pow(uVal, 1.0/mm);
	else  return 0.0;
    }
    
    else if (uVal >= u[M]) return linearInterp(uVal, u[M-1], u[M], H[M-1], H[M]);
    else
    {
	for (i=2; i<=M; ++i)
	{
	    if (uVal <= u[i])
	    {
		if (equal(uVal,u[i]) && equal(u[i-1],u[i]))
		   return 0.5 * (H[i-1] + H[i]);
		
		else return linearInterp(uVal, u[i-1], u[i], H[i-1], H[i]);
	    }
	}
    }

    cout << "\n*** NonLinearity::HValue(...): uVal " << uVal << " not located\n";
    cout.flush(); abort();
    return 0.0;
}
//-------------------------------------------------------------------------



Casting:: Casting() : NonLinearity()
{
    kappa0 = 7.0/20000.0;     	Cmd.get("kappa0", &kappa0);
    kappa1 = 7.0/4000.0;     	Cmd.get("kappa1", &kappa1);

    M=8;
    u.resize(M);
    H.resize(M);
/*
    u[1]= -0.315;   	       H[1]= 3.27;
    u[2]= -0.0105;   	       H[2]= 6.94;
    u[3]= -0.0105;   	       H[3]= 7.15;
    u[4]= -0.007;   	       H[4]= 7.30;
    u[5]= -0.007;   	       H[5]= 7.64;
    u[6]=  0.0;   	       H[6]= 7.92;
    u[7]=  0.0;   	       H[7]= 9.44;
    u[8]=  0.175;   	       H[8]= 9.84;
*/

    u[1]= -0.45;   	       H[1]= 3.27;
    u[2]= -0.015;   	       H[2]= 6.94;
    u[3]= -0.015;   	       H[3]= 7.15;
    u[4]= -0.01;   	       H[4]= 7.30;
    u[5]= -0.01;   	       H[5]= 7.64;
    u[6]=  0.0;   	       H[6]= 7.92;
    u[7]=  0.0;   	       H[7]= 9.44;
    u[8]=  0.175;   	       H[8]= 9.84;

}
//-------------------------------------------------------------------------

Stefan:: Stefan() : NonLinearity()
{
    kappa0 = 1.0;     	Cmd.get("kappa0", &kappa0);
    kappa1 = 1.0;     	Cmd.get("kappa1", &kappa1);
    c0 = 1.0;     	Cmd.get("c0", &c0);
    c1 = 1.0;     	Cmd.get("c1", &c1);
    s0 = 0.0;     	Cmd.get("s0", &s0);
    s1 = 0.0;     	Cmd.get("s1", &s1);
    theta1 = 0.0;     	Cmd.get("theta1", &theta1);


    Real a0 = c0 / kappa0; 
    Real a1 = c1 / kappa1;


    M=4;
    u.resize(M);
    H.resize(M);

    u[1]= theta1 - 1.0;   	H[1]= -(a0+s0);
    u[2]= theta1;   		H[2]= -s0;
    u[3]= theta1;   		H[3]=  s1;
    u[4]= theta1 + 1.0;   	H[4]=  a1+s1;

    cout << "\nStefan constructed\n";
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

Obstacle:: Obstacle() : NonLinearity()
{ 
    M=2;
    u.resize(M);
    H.resize(M);

    u[1]=-1.0;   H[1]= 0.0;
    u[2]= 0.0;   H[2]= 0.0; 
 
    cout << "\nObstacle constructed\n";
}

Real Obstacle:: lowerObstacle(Vector<Real>& /*xvec*/) const { return -MachMax; }

Real Obstacle:: upperObstacle(Vector<Real>& xvec) const
{
    Real x = xvec[1];
    Real y = xvec[2];

    if (y < 0.5)
    {
	if (x < 0.5) return (x<y)? x:y;
	else return ((1-x)<y)? (1-x):y;
    }
    else
    {
	if (x < 0.5) return (x<1.-y)? x:(1.-y);
	else return ((1-x)< (1.-y) )? (1.-x):(1.-y);
    }
    //return 1000.;
}
//-------------------------------------------------------------------------

Real Obstacle:: HValue(int /*node*/, Real /*uVal*/) const
{
    cout << "\n*** Obstacle::HValue(...): H makes no sense\n";
    cout.flush(); abort();
    return 0.0;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

PorousMedia:: PorousMedia() : NonLinearity(), critic(1)
{
    int i;
    Real hM, rightEnd=1.;
    mm = 2.;
    
    M=100;
    u.resize(M);
    H.resize(M);
    

    hM = pow(rightEnd,1./mm)/M;

    for (i=1; i<=M; i++)
    {
	u[i] = pow(i*hM,mm);
	H[i] = pow(u[i],1./mm);
    }
    
    u[1] = u[2]/2.;
    H[1] = pow(u[1],1./mm);
    
    cout << "\nPorousMedia constructed\n";
}
//-------------------------------------------------------------------------

Real PorousMedia:: lowerObstacle(Vector<Real>& /*xvec*/)  const { return 0.0; }

//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

void PorousMedia:: update(const Interface& interface) 
{ 
    int i;

    NonLinearity::update(interface); 

    critic.resize(interface.Dim());
    FORALL(critic,i) critic[i] = True;	   // to include dirichlet nodes
}
//-------------------------------------------------------------------------


void Obstacle:: defectCorrection(Real& e, Real& r, Real oldSol, 
				 Real diag, Real& a, Real& uppDefO, 
				 Real& lowDefO, Bool& critical, 
				 Real uppO, Real lowO, Real /*arWeight*/) const
{
    Real newSol = oldSol + r/diag;

    critical = False;

    if (newSol <= lowO)
    {
	newSol = lowO;
	lowDefO = 0.0;
	uppDefO = 0.0;
	critical = True;
    }
    else if (newSol >= uppO)
    {
	newSol = uppO;
	uppDefO = 0.0;
	lowDefO = 0.0;
	critical = True;
    }
    else
    {
	lowDefO = lowO-newSol;
	uppDefO = uppO-newSol;
    }

    a = 0.0;
    r = 0.0;
    e = newSol - oldSol;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void NonLinearity:: defectCorrection(int node, Real& e, Real& r, Real oldSol, 
				     Real diag) 
{
    Real a, uppDefO, lowDefO;		// dummies
    Bool critical = True;

    defectCorrection(e, r, oldSol, diag, a, uppDefO, lowDefO, critical,
		     uppObstacle[node], lowObstacle[node], areaWeight[node]);

    if (critical) setCritical(node);
    else 	  setUnCritical(node);
}
//-------------------------------------------------------------------------


void NonLinearity:: defectCorrection(int node, Real& e, Real& r, Real oldSol, 
				     Real diag, Real& a, Real& uppDefO, 
				     Real& lowDefO, Bool& critical)
{
    defectCorrection(e, r, oldSol, diag, a, uppDefO, lowDefO, critical,  
		     uppObstacle[node], lowObstacle[node], areaWeight[node]);

    if (critical) setCritical(node);
    else 	  setUnCritical(node);
}
//-------------------------------------------------------------------------


//  --  apply nonlinear contributions: 
//	-> defect correction e, residual correction r; 
//	for multi-grid: correction a on diagonal, obstacles, and critical flag
//


void NonLinearity:: defectCorrection(Real& e, Real& r, Real oldSol, 
				     Real diag, Real& a, Real& uppDefO, 
				     Real& lowDefO, Bool& critical, 
				     Real uppO, Real lowO, Real arWeight) const
{
    int  k;
    Real  newSol, b, u0, u1, H0, H1, hSquare, s=0.0;

    critical = False;

    a = 0.0;

    b = r + diag*oldSol;
    hSquare = arWeight;
	

    H0 = H[1]; u0 = u[1];
    H1 = H[2]; u1 = u[2];
	
    if (b < H1*hSquare+diag*u1)			// the first interval
    {
	a = (H1-H0)/(u1-u0);
	
	if (lowO > u0) cout << "Help!!!\n";

	u0 = lowO;
	H0 = a*(u0-u1) + H1;
	    
	if (b < H0*hSquare + diag*u0)
	{
	    a=0.; 
	    s=0.;
	    newSol = u0;
	    
	    lowDefO = 0.0;
	    uppDefO = 0.0;
	    critical = True;
	}
	else 
	{
	    a = a*hSquare; s = H1*hSquare-a*u1;
	    newSol = (b-s)/(diag+a);
	    
	    lowDefO = u0-newSol;
	    uppDefO = u1-newSol;
	}
	    
	e = newSol - oldSol;
    }
    else
    {
	for (k=2; k < M-1; k++)
	{
	    H0 = H[k];   u0 = u[k];
	    H1 = H[k+1]; u1 = u[k+1];
	    
	    if ((H0*hSquare + diag*u0 <= b)&&(b < H1*hSquare+diag*u1) )
	    {
		if (equal(u0,u1))
		{
		    a=0.; s=0.;
		    newSol = u0;
		    
		    lowDefO = 0.0;
		    uppDefO = 0.0;
		    critical = True;
		}
		else 
		{
		    a = (H0-H1)*hSquare/(u0-u1); s = H0*hSquare-a*u0;
		    newSol = (b-s)/(diag+a);
		    
		    lowDefO = u0 - newSol;
		    uppDefO = u1 - newSol;
		}
		    
		e = newSol - oldSol;
		    
		//if  (a < theta)
		//{
		//lowDefO = 0.0;
		//uppDefO = 0.0;
		//critical = True;
		//}
		break; 
	    } 
	}
    }
	
    H0 = H[M-1];   u0 = u[M-1];
    H1 = H[M];     u1 = u[M];
	

    if (H0*hSquare+diag*u0 <= b)	// the last interval
    {
	a = (H1-H0)/(u1-u0);
	
	if (uppO < u1) cout << "MayDay\n";

	u1 = uppO;
	H1 = a*(u1 - u0) + H0;

	if (H1*hSquare+diag*u1 <= b)
	{
	    a=0.; s=0.;
	    newSol = u1;
	    
	    lowDefO = 0.0;
	    uppDefO = 0.0;
	    critical = True;
	}
	else 
	{
	    a = a*hSquare; s = H0*hSquare-a*u0;
	    newSol = (b-s)/(diag+a);
	    
	    lowDefO = u0-newSol;
	    uppDefO = u1-newSol;
	}
	
	e = newSol - oldSol;
    }

    r = s + a*oldSol;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


NLTest:: NLTest() : NonLinearity()
{ 
    M=2;
    u.resize(M);
    H.resize(M);

    u[1]=-1.0;   H[1]= 0.0;
    u[2]= 0.0;   H[2]= 0.0; 
 
    m = 1000.0;
    cout << "\nNLTest constructed\n";
    Cmd.get("MValue", &m);
}
//-------------------------------------------------------------------------


void NLTest:: updateObstacles(Vector<Real>& uppO, Vector<Real>& lowO, 
			      const Interface& interface) const
{
    int i;

    uppO.resize(interface.Dim());
    lowO.resize(interface.Dim());

    FORALL(uppO,i) { lowO[i] = -10000.0;  uppO[i] = 10000.0; }
}
//-------------------------------------------------------------------------


Real NLTest:: HValue(int /*node*/, Real uVal) const { return HValue(uVal); }

Real NLTest:: HValue(Real uVal) const
{
    return m * uVal * uVal;
}
//-------------------------------------------------------------------------


void NLTest:: defectCorrection(Real& e, Real& r, Real oldSol, 
			       Real diag, Real& a, Real& uppDefO, 
			       Real& lowDefO, Bool& critical, 
			       Real uppO, Real lowO, Real arWeight) const
{
    Real newSol;
    Real dH = m*arWeight * oldSol;
    Real H0 = dH * oldSol;
    dH *= 2.0;

    newSol = oldSol + (r - H0)/(diag + dH);

    critical = False;

    if (newSol <= lowO)
    {
	newSol = lowO;
	lowDefO = 0.0;
	uppDefO = 0.0;
	critical = True;
    }
    else if (newSol >= uppO)
    {
	newSol = uppO;
	uppDefO = 0.0;
	lowDefO = 0.0;
	critical = True;
    }
    else
    {
	lowDefO = lowO-newSol;
	uppDefO = uppO-newSol;
    }

    e = newSol - oldSol;
    a = dH;
    r = H0;
}
//-------------------------------------------------------------------------