/*
 $Id: userStatic.cc,v 1.3 1996/11/19 09:47:01 bzferdma Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "problem.h"
#include "materialsA.h"
#include "dirichletA.h"

#include "physics.h"
#include "utils.h"

#include "numerics.h"

#include "cmdpars.h"
extern CmdPars Cmd;

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

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

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

Real UserStaticMaterial:: S1d(Real x)
{
    Real s;

    s = 1.0 - 4.0*x + x*x*(x-1.0);
    return s;
}
    
Real UserStaticMaterial:: S2d(Real x, Real y)
{
    Real s;

    s =  x*y*(x+y) - x*x - y*y;
    return s;
}

Real UserStaticMaterial:: S3d(Real x, Real y, Real z)
{
    Real s;

    s = (x*y*z)*(x+y+z) - y*z;
    return s;
}
//-------------------------------------------------------------------------


Real UserStaticMaterial:: M(int /*type*/, Vector<Real>* x, int /*i*/, int /*j*/)
{
    Vector<Real>& xx = *x;

    switch(spaceDim)
    {
      case 1: return M1d(xx[1]);
      case 2: return M2d(xx[1], xx[2]);
      case 3: return M3d(xx[1], xx[2], xx[3]);
      default: cout.flush(); abort(); return 0;
    }
}

Real UserStaticMaterial:: M1d(Real x)
{
    Real q;

    q = x;
    return q;
}
    
Real UserStaticMaterial:: M2d(Real x, Real y)
{
    Real q;

    q = x+y;
    return q;
}

Real UserStaticMaterial:: M3d(Real x, Real y, Real z)
{
    Real q;

    q = x+y+z;
    return q;
}
//-------------------------------------------------------------------------



void UserStaticMaterial:: setMaterialType()
{
    //cout << "In UserStaticMaterial::setMaterialType: isotropic" << endl;
    matType = isotropic;
}

Real UserStaticMaterial:: E(int /*type*/, Vector<Real>* x, int /*i*/, int /*j*/)
{
    Vector<Real>& xx = *x;

    switch(spaceDim)
    {
      case 1: return E1d(xx[1]);
      case 2: return E2d(xx[1], xx[2]);
      case 3: return E3d(xx[1], xx[2], xx[3]);
      default: cout.flush(); abort(); return 0;
    }
}

Real UserStaticMaterial:: E1d(Real x)
{
    Real e;
    
    e = x;
    return e;
}
    
Real UserStaticMaterial:: E2d(Real x, Real y)
{
    Real e;
    
    e = x*y;
    return e;
}

Real UserStaticMaterial:: E3d(Real x, Real /*y*/, Real /*z*/)
{
    Real e;
    
    e = x;
    return e;
}
//-------------------------------------------------------------------------


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

    switch(spaceDim)
    {
      case 1: return 1.0; // xx[1]*(2.0*xx[1] - 1.0);
      case 2: if (type==1) return xx[1]*xx[2]*xx[2];
              else if (type==2) return 2.0*xx[1]*xx[1]*xx[2];
      case 3: return xx[1]*xx[2]*xx[3];
      default: cout.flush(); abort(); return 0;
    }
}

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

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

    switch(spaceDim)
    {
      case 1: return xx[1]*(2.0*xx[1] - 1.0);
      case 2: return xx[1];
      case 3: return xx[1]*xx[2]*xx[3];
      default: cout.flush(); abort(); return 0;
    }
}
//-------------------------------------------------------------------------
Real UserStaticMaterial:: Inner(int /*type*/, Vector<Real>* x, Real /*time*/)
{
    Vector<Real>& xx = *x;

    switch(spaceDim)
    {
      case 2: 
	if (!eTerm[1])  unknownID(1);
	if (!eTerm[2])  unknownID(2);
	return (eVal[1]-eVal[2])*(xx[1] + xx[2]); 
      case 3: 	
	if (!eTerm[1])  unknownID(1);
	if (!eTerm[2])  unknownID(2);
	return (eVal[1]-eVal[2])*(xx[1] + xx[2] + xx[3]); 
      default: cout.flush(); abort(); return 0;
    }
}

Bool  UserStaticMaterial:: trueSolKnown()
{
  return True;
}; 

Num UserStaticMaterial:: trueSolInPoint(const Vector<Real>& x, const Real /*time*/)
{
  switch(spaceDim)
    {
      case 1:          // user-static-1d.cmd
              {
		return x[1]*(x[1]-1.0);                             
	      }
      case 2:         // user-static-2d.cmd
	      { 
		return x[1]*x[2];
	      }
      case 3:         // user-static-3d.cmd
	      {
		return x[1]*x[2]*x[3];
	      }
      default: cout.flush(); abort(); return 0;
    }
}

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

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

    int spaceDim = x.h;
    switch(spaceDim)
    {
      case 1: Values[node] = userBC1d(x[1]); return;
      case 2: Values[node] = userBC2d(x[1], x[2]); return;
      case 3: Values[node] = userBC3d(x[1], x[2], x[3]); return;
      default: cout.flush(); abort(); return;
    }
}

Real UserDirichletBCs:: userBC1d(Real x)
{
    Real sum;
    
    sum = x*(x-1.0);
    return sum;
}

Real UserDirichletBCs:: userBC2d(Real x, Real y)
{
    Real sum;
    
    sum = x * y;
    return sum;
}


Real UserDirichletBCs:: userBC3d(Real x, Real y, Real z)
{
    Real sum;
    
    sum = x * y * z;
    return sum;
}

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

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



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

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

Real CylindricCoord:: S1d(Real x)
{
    Real s;

    s = 1.0 - 4.0*x + x*x*(x-1.0);
    return s;
}
    
Real CylindricCoord:: S2d(Real x, Real y)
{
    Real s;

    // s =  -(4.0 + x*x)*x*exp(-y);
    s =  -x*exp(-y);
    return s;
}

Real CylindricCoord:: S3d(Real x, Real y, Real z)
{
    Real s;

    s = -(4.0 + x*x + y*y)*exp(-z);
    return s;
}
//-------------------------------------------------------------------------


Real CylindricCoord:: M(int /*type*/, Vector<Real>* x, int /*i*/, int /*j*/)
{
    Vector<Real>& xx = *x;

    switch(spaceDim)
    {
      case 1: return M1d(xx[1]);
      case 2: return M2d(xx[1], xx[2]);
      case 3: return M3d(xx[1], xx[2], xx[3]);
      default: cout.flush(); abort(); return 0;
    }
}

Real CylindricCoord:: M1d(Real x)
{
    Real q;

    q = x;
    return q;
}
    
Real CylindricCoord:: M2d(Real /*x*/, Real /*y*/)
{  
    Real q;

    q = 0.0;
    return q;
}

Real CylindricCoord :: M3d(Real /*x*/, Real /*y*/, Real /*z*/)
{
    Real q;

    q = 0.0;
    return q;
}
//-------------------------------------------------------------------------



void CylindricCoord:: setMaterialType()
{
    //cout << "In CylindricCoord::setMaterialType: isotropic" << endl;
    matType = isotropic;
}

Real CylindricCoord:: E(int /*type*/, Vector<Real>* x, int /*i*/, int /*j*/)
{
    Vector<Real>& xx = *x;

    switch(spaceDim)
    {
      case 1: return E1d(xx[1]);
      case 2: return E2d(xx[1], xx[2]);
      case 3: return E3d(xx[1], xx[2], xx[3]);
      default: cout.flush(); abort(); return 0;
    }
}

Real CylindricCoord:: E1d(Real x)
{
    Real e;
    
    e = x;
    return e;
}
    
Real CylindricCoord:: E2d(Real x, Real /*y*/)
{
    Real e;
    
    e = x;
    return e;
}

Real CylindricCoord:: E3d(Real /*x*/, Real /*y*/, Real /*z*/)
{
    Real e;
    
    e = 1.0;
    return e;
}
//-------------------------------------------------------------------------


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

    switch(spaceDim)
    {
      case 1: return 0.0;
      case 2: return 0.0;
      case 3: return xx[1]*xx[2]*xx[3];
      default: cout.flush(); abort(); return 0;
    }
}

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

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

    switch(spaceDim)
    {
      case 1: return xx[1]*(2.0*xx[1] - 1.0);
      case 2: return xx[1];
      case 3: return xx[1]*xx[2]*xx[3];
      default: cout.flush(); abort(); return 0;
    }
}

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

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

    int spaceDim = x.h;
    switch(spaceDim)
    {
      case 1: Values[node] = userBC1d(x[1]); return;
      //case 2: Values[node] = userBC2d(x[1], x[2]); return;
      case 2: if (id==1)
              { if (x[2]<=-0.9999999) Values[node] = 300.0*(1.0-x[1])+50.0;
	        else Values[node] = 50.0;
		return;
	      }
              else if (id==2){ Values[node] = 200.0; return;}
      case 3: Values[node] = userBC3d(x[1], x[2], x[3]); return;
      default: cout.flush(); abort(); return;
    }
}

Real CylDirichletBCs:: userBC1d(Real x)
{
    Real sum;
    
    sum = x*(x-1.0);
    return sum;
}

Real CylDirichletBCs:: userBC2d(Real x, Real y)
{
    Real sum;
    
    sum = x*x*exp(-y);
    return sum;
}


Real CylDirichletBCs:: userBC3d(Real x, Real y, Real z)
{
    Real sum;
    
    sum = (x*x + y*y)*exp(-z);
    return sum;
}

Bool CylindricCoord:: trueSolKnown()
{
  switch(spaceDim)
    {
      case 1:  return False;
      case 2:  return False;
      case 3:  return True;
      default: cout.flush(); abort(); return 0;
    }
}; 

Num CylindricCoord:: trueSolInPoint(const Vector<Real>& x, const Real /*time*/)
{
  switch(spaceDim)
    {
      case 1:   return 0.0;                             
      case 2:   return 0.0;
      case 3:   return (x[1]*x[1] + x[2]*x[2])*exp(-x[3]);
      default: cout.flush(); abort(); return 0;
    }
}


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