/*
 $Id: materialstr.cc,v 1.4 1996/11/19 10:02:53 bzferdma Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "materialstr.h"

#include "cmdpars.h"
extern CmdPars Cmd;


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

UserTransMaterial:: UserTransMaterial(const char* fileName, int spaceDim0)
	: DefaultMaterial(fileName, spaceDim0)
{ 
    setMaterialType();
}

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


StepSource:: StepSource(const char* fileName, int spaceDim0)
	: DefaultMaterial(fileName, spaceDim0)
{ 
    setMaterialType();
}


Real StepSource:: S(int type, Vector<Real>* /*x*/, Real time)
{
    if (type == 1 && time <= 10.0) return 1.0;
    else   return 0.0;
}
//-------------------------------------------------------------------------


TransPeakSource:: TransPeakSource(const char* fileName, int spaceDim0)

	:  PeakSource(fileName, spaceDim0)
{ 
    setMaterialType();
}


Real TransPeakSource:: S(int /*type*/, Vector<Real>* x, Real time)
{
    Vector<Real>& xx = *x;
    Real val;

    switch(spaceDim)
    {
      case 1: val = S1d(xx[1],time); return val;
      case 2: val = S2d(xx[1], xx[2],time);  return val;
      case 3: val = S3d(xx[1], xx[2], xx[3],time);  return val;
      default: cout.flush(); abort(); return 0;
    }
}


Real TransPeakSource:: S1d(Real x, Real time)
{
    Real xm, t = time, x0;
    Real f, fx, fxx, ft;
    Real e, ex, exx;
    Real uxx, p=100.0;
    Real pi = 3.141592653589;

    x0 = 0.5 + 0.375*sin(2.0*pi*t/20.0);
    xm = x-x0;
    
    e   = exp(-p*(xm*xm));
    ex  = -2.0*p*xm*e;
    exx = -2.0*p*e + (-2.0*p*xm)*ex;
    
    f   = x*(x-1.);
    fx  = (2.0*x-1.0);
    fxx = 2.0;
    
    ft = f*ex*(pi/10.0*cos(pi*t/10.0));
    uxx = ft - ( fxx*e + 2.0*fx*ex + f*exx );
     

    return uxx;
}
  
Real TransPeakSource:: S2d(Real x, Real y, Real time)
{
    Real xm, ym;
    Real f, fx, fy, fxx, fyy, ft;
    Real e, ex, exx, ey, eyy;
    Real uxx, uyy, p=100.0;
 
    xm=x-0.5;
    ym=y-0.5;

    Real r = 0.375, t = time;
    Real pi = 3.1414;
    xm = 0.5 + r * cos(2.0*pi*t/20.0);
    ym = 0.5 + r * sin(2.0*pi*t/20.0);
    xm = x-xm;
    ym = y-ym;

    e   = exp(-p*(xm*xm+ym*ym));
    ex  = -2.0*p*xm*e;
    exx = -2.0*p*e + (-2.0*p*xm)*ex;
    ey  = -2.0*p*ym*e;
    eyy = -2.0*p*e + (-2.0*p*ym)*ey;

    f   = x*(x-1.) * y*(y-1.);
    fx  = (2.0*x-1.0)*y*(y-1.);
    fxx = 2.0*y*(y-1.);
    fy  = (2.0*y-1.0)*x*(x-1.);
    fyy = 2.0*x*(x-1.);

    uxx = fxx*e + 2.0*fx*ex + f*exx;
    uyy = fyy*e + 2.0*fy*ey + f*eyy;

    ft = -2.0/10.0*p*e*f*r*pi*( - xm*sin(pi*t/10.0) + ym*cos(pi*t/10.0) );
    return (ft - uxx - uyy);
}

/*   constant peak
Real TransPeakSource:: S2d(Real x, Real y, Real time)
{
    Real xm, ym;
    Real f, fx, fy, fxx, fyy;
    Real e, ex, exx, ey, eyy;
    Real uxx, uyy, p=100.0;
 
    xm=x-0.5;
    ym=y-0.5;

    e   = exp(-p*(xm*xm+ym*ym));
    ex  = -2.0*p*xm*e;
    exx = -2.0*p*e + (-2.0*p*xm)*ex;
    ey  = -2.0*p*ym*e;
    eyy = -2.0*p*e + (-2.0*p*ym)*ey;

    f   = x*(x-1.) * y*(y-1.);
    fx  = (2.0*x-1.0)*y*(y-1.);
    fxx = 2.0*y*(y-1.);
    fy  = (2.0*y-1.0)*x*(x-1.);
    fyy = 2.0*x*(x-1.);

    uxx = fxx*e + 2.0*fx*ex + f*exx;
    uyy = fyy*e + 2.0*fy*ey + f*eyy;

    return (uxx + uyy);
}
*/

Real TransPeakSource:: S3d(Real x, Real y, Real z, Real /*time*/)
{
   Real a, b, c, p;
   Real xm1,ym1,zm1,xm1x,ym1y,zm1z,xma,ymb,zmc;
   Real sum,divisor;

   a = 0.5;
   b = 0.5;
   c = 0.5;
   p = 150.0;

   xma  = x-a;
   ymb  = y-b;
   zmc  = z-c;

   xm1  = x-1.0;
   ym1  = y-1.0;
   zm1  = z-1.0;
   xm1x = x*xm1;
   ym1y = y*ym1;
   zm1z = z*zm1;

   divisor = exp(p*(xma*xma + ymb*ymb + zmc*zmc));

   sum =       2.0*(xm1x*ym1y+xm1x*zm1z+ym1y*zm1z);
   sum = sum - 6.0*p*(xm1x*ym1y*zm1z);
   sum = sum - 4.0*p*(xm1*xma*ym1y*zm1z + x*xma*ym1y*zm1z);
   sum = sum + 4.0*p*p*(xm1x*xma*xma*ym1y*zm1z);
   sum = sum - 4.0*p*(xm1x*ym1*ymb*zm1z + xm1x*y*ymb*zm1z);
   sum = sum + 4.0*p*p*(xm1x*ym1y*ymb*ymb*zm1z);
   sum = sum - 4.0*p*(xm1x*ym1y*zm1*zmc + xm1x*ym1y*z*zmc);
   sum = sum + 4.0*p*p*(xm1x*ym1y*zm1z*zmc*zmc);

    return -100.0*sum/divisor;
}
//-------------------------------------------------------------------------

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

CastingMaterial:: CastingMaterial(const char* fileName, int spaceDim0)

	: DefaultMaterial(fileName, spaceDim0)
{
   qMold  = 1.0;   Cmd.get("qMold",&qMold);
   qSpray = 0.5;  Cmd.get("qSpray",&qSpray);

   setMaterialType();
}



Real CastingMaterial:: Neumann(int /*type*/, Vector<Real>* /*x*/, Real time)
{
    switch(spaceDim)
    {
      case 1:  cout.flush(); abort(); return 0.0;
      case 2:  if ( time < 12.5  )    return -qMold*0.142;
               else                   return -0.1403*qSpray; // - 0.0328;
      case 3:  cout.flush(); abort(); return 0;
      default: cout.flush(); abort(); return 0;
    }
}

Real CastingMaterial:: Cauchy(int /*type*/, Vector<Real>* /*x*/, Real time)
{
    switch(spaceDim)
    {
      case 1:  cout.flush(); abort(); return 0;
      case 2:  if ( time < 12.5  ) return qMold * 2.0/10.0; // U <= 0 recommended!
               else                return qSpray*2.0/10.0;
               // if ( (time > 12.5) && (time < 121) ) return 1.38*2.0/7.0;
      case 3:  cout.flush(); abort(); return 0;
      default: cout.flush(); abort(); return 0;
    }
}

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


StefanSource:: StefanSource(const char* fileName, int spaceDim0)

	:  VarSourceMaterial(fileName, spaceDim0)
{ 
    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("StefanTheta1", &theta1);

    setMaterialType();
}
//-------------------------------------------------------------------------


Real StefanSource:: S(int /*type*/, Vector<Real>* x, Real time)
{
    Vector<Real>& xx = *x;

    switch(spaceDim)
    {
      case 1: cout.flush(); abort(); return 0;;
      case 2: return S2d(xx[1], xx[2], time);
      case 3: cout.flush(); abort(); return 0;
      default: cout.flush(); abort(); return 0;
    }
//    return 0;
}
//-------------------------------------------------------------------------


Real StefanSource:: S2d(Real x, Real y, Real t)
{
    Real f=0.0, theta;
    
    theta = (x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)-exp(-4.*t)/4.;
    
    if (theta <  theta1) f = c0    *exp(-4.*t)    - 4*kappa0;
    if (theta == theta1) f =(c0+c1)*exp(-4.*t)/2. - 2*(kappa0+kappa1);
    if (theta >  theta1) f = c1    *exp(-4.*t)    - 4*kappa1;
    
    return f;
}