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

#include "problem3.h"
#include "physics.h"

#include "triang3tr.h"
#include "elements3.h"
#include "elements3mc.h"
#include "intB.h"

#include "adapt.h"

#include "dirichletA.h"
#include "materialsA.h"
#include "linsystem.h"
#include "precond.h"

#include "cmdpars.h"
extern CmdPars Cmd;

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

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

    if (Cmd.isSet("Problem", "StaticHeatConduction"))    
    	  problem = new StaticHeatConduction3(inFile);
    else if (Cmd.isSet("Problem", "QuadStaticHeatConduction"))    
    	  problem = new QuadStaticHeatConduction3(inFile);

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

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

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


StaticHeatConduction3:: StaticHeatConduction3(char* inFile) 

 	: Problem(inFile,sym,3), StaticHeatConduction(), mesh(0)
{

    dirichletBCs = newDirichletBCs();
    newMaterial();

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

StaticHeatConduction3:: StaticHeatConduction3() : mesh(0) { }

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

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

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

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

void StaticHeatConduction3:: 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 StaticHeatConduction3:: newMesh()
{
    delete mesh;

    if (Cmd.isSet("errorEstimator","EFdly"))  mesh = new EFMESH3(fileName);
    else  mesh = new MESH3(fileName);
}
//-------------------------------------------------------------------------


void StaticHeatConduction3:: newErrorEstimator()
{
    Element*   elementDLY;
    Interface* interfaceDLY;

    delete errorEstimator;


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

    else if (Cmd.isSet("errorEstimator","dly") ||
	     Cmd.isSet("errorEstimator","moddly"))
    {
	if (Cmd.isSet("errorEstimator","dly")) 
	     elementDLY = new HQuadTetra(material);
	else elementDLY = new HQuadTetraFast(material);

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

    else if (Cmd.isSet("errorEstimator","EFdly")) 
    {
	elementDLY   = new EFTetraFast(material);
	interfaceDLY = new EFTetraSG(mesh, elementDLY, dirichletBCs, 
				     Ab, precond, spaceDim, nComp);
	errorEstimator = new EFDLY(elementDLY, interfaceDLY);
    }

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


QuadStaticHeatConduction3:: QuadStaticHeatConduction3(char* inFile) 

 	: Problem(inFile,sym,3), StaticHeatConduction3()
{
    dirichletBCs = newDirichletBCs();
    newMaterial();
 
    element = new LQuadTetra(material);
}
//-------------------------------------------------------------------------


void QuadStaticHeatConduction3:: 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*** QuadStaticHeatConduction3: only SG Precond.!\n";
	cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------


void QuadStaticHeatConduction3:: newErrorEstimator()
{
    delete errorEstimator;

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


TransientHeatConduction3:: TransientHeatConduction3(char* inFile) 

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

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

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

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

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

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

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


void TransientHeatConduction3:: 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 TransientHeatConduction3:: newErrorEstimator()
{
    Element*   elementDLY;
    Interface* interfaceDLY;

    delete errorEstimator;


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

    else if (Cmd.isSet("errorEstimator","dly") ||
	     Cmd.isSet("errorEstimator","moddly"))
    {
	if (Cmd.isSet("errorEstimator","dly")) 
	     elementDLY = new HQuadTetra(material);
	else elementDLY = new HQuadTetraFast(material);

	interfaceDLY = new HQuadElemSG(mesh, elementDLY, dirichletBCs,
				       Ab, precond, spaceDim,nComp);
	
	errorEstimator = new DLY(elementDLY, interfaceDLY);
	
    }
    // EFTetraFast may not be used here, as there is no adequate 
    //   lagrange element for the lumping procedure

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


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


MCStaticHeatConduction3:: MCStaticHeatConduction3(char* inFile) 

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

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


void MCStaticHeatConduction3:: 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*** MCStaticHeatConduction3:: interface not available\n";
	cout.flush(); abort();
	break;
    }
}
//-------------------------------------------------------------------------


void MCStaticHeatConduction3:: newErrorEstimator()
{
    delete errorEstimator;

    if (Cmd.isSet("errorEstimator","none")) 
			errorEstimator = new ErrorEstimator();
    else
    {
	if (Cmd.isSet("errorEstimator","dly"))
	{
	    Element*   elementDLY   = new HQuadMCTetra(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();
	}
    }
}