/*
$Id: intA.cc,v 1.3 1996/10/11 15:15:29 bzferdma Exp $
(C)opyright 1996 by Konrad-Zuse-Center, Berlin
All rights reserved.
Part of the Kaskade distribution
*/
#include "intA.h"
#include "utils.h"
#include "familyA.h"
#include "nodeco.h"
#include "dirichlet.h"
#include "triang.h"
#include "triangA.h"
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
// -- store node Numbers to vector node for change local <-> global:
void LinElemInt:: getGlobalNodes(const PATCH* t, Vector<int>& node) const
{
int i, n, nGlob, nLoc = 1;
const Vector<PT*>& pt = t->getPoints();
FORALL(pt,i)
{
nGlob = pt[i]->getNode();
if (nGlob) FORNCOMP(n) node[nLoc++] = nGlob++;
else FORNCOMP(n) node[nLoc++] = 0;
}
//node[1] = e->p1->getNode();
//node[2] = e->p2->getNode(); ...
}
//-------------------------------------------------------------------------
void LinElemInt:: getGlobalMLNodes(const PATCH* t, Vector<int>& node, int depth)
const
{
int i, n, nGlob, nLoc = 1;
const Vector<PT*>& pt = t->getPoints();
FORALL(pt,i)
{
nGlob = pt[i]->MLNode(depth);
if (nGlob) FORNCOMP(n) node[nLoc++] = nGlob++;
else FORNCOMP(n) node[nLoc++] = 0;
}
//node[1] = e->p1->MLNode(depth);
//node[2] = e->p2->MLNode(depth); ...
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void LQuadElemInt:: getGlobalNodes(const PATCH* t, Vector<int>& node) const
{
int i, n, nGlob, nLoc = 1;
const Vector<PT*>& pt = t->getPoints();
FORALL(pt,i)
{
nGlob = pt[i]->getNode();
if (nGlob) FORNCOMP(n) node[nLoc++] = nGlob++;
else FORNCOMP(n) node[nLoc++] = 0;
}
const Vector<EDG*>& edg = t->getEdges();
FORALL(edg,i)
{
nGlob = edg[i]->getNode();
if (nGlob) FORNCOMP(n) node[nLoc++] = nGlob++;
else FORNCOMP(n) node[nLoc++] = 0;
}
}
//-------------------------------------------------------------------------
void LQuadElemInt:: getGlobalMLNodes(const PATCH* t, Vector<int>& node,
int depth) const
{
int i, n, nGlob, nLoc = 1;
const Vector<PT*>& pt = t->getPoints();
FORALL(pt,i)
{
nGlob = pt[i]->MLNode(depth);
FORNCOMP(n) node[nLoc++] = nGlob++;
}
const Vector<EDG*>& edg = t->getEdges();
FORALL(edg,i)
{
nGlob = edg[i]->getNode();
FORNCOMP(n) node[nLoc++] = nGlob++;
}
// node[1] = e->p1->MLNode(depth);
// node[2] = e->p2->MLNode(depth);
// node[3] = e->getNode(); ...
}
//-------------------------------------------------------------------------
void EFTetraInt:: getGlobalNodes(const PATCH* t, Vector<int>& node) const
{
int i, n, nGlob, nLoc = 1;
const Vector<PT*>& pt = t->getPoints();
FORALL(pt,i)
{
nGlob = pt[i]->getNode();
if (nGlob) FORNCOMP(n) node[nLoc++] = nGlob++;
else FORNCOMP(n) node[nLoc++] = 0;
}
const Vector<EDG*>& edg = t->getEdges();
FORALL(edg,i)
{
nGlob = edg[i]->getNode();
if (nGlob) FORNCOMP(n) node[nLoc++] = nGlob++;
else FORNCOMP(n) node[nLoc++] = 0;
}
Vector<int> faceNodes(4);
t->getFaceNodes(faceNodes);
FORALL(faceNodes,i)
{
nGlob = faceNodes[i];
if (nGlob) FORNCOMP(n) node[nLoc++] = nGlob++;
else FORNCOMP(n) node[nLoc++] = 0;
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void LinElemInt:: setNodeNumbers()
{
int node;
PATCH* p;
mesh->resetPtIter();
while ((p=mesh->ptIterAll())) p->unMark();
mesh->resetPtIter();
if (noOfNodes.high() < 0)
{
node = 1;
while ((p=mesh->ptIterAll()))
{
p->setNode(node);
node += nComp;
}
}
else
{
node = noOfNodes.Top() + 1;
while ((p=mesh->ptIterAll()))
if (p->getNode() == 0)
{
p->setNode(node);
p->setMark(); // new nodes are marked
node += nComp;
}
}
node -= 1;
noOfNodes.push(node);
}
//-------------------------------------------------------------------------
void LQuadElemInt:: setNodeNumbers()
{
int node;
PATCH* p;
mesh->resetPtIter();
while ((p=mesh->ptIterAll())) p->unMark();
mesh->resetPtIter();
mesh->resetEdgIter();
if (noOfNodes.high() < 0)
{
node = 1;
while ((p=mesh->ptIterAll()))
{
p->setNode(node);
p->setMark();
node += nComp;
}
node -= 1;
pointNodes.push(node);
}
else
{
node = pointNodes.Top() + 1;
while ((p=mesh->ptIterAll()))
if (p->getNode() == 0) // new node
{
p->setNode(node);
p->setMark();
node += nComp;
}
node -= 1;
pointNodes.push(node);
}
node += 1;
while ((p=mesh->edgIterAll()))
{
p->setNode(node);
node += nComp;
}
node -= 1;
noOfNodes.push(node);
}
//-------------------------------------------------------------------------
void HQuadElemInt:: setNodeNumbers()
{
int node;
PATCH* p;
mesh->resetPtIter();
while ((p=mesh->ptIterAll())) p->unMark();
mesh->resetPtIter();
mesh->resetEdgIter();
if (noOfNodes.high() < 0)
{
node = 1;
while ((p=mesh->ptIterAll()))
{
p->setNode(node);
p->setMark();
node += nComp;
}
node -= 1;
pointNodes.push(node);
}
else
{
node = pointNodes.Top() + 1;
while ((p=mesh->ptIterAll()))
if (p->getNode() == 0) // new node
{
p->setNode(node);
p->setMark();
node += nComp;
}
node -= 1;
pointNodes.push(node);
}
node += 1;
if (dirichletBCs->constant())
{
while ((p=mesh->edgIterAll()))
if (p->isDirichlet()) p->setNode(0);
else
{
p->setNode(node);
node += nComp;
}
}
else
{
while ((p=mesh->edgIterAll()))
{
p->setNode(node);
node += nComp;
}
}
node -= 1;
noOfNodes.push(node);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void LinElemInt:: interpolateSolution(Vector<Num>& u) const
{
PATCH* ed, *pm;
int n, node1, node2, pmNode;
mesh->resetEdgIter();
while ((ed=mesh->edgIterAllOfHistory()))
{
if (ed->refined())
{
pm = ed->getMidPoint();
if (pm->marked()) // new node here
{
const Vector<PT*>& p = ed->getPoints();
pmNode = pm->getNode();
node1 = p[1]->getNode();
node2 = p[2]->getNode();
FORNCOMP(n) u[pmNode++] = 0.5*(u[node1++] + u[node2++]);
}
}
}
}
//-------------------------------------------------------------------------
void LQuadElemInt:: interpolateSolution(Vector<Num>& u) const
{
int i, n;
PATCH* ed, *pm;
int edNode, node1, node2, node;
Vector<Num> uOld(u.high());
FORALL(u,i) uOld[i] = u[i];
mesh->resetEdgIter();
while ((ed=mesh->edgIterAllOfHistory())) // set the new point nodes
{
if (ed->refined())
{
pm = ed->getMidPoint();
if (pm->marked())
{
node = pm->getNode();
edNode = ed->getNode();
FORNCOMP(n) u[node++] = uOld[edNode++];
}
}
}
mesh->resetEdgIter();
while ((ed=mesh->edgIterAll())) // the edge nodes
{
const Vector<PT*>& p = ed->getPoints();
edNode = ed->getNode();
node1 = p[1]->getNode();
node2 = p[2]->getNode();
FORNCOMP(n) u[edNode++] = 0.5*(u[node1++] + u[node2++]);
}
}
//-------------------------------------------------------------------------
void HQuadElemInt:: interpolateSolution(Vector<Num>& u) const
{
PATCH* ed, *pm;
int i, n, node, node1, node2, edNode;
Vector<Num> uOld(u.high());
FORALL(u,i) uOld[i] = u[i];
mesh->resetEdgIter();
while ((ed=mesh->edgIterAllOfHistory())) // the new point nodes
{
if (ed->refined())
{
pm = ed->getMidPoint();
if (pm->marked()) // new point node here
{
edNode = ed->getNode();
node = pm->getNode();
if (edNode) FORNCOMP(n) u[node++] = uOld[edNode++];
else FORNCOMP(n) u[node++] = 0.0;
const Vector<PT*>& p = ed->getPoints();
node = pm->getNode();
node1 = p[1]->getNode();
node2 = p[2]->getNode();
FORNCOMP(n) u[node++] += 0.5 * (uOld[node1++] + uOld[node2++]);
}
}
}
mesh->resetEdgIter();
while ((ed=mesh->edgIterAll())) // the edge nodes
{
edNode=ed->getNode();
if (edNode) FORNCOMP(n) u[edNode++] = 0.0;
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void LinElemInt:: updateDirichletBCs(Real time) const
{
int n, node;
PATCH* p;
Vector<Real> x(spaceDim);
const Bool constDBCs = dirichletBCs->constant();
dirichletBCs->extend(noOfNodes.Top());
mesh->resetPtIter();
while ((p=mesh->ptIterAll()))
{
FORNCOMP(n)
{
if (dirichletBCs->isDirichlet(p->isDirichlet(), p->Class(), n, time))
{
if (!constDBCs) p->getNodeCoord(x);
node = p->getNode() + n-1;
dirichletBCs->setBC(node, p->Class(), x, n, time);
}
}
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void LQuadElemInt:: updateDirichletBCs(Real time) const
{
int n, node;
PATCH* p;
Vector<Real> x(spaceDim);
const Bool constDBCs = dirichletBCs->constant();
dirichletBCs->extend(noOfNodes.Top());
mesh->resetPtIter();
while ((p=mesh->ptIterAll()))
{
FORNCOMP(n)
{
if (dirichletBCs->isDirichlet(p->isDirichlet(), p->Class(), n, time))
{
if (!constDBCs) p->getNodeCoord(x);
node = p->getNode() + n-1;
dirichletBCs->setBC(node, p->Class(), x, n, time);
}
}
}
mesh->resetEdgIter();
while ((p=mesh->edgIterAll()))
{
FORNCOMP(n)
{
if (dirichletBCs->isDirichlet(p->isDirichlet(), p->Class(), n, time))
{
if (!constDBCs) p->getNodeCoord(x);
node = p->getNode() + n-1;
dirichletBCs->setBC(node, p->Class(), x, n, time);
}
}
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void HQuadElemInt:: updateDirichletBCs(Real time) const
{
int n, node;
PATCH* p;
Real val, valP1, valP2;
Vector<Real> x(spaceDim);
const Bool constDBCs = dirichletBCs->constant();
dirichletBCs->extend(noOfNodes.Top());
mesh->resetPtIter();
while ((p=mesh->ptIterAll()))
{
FORNCOMP(n)
{
if (dirichletBCs->isDirichlet(p->isDirichlet(), p->Class(), n, time))
{
if (!constDBCs) p->getNodeCoord(x);
node = p->getNode() + n-1;
dirichletBCs->setBC(node, p->Class(), x, n, time);
}
}
}
mesh->resetEdgIter();
while ((p=mesh->edgIterAll()))
{
FORNCOMP(n)
{
if (dirichletBCs->isDirichlet(p->isDirichlet(), p->Class(), n, time))
{
if (!constDBCs) p->getNodeCoord(x);
node = p->getNode() + n-1;
dirichletBCs->setBC(node, p->Class(), x, n, time);
const Vector<PT*>& pt = p->getPoints();
valP1 = dirichletBCs->value(pt[1]->getNode() + n-1);
valP2 = dirichletBCs->value(pt[2]->getNode() + n-1);
val = dirichletBCs->value(node);
val -= 0.5 * (valP1 + valP2); // hier. basis
dirichletBCs->setValue(val, node);
}
}
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void LinElemInt:: getNodeCoordinates(NodeCoordinates& nc) const
{
PATCH* p;
Vector<Real> x(spaceDim);
mesh->resetPtIter();
while ((p=mesh->ptIterAll()))
{
p->getNodeCoord(x);
nc.setCoordinate(p->getNode(), x);
}
}
//-------------------------------------------------------------------------
void LQuadElemInt:: getNodeCoordinates(NodeCoordinates& nc) const
{
PATCH* p;
Vector<Real> x(spaceDim);
mesh->resetPtIter();
mesh->resetEdgIter();
while ((p=mesh->ptIterAll()))
{
p->getNodeCoord(x);
nc.setCoordinate(p->getNode(), x);
}
while ((p=mesh->edgIterAll()))
{
if (p->getNode())
{
p->getNodeCoord(x);
nc.setCoordinate(p->getNode(), x);
}
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void LinElemInt:: setMLNodeNumbers(int maxDepth, int maxDepthM1, int targetDepth)
{
int depth;
PATCH* p;
if (maxDepth > maxDepthM1) noOfNodes.push(0);
mesh->resetPtIter();
if (maxDepth == 0)
{
while ((p=mesh->ptIterAll())) p->setNode(0);
mesh->resetPtIter();
}
while ((p=mesh->ptIterAll()))
{
PT* pt = p->castToPT();
if (pt->getNode() == 0) // a new node
{
depth = pt->Depth();
mesh->newNodeStack(pt, depth, depth+targetDepth);
pt->pushNode(noOfNodes[depth] + 1);
noOfNodes[depth] += nComp;
pt->setMark(); // new nodes are marked
}
else pt->unMark();
for (depth=pt->maxDepth()+1; depth<=maxDepth; ++depth)
{
pt->pushNode(noOfNodes[depth] + 1);
noOfNodes[depth] += nComp;
}
p->setNode(p->MLNode(maxDepth));
}
}
//-------------------------------------------------------------------------
void LinElemInt:: interpolateMLSolution(Vector<Num>& u, int maxDepth,
int maxDepthM1) const
{
int i, n, node, node0;
PATCH* p;
Vector<Num> uPrev(u.high());
FORALL(u,i) uPrev[i] = u[i];
u.resize(noOfNodes.Top());
FORALL(u,i) u[i] = 0.0;
if (maxDepth > maxDepthM1) // a new level was created
{
mesh->resetPtIter();
while ((p=mesh->ptIterAll()))
if (!p->marked()) // 'old' node
{
node = p->getNode();
node0 = p->MLNode(maxDepthM1);
FORNCOMP(n) u[node++] = uPrev[node0++];
}
}
else
{
mesh->resetPtIter();
while ((p=mesh->ptIterAll()))
if (!p->marked()) // 'old' node
{
node = p->getNode();
FORNCOMP(n) { u[node] = uPrev[node]; ++node; }
}
}
interpolateSolution(u); // handle new nodes
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void LQuadElemInt:: setHighOrderNodes(int* edNode)
{
PATCH* patch = 0;
mesh->resetEdgIter();
*edNode += 1;
while ((patch=mesh->edgIterAll()))
{
patch->setNode(*edNode);
*edNode += nComp;
}
*edNode -= 1;
noOfNodes.push(*edNode);
}
//-------------------------------------------------------------------------
void HQuadElemInt:: setHighOrderNodes(int* edNode)
{
PATCH* patch = 0;
mesh->resetEdgIter();
*edNode += 1;
while ((patch=mesh->edgIterAll()))
{
if (patch->isDirichlet()) patch->setNode(0);
else
{
patch->setNode(*edNode);
*edNode += nComp;
}
}
*edNode -= 1;
noOfNodes.push(*edNode);
}
//-------------------------------------------------------------------------
void EFTetraInt:: setHighOrderNodes(int* edNode)
{
int i;
Vector<PATCH*> neighb(4);
PATCH* patch = 0;
mesh->resetEdgIter();
*edNode += 1;
while ((patch=mesh->edgIterAll()))
{
if (patch->isDirichlet()) patch->setNode(0);
else
{
patch->setNode(*edNode);
*edNode += nComp;
}
}
mesh->resetElemIter(); // first reset faceNodes
while ((patch=mesh->elemIterAll()))
{
for (i=1; i<=4; ++i)
{
int& node = patch->getFaceNode(i);
node = 0;
}
}
mesh->resetElemIter();
while ((patch=mesh->elemIterAll()))
{
patch->getNeighbours(neighb);
const Vector<PATCH*>& faces = patch->getFaces();
FORALL(faces,i)
{
if (faces[i])
{
int& node = patch->getFaceNode(i);
if (!faces[i]->isDirichlet())
{
node = *edNode;
*edNode += nComp;
}
}
else if (neighb[i])
{
int& node = patch->getFaceNode(i);
if (node == 0)
{
int& neighbNode = neighb[i]->getMyFaceNode(patch);
if (neighbNode)
{
cout << "\n*** EFTetraInt:: setHighOrderNodes: "
<< "faceNode in Neighbour already set\n";
cout.flush(); abort();
}
node = *edNode;
neighbNode = *edNode;
*edNode += nComp;
}
}
}
}
*edNode -= 1;
noOfNodes.push(*edNode);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
// -- solution vector to quad. nodal basis: uN = P * uH
void LQuadElemInt:: solToNB(Vector<Num>& u)
{
int node, n, n1, n2;
PATCH* edg = 0;
mesh->resetEdgIter();
while ((edg=mesh->edgIterAll()))
{
if ((node = edg->getNode()))
{
const Vector<PT*>& p = edg->getPoints();
n1 = p[1]->getNode();
n2 = p[2]->getNode();
FORNCOMP(n) u[node++] += 0.5*(u[n1++] + u[n2++]);
}
}
}
//-------------------------------------------------------------------------
// -- solution vector to quad. hierarchical basis: uH = P**(-1) * uN
void LQuadElemInt:: solToHB(Vector<Num>& u)
{
int node, n, n1, n2;
PATCH* edg = 0;
mesh->resetEdgIter();
while ((edg=mesh->edgIterAll()))
{
if ((node = edg->getNode()))
{
const Vector<PT*>& p = edg->getPoints();
n1 = p[1]->getNode();
n2 = p[2]->getNode();
FORNCOMP(n) u[node++] -= 0.5*(u[n1++] + u[n2++]);
}
}
}
//-------------------------------------------------------------------------
// -- right-hand-side to quad. nodal basis: bN = P**(-T) * bH
void LQuadElemInt:: rhsToNB(Vector<Num>& b)
{
int n, node, nn, n1;
PATCH* edg = 0;
mesh->resetEdgIter();
while ((edg=mesh->edgIterAll()))
{
if ((node = edg->getNode()))
{
const Vector<PT*>& p = edg->getPoints();
if (!p[1]->isDirichlet())
{
n1 = p[1]->getNode();
nn = node;
FORNCOMP(n) b[n1++] -= 0.5 * b[nn++];
}
if (!p[2]->isDirichlet())
{
n1 = p[2]->getNode();
nn = node;
FORNCOMP(n) b[n1++] -= 0.5 * b[nn++];
}
}
}
}
//-------------------------------------------------------------------------
// -- right-hand-side to quad. hierarchical basis: bH = P**T * bN
void LQuadElemInt:: rhsToHB(Vector<Num>& b)
{
int n, node, nn, n1;
PATCH* edg = 0;
mesh->resetEdgIter();
while ((edg=mesh->edgIterAll()))
{
if ((node = edg->getNode()))
{
const Vector<PT*>& p = edg->getPoints();
if (!p[1]->isDirichlet())
{
n1 = p[1]->getNode();
nn = node;
FORNCOMP(n) b[n1++] += 0.5 * b[nn++];
}
if (!p[2]->isDirichlet())
{
n1 = p[2]->getNode();
nn = node;
FORNCOMP(n) b[n1++] += 0.5 * b[nn++];
}
}
}
}
//-------------------------------------------------------------------------
void LQuadElemInt:: setHBGeneration(Generation& gen)
{
int node;
PATCH* pt=0, *edg=0;
mesh->resetPtIter();
mesh->resetEdgIter();
while ((pt=mesh->ptIterAll()))
{
PointSon* son = new PointSon;
node = pt->getNode();
son->father = node;
gen.insertSon(node, son);
}
mesh->resetEdgIter();
while ((edg=mesh->edgIterAll()))
{
if ((node = edg->getNode()))
{
EdgeSon* eSon = new EdgeSon;
const Vector<PT*>& p = edg->getPoints();
eSon->father1 = p[1]->getNode();
eSon->father2 = p[2]->getNode();
gen.insertSon(node, eSon);
}
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void LinElemInt:: updateMGFamilyTree()
{
PATCH* ed, *pm;
int son;
LinEdgeGeneration* newGen = new LinEdgeGeneration(noOfNodes.Prev()+1,
noOfNodes.Top());
familyTree->setGeneration(newGen);
mesh->resetEdgIter();
while ((ed=mesh->edgIterAllOfHistory()))
{
if (ed->refined())
{
pm = ed->getMidPoint();
if (pm->marked()) // new node here
{
son = pm->getNode();
const Vector<PT*>& p = ed->getPoints();
newGen->setFathers(p[1]->getNode(), p[2]->getNode(), son);
}
}
}
}
//-------------------------------------------------------------------------
void LinElemInt:: updateMLFamilyTree(int maxDepth, int maxDepthM1)
{
int depth, nodeNo, father1, father2;
PATCH* ed, *pm, *pt;
for (depth=1; depth<=maxDepthM1; ++depth)
familyTree->extendGeneration(depth, noOfNodes[depth]);
if (maxDepth > maxDepthM1) // propagate 'old' nodes
{
Generation* gen;
if (nComp == 1) gen = new Generation(noOfNodes[maxDepth]);
else gen = new MCGeneration(noOfNodes[maxDepth], nComp);
familyTree->setGeneration(gen);
mesh->resetPtIter();
while ((pt=mesh->ptIterAll()))
{
if (!pt->marked()) // 'old' node
{
nodeNo = pt->MLNode(maxDepth);
father1 = pt->MLNode(maxDepth-1);
familyTree->newPointSon(nodeNo, father1, maxDepth);
}
}
}
// -- care for the new nodes (loop over edges to get the fathers):
mesh->resetEdgIter();
while ((ed=mesh->edgIterAllOfHistory()))
{
if (ed->refined())
{
pm = ed->getMidPoint();
if (pm->marked()) // a new node
{
depth = pm->Depth();
nodeNo = pm->MLNode(depth);
const Vector<PT*>& p = ed->getPoints();
father1 = p[1]->MLNode(depth-1);
father2 = p[2]->MLNode(depth-1);
familyTree->newEdgeSon(nodeNo, father1, father2, depth);
// -- set the remaining levels:
for (depth=pm->Depth()+1; depth<=maxDepth; ++depth)
{
nodeNo = pm->MLNode(depth);
father1 = pm->MLNode(depth-1);
familyTree->newPointSon(nodeNo, father1, depth);
}
}
}
}
}
syntax highlighted by Code2HTML, v. 0.9.1