/*
 $Id: problem2.cc,v 1.4 1996/11/20 10:01:26 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "problem2.h"

#include "physics.h"

#include "triang2tr.h"
#include "elements2.h"
#include "elements2mc.h"
#include "intB.h"

#include "adaptnl.h"

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

#include "linsystem.h"
#include "precond.h"

//#include "structure.h"
//#include "boxes.h"

#include "cmdpars.h"
extern CmdPars Cmd;
extern Problem* createProblem2(char* inFile, int /*spaceDim*/);


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

Problem* createProblem2(char* inFile, int /*spaceDim*/)
{
    Problem* problem = 0;

    if (Cmd.isSet("Problem", "StaticHeatConduction"))    
	problem = new StaticHeatConduction2(inFile);
    if (Cmd.isSet("Problem", "QuadStaticHeatConduction"))    
	problem = new QuadStaticHeatConduction2(inFile);

    else if (Cmd.isSet("Problem", "SkinEquation")) 
    	problem = new SkinEquation(inFile);

    else if (Cmd.isSet("Problem", "TransientHeatConduction"))    
	  problem = new TransientHeatConduction2(inFile);

    else if (Cmd.isSet("Problem", "Obstacle"))
 	  problem = new StaticNonLinearProblem2(inFile);

    else if (Cmd.isSet("Problem", "Stefan") ||   
	     Cmd.isSet("Problem", "Casting") || 
	     Cmd.isSet("Problem", "PorousMedia"))   
        problem = new TransientNonLinearProblem2(inFile);

    else if (Cmd.isSet("Problem", "MCStaticHeat"))    
	  problem = new MCStaticHeatConduction2(inFile);
    else if (Cmd.isSet("Problem", "MCTransientHeatConduction"))    
	  problem = new MCTransientHeatConduction2(inFile);

    return problem;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


StaticHeatConduction2:: StaticHeatConduction2(char* inFile) 

 	: Problem(inFile,sym,2), StaticHeatConduction(), mesh(0)
{
    dirichletBCs = newDirichletBCs();
    newMaterial();
 
    element = new Triangle(material);
}
//-------------------------------------------------------------------------

StaticHeatConduction2:: StaticHeatConduction2() : mesh(0) { }

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

StaticHeatConduction2:: ~StaticHeatConduction2()
{
    if (mesh) { mesh->Info(); delete mesh; }
}
//-------------------------------------------------------------------------

MESH* StaticHeatConduction2:: Mesh() const { return mesh; }

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

void StaticHeatConduction2:: newMesh()
{
    delete mesh;
    if (biSection) mesh = new MESH2BISECT(fileName);
    else	   mesh = new MESH2(fileName);
}
//-------------------------------------------------------------------------

void StaticHeatConduction2:: newInterface()
{
    delete interface;

    if (mesh 	== 0) missingObject("newInterface","mesh");
    if (element == 0) missingObject("newInterface","element");
    if (Ab 	== 0) missingObject("newInterface","Ab");
    if (precond == 0) missingObject("newInterface","precond");
    if (dirichletBCs == 0) missingObject("newInterface","dirichletBCs");


    switch(precond->mode())
    {
      case singleGrid:
	interface = new LinElemSG(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp);
	break;

      case multiGrid:
	interface = new LinElemMG(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp);
	break;

      case multiLevel:
	interface = new LinElemML(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp); 
	break;
    }
}
//-------------------------------------------------------------------------


void StaticHeatConduction2:: newErrorEstimator()
{
    delete errorEstimator;

    if (Cmd.isSet("errorEstimator","none")) 
			errorEstimator = new ErrorEstimator();

    else if (Cmd.isSet("errorEstimator","dly") ||
	     Cmd.isSet("errorEstimator","moddly"))
    {
	Element* elementDLY;

	if (Cmd.isSet("errorEstimator","dly")) 
	     elementDLY = new HQuadTriangle(material);
	else elementDLY = new HQuadTriangleFast(material);

	Interface* interfaceDLY = new HQuadElemSG(mesh, elementDLY, 
						  dirichletBCs, Ab, precond,
						  spaceDim,nComp);
	errorEstimator = new DLY(elementDLY, interfaceDLY);

    }
    else  
    {
	cout << "\n*** Problem:: no valid error estimator specified\n"; 
	cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


QuadStaticHeatConduction2:: QuadStaticHeatConduction2(char* inFile) 

 	: Problem(inFile,sym,2), StaticHeatConduction2()
{
    dirichletBCs = newDirichletBCs();
    newMaterial();
 
    element = new LQuadTriangle(material);
}
//-------------------------------------------------------------------------


void QuadStaticHeatConduction2:: newInterface()
{
    delete interface;

    if (mesh 	== 0) missingObject("newInterface","mesh");
    if (element == 0) missingObject("newInterface","element");
    if (Ab 	== 0) missingObject("newInterface","Ab");
    if (precond == 0) missingObject("newInterface","precond");
    if (dirichletBCs == 0) missingObject("newInterface","dirichletBCs");


    switch(precond->mode())
    {
      case singleGrid:
	interface = new LQuadElemSG(mesh,element,dirichletBCs, Ab,precond,
				    spaceDim,nComp);
	break;

      default: cout << "\n*** LQuadStaticHeatConduction2: only SG Precond.!\n";
	cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------


void QuadStaticHeatConduction2:: newErrorEstimator()
{
    delete errorEstimator;

    //if (Cmd.isSet("errorEstimator","none")) 
    errorEstimator = new ErrorEstimator();
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


SkinEquation:: SkinEquation(char* inFile) : Problem(inFile,sym,2), mesh(0)
{
    if(!Cmd.get("frequency", &omega)) 
	{ cout << "\n*** SkinEquation: no frequency given\n";  cout.flush(); abort(); }
    omega *= 2.*Pi;

    dirichletBCs = newDirichletBCs();
    newMaterial();

    element = new Triangle(material);
}
//-------------------------------------------------------------------------

SkinEquation:: ~SkinEquation()
{
    if (mesh) { mesh->Info();  delete mesh; }
}
//-------------------------------------------------------------------------

void SkinEquation:: newMaterial()
{
    delete material;
    material = new DefaultMaterial(fileName, spaceDim);
}
//-------------------------------------------------------------------------

DirichletBCs* SkinEquation:: newDirichletBCs()
{
    return new ConstDirichletBCs(fileName);
}
//-------------------------------------------------------------------------

MESH* SkinEquation:: Mesh() const { return mesh; }

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

void SkinEquation:: newMesh()
{
    delete mesh;
    mesh = new MESH2(fileName);

    //int nSteps = 0; Cmd.get("stepsResolve",&nSteps);
    //if (nSteps > 0) mesh->Resolve(0.0, nSteps);
}
//-------------------------------------------------------------------------

void SkinEquation:: newInterface()
{
    delete interface;

    if (mesh 	== 0) missingObject("newInterface","mesh");
    if (element == 0) missingObject("newInterface","element");
    if (Ab 	== 0) missingObject("newInterface","Ab");
    if (precond == 0) missingObject("newInterface","precond");
    if (dirichletBCs == 0) missingObject("newInterface","dirichletBCs");


    switch(precond->mode())
    {
      case singleGrid:
	interface = new LinElemSG(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp);
	break;
	
      case multiGrid:
	interface = new LinElemMG(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp);
	break;

      case multiLevel:
	interface = new LinElemML(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp); 
	break;
    }
}
//-------------------------------------------------------------------------


void SkinEquation:: newErrorEstimator()
{
    delete errorEstimator;

    if (Cmd.isSet("errorEstimator","none")) 
			errorEstimator = new ErrorEstimator();

    else if (Cmd.isSet("errorEstimator","dly") ||
	     Cmd.isSet("errorEstimator","moddly") )
    {
	Element* elementDLY;

	if (Cmd.isSet("errorEstimator","dly")) 
	     elementDLY = new HQuadTriangle(material);
	else elementDLY = new HQuadTriangleFast(material);

	Interface* interfaceDLY = new HQuadElemSG(mesh, elementDLY, 
						  dirichletBCs, Ab, precond,
						  spaceDim,nComp);
	errorEstimator = new DLY(elementDLY, interfaceDLY);
    }
    else  
    {
	cout << "\n*** Problem:: no valid error estimator specified\n"; 
	cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------


void SkinEquation:: 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;

    int dim = elem.NoOfNodes();
    Matrix<Real> AReal(dim, dim);
    Vector<Real> bReal(dim);

    elem.initAb(AReal,bReal);

    if (!elem.assembleEllip(t, AReal, Jac, pattern)) 
      missingMaterialTerm("elliptic", t.Class());

    for (i=1; i<=dim; ++i)
	for (k=1; k<=i; ++k)  A(i,k) += num(AReal(i,k),0.0);


    if (elem.assembleSource(t, bReal, Jac, pattern))
      FORALL(b,i) b[i] = num(bReal[i],0.0);

    elem.initAb(AReal);
    
    if (elem.assembleMass(t, AReal, Jac, pattern))   // conductor element
    {
	FORALL_ROWS(A,i)
	  for (k=1; k<=i; ++k)  A(i,k) += num(0.0, -omega*AReal(i,k));
    }
}
//-------------------------------------------------------------------------


void SkinEquation:: compEnergy(Vector<Num>& v, Real* eNorm, Num* fct)
{
    Ab->compEnergy(v, eNorm, fct);
}
//-------------------------------------------------------------------------


TransientHeatConduction2:: TransientHeatConduction2(char* inFile) 

       	: Problem(inFile,sym,2), TransientHeatConduction(),
	  mesh(0), prevMesh(0)
{
    dirichletBCs = newDirichletBCs();
    newMaterial();

    element = new Triangle(material);
}
//-------------------------------------------------------------------------

TransientHeatConduction2:: ~TransientHeatConduction2()
{
    if (mesh) { mesh->Info(); delete mesh; }
    delete prevMesh;
}
//-------------------------------------------------------------------------

MESH* TransientHeatConduction2:: Mesh()     const { return mesh; }
MESH* TransientHeatConduction2:: PrevMesh() const { return prevMesh; }

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

void TransientHeatConduction2:: newMesh()
{
    delete mesh;
    mesh = new MESH2Trans(fileName,prevMesh);
}
//-------------------------------------------------------------------------

void TransientHeatConduction2:: shiftMesh()
{
    delete prevMesh;
    prevMesh = mesh;
    mesh = 0;
}
//-------------------------------------------------------------------------


void TransientHeatConduction2:: newInterface()
{
    delete interface;

    if (mesh 	== 0) missingObject("newInterface","mesh");
    if (element == 0) missingObject("newInterface","element");
    if (Ab 	== 0) missingObject("newInterface","Ab");
    if (precond == 0) missingObject("newInterface","precond");
    if (dirichletBCs == 0) missingObject("newInterface","dirichletBCs");

    switch(precond->mode())
    {
      case singleGrid:
	interface = new LinElemSG(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp);
	break;
	
      case multiGrid:
	interface = new LinElemMG(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp); 
	break;
	
      case multiLevel:
	interface = new LinElemML(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp); 
	break;
    }
}
//-------------------------------------------------------------------------

void TransientHeatConduction2:: newErrorEstimator()
{
    delete errorEstimator;

    if (Cmd.isSet("errorEstimator","none")) 
			errorEstimator = new ErrorEstimator();
    else if (Cmd.isSet("errorEstimator","dly") ||
	     Cmd.isSet("errorEstimator","moddly"))
    {
	Element* elementDLY;
	
	if (Cmd.isSet("errorEstimator","dly")) 
	     elementDLY = new HQuadTriangle(material);
	else elementDLY = new HQuadTriangleFast(material);
	
	Interface* interfaceDLY = new HQuadElemSG(mesh, elementDLY, 
						  dirichletBCs, Ab, precond,
						  spaceDim,nComp);
	errorEstimator = new DLY(elementDLY, interfaceDLY);
    }	
    else  
    {
	cout << "\n*** Problem:: no valid error estimator specified\n"; 
	cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


StaticNonLinearProblem2:: StaticNonLinearProblem2(char* inFile) 

     	: Problem(inFile,sym,2),  StaticNonLinearProblem(), 
	  mesh(0)
{
    dirichletBCs = newDirichletBCs();
    newMaterial();

    element = new Triangle(material);
 }
//-------------------------------------------------------------------------

StaticNonLinearProblem2:: ~StaticNonLinearProblem2()
{
    if (mesh) { mesh->Info(); delete mesh; }
}
//-------------------------------------------------------------------------

MESH* StaticNonLinearProblem2:: Mesh() const { return mesh; }

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

void StaticNonLinearProblem2:: newMesh()
{
    delete mesh;
    if (biSection) mesh = new MESH2BISECT(fileName);
    else	   mesh = new MESH2(fileName);
}
//-------------------------------------------------------------------------

void StaticNonLinearProblem2:: newInterface()
{
    delete interface;

    if (mesh 	== 0) missingObject("newInterface","mesh");
    if (element == 0) missingObject("newInterface","element");
    if (Ab 	== 0) missingObject("newInterface","Ab");
    if (precond == 0) missingObject("newInterface","precond");
    if (dirichletBCs == 0) missingObject("newInterface","dirichletBCs");

    switch(precond->mode())
    {
      case singleGrid:
	interface = new LinElemSG(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp);
	break;

      case multiGrid:
	interface = new LinElemMG(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp);
	break;

      case multiLevel:
	interface = new LinElemML(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp); 
	break;
    }
}
//-------------------------------------------------------------------------

void StaticNonLinearProblem2:: newErrorEstimator()
{
    delete errorEstimator;

    if (Cmd.isSet("errorEstimator","none")) 
			errorEstimator = new ErrorEstimator();

    else if (Cmd.isSet("errorEstimator","rk1")) 
    {
	Element* elementDLY, *elementDLYLagrange;
	elementDLY 	   = new HQuadTriangle(material);
	elementDLYLagrange = new LQuadTriangle(material);
	Interface* interfaceDLY = new HQuadElemSG(mesh, elementDLY, 
						  dirichletBCs, Ab, precond,
						  spaceDim,nComp);
	errorEstimator = new RK1(elementDLY, elementDLYLagrange, interfaceDLY);
    }

    else if (Cmd.isSet("errorEstim","rk"))
    {
	Element* elementDLY, *elementDLYLagrange;
	elementDLY 	   = 0;
	elementDLYLagrange = new LQuadTriangle(material);
        dirichletBCsDLY    = newDirichletBCs();
	Interface* interfaceDLY = new LQuadElemSG(mesh, elementDLYLagrange, 
						  dirichletBCsDLY, Ab, precond,
						  spaceDim,nComp);
	if      (Cmd.isSet("errorEstimator","rkA"))
	 errorEstimator = new RKA(elementDLY, elementDLYLagrange, interfaceDLY);
	else if (Cmd.isSet("errorEstimator","rk2"))
	 errorEstimator = new RK2(elementDLY, elementDLYLagrange, interfaceDLY);
    }
	
    else  
    {
	cout << "\n*** Problem:: no valid error estimator specified\n"; 
	cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


TransientNonLinearProblem2:: TransientNonLinearProblem2(char* inFile) 

     	: Problem(inFile,sym,2), TransientNonLinearProblem(), 	  
	  mesh(0), prevMesh(0)
{
    dirichletBCs = newDirichletBCs();
    newMaterial();

    element = new Triangle(material);
}
//-------------------------------------------------------------------------

TransientNonLinearProblem2:: ~TransientNonLinearProblem2()
{
    if (mesh) { mesh->Info(); delete mesh; }
    delete prevMesh;
}
//-------------------------------------------------------------------------

MESH* TransientNonLinearProblem2:: Mesh()     const { return mesh; }
MESH* TransientNonLinearProblem2:: PrevMesh() const { return prevMesh; }

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

void TransientNonLinearProblem2:: newMesh()
{
    delete mesh;
    mesh = new MESH2Trans(fileName,prevMesh);
}
//-------------------------------------------------------------------------

void TransientNonLinearProblem2:: shiftMesh()
{
    delete prevMesh;
    prevMesh = mesh;
    mesh = 0;
}
//-------------------------------------------------------------------------

void TransientNonLinearProblem2:: newInterface()
{
    delete interface;

    if (mesh 	== 0) missingObject("newInterface","mesh");
    if (element == 0) missingObject("newInterface","element");
    if (Ab 	== 0) missingObject("newInterface","Ab");
    if (precond == 0) missingObject("newInterface","precond");
    if (dirichletBCs == 0) missingObject("newInterface","dirichletBCs");

    switch(precond->mode())
    {
      case singleGrid:
	interface = new LinElemSG(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp);
	break;

      case multiGrid:
	interface = new LinElemMG(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp);
	break;

      case multiLevel:
	interface = new LinElemML(mesh,element,dirichletBCs, Ab,precond,
				  spaceDim,nComp); 
	break;
    }
}
//-------------------------------------------------------------------------


void TransientNonLinearProblem2:: newErrorEstimator()
{
    delete errorEstimator;

    if (Cmd.isSet("errorEstimator","none")) 
			errorEstimator = new ErrorEstimator();

    else if (Cmd.isSet("errorEstimator","rk1")) 
    {
	Element* elementDLY, *elementDLYLagrange;
	elementDLY 	   = new HQuadTriangle(material);
	elementDLYLagrange = new LQuadTriangle(material);
	Interface* interfaceDLY = new HQuadElemSG(mesh, elementDLY, 
						  dirichletBCs, Ab, precond,
						  spaceDim,nComp);
	errorEstimator = new RK1(elementDLY, elementDLYLagrange, interfaceDLY);
    }
	
    else if (Cmd.isSet("errorEstim","rk"))
    {
	Element* elementDLY, *elementDLYLagrange;
	elementDLY 	   = 0;
	elementDLYLagrange = new LQuadTriangle(material);
        dirichletBCsDLY    = newDirichletBCs();
	Interface* interfaceDLY = new LQuadElemSG(mesh, elementDLYLagrange, 
						  dirichletBCsDLY, Ab, precond,
						  spaceDim, nComp);
	if      (Cmd.isSet("errorEstimator","rkA"))
	 errorEstimator = new RKA(elementDLY, elementDLYLagrange, interfaceDLY);
	else if (Cmd.isSet("errorEstimator","rk2"))
	 errorEstimator = new RK2(elementDLY, elementDLYLagrange, interfaceDLY);
    }

    else  
    {
	cout << "\n*** Problem:: no valid error estimator specified\n"; 
	cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


// --  	  test example for multi-component field:  nComp = 3


MCStaticHeatConduction2:: MCStaticHeatConduction2(char* inFile) 

 	: Problem(inFile,sym,2,3), StaticHeatConduction2() 
{
    dirichletBCs = newDirichletBCs();
    newMaterial();

    element = new MCTriangle(material, nComp);
}
//-------------------------------------------------------------------------


void MCStaticHeatConduction2:: newInterface()
{
    delete interface;

    if (mesh 	== 0) missingObject("newInterface","mesh");
    if (element == 0) missingObject("newInterface","element");
    if (Ab 	== 0) missingObject("newInterface","Ab");
    if (precond == 0) missingObject("newInterface","precond");
    if (dirichletBCs == 0) missingObject("newInterface","dirichletBCs");


    switch(precond->mode())
    {
      case singleGrid:
	interface = new LinElemSG(mesh, element, dirichletBCs, Ab, precond,
				  spaceDim, nComp);
	break;
	
      case multiLevel:
	interface = new LinElemML(mesh, element, dirichletBCs, Ab, precond, 
				  spaceDim, nComp); 
	break;

      default:
	cout << "\n*** MCStaticHeatConduction2:: interface not available\n";
	cout.flush(); abort();
	break;
    }
}
//-------------------------------------------------------------------------


void MCStaticHeatConduction2:: newErrorEstimator()
{
    delete errorEstimator;

    if (Cmd.isSet("errorEstimator","none")) 
			errorEstimator = new ErrorEstimator();
    else
    {
	if (Cmd.isSet("errorEstimator","dly"))
	{
	    Element*   elementDLY   = new HQuadMCTriangle(material,nComp);
	    Interface* interfaceDLY = new HQuadElemSG(mesh,elementDLY,
						      dirichletBCs,Ab,precond,
						      spaceDim,nComp);
	    errorEstimator = new DLY(elementDLY, interfaceDLY);
	}
	else  
	{
	    cout << "\n*** Problem:: no valid error estimator specified\n"; 
	    cout.flush(); abort();
	}
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


// --  test example for multi-component field:  nComp = 3


MCTransientHeatConduction2:: MCTransientHeatConduction2(char* inFile) 

 	: Problem(inFile,sym,2,3), TransientHeatConduction2()
{
    dirichletBCs = newDirichletBCs();
    newMaterial();

    element = new MCTriangle(material, nComp);
}
//-------------------------------------------------------------------------


void MCTransientHeatConduction2:: newInterface()
{
    delete interface;

    if (mesh 	== 0) missingObject("newInterface","mesh");
    if (element == 0) missingObject("newInterface","element");
    if (Ab 	== 0) missingObject("newInterface","Ab");
    if (precond == 0) missingObject("newInterface","precond");
    if (dirichletBCs == 0) missingObject("newInterface","dirichletBCs");

    switch(precond->mode())
    {
      case singleGrid:
	interface = new LinElemSG(mesh,element,dirichletBCs,Ab,precond,
				  spaceDim,nComp);
	break;
	
      case multiLevel:
	interface = new LinElemML(mesh,element,dirichletBCs,Ab,precond,
				  spaceDim,nComp); 
	break;

     default:
	cout << "\n*** Interface not available\n"; cout.flush(); abort(); break;
    }
}
//-------------------------------------------------------------------------


void MCTransientHeatConduction2:: newErrorEstimator()
{
    delete errorEstimator;

    if (Cmd.isSet("errorEstimator","none")) 
			errorEstimator = new ErrorEstimator();
    else
    {
	if (Cmd.isSet("errorEstimator","dly"))
	{
	    Element*   elementDLY   = new HQuadMCTriangle(material,nComp);
	    Interface* interfaceDLY = new HQuadElemSG(mesh,elementDLY,
						      dirichletBCs,Ab,precond,
						      spaceDim,nComp);
	    errorEstimator = new DLY(elementDLY, interfaceDLY);
	}
	else  
	{
	    cout << "\n*** Problem:: no valid error estimator specified\n"; 
	    cout.flush(); abort();
	}
    }
}