/*
$Id: triang.cc,v 1.4 1996/10/11 15:15:34 bzferdma Exp $
(C)opyright 1996 by Konrad-Zuse-Center, Berlin
All rights reserved.
Part of the Kaskade distribution
*/
#include "triang.h"
#include "triangA.h"
#include "utils.h"
#include "numerics.h"
#include "cmdpars.h"
extern CmdPars Cmd;
//-------------------------------------------------------------------------
Jacobian:: Jacobian() : J(3,3), Jinv(3,3) { }
//-------------------------------------------------------------------------
void PATCH:: notImplemented(const char* s) const
{
cout << "\n*** PATCH:: " << s << " not implemented\n"; cout.flush(); abort();
}
//-------------------------------------------------------------------------
int PATCH:: MLNode(int /*depth*/) const { notImplemented("MLNode"); return 0; }
PATCH* PATCH:: getMidPoint() const { notImplemented("getMidPoint"); return 0; }
static const Vector<PT*> pDummyForCompiler(1);
static const Vector<EDG*> eDummyForCompiler(1);
const Vector<PT*>& PATCH:: getPoints() const { notImplemented("getPoints()");
return pDummyForCompiler; }
const Vector<EDG*>& PATCH:: getEdges() const { notImplemented("getEdges()");
return eDummyForCompiler; }
static const Vector<PATCH*> fDummyForCompiler(1);
const Vector<PATCH*>& PATCH:: getFaces() const { notImplemented("getFaces()");
return fDummyForCompiler; }
void PATCH:: getAllFaces (Vector<PATCH*>& /*f*/) const { notImplemented("getAllFaces"); }
void PATCH:: returnFaces (Vector<PATCH*>& /*f*/) const { return; }
void PATCH:: getNeighbours(Vector<PATCH*>& /*neighbours*/) const
{ notImplemented("getNeighbours");}
void PATCH:: getEdgeOrientation(Vector<Bool>& /*orient*/) const
{ notImplemented("getEdgeOrientation");}
int PATCH:: noOfPoints() const { notImplemented("noOfPoints"); return 0; }
int PATCH:: noOfEdges() const { notImplemented("noOfEdges"); return 0; }
int PATCH:: noOfFaces() const { notImplemented("noOfFaces"); return 0; }
int PATCH:: noOfNeighbours() const
{ notImplemented("noOfNeighboursFaces"); return 0;}
float PATCH:: getError() const { notImplemented("getError"); return 0; }
void PATCH:: setError(float /*x*/) { notImplemented("setError"); }
//-------------------------------------------------------------------------
void PATCH:: getFaceNodes(Vector<int>& /*nodes*/) const
{
notImplemented("getFaceNodes");
}
//-------------------------------------------------------------------------
int& PATCH:: getFaceNode(int /*face*/) // for EFTET3
{
static int dummy;
notImplemented("getFaceNode");
return dummy;
}
int& PATCH:: getMyFaceNode(PATCH* /*patch*/) // for EFTET3
{
static int dummy;
notImplemented("getMyFaceNode");
return dummy;
}
//-------------------------------------------------------------------------
void PATCH:: realCoordinates(const Vector<Real>& /*u*/, Vector<Real>& /*x*/) const
{
notImplemented("realCoordinates(p1)");
}
void PATCH:: unitCoordinates(const Vector<Real>& /*u*/, Vector<Real>& /*x*/) const
{
notImplemented("unitCoordinates(p1)");
}
//-------------------------------------------------------------------------
void PATCH:: centerOfGravity(Vector<Real>& /*xG*/) const
{
notImplemented("centerOfGravity");
}
//-------------------------------------------------------------------------
void PATCH:: getNodeCoord(Vector<Real>& /*x*/) const
{
notImplemented("getCoordinates(p1)");
}
//-------------------------------------------------------------------------
void PATCH:: getCoordinates(Vector<Real>& /*p1*/) const
{
notImplemented("getCoordinates(p1)");
}
void PATCH:: getCoordinates(Vector<Real>& /*p1*/, Vector<Real>& /*p2*/) const
{
notImplemented("getCoordinates(p1,p2)");
}
void PATCH:: getCoordinates(Vector<Real>& /*p1*/, Vector<Real>& /*p2*/,
Vector<Real>& /*p3*/) const
{
notImplemented("getCoordinates(p1,p2,p3)");
}
void PATCH:: getCoordinates(Vector<Real>& /*p1*/, Vector<Real>& /*p2*/, Vector<Real>& /*p3*/,
Vector<Real>& /*p4*/) const
{
notImplemented("getCoordinates(p1,p2,p3,p4)");
}
//-------------------------------------------------------------------------
void PATCH:: compJ(Jacobian& /*j*/, Bool /*checkDetJ*/) const { notImplemented("compJ"); }
void PATCH:: compJinv(Jacobian& /*j*/) const { notImplemented("compJinv"); }
//-------------------------------------------------------------------------
Real PATCH:: volume() const { notImplemented("volume"); return 0; }
//-------------------------------------------------------------------------
Bool PATCH:: onBoundary() const { notImplemented("onBoundary"); return 0; }
Bool PATCH:: isDirichlet()const { notImplemented("isDirichlet"); return 0; }
Bool PATCH:: isNeumann() const { notImplemented("isNeumann"); return 0; }
Bool PATCH:: isCauchy() const { notImplemented("isCauchy"); return 0; }
Bool PATCH:: isGreen() const { notImplemented("isGreen"); return 0; }
Bool PATCH:: isVolatile() const { notImplemented("isVolatile"); return 0; }
PATCH* PATCH:: Father() const { notImplemented("Father"); return 0; }
PATCH* PATCH:: FirstSon() const { notImplemented("FirstSon"); return 0; }
PATCH* PATCH:: Partner() const { notImplemented("Partner"); return 0; }
void PATCH:: setPartner(PATCH* /*p*/) { notImplemented("setPartner"); }
int PATCH:: NoOfSons() const { notImplemented("NoOfSons"); return 0; }
//-------------------------------------------------------------------------
Real PATCH:: hMax() const
{
int i;
Real h, maxH = 0.0;
const Vector<EDG*>& edges = getEdges();
FORALL(edges,i)
{
h = edges[i]->hMax();
maxH = Max(h,maxH);
}
return maxH;
}
//-------------------------------------------------------------------------
void PATCH:: setMark() { notImplemented("setMark"); }
void PATCH:: unMark() { notImplemented("unMark"); }
Bool PATCH:: marked() const { notImplemented("marked"); return 0; }
//-------------------------------------------------------------------------
const PATCH* PATCH:: findPatch(const Vector<Real>& /*x*/, Vector<Real>& /*xUnit*/) const
{
notImplemented("findPatch"); return 0;
}
Bool PATCH:: inPatch(const Vector<Real>& /*x*/, Vector<Real>& /*xUnit*/) const
{
notImplemented("inPatch"); return 0;
}
PATCH* PATCH:: findPartner(const Vector<Real>& /*xG*/, int /*maxDepth*/)
{
notImplemented("findPartner"); return 0;
}
//-------------------------------------------------------------------------
void PATCH:: readBCValues(char* boundP, char* classA, BufferedParser& parser)
{
char BCType;
parser.getValue(&BCType);
BCType = toupper(BCType);
//if (BCType == 'I') { *boundP = INTERIOR; *classA = 0; }
// else
if (BCType == 'D') *boundP = DIRICHLET;
else if (BCType == 'N') *boundP = NEUMANN;
else if (BCType == 'C') *boundP = CAUCHY;
else {
cout << "\n*** readBC: Unknown Boundary Type " << BCType << "\n";
cout.flush(); abort();
}
int cl;
parser.getValue(&cl);
if (cl <= 0)
{
cout << "\n*** Boundary condition class must be > 0 (read: "
<< cl << ")\n";
cout.flush(); abort();
}
*classA = cl; // int -> char
}
void PATCH:: readBCValues(char* boundP, char* classA, char* isoP, BufferedParser& parser)
{
char BCType;
parser.getValue(&BCType);
BCType = toupper(BCType);
//if (BCType == 'I') { *boundP = INTERIOR; *classA = 0; }
// else
if (BCType == 'D') *boundP = DIRICHLET;
else if (BCType == 'N') *boundP = NEUMANN;
else if (BCType == 'C') *boundP = CAUCHY;
else {
cout << "\n*** readBC: Unknown Boundary Type " << BCType << "\n";
cout.flush(); abort();
}
int cl;
parser.getValue(&cl);
if (cl <= 0)
{
cout << "\n*** Boundary condition class must be > 0 (read: "
<< cl << ")\n";
cout.flush(); abort();
}
*classA = cl; // int -> char
char is;
if(parser.getSpecValue(&is) == parser.ok)
{
is = toupper(is);
if (is == 'S') *isoP = CIRCULAR;
else if (is == 'C') *isoP = CYLINDRICAL;
else *isoP = PLANE;
}
else *isoP = PLANE;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
//-------------------------- base class MESH ----------------------------
MESH:: MESH(const char* inFileName) : varAllocator(ElementsInBlock*100),
elemIter(0), ptIter(0), edgIter(0)
{
fileName = new char[strlen(inFileName)+1];
strcpy(fileName, inFileName);
maxDepth = 0;
refStep = 0;
innerBoundary = False; if (Cmd.isTrue("innerBoundary")) innerBoundary=True;
infoRefinement = 0; Cmd.get("infoRefinement", &infoRefinement);
timeRefinement = 0; Cmd.get("timeRefinement", &timeRefinement);
accTimeRefinement = 0; Cmd.get("accTimeRefinement", &accTimeRefinement);
accTime = 0.0;
}
//-------------------------------------------------------------------------
MESH:: ~MESH() { delete fileName; }
//-------------------------------------------------------------------------
void MESH:: notImplemented(const char* s) const
{
cout << "\n*** MESH::" << s << " not implemented\n"; cout.flush(); abort();
}
//-------------------------------------------------------------------------
void MESH:: Flatten() { notImplemented("Flatten"); }
//-------------------------------------------------------------------------
void MESH:: writeTriangulation(const char* /*outFileText*/,const Vector<Real>& /*u*/)
const
{
notImplemented("writeTriangulation");
}
//-------------------------------------------------------------------------
void MESH:: missingItem(const char* item, const char* fileName0)
{
cout << "\n*** File " << fileName0 << " without " << item << "\n";
cout.flush(); abort();
}
//-------------------------------------------------------------------------
void MESH:: pointIndexError(int no)
{
cout << "\n*** readTriangulation: point index too high: " << no << "\n";
cout.flush(); abort();
}
//-------------------------------------------------------------------------
void MESH:: checkMaxIndex(int maxIndex)
{
if (maxIndex < 0) {
cout << "\n*** readTriangulation: maxIndex must be >= 0 (read: "
<< maxIndex << ")\n";
cout.flush(); abort();
}
}
//-------------------------------------------------------------------------
void MESH:: boundaryWarning(int pNo1, int pNo2, int pNo3)
{
cout << "\n*** Warning --- readTriangulation: boundary element specified"
<< "\n\tthat is not on exterior boundary; points: " << pNo1;
if (pNo2) cout << pNo2;
if (pNo3) cout << pNo3;
cout << "\n";
}
//-------------------------------------------------------------------------
void MESH:: timeInfo(Timer& timer, Timer& accTimer)
{
if (timeRefinement) { cout << "\n\tGrid Refinement: "; timer.cpu(); }
if (accTimeRefinement)
{
accTime += accTimer.cpu(False);
cout << Form("\tAccumulated time: %1.2f sec.\n", accTime);
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void MESH:: resetEdgIter()
{
if (edgIter != 0)
{
cout << "\n*** resetEdgIter: edgIter != 0"
<< " \t (previous operation not completed?)\n";
cout.flush(); abort();
}
edgIterLevel = 0;
}
//-------------------------------------------------------------------------
void MESH:: resetPtIter()
{
if (ptIter != 0)
{
cout << "\n*** resetPtIter: ptIter != 0"
<< " \t (previous operation not completed?)\n";
cout.flush(); abort();
}
ptIterLevel = 0;
}
//-------------------------------------------------------------------------
const PATCH* MESH:: findPatch(const Vector<Real>& /*x*/, Vector<Real>& /*xUnit*/,
const PATCH* /*newPatch*/) const
{
notImplemented("findPatch"); return 0;
}
//-------------------------------------------------------------------------
void MESH:: patchNotFound() const
{
cout << "\n*** MESH: patch not found\n"; cout.flush(); abort();
}
//-------------------------------------------------------------------------
// -- set the partners of depth 0
void MESH:: setInitialPartners(MESH* prevMesh)
{
PATCH* patch, *partner;
if (prevMesh==0) return;
prevMesh->resetElemIter();
resetElemIter();
while ((patch=elemIterAll()))
{
partner = prevMesh->elemIterAllOfHistory();
patch->setPartner(partner);
}
// -- dummy operation on rest of elements:
while ((partner=prevMesh->elemIterAllOfHistory())) { }
}
//-------------------------------------------------------------------------
void MESH:: setPartners()
{
Timer timer;
int k, noOfSons;
PATCH* fPartner, *p, *son, *father, *fPartnerSon;
const int lastStep = True;
resetElemIter(lastStep);
while ((p = elemIterAllOfHistory()))
{
if (p->Partner() == 0)
{
father = p->Father();
fPartner = father->Partner();
son = p;
noOfSons = father->NoOfSons();
if (father->isGreen() || fPartner->isGreen() ||
!fPartner->refined() || fPartner->Depth() < father->Depth())
{
for (k=1; k<=noOfSons; ++k)
{
son->setPartner(fPartner);
son = son->Next();
}
}
else
{
fPartnerSon = fPartner->FirstSon();
for (k=1; k<=noOfSons; ++k)
{
son->setPartner(fPartnerSon);
son = son->Next();
fPartnerSon = fPartnerSon->Next();
}
}
}
}
if (Cmd.isTrue("TimeTransport"))
{ cout << "\n\tFind Partners: "; timer.cpu(); }
}
//-------------------------------------------------------------------------
void MESH:: newNodeStack(PT* p, int depth, int targetDepth)
{
static const int nodeStackSize = sizeof(NodeStack);
if (depth >= targetDepth)
{
cout << "\n***MESH:: newNodeStack: depth >= targetDepth ("
<< depth << ">=" << targetDepth << ")\n";
cout.flush(); abort();
}
p->nodeStack = (NodeStack*) varAllocator.Get(nodeStackSize);
p->nodeStack->init(&varAllocator, depth, targetDepth);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void NodeStack:: init(VarSizeAllocator* varAlloc0, int depth, int targetDepth)
{
static const int intSize = sizeof(int);
varAlloc = varAlloc0;
l = depth;
h = depth-1;
top = targetDepth;
v = (int*) varAlloc->Get(intSize*(top-l+1));
v -= l;
}
//-------------------------------------------------------------------------
void NodeStack:: push(int node)
{
if (h==top) extend();
v[++h] = node;
}
//-------------------------------------------------------------------------
void NodeStack:: extend()
{
static const int intSize = sizeof(int);
top += top;
int* vnew = (int*) varAlloc->Get(intSize*(top-l+1));
vnew -= l;
for (int i=l; i<=h; ++i) vnew[i] = v[i];
// delete v; was allocated from varAlloc
v = vnew;
}
//-------------------------------------------------------------------------
void NodeStack:: checkBounds(int i) const
{
if (i<l || i>h)
{
cout << "\n*** Vector index " << i << " out of range\n";
cout.flush(); abort();
}
}
//-------------------------------------------------------------------------
void NodeStack:: print() const
{
int i;
for (i=l; i<=h; ++i) cout << v[i] << " ";
cout << "\n";
}
syntax highlighted by Code2HTML, v. 0.9.1