/*
 $Id: dirichletA.cc,v 1.3 1996/10/18 10:36:19 bzferdma Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "dirichletA.h"

#include "physics.h"
#include "mzibutil.h"

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

#include "cmdpars.h"
extern CmdPars Cmd;

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


ConstDirichletBCs:: ConstDirichletBCs(char* fileName) : DirichletBCs()
{ 
    readValues(fileName);
}
//-------------------------------------------------------------------------

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

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


void ConstDirichletBCs:: readValues(const char* fileName)
{
    int   i, classA, maxClass = 0, minClass = 0;
    float val;
    char  s[MaxLineLength], BCType;

    if (Cmd.get("matFile",s)) 
    { 

     if (!strchr(s,'.')) strcat(s,".mat");
    }
    else
    {
	strcpy(s,fileName);
	if (strchr(s,'.')) strcpy(strchr(s,'.'),".mat");
	else strcat(s,".mat");
    }

    BufferedParser parser(s, "end", CommentFlag);


    if (parser.findWord("Bound") != parser.ok) 
    {
	cout <<"\n*** No boundary conditions encountered in " << s << "\n";
	cout.flush(); abort();
    }
    long position = parser.ftell();


    // --		   scan the entries:


    while (parser.getValue(&classA) == parser.ok)
    {
	if (classA <= 0) 
	{
	    cout << "\n*** Class DirichletBCs: identifiers"
	         << " must be > 0 (read: " << classA << ")\n";
	    cout.flush(); abort();
	}

	parser.getValue(&BCType);
	BCType = toupper(BCType);
	if (BCType == 'D')  if (classA > maxClass) maxClass = classA;

	parser.getValue(&val);
	if (BCType == 'C') parser.getValue(&val);
    }
/*
    if (maxClass == 0)
    {
	cout << "\n*** No Dirichlet boundary conditions encountered in " 
	     << s << "\n";
	cout.flush(); abort();
    }
*/

    // --		   ... allocate and read:

    inputValues.resize(minClass, maxClass);
    FORALL(inputValues,i) inputValues[i] = -9999;

    parser.fseek(position);

    while (parser.getValue(&classA) == parser.ok)
    {
	parser.getValue(&BCType);
	BCType = toupper(BCType);
	parser.getValue(&val);
	if (BCType == 'D') inputValues[classA] = val;
	if (BCType == 'C') parser.getValue(&val);
    }


    		    // ... read factor and apply it:

    parser.rewind();

    float dFact = 1.0, fact;

    if (!(parser.findWord("Fac") == parser.ok))  return;
    
    while (parser.getString(s) == parser.ok)
    {
	parser.getValue(&fact);
	strToLower(s);
	if (s[0]=='d') dFact = fact;
    }
    FORALL(inputValues,i) inputValues[i] *= dFact;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


VarDirichletBCs:: VarDirichletBCs () : DirichletBCs()  { }

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


LinDirichletBCs:: LinDirichletBCs() : VarDirichletBCs() { }

void LinDirichletBCs:: setBC(int node, int id, Vector<Real>& x, int /*comp*/, 
			     Real /*time*/)
{ 
    int  i;
    Real sum = 0.0;
    FORALL(x,i)  sum += x[i];

    Id[node] = id;
    Values[node] = sum;
}
//-------------------------------------------------------------------------


QuadDirichletBCs:: QuadDirichletBCs() : VarDirichletBCs() { }

void QuadDirichletBCs:: setBC(int node, int id, Vector<Real>& x, int /*comp*/, 
			      Real /*time*/)
{ 
    int  i;
    Real sum = 0.0;
    FORALL(x,i)  sum += x[i]*x[i];

    Id[node] = id;
    Values[node] = sum;
}
//-------------------------------------------------------------------------


UserDirichletBCs:: UserDirichletBCs() : VarDirichletBCs() { }

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


CylDirichletBCs:: CylDirichletBCs() : UserDirichletBCs() { }

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

RootOfRBCs:: RootOfRBCs() : VarDirichletBCs() { }

void RootOfRBCs:: setBC(int node, int id, Vector<Real>& x, int /*comp*/, 
			      Real /*time*/)
{ 
    int  i;
    Real sum = 0.0, r;
    FORALL(x,i)  sum += x[i]*x[i];
    r = sqrt(sum);

    Values[node]=sqrt(r);
    Id[node] = id;

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

LayerBCs:: LayerBCs() : VarDirichletBCs() { }

void LayerBCs:: setBC(int node, int id, Vector<Real>& x, int /*comp*/, 
			      Real /*time*/)
{ 
    
    Values[node]= (cosh(10.0*x[1]) + cosh(10.0*x[2])) / (2.0*cosh(10.0));
    Id[node] = id;

}
//-------------------------------------------------------------------------
static Real arccos(Real x)
{
    Real y;

    y=(fabs(x)<1.0e-10)?REALPI2:atan(sqrt(1.0-x*x)/x);
    if (x<=-1.0e-10) y = REALPI+y;
    return y;
}
 
SlitBCs:: SlitBCs() : VarDirichletBCs() { }

void SlitBCs:: setBC(int node, int id, Vector<Real>& x, int /*comp*/, 
			      Real /*time*/)
{ 
    Real value, phi, r = sqrt(x[1]*x[1] + x[2]*x[2]);

    if (id==2) value = sqrt(sqrt(r));
    else
    {
	phi = (r>0.0)?arccos(x[1]/r):0.0;
	if (x[2]<0.0) phi = 2*REALPI-phi;
    
	value = sqrt(sqrt(r)) *sin(phi/4);
    }

    Values[node]= value;
    Id[node] = id;

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