/*
$Id: triang1.cc,v 1.5 1996/11/06 12:26:19 roitzsch Exp $
(C)opyright 1996 by Konrad-Zuse-Center, Berlin
All rights reserved.
Part of the Kaskade distribution
*/
#include "triang1.h"
#include "utils.h"
#include "numerics.h"
#include "mzibutil.h"
#if NOKASKPLOT == 0
#include "feplot1.h"
#endif
#include "cmdpars.h"
extern CmdPars Cmd;
Vector<PT*> EDG1::p(2);
Vector<EDG*> EDG1::e(1);
Vector<PATCH*> EDG1::f(2);
MESH1:: MESH1(const char* inFileName, Bool readMesh)
: MESH(inFileName), ptList(0,25), edgList(0,25),
ptAlloc(ElementsInBlock), edgAlloc(ElementsInBlock)
{
noOfPoints = 0;
noOfEdges = 0;
if (readMesh) readTriangulation(fileName);
}
//-------------------------------------------------------------------------
MESH1:: ~MESH1()
{
int i;
FORALL( ptList,i) delete ptList[i];
FORALL(edgList,i) delete edgList[i];
}
//-------------------------------------------------------------------------
FEPlot* MESH1:: newFEPlot(int plotType, char* caption, float size)
{
#if NOKASKPLOT == 0
return new FEPlotMESH1(this, plotType, caption, size);
#else
return 0;
#endif
}
//-------------------------------------------------------------------------
void MESH1:: readTriangulation(const char* fileName0)
{
int i, maxIndex;
DList<PT1>* ptL = new DList<PT1>; ptList.push(ptL);
DList<EDG1>* edgL = new DList<EDG1>; edgList.push(edgL);
BufferedParser parser(fileName0, "end", CommentFlag);
if (parser.findWord("Poin") != parser.ok) missingItem("points",fileName0);
if (parser.findWord("max") != parser.ok) missingItem("maxIndex",fileName0);
parser.getValue(&maxIndex);
checkMaxIndex(maxIndex);
Vector<PT1*> ptAdr(0,maxIndex); // for point addresses
FORALL(ptAdr,i) ptAdr[i] = 0;
readPoints(parser, ptAdr);
if (parser.findWord("Elem") != parser.ok) missingItem("triangles",fileName0);
readElements(parser, ptAdr);
if (parser.findWord("Boun") != parser.ok) missingItem("boundary",fileName0);
readBoundary(parser, ptAdr);
EDG1* e;
DListIter<EDG1> iter(edges());
while ((e=iter.all())) e->setBoundary();
if (infoRefinement) Info();
if (Cmd.isTrue("printMesh")) cout << *this;
}
//-------------------------------------------------------------------------
void MESH1:: readPoints(BufferedParser& parser, Vector<PT1*>& ptAdr)
{
int ptNo, count=0;
PT1* p;
while (parser.getValue(&ptNo) == parser.ok)
{
if (ptNo > ptAdr.high()) pointIndexError(ptNo);
p = ptAlloc.Get(); //new PT1;
ptAdr[ptNo] = p;
parser.getValue(&p->x);
p->node = ++count; // used in Find...
ptList[0]->add(p);
++noOfPoints;
}
}
//-------------------------------------------------------------------------
void MESH1:: readElements(BufferedParser& parser, Vector<PT1*>& ptAdr)
{
int pNo1, pNo2, classA;
EDG1 *e;
while (parser.getValue(&pNo1) == parser.ok)
{
parser.getValue(&pNo2);
e = getEDG1(); // new EDG1;
e->p1 = ptAdr[pNo1];
e->p2 = ptAdr[pNo2];
setRightHand(e);
parser.getValue(&classA);
e->classA = classA;
SetEP(e);
edgList[0]->add(e);
++noOfEdges;
}
}
//-------------------------------------------------------------------------
void MESH1:: writeTriangulation(const char* outFileText, const Vector<Real>& u)
const
{
int i;
int node;
EDG1 *ed;
PT1 *p;
FILE *outFile;
Vector<int> PNodes(noOfPoints);
// saving the node numbers for later reuse
i = 0;
DListIter<PT1> PIter(points());
while ((p=PIter.all())) PNodes[++i] = p->node;
char* outText = new char[strlen(outFileText)+10];
strcpy(outText, outFileText);
strcpy(strchr(outText,'.'),".zib.geo");
outFile = fopen(outText, "w");
if (outFile==0)
{
cout << "\n*** writeTriangulation: could not open file " << outText
<< "\n";
cout.flush(); abort();
}
fprintf(outFile, "# %s\n\n", fileName);
fprintf(outFile, "Points: maxIndex %d\n", noOfPoints);
node = 0;
DListIter<PT1> iterPts(points());
while ((p=iterPts.all())) {
p->node = node+1;
fprintf(outFile, " %d %e\n", p->node, p->x);
node++;
}
fprintf(outFile, "END\n");
fprintf(outFile, "\nElements:\n");
DListIter<EDG1> iterEds(edges());
while ((ed=iterEds.all())) {
fprintf(outFile, "%d %d %d\n",
(ed->p1)->node, (ed->p2)->node, ed->classA );
}
fprintf(outFile, "END\n");
fprintf(outFile, "\nBoundary: \n");
iterPts.reset();
while ((p=iterPts.all())) {
if (!p->onBoundary()) continue;
if (p->isDirichlet())
fprintf(outFile, "%d d %d\n",
p->node, p->classA );
else if (p->isNeumann())
fprintf(outFile, "%d n %d\n",
p->node, p->classA );
else if (p->isCauchy())
fprintf(outFile, "%d c %d\n",
p->node, p->classA );
}
fprintf(outFile, "END\n");
// output of values at the nodes located in the points
fprintf(outFile, "\nValues in nodes: \n");
i = 0;
iterPts.reset();
while ((p=iterPts.all())) fprintf(outFile, "%e\n", u[PNodes[++i]]);
fprintf(outFile, "END\n");
cout << "\n WriteZIB-Data ok.; written to " << outText << "\n";
delete outText;
fclose(outFile);
// restore of node numbers
i = 0;
PIter.reset();
while ((p=PIter.all())) p->node = PNodes[++i];
}
//-------------------------------------------------------------------------
void MESH1:: setRightHand(EDG1* e)
{
if (e->p2->x < e->p1->x)
{
PT1* p = e->p2;
e->p2 = e->p1;
e->p1 = p;
}
}
//-------------------------------------------------------------------------
void MESH1:: readBoundary(BufferedParser& parser, Vector<PT1*>& ptAdr)
{
int pNo;
PT1 *p;
while (parser.getValue(&pNo) == parser.ok)
{
p = ptAdr[pNo];
if (!p->onBoundary()) boundaryWarning(pNo);
p->readBC(parser);
}
}
//-------------------------------------------------------------------------
void MESH1:: Refine()
{
EDG1* e;
int noPt[2], noEdg[2];
Timer timer, accTimer;
noPt[0] = noOfPoints;
noEdg[0] = noOfEdges;
DList<EDG1>* edgL = new DList<EDG1>;
edgList.push(edgL);
DListIter<EDG1> iter(edges());
while ((e=iter.all())) { if (e->mark) RefineEdge(e); }
noPt[1] = noOfPoints;
noEdg[1] = noOfEdges;
++refStep;
if (infoRefinement) refinementInfo (noPt, noEdg);
timeInfo(timer, accTimer);
if (Cmd.isTrue("printMesh")) cout << *this;
}
//-------------------------------------------------------------------------
void MESH1:: Flatten()
{
int noReleasedEdges=0;
int noPoints=0, noEdges=0;
int k;
EDG1* ed;
PT1* p;
DList<EDG1>* edgL = new DList<EDG1>;
DList<PT1>* ptL = new DList<PT1>;
Stack<EDG1*> edgToRetain(25), edgToRelease(25);
Stack<PT1*> ptToRetain(25);
DListIter<EDG1> iterEdg(edges());
DListIter<PT1> iterPt(points());
while ((ed=iterEdg.allOfHistory()))
{
if (!ed->refined())
{
edgToRetain.push(ed);
noEdges++;
}
else
{
edgToRelease.push(ed);
noReleasedEdges++;
}
}
while ((p=iterPt.allOfHistory()))
{
ptToRetain.push(p);
noPoints++;
}
while (!edgList.empty()) { delete edgList.pop(); }
edgList.push(edgL);
while (!ptList.empty()) { delete ptList.pop(); }
ptList.push(ptL);
for (k=ptToRetain.l; k<=ptToRetain.h; k++)
{
p = ptToRetain[k];
ptL->add(p);
p->depth = 0;
}
for (k=edgToRetain.l; k<=edgToRetain.h; k++)
{
ed = edgToRetain[k];
ed->father = nil;
ed->depth = 0;
edgL->add(ed);
}
while (!edgToRelease.empty())
{
ed = edgToRelease.pop();
edgAlloc.Return(ed);
}
refStep = 0;
maxDepth = 0;
ptIterLevel = 0;
iterLevel = 0;
noOfPoints = noPoints;
noOfEdges = noEdges;
if (infoRefinement)
{
cout << "\n Flatten: " << noReleasedEdges << " edges released\n";
cout << " " << noPoints << " points, " <<
noEdges << " edges retained\n";
}
return;
}
//-------------------------------------------------------------------------
Bool MESH1:: RefineEdge(EDG1 *ed)
{
EDG1 *e1, *e2;
PT1 *pm;
e1 = getEDG1();
e2 = getEDG1();
++noOfEdges;
pm = ptAlloc.Get(); // new PT1; assert(pm != 0);
++noOfPoints;
ed->pm = pm;
pm->depth = ed->depth + 1;
pm->x = 0.5*(ed->p1->x + ed->p2->x);
e1->p1 = ed->p1;
e1->p2 = pm;
e2->p1 = pm;
e2->p2 = ed->p2;
e1->father = e2->father = ed;
e1->type = e2->type = ed->type;
e1->classA = e2->classA = ed->classA;
e1->depth = e2->depth = ed->depth + 1;
ed->firstSon = e1;
DelEP(ed);
SetEP(e1);
SetEP(e2);
edgList[edgList.h]->add(e1);
edgList[edgList.h]->add(e2);
if (pm->depth > ptList.h)
{
DList<PT1>* ptL = new DList<PT1>;
ptList.push(ptL);
}
ptList[pm->depth]->add(pm);
if (e1->depth > maxDepth) maxDepth = e1->depth;
return True;
}
//-------------------------------------------------------------------------
// -- SetEP(e): set edge e on its points
void MESH1:: SetEP(EDG1 *e)
{
PT1 *p1, *p2;;
int k = 0;
p1 = e->p1;
if (p1->e1 == 0) { p1->e1 = e; ++k; }
else if (p1->e2 == 0) { p1->e2 = e; ++k; }
p2 = e->p2;
if (p2->e1 == 0) { p2->e1 = e; ++k; }
else if (p2->e2 == 0) { p2->e2 = e; ++k; }
if (k != 2)
{
cout << "\n*** SetEP: Edge with impossible state, k = " << k << "\n";
printf("\n e->p1: %p",e->p1);
printf("\n e->p2: %p",e->p2);
e->print(cout);
cout.flush(); abort();
}
if (p1->onBoundary() || p2->onBoundary()) e->boundP = BOUNDARY;
}
//-------------------------------------------------------------------------
// -- DelEP(e): delete edge e on its points
void MESH1:: DelEP(EDG1 *e)
{
PT1* p;
int k = 0;
p = e->p1;
if (p->e1 == e) { p->e1 = 0; ++k; }
else if (p->e2 == e) { p->e2 = 0; ++k; }
p = e->p2;
if (p->e1 == e) { p->e1 = 0; ++k; }
else if (p->e2 == e) { p->e2 = 0; ++k; }
if (k != 2)
{
cout << "\n*** DelEP: impossible state, k = " << k << ")\n";
cout.flush(); abort();
}
}
//-------------------------------------------------------------------------
void MESH1:: Info()
{
cout << "\n Triangulation: points edges "
<< " Depth=" << maxDepth << " Steps=" << refStep
<< Form("\n%15s%10d%10d", " ", noOfPoints, noOfEdges)
<< "\n";
}
//-------------------------------------------------------------------------
int MESH1:: MemSpace()
{
int space = varAllocator.MemSpace();
space += edgAlloc.MemSpace();
space += ptAlloc.MemSpace();
return space;
}
//-------------------------------------------------------------------------
void MESH1:: refinementInfo(int* noPt, int* noEdg)
{
int blocks = varAllocator.NoOfBlocks();
blocks += edgAlloc.NoOfBlocks();
blocks += ptAlloc.NoOfBlocks();
int space = MemSpace();
/*
FORALL (edgList,depth)
{
space += edgList[depth]->noOfElements * sizeof(EDG1);
if (depth <= ptList.h)
space += ptList[depth]->noOfElements * sizeof(PT1);
}
*/
cout << "\n Refine: points edges "
<< " Depth=" << maxDepth << " Steps=" << refStep
<< Form("\n%15s%10d%10d", "previous", noPt[0], noEdg[0]);
if (infoRefinement > 1)
{
cout << Form("\n%15s%10d%10d", "generated",
noPt[1]-noPt[0], noEdg[1]-noEdg[0]);
}
cout << Form("\n%15s%10d%10d Space : %1.6f MB","total",
noOfPoints, noOfEdges, float(space)/1e6);
cout << Form("\n%15s%10.2f%10.2f Blocks: %1d\n","ratio",
float(noOfPoints)/noPt[0], float(noOfEdges)/noEdg[0], blocks);
if (infoRefinement > 2)
{
int npDepth, step, np=0, ne=0;
cout << "\n Grid Statistics:"
<< Form("\n%15s%10s%10s%10s", "depth/step", "points", "edges",
"(total)");
FORALL (edgList,step)
{
npDepth = (step<=ptList.h)? ptList[step]->noOfElements : 0;
cout << Form("\n%15d%10d%10d", step,
npDepth, edgList[step]->noOfElements);
np += npDepth;
ne += edgList[step]->noOfElements;
}
cout << Form("\n%15s%10d%10d", "total sum", np, ne);
cout << "\n";
}
}
//-------------------------------------------------------------------------
ostream& operator << (ostream& os, const MESH1& t)
{
PT1* p;
EDG1* e;
os << "\nTriangulation";
if (t.fileName) os << " " << t.fileName;
os <<":\n";
DListIter<PT1> iterPT (t.points());
DListIter<EDG1> iterEDG(t.edges());
while ((p = iterPT.all())) p->print(cout);
while ((e = iterEDG.all())) e->print(cout);
os << "\nRefinement step " << t.refStep;
os << " -- points: " << t.noOfPoints;
os << " edges: " << t.noOfEdges;
os << "\n";
return os;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void PT1:: print(ostream& os) const
{
os << "Point " << node << " BC " << (int) boundP << " class "
<< (int) classA << " mark " << (int) mark
<< " x: " << x << "\n";
}
void PT1:: reset()
{
boundP = INTERIOR;
mark = classA = depth = NUL;
node = 0;
e1 = e2 = 0;
}
//-------------------------------------------------------------------------
void PT1:: getNodeCoord(Vector<Real>& v) const { v[1] = x; }
//-------------------------------------------------------------------------
void PT1:: getCoordinates(Vector<Real>& v) const { v[1] = x; }
//-------------------------------------------------------------------------
void PT1:: readBC(BufferedParser& parser)
{
readBCValues(&boundP, &classA, parser);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void EDG1:: compJ(Jacobian& j, Bool /*checkDetJ*/) const
{
j.J(1,1) = p2->x - p1->x;
j.detJ = j.J(1,1);
}
void EDG1:: compJinv(Jacobian& j) const
{
compJ(j);
j.Jinv(1,1) = 1./j.J(1,1);
}
//-------------------------------------------------------------------------
void EDG1:: print(ostream& os) const
{
os << "Edge with points " << p1->node << " " << p2->node
<< " type " << (int) type
<< " class " << (int) classA << " depth " << (int) depth << "\n";
}
void EDG1:: reset()
{
mark = classA = type = depth = NUL;
boundP = INTERIOR;
p1 = p2 = pm = 0;
node = 0;
father = firstSon = 0;
error = 0.0;
}
void EDG1:: setBoundary()
{
if (p1->onBoundary() || p2->onBoundary()) boundP = BOUNDARY;
}
//-------------------------------------------------------------------------
void EDG1:: getNeighbours(Vector<PATCH*>& neighbours) const
{
if (p1->e1 != this) neighbours[1] = p1->e1;
else neighbours[1] = p1->e2;
if (p2->e1 != this) neighbours[2] = p2->e1;
else neighbours[2] = p2->e2;
}
//-------------------------------------------------------------------------
// -- point at unit coordinate -> point at x
void EDG1:: realCoordinates(const Vector<Real>& u, Vector<Real>& x) const
{
Real L1 = 1.0 - u[1];
Real L2 = u[1];
x[1] = L1*p1->x + L2*p2->x;
}
//-------------------------------------------------------------------------
// -- point at x -> unit coordinates
void EDG1:: unitCoordinates(const Vector<Real>& x, Vector<Real>& u) const
{
u[1] = (x[1] - p1->x) / (p2->x - p1->x);
}
//-------------------------------------------------------------------------
void EDG1:: getNodeCoord(Vector<Real>& v) const
{
v[1] = 0.5*(p1->x + p2->x);
}
//-------------------------------------------------------------------------
void EDG1:: getCoordinates(Vector<Real>& cP1, Vector<Real>& cP2) const
{
p1->getCoordinates(cP1);
p2->getCoordinates(cP2);
}
void EDG1:: compNormalVector(Vector<Real>& /*n*/, Real* /*length*/) const
{
notImplemented("compNormalVector");
}
//-------------------------------------------------------------------------
Real EDG1:: hMax() const { return Abs(p2->x - p1->x); }
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
void MESH1:: resetElemIter(Bool lastStep)
{
if (elemIter != 0)
{
cout << "\n*** resetElemIter: elemIter != 0"
<< " \t (previous operation not completed?)\n";
cout.flush(); abort();
}
if (lastStep) iterLevel = edgList.high();
else iterLevel = 0;
}
//-------------------------------------------------------------------------
PATCH* MESH1:: elemIterAll()
{
if (elemIter != 0) // iterate rest of list
{
for (elemIter=elemIter->Next(); elemIter; elemIter=elemIter->Next())
if (!elemIter->refined()) return elemIter;
if (++iterLevel > edgList.high()) return 0;
}
while (1) // iterate new lists on stack
{
for (elemIter=edgList[iterLevel]->first; elemIter;
elemIter=elemIter->Next())
if (!elemIter->refined()) return elemIter;
if (++iterLevel > edgList.high()) return 0;
}
}
//-------------------------------------------------------------------------
PATCH* MESH1:: elemIterAllOfHistory()
{
if (elemIter != 0) // iterate rest of list
{
for (elemIter=elemIter->Next(); elemIter; elemIter=elemIter->Next())
return elemIter;
if (++iterLevel > edgList.high()) return 0;
}
while (1) // iterate new lists on stack
{
for (elemIter=edgList[iterLevel]->first; elemIter;
elemIter=elemIter->Next()) return elemIter;
if (++iterLevel > edgList.high()) return 0;
}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
PATCH* MESH1:: edgIterAll()
{
if (edgIter != 0) // iterate rest of list
{
for (edgIter=edgIter->Next(); edgIter; edgIter=edgIter->Next())
if (!edgIter->refined()) return edgIter;
if (++edgIterLevel > edgList.high()) return 0;
}
while (1) // iterate new lists on stack
{
for (edgIter=edgList[edgIterLevel]->first; edgIter;
edgIter=edgIter->Next())
if (!edgIter->refined()) return edgIter;
if (++edgIterLevel > edgList.high()) return 0;
}
}
//-------------------------------------------------------------------------
PATCH* MESH1:: edgIterAllOfHistory()
{
if (edgIter != 0) // iterate rest of list
{
for (edgIter=edgIter->Next(); edgIter; edgIter=edgIter->Next())
return edgIter;
if (++edgIterLevel > edgList.high()) return 0;
}
while (1) // iterate new lists on stack
{
for (edgIter=edgList[edgIterLevel]->first; edgIter;
edgIter=edgIter->Next()) return edgIter;
if (++edgIterLevel > edgList.high()) return 0;
}
}
//-------------------------------------------------------------------------
PATCH* MESH1:: ptIterAll()
{
if (ptIter != 0) // iterate rest of list
{
for (ptIter=ptIter->Next(); ptIter; ptIter=ptIter->Next())
if (!ptIter->refined()) return ptIter;
if (++ptIterLevel > ptList.high()) return 0;
}
while (1) // iterate new lists on stack
{
for (ptIter=ptList[ptIterLevel]->first; ptIter;
ptIter=ptIter->Next())
if (!ptIter->refined()) return ptIter;
if (++ptIterLevel > ptList.high()) return 0;
}
}
syntax highlighted by Code2HTML, v. 0.9.1