/*
 $Id: problemnl.cc,v 1.3 1996/10/08 07:53:02 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "problemnl.h"

#include "physics.h"

#include "nonlin.h"
#include "int.h"
#include "triang.h"
#include "elements.h"
#include "dirichlettr.h"
#include "materialstr.h"

#include "precondnl.h"

#include "cmdpars.h"
extern CmdPars Cmd;


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


NonLinearProblem:: NonLinearProblem() : nonLinearity(0)
{ 
    
    if      (Cmd.isSet("Problem","Obstacle"))  	nonLinearity = new Obstacle(); 
    else if (Cmd.isSet("Problem","Stefan")) 	nonLinearity = new Stefan();
    else if (Cmd.isSet("Problem","Casting")) 	nonLinearity = new Casting();
    else if (Cmd.isSet("Problem","PorousMedia"))nonLinearity = new PorousMedia();
    if (nonLinearity == 0) 
    {
	cout << "\n*** NonLinearProblem: no NonLinearity specified\n";
	cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------

NonLinearProblem:: ~NonLinearProblem() { delete nonLinearity; }

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


DirichletBCs* NonLinearProblem:: newDirichletBCs()
{
    if (Cmd.isSet("DirichletBCs","ConstDirichlet"))
         return new ConstDirichletBCs(fileName);

    else if (Cmd.isSet("DirichletBCs","StefanDirichlet"))
         return new StefanDirichletBCs();

    else if (Cmd.isSet("DirichletBCs","CastingDirichlet"))
         return new CastingDirichletBCs();

    else if (Cmd.isSet("DirichletBCs","PorousMediaDirichlet"))
         return new PorousMediaDirichletBCs();

    else MissingParameter("DirichletBCs");
    return 0;
}
//-------------------------------------------------------------------------

void NonLinearProblem:: newMaterial()
{
    delete material;

    if (Cmd.isSet("Material","DefaultMaterial"))
    	 material = new DefaultMaterial(fileName, spaceDim);

    else if (Cmd.isSet("Material","CastingMaterial"))
    	 material = new CastingMaterial(fileName, spaceDim);

    else if (Cmd.isSet("Material","StefanSource"))
    	 material = new StefanSource(fileName, spaceDim);

    else MissingParameter("Material");
}
//-------------------------------------------------------------------------


void NonLinearProblem:: newPreconditioner()
{
    delete precond;
    if      (Cmd.isSet("precond","NonLinSGGS")) 
    				precond = new NonLinearSGGS(*nonLinearity);
    else if (Cmd.isSet("precond","NonLinMLGS")) 
    				precond = new NonLinearMLGS(*nonLinearity);
    else if (Cmd.isSet("precond","TrcNonLinMLGS")) 
    				precond = new TrcNonLinearMLGS(*nonLinearity);
    else MissingParameter("Non-Linear Preconditioner");
}    
//-------------------------------------------------------------------------


void NonLinearProblem:: assembleGlobal(Real time0)
{
    nonLinearity->update(*interface);
    Problem::assembleGlobal(time0);
};
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void StaticNonLinearProblem:: assemble(const Element& elem, const PATCH& t, 
				       const Jacobian& Jac, Matrix<Num>& A, 
				       Vector<Num>& b, 
				       const Matrix<Bool>* pattern, 
				       Bool errorEstimatorCall)
{
    int i;

    if (!elem.assembleEllip(t,A,Jac,pattern))  missingMaterialTerm("elliptic",
								   t.Class());
    elem.assembleMass(t,A,Jac,pattern);
    elem.assembleSource(t,b,Jac,pattern);

    if (t.onBoundary())
    {
	elem.assembleCauchyBCs(t,A,b,pattern);
	elem.assembleNeumannBCs(t,b,pattern);
    }

    if (Mesh()->innerBoundary) elem.assembleInnerBCs(t,b,pattern);
    //if (Cmd.isTrue("innerBoundary")) elem.assembleInnerBCs(t,b,pattern);

    // --    assemble area weights for non-linear contributions:

    if (!errorEstimatorCall)		    // if not called by error estimator
    {
	const int dim = elem.NoOfNodes();
	Vector<Real> LNorm(dim);
	Vector<int> nodes(dim);
	
	FORALL(LNorm,i) LNorm[i] = 0.0;

	elem.assembleLNorm(t,LNorm,Jac,pattern);
	interface->getGlobalNodes(&t,nodes);
	
	nonLinearity->addAreaWeight(LNorm,nodes);
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void TransientNonLinearProblem:: assembleGlobal(Real time0)
{
    uPrevOnNewMesh.resize(1,interface->Dim());
    solutionToNewMesh(uPrevOnNewMesh);
    nonLinearity->update(*interface);

    Problem::assembleGlobal(time0);
};
//-------------------------------------------------------------------------


void TransientNonLinearProblem:: assemble(const Element& elem, const PATCH& t, 
					  const Jacobian& Jac, Matrix<Num>& A, 
					  Vector<Num>& b, 
					  const Matrix<Bool>* pattern,
					  Bool errorEstimatorCall)
{
    int i, k, node;

    const int dim = elem.NoOfNodes();

    Vector<Num> s(dim);


    // --		the stiffness matrix:

    if (!elem.assembleEllip(t,A,Jac,pattern)) missingMaterialTerm("elliptic",
								  t.Class());
    for (i=1; i<=dim; ++i)
    for (k=1; k<=i;   ++k)  A(i,k) = tau*A(i,k);


    // --      source term and boundary conditions:


    for (i=1; i<=dim; ++i) s[i] = 0.0;

    Bool update;
    update = elem.assembleSource(t,s,Jac,pattern,time);

    if (t.onBoundary()) 
	update = update || elem.assembleNeumannBCs(t,s,pattern,time);

    if (update) FORALL(b,i) b[i] += tau*s[i];


    if (t.onBoundary()) 
    {
	Matrix<Num> M(dim,dim);	 

	for (i=1; i<=dim; ++i) 
	{
	    s[i] = 0.0;
	    for (k=1; k<=i; ++k) M(i,k) = 0.0;
	}
	if (elem.assembleCauchyBCs(t,M,s,pattern,time))
	{
	    for (i=1; i<=dim; ++i) 
	    {
		b[i] += tau*s[i];
		for (k=1; k<=i; ++k)  A(i,k) += tau*M(i,k);
	    }
	}
    }

    if (!errorEstimatorCall)
    {
	Vector<int>  nodes(dim);
	Vector<Real> LNorm(dim);
	

	// --       area weights for non-linear contributions:
	
	FORALL(LNorm,i) LNorm[i] = 0.0;
	
	elem.assembleLNorm(t,LNorm,Jac,pattern);
	interface->getGlobalNodes(&t,nodes);
	
	nonLinearity->addAreaWeight(LNorm,nodes);


	// -- the contribution of the previous step to the right-hand-side:
    
	for (i=1; i<=dim; ++i)  
	{
	    node = nodes[i];
	    b[i] += LNorm[i] * nonLinearity->HValue(node,uPrevOnNewMesh[node]);
	}
    }
}
//-------------------------------------------------------------------------


Real TransientNonLinearProblem:: compUNorm(Vector<Num>& v, Vector<Num>* /*MDiag*/)
{
    int i;
    
    Real minU = ::real(v[1]);

    FORALL(v,i) { if (::real(v[1]) < minU)  minU = ::real(v[1]); }

    FORALL(v,i) v[i] -= minU;
		  
    Real UNorm = 0.0;

    UNorm = 2.0*EllipticNorm(v);
		  
    FORALL(v,i) v[i] += minU;
		  
    return sqrt(UNorm);
}