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

#include "dirichlettr.h"

#include "physics.h"

#include "utils.h"
#include "numerics.h"

#include "cmdpars.h"
extern CmdPars Cmd;

//-------------------------------------------------------------------------
ConstTransDirichlet:: ConstTransDirichlet(char* fileName)
                    : ConstDirichletBCs(fileName)
{ 
}

Real ConstTransDirichlet:: initialValue(Vector<Real>& /*x*/, Real /*time*/, int /*comp*/) const
{
    return 0.0;
}
//-------------------------------------------------------------------------

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

UserTransDirichlet:: UserTransDirichlet()
{ 
}

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

StepDirichlet:: StepDirichlet()
{ 
    xStep = 0.5;   	Cmd.get("xStep", &xStep);
}
//-------------------------------------------------------------------------

void StepDirichlet:: setBC(int node, int id, Vector<Real>& x, int /*comp*/, 
			   Real /*time*/)
{
    Id[node] = id;
    Values[node] = value(x);
}
//-------------------------------------------------------------------------
 
Real StepDirichlet:: value(Vector<Real>& x) const
{
    if (x[1] <= xStep) return 1.0;
    else  return 0.0;
}
//-------------------------------------------------------------------------

Real StepDirichlet:: initialValue(Vector<Real>& x, Real /*time*/, int /*comp*/) const
{
    return value(x);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

void JumpDirichlet:: setBC(int node, int id, Vector<Real>& x, int /*comp*/, 
			   Real /*time*/)
{
    Id[node] = id;
    Values[node] = value(x);
}
//-------------------------------------------------------------------------

Real JumpDirichlet:: value(Vector<Real>& /*x*/) const { return 0.0; }

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

Real JumpDirichlet:: initialValue(Vector<Real>& x, Real /*time*/, int /*comp*/) const
{
    if (x[1] <= xStep) return -1.0;
    else  return 1.0;
}
//-------------------------------------------------------------------------

Real TransDirichlet:: value(Vector<Real>& x) const
{
    Real val=1.0, p=100.0, x0=0.5, xm;
    int i;

    FORALL(x,i)
    {
    	xm = x[i] - x0;
    	val *= x[i]*(x[i]-1.0)*exp(-p*xm*xm);
    }

    return val;
}

void TransDirichlet:: setBC(int node, int id, Vector<Real>& x, int /*comp*/, 
			    Real /*time*/)
{
    Id[node] = id;
    Values[node] = value(x);
}
//-------------------------------------------------------------------------
Real TransDirichlet:: initialValue(Vector<Real>& x, Real /*time*/, int /*comp*/) const
{
    return value(x);
}
//-------------------------------------------------------------------------

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

StefanDirichletBCs:: StefanDirichletBCs() : VarDirichletBCs()
{ 
    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);

    theta1 = 0.0;     	Cmd.get("StefanTheta1", &theta1);
}
//-------------------------------------------------------------------------

void StefanDirichletBCs:: setBC(int node, int id, Vector<Real>& x, int /*comp*/, 
				Real time)
{
    Id[node] = id;
    Values[node] = value(x, time);
}
//-------------------------------------------------------------------------

Real StefanDirichletBCs:: initialValue(Vector<Real>& x, Real t, int /*comp*/) const
{
    return value(x,t);
}
//-------------------------------------------------------------------------

Real StefanDirichletBCs:: value(Vector<Real>& x, Real t) const
{
    Real theta = sqr(x[1]-0.5) + sqr(x[2]-0.5) - exp(-4.*t)/4.;

    if (theta > theta1) theta *= kappa1;
    else		theta *= kappa0;

    return theta;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


CastingDirichletBCs:: CastingDirichletBCs() : VarDirichletBCs()
{ 
    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);

    theta1 = 0.0;     	Cmd.get("StefanTheta1", &theta1);
}
//-------------------------------------------------------------------------

void CastingDirichletBCs:: setBC(int node, int id, Vector<Real>& x, int /*comp*/, 
				Real time)
{
    Id[node] = id;
    Values[node] = value(x, time);
}
//-------------------------------------------------------------------------

Real CastingDirichletBCs:: initialValue(Vector<Real>& /*x*/, Real /*t*/, int /*comp*/) const
{
    return 0.0000001;
    // return 0.0;
}
//-------------------------------------------------------------------------

Real CastingDirichletBCs:: value(Vector<Real>& /*x*/, Real /*t*/) const
{
    
    return -0.3;
}
//-------------------------------------------------------------------------

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

PorousMediaDirichletBCs:: PorousMediaDirichletBCs() : VarDirichletBCs() { }

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

void PorousMediaDirichletBCs:: setBC(int node, int id, Vector<Real>& /*x*/, 
				     int /*comp*/, Real /*time*/)
{
    Id[node] = id;
    Values[node] = 0.0;
}
//-------------------------------------------------------------------------


Real PorousMediaDirichletBCs:: initialValue(Vector<Real>& x, Real /*t*/, int /*comp*/)
									  const 
{
    Real r=0.0, phi=0.0, val=0.0;

    r = x[1]*x[1]+x[2]*x[2];
    
    if (x[1] == 0.0) phi = 4.*atan(1.)/2.;
    else             phi = atan(x[2]/x[1]);
    
    val = 0.4 - r*(1.+0.5*sin(14*phi));

    if (val < 0.0) val = 0.0;
    return val;
}