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

#include "intB.h"

#include "utils.h"
#include "familyA.h"
#include "connect.h"
#include "nodeco.h"
#include "dirichlet.h"

#include "triang.h"
#include "elements.h"

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

#include "sysmatsp.h"
#include "sysmatml.h"
#include "sysmatbl.h"

#include "cmdpars.h"
extern CmdPars Cmd;


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


void SGInt:: updateLevel0(Vector<Num>& u)
{
    int i;
    Timer timer, accTimer;

    setNodeNumbers();

    u.resize(1,noOfNodes.Top());
    FORALL(u,i) u[i] = 0.0;

    if (nComp == 1)
    {
	ConnectionPattern cPattern(*mesh, *element, *this, noOfNodes.Top());
	A = new SparseMatrix(cPattern, Ab->symmetry, spaceDim);
    }
    else A = new MLBlockMatrix(Ab->symmetry, spaceDim, noOfNodes.Top(), nComp); 

    Ab->update(A);

    timeInfo(timer, accTimer);
}
//-------------------------------------------------------------------------


void SGInt:: refine(Vector<Num>& u)
{
    mesh->Refine();

    Timer timer, accTimer;

    setNodeNumbers();

    u.extendAndCopy(noOfNodes.Top());
    interpolateSolution(u);

    if (nComp == 1)
    {
	ConnectionPattern cPattern(*mesh, *element, *this, noOfNodes.Top());
	A = new SparseMatrix(cPattern, Ab->symmetry, spaceDim);
    }
    else A = new MLBlockMatrix(Ab->symmetry, spaceDim, noOfNodes.Top(), nComp); 

    Ab->update(A);

    timeInfo(timer, accTimer);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


MGInt:: MGInt() { familyTree = new FamilyTree(1,DefaultStackSize); }

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


void MGInt:: refine(Vector<Num>& u)
{
    int i;

    mesh->Refine();

    Timer timer, accTimer;

    setNodeNumbers();

    updateMGFamilyTree();

    Vector<Num> uPrev(u.high());
    FORALL(u,i) uPrev[i] = u[i];
    u.resize(noOfNodes.Top());
    FORALL(u,i) u[i] = 0.0;

    familyTree->prolong(uPrev,u,noOfNodes.high());


    ConnectionPattern cPattern(*mesh, *element, *this, noOfNodes.Top());
    
    A = new SparseMatrix(cPattern, Ab->symmetry, spaceDim);

    Ab->update(A);

    timeInfo(timer, accTimer);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


MLInt:: MLInt() : Al(0,targetDepth),  maxDepth(-1), maxDepthM1(-2)
{
    familyTree = new FamilyTree(1,DefaultStackSize);

    localSmooth = True;   Cmd.get("localSmooth",&localSmooth);
    localTop    = True;   Cmd.get("localTop",   &localTop);
    localExtend = 1;      Cmd.get("localExtend",&localExtend);

    int baseLevel = 0;	  Cmd.get("BaseLevel", &baseLevel);

    if (baseLevel > 0 && localSmooth)
    {
	cout << "\n** Multi-level-options inconsistent: baseLevel > 0\n"
	     << "\t and localSmooth == True (set to False)\n";
	localSmooth = False;
    }
}
//-------------------------------------------------------------------------


void MLInt:: updateLevel0(Vector<Num>& u)	
{
    int i;
    Timer timer, accTimer;
    MLMatrix* AML;

    maxDepth = 0;

    setMLNodeNumbers(maxDepth, maxDepthM1, targetDepth);

    u.resize(1,noOfNodes[0]);
    FORALL(u,i) u[i] = 0.0;

    if (nComp == 1)  
         AML = new MLSparseMatrix(Ab->symmetry, spaceDim, noOfNodes.Top()); 
    else AML = new MLBlockMatrix(Ab->symmetry, spaceDim, noOfNodes.Top(), nComp);


    Al.push(AML);
    A = AML;

    Ab->update(A);

    timeInfo(timer, accTimer);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void MLInt:: refine(Vector<Num>& u)
{
    int depth;
    MLMatrix* AL;

    if (spaceDim > 1) removeGreenNodes();		// from diagonal !

    mesh->Refine();

    Timer timer, accTimer;

    maxDepthM1 = maxDepth;
    maxDepth   = mesh->MaxDepth();

    setMLNodeNumbers(maxDepth, maxDepthM1, targetDepth);
    updateMLFamilyTree(maxDepth, maxDepthM1);  	// grid-markers are used here
    interpolateMLSolution(u, maxDepth, maxDepthM1);

    if (maxDepth > maxDepthM1)  		// create new matrix
    {
	if (nComp == 1)	    
	{
	    if (localSmooth) 
	           AL = new LocalMLSparseMatrix(Ab->symmetry, spaceDim);
	    else   AL = new MLSparseMatrix(Ab->symmetry, spaceDim);
	}
	else
	{
	    if (localSmooth) AL = new LocalMLBlockMatrix(Ab->symmetry, spaceDim, 
							 noOfNodes.Top(),nComp);
	    else             AL = new MLBlockMatrix(Ab->symmetry, spaceDim, 
						    noOfNodes.Top(),nComp);
	}
	Al.push(AL);
    }

    for (depth=1; depth<=maxDepth; ++depth)     // depth 0 matrix never extended
      Al[depth]->extend(noOfNodes[depth]);


    setDirichletFlags();   
    if (localSmooth) setSmoothingFlags();
 
    A = Al[maxDepth];				// for preconditioner-update
    Ab->update(A);

    
    timeInfo(timer, accTimer);
    if (Cmd.isTrue("printInterFace")) print();
}
//-------------------------------------------------------------------------


void MLInt:: setDirichletFlags()
{
    int minD, depth, n, node;
    PATCH* p;

    //	--  the top-level matrix is handled in setDirichletBCs

    mesh->resetPtIter();
    while ((p=mesh->ptIterAll()))
    {
	if (p->isDirichlet())
	{
	    minD = Max(1,p->Depth());

	    for (depth=minD; depth<=maxDepth-1; ++depth)
	    {
		node = p->MLNode(depth);
		FORNCOMP(n) Al[depth]->setDirichlet(node++);
	    }
	}
    }
}
//-------------------------------------------------------------------------


void MLInt:: setSmoothingFlags()
{
    int  i, maxD, depth;
    PATCH* p;

    if (localTop) maxD = maxDepth;
    else	
    {
	maxD = maxDepth-1;
	for (i=1; i<=Al[maxDepth]->Dim(); ++i) Al[maxDepth]->Smooth(i);
    }
    
    for (depth=1; depth<=maxD; ++depth)	
 	for (i=1; i<=Al[depth]->Dim(); ++i)  Al[depth]->dontSmooth(i);
    
    mesh->resetPtIter();
    while ((p=mesh->ptIterAll()))
    {
	depth = p->Depth();         // point is new at the depth it was created
	
	if (depth>0 && depth<=maxD) Al[depth]->Smooth(p->MLNode(depth));
    }

    if (localExtend) extendLocalSmoothingPattern(); 
}
//-------------------------------------------------------------------------


void MLInt:: extendLocalSmoothingPattern()
{
    int i, n, depth, maxD;
    PATCH* patch;
    Bool smoothElement;

    Vector<int> globalNodes(element->NoOfNodes());
  
    if (localTop) maxD = maxDepth;
    else	  maxD = maxDepth-1;		// smooth all nodes on top Level

    for (depth=1; depth<=maxD; ++depth) Al[depth]->unMarkNodes();
	

    for (n=localExtend; n > 0; --n)
    {  
	mesh->resetElemIter();

	while ((patch=mesh->elemIterAllOfHistory()))
	{
	    depth = patch->Depth();

	    if (depth > 0 && depth <= maxD)
	    {
		getGlobalMLNodes(patch, globalNodes, depth);
  
		smoothElement = False;
		FORALL(globalNodes,i) 		// a smoothed node in element ?
		{
		    if (globalNodes[i])
		      if (Al[depth]->isSmoothed(globalNodes[i])) 
		      { 
			  smoothElement = True; 
			  break; 
		      }
		}
		if (smoothElement)
		{
		    FORALL(globalNodes,i) 
		    { 
			if (globalNodes[i])
			  if (!Al[depth]->isSmoothed(globalNodes[i]))
			    		     Al[depth]->setMark(globalNodes[i]);
		    }
		}
	    }
	}

	for (depth=1; depth<=maxD; ++depth)
	{
	    for (i=1; i<=Al[depth]->Dim(); ++i) 
		if (Al[depth]->isMarked(i)) Al[depth]->Smooth(i);
	}
    }
}  
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


LinElemSG:: LinElemSG(MESH* mesh0, const Element* element0, 
		      DirichletBCs* dirichletBCs0, LinSystem* Ab0, 
		      Preconditioner* precond0, int spaceDim0, int nComp0)

	: Interface(mesh0, element0, dirichletBCs0, Ab0, precond0, spaceDim0, 
		    nComp0), LinElemInt(), SGInt()
{ }
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


LinElemMG:: LinElemMG(MESH* mesh0, const Element* element0,
		      DirichletBCs* dirichletBCs0, LinSystem* Ab0, 
		      Preconditioner* precond0, int spaceDim0, int nComp0)

	: Interface(mesh0, element0, dirichletBCs0, Ab0, precond0, spaceDim0),
	  LinElemInt(), MGInt()
{ 
    if (nComp0 > 1) 
    {
	cout << "\n*** LinElemMG: nComp must be == 1\n";  cout.flush(); abort();
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


LinElemML:: LinElemML(MESH* mesh0, const Element* element0,
		      DirichletBCs* dirichletBCs0, LinSystem* Ab0, 
		      Preconditioner* precond0, int spaceDim0, int nComp0)

	: Interface(mesh0, element0, dirichletBCs0, Ab0, precond0, spaceDim0,
		    nComp0), LinElemInt(), MLInt()
{ }
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


LQuadElemSG:: LQuadElemSG(MESH* mesh0, const Element* element0,
			  DirichletBCs* dirichletBCs0, LinSystem* Ab0, 
			  Preconditioner* precond0, int spaceDim0, int nComp0)

	: Interface(mesh0, element0, dirichletBCs0, Ab0, precond0, spaceDim0,
		    nComp0), LQuadElemInt(), SGInt()
{ }
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


HQuadElemSG:: HQuadElemSG(MESH* mesh0, const Element* element0,
			  DirichletBCs* dirichletBCs0, LinSystem* Ab0, 
			  Preconditioner* precond0, int spaceDim0, int nComp0)

	: Interface(mesh0, element0, dirichletBCs0, Ab0, precond0, spaceDim0,
		    nComp0), LQuadElemInt(), SGInt()
{ }
//-------------------------------------------------------------------------


EFTetraSG:: EFTetraSG(MESH* mesh0, const Element* element0,
		      DirichletBCs* dirichletBCs0, LinSystem* Ab0, 
		      Preconditioner* precond0, int spaceDim0, int nComp0)

	: Interface(mesh0, element0, dirichletBCs0, Ab0, precond0, spaceDim0,
		    nComp0), EFTetraInt(), SGInt()
{ }