/*
 $Id: familyA.cc,v 1.2 1996/11/20 09:57:25 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "familyA.h"
#include "utils.h"

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

StaticAllocator<PointSon> PointSon:: alloc(ElementsInBlock);
StaticAllocator<EdgeSon>  EdgeSon::  alloc(ElementsInBlock);

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


MCGeneration:: MCGeneration(int noOfNodes, int nComp0) 

		: Generation(noOfNodes/nComp0), nComp(nComp0)
{ }
//-------------------------------------------------------------------------

void MCGeneration:: extend(int noOfNodes) 
{
    Generation::extend(noOfNodes/nComp);
}
//-------------------------------------------------------------------------

void MCGeneration:: insertSon(int node, Son* newSon) 
{ 
    int block = MCNode.block(node, nComp);
    son[block] = newSon; 
}
//-------------------------------------------------------------------------

Son* MCGeneration:: getSon(int block) const   
{ 
    //block = MCNode.baseNode(node, nComp);
    return son[block]; 
}
//-------------------------------------------------------------------------


void MCGeneration:: prolong(const Vector<Num>& el, Vector<Num>& eh) const
{
    int i, node = 1;
    FORALL(son,i) 
    {
	son[i]->prolong(el, eh, node, nComp);
	node += nComp;
    }
}
//-------------------------------------------------------------------------


void MCGeneration:: restr(Vector<Num>& rh, Vector<Num>& rl) const
{
    int i, node = 1;
    FORALL(rl,i) rl[i] = 0.0;
    FORALL(son,i) 
    {
	son[i]->restr(rh, rl, node, nComp);
	node += nComp;
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void MCGeneration:: solToNB(Vector<Num>& e) const
{
    int i, node = 1;
    FORALL(son,i) 
    {
	son[i]->prolong(e, e, node, nComp);
	node += nComp;
    }
}
//-------------------------------------------------------------------------


void MCGeneration:: rhsToHB(Vector<Num>& r) const
{
    int i, node = 1;
    FORALL(son,i) 
    {
	son[i]->restr(r, r, node, nComp);	
	node += nComp;
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


LinEdgeGeneration:: LinEdgeGeneration(int low, int high) : Generation(1)
{ 
    fathers = new Matrix<int>(1,2,low,high); 
}

LinEdgeGeneration:: ~LinEdgeGeneration() { delete fathers; }

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

void LinEdgeGeneration:: prolong(const Vector<Num>& el, Vector<Num>& eh) const
{
    int i, n, fatherNode;
    const int low = fathers->cl;

    for (i=1; i<low; ++i)  eh[i] += el[i];	// point nodes 

    FORALL_COLUMNS(*fathers,i) 		
    {
	for (n=1; n<=2; ++n)
	{
	    fatherNode = (*fathers)(n,i);
	    if (fatherNode > low) eh[i] += 0.5 * eh[fatherNode];
	    else	          eh[i] += 0.5 * el[fatherNode];
	}
    }
    // the case (fatherNode > low) may occur with refinement by bisection!
}
//-------------------------------------------------------------------------

void LinEdgeGeneration:: restr(Vector<Num>& rh, Vector<Num>& rl) const
{
    int i, n, fatherNode;
    const int low = fathers->cl;

    for (i=1; i<low; ++i)  rl[i] = rh[i];

    for (i=fathers->ch; i>=low; --i)
    {
	for (n=1; n<=2; ++n)
	{
	    fatherNode = (*fathers)(n,i);
	    if (fatherNode > low) rh[fatherNode] += 0.5 * rh[i];
	    else	          rl[fatherNode] += 0.5 * rh[i];
	}
    }
    // the case (fatherNode > low) may occur with refinement by bisection!
}
//-------------------------------------------------------------------------

	  // Hierarchical basis transfer operations P and P**T:

void LinEdgeGeneration:: solToNB(Vector<Num>& e) const
{
    int i;
    FORALL_COLUMNS(*fathers,i) 
      e[i] += 0.5* (e[(*fathers)(1,i)] + e[(*fathers)(2,i)]);
}
//-------------------------------------------------------------------------

void LinEdgeGeneration:: rhsToHB(Vector<Num>& r) const
{
    int i;
    for (i=fathers->ch; i>=fathers->cl; --i)
    {
	    r[(*fathers)(1,i)] += 0.5*r[i];
	    r[(*fathers)(2,i)] += 0.5*r[i];
    }
}
//-------------------------------------------------------------------------

void LinEdgeGeneration:: setFathers(int father1, int father2, int node) 
{ 
    (*fathers)(1,node) = father1;
    (*fathers)(2,node) = father2;
}
//-------------------------------------------------------------------------

void LinEdgeGeneration:: print() const 
{
    int i;
    FORALL_COLUMNS(*fathers,i)
    {
	cout << "  node " << i << "  fathers: " 
	<< (*fathers)(1,i) << " " << (*fathers)(2,i) << "\n";
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void PointSon:: prolong(const Vector<Num>& el, Vector<Num>& eh, int node, 
			int nComp) const
{
    int n, f1 = father;
    FORNCOMP(n) eh[node++] += el[f1++];
}
//-------------------------------------------------------------------------

void PointSon:: restr(Vector<Num>& rh, Vector<Num>& rl, int node, 
			 int nComp) const
{
    int n, f1 = father;
    FORNCOMP(n) rl[f1++] += rh[node++];
}
//-------------------------------------------------------------------------

void PointSon:: print() const { cout << "   father: " << father << "\n"; }

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


void EdgeSon:: prolong(const Vector<Num>& el, Vector<Num>& eh, int node, 
		       int nComp) const
{
    int n, f1 = father1, f2 = father2;

    FORNCOMP(n)  eh[node++] += 0.5*(el[f1++] + el[f2++]); 
}
//-------------------------------------------------------------------------


void EdgeSon:: restr(Vector<Num>& rh, Vector<Num>& rl, int node, 
			int nComp) const
{
    int n, f1 = father1, f2 = father2;

    FORNCOMP(n)  
    { 
	rl[f1++] += 0.5 * rh[node];
	rl[f2++] += 0.5 * rh[node];
	++node;
    }
}
//-------------------------------------------------------------------------


void EdgeSon:: print() const
{
    cout << "   fathers: " << father1 << " " << father2 << "\n";
}
//-------------------------------------------------------------------------

int EdgeSon:: getFather(int no)  const
{
    if      (no == 1) return father1;
    else if (no == 2) return father2;
    else 
    {
	cout << "\n***  EdgeSon::Father(int no): wrong number no ("
	     << no << ") given\n";
	cout.flush(); abort();
    }
    return 0;
}
//-------------------------------------------------------------------------