/*
$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()
{ }
syntax highlighted by Code2HTML, v. 0.9.1